unified_real/
lib.rs

1use computable_real::Real as ComputableReal;
2use num::traits::Inv;
3use num::{BigInt, BigRational, One, Signed, ToPrimitive, Zero};
4use std::convert::From;
5use std::str::FromStr;
6
7pub mod property;
8pub use property::{Kind, Property};
9
10pub mod real;
11pub use real::Real;
12
13pub mod traits;
14
15const INIT_CMP_TOLERANCE: i32 = -100;
16const REL_CMP_TOLERANCE: i32 = -1000;
17const CMP_TOLERANCE: i32 = -3500;
18const ZERO_CMP_TOLERANCE: i32 = -5000;
19
20const EXTRACT_SQUARE_MAX: u64 = 43;
21const EXTRACT_SQUARE_MAX_LEN: u64 = 5000;
22
23impl From<Real> for ComputableReal {
24    fn from(value: Real) -> Self {
25        if value.rat.is_zero() {
26            ComputableReal::zero()
27        } else if value.rat.is_one() {
28            value.cr
29        } else {
30            value.cr * value.rat.into()
31        }
32    }
33}
34
35fn exact_log(n: &BigInt, base: u32) -> Option<i64> {
36    if !n.is_positive() {
37        return None;
38    }
39    if n.is_one() {
40        return Some(0);
41    }
42    if base == 1 {
43        return None;
44    }
45
46    let double = n.to_f64().unwrap_or(f64::INFINITY);
47    let approx: f64 = double.log(base as f64);
48    if double.is_infinite() {
49        // Rule out some simple cases
50        if (base % 2 != 0 && !n.bit(0))
51            || (base % 3 != 0 && n % 3 == BigInt::ZERO)
52            || (base % 5 != 0 && n % 5 == BigInt::ZERO)
53        {
54            return None;
55        }
56    } else if (approx - approx.round()).abs() > 1e-6 {
57        // If the floating point log is not close to an integer, then n dooesn't
58        //  have an exact log.
59        return None;
60    }
61
62    let mut res = 0;
63    let mut powers: Vec<BigInt> = vec![BigInt::from(base)];
64    let mut n_reduced = n.clone();
65
66    // We build a power table, where powers[i] = base^(2^i)
67    loop {
68        let last = powers.last().unwrap();
69        let next = last * last;
70        if next.bits() > n_reduced.bits() {
71            break;
72        }
73        let q = &n_reduced / &next;
74        let r = n_reduced % &next;
75        if !r.is_zero() {
76            return None;
77        }
78
79        res += 1 << powers.len();
80        powers.push(next);
81        n_reduced = q;
82    }
83
84    // Using the power table, we divide n_reduced down to 1
85    for (i, power) in powers.iter().rev().enumerate() {
86        if power.bits() <= n_reduced.bits() {
87            let q = &n_reduced / power;
88            let r = n_reduced % power;
89            if !r.is_zero() {
90                return None;
91            }
92            res += 1 << i;
93            n_reduced = q;
94        }
95    }
96
97    if n_reduced.is_one() {
98        Some(res)
99    } else {
100        None
101    }
102}
103
104/// Returns true if the trigonometric function of the given argument can be
105/// simplified, or reduced to the (0, 1/2) interval.
106fn reducible_trig_arg(n: &BigRational) -> bool {
107    !n.is_positive()
108        || n >= &BigRational::new(1.into(), 2.into())
109        || n == &BigRational::new(1.into(), 3.into())
110        || n == &BigRational::new(1.into(), 4.into())
111        || n == &BigRational::new(1.into(), 6.into())
112}
113
114// TODO: comment
115fn irreducible_sqrt_arg(n: &BigRational) -> bool {
116    n.numer().bits() <= 30
117        && n.numer().abs() <= EXTRACT_SQUARE_MAX.into()
118        && n.denom().bits() <= 30
119        && n.denom().abs() <= EXTRACT_SQUARE_MAX.into()
120}
121
122fn nth_root(num: &BigInt, n: i32) -> Option<BigInt> {
123    if num.is_negative() {
124        if n & 1 == 0 {
125            return None;
126        } else {
127            return nth_root(&-num, n).map(|r| -r);
128        }
129    }
130    if num.is_zero() {
131        return Some(BigInt::ZERO);
132    }
133
134    let mut cr = ComputableReal::int(num.clone());
135    cr = if n == 2 {
136        cr.sqrt()
137    } else if n == 4 {
138        cr.sqrt().sqrt()
139    } else {
140        (cr.ln() / ComputableReal::int(n.into())).exp()
141    };
142
143    let scale = -10;
144    let appr = cr.appr(scale).value;
145    let mask: BigInt = ((1 << -scale) - 1).into();
146    let bits = &appr & &mask;
147
148    if !bits.is_zero() && bits != mask {
149        return None;
150    }
151
152    let res = if bits.is_zero() {
153        appr >> -scale
154    } else {
155        (appr + 1) >> -scale
156    };
157
158    if res.pow(n as u32) == *num {
159        Some(res)
160    } else {
161        None
162    }
163}
164
165/// Returns i1 and i2 such that n = i1^2 * i2.
166fn int_extract_square(n: &BigInt) -> (BigInt, BigInt) {
167    const PRIMES: [u64; 6] = [2, 3, 5, 7, 11, 13];
168
169    let mut square = BigInt::one();
170    let mut n = n.clone();
171    if n.bits() > EXTRACT_SQUARE_MAX {
172        return (square, n);
173    }
174
175    for &prime in &PRIMES {
176        if n.is_one() {
177            break;
178        }
179        loop {
180            let q = &n / prime;
181            let r = &n % prime;
182            if r.is_zero() {
183                n = q;
184                square *= prime;
185            } else {
186                break;
187            }
188        }
189    }
190
191    for i in 1..=10 {
192        let i = BigInt::from(i);
193        let q = &n / &i;
194        let r = &n % &i;
195        if r.is_zero() {
196            if let Some(root) = nth_root(&q, 2) {
197                n = i;
198                square *= root;
199                break;
200            }
201        }
202    }
203
204    (square, n)
205}
206
207/// Returns r1 and r2 such that n = r1^2 * r2.
208pub fn extract_square(n: &BigRational) -> (BigRational, BigRational) {
209    if n.is_zero() {
210        return (BigRational::zero(), BigRational::one());
211    }
212    let (num_square, mut num_n) = int_extract_square(&n.numer().abs());
213    let (den_square, den_n) = int_extract_square(&n.denom().abs());
214    if n.is_negative() {
215        num_n = -num_n;
216    }
217    (
218        BigRational::new(num_square, den_square),
219        BigRational::new(num_n, den_n),
220    )
221}
222
223/// Returns a non-zero r such that a=b^r, or None if no such r exists.
224/// Effectively, we find a common power of a and b, and return the exponents
225/// as a rational number.
226fn common_power(a: &BigInt, b: &BigInt) -> Option<BigRational> {
227    const MAX_LENGTH: u64 = 200;
228    if !a.is_positive() || !b.is_positive() {
229        return None;
230    }
231    if a == b {
232        Some(BigRational::one())
233    } else if a < b {
234        common_power(b, a).map(|r| r.inv())
235    } else if a == &BigInt::one() || b == &BigInt::one() {
236        None
237    } else if a.bits() > MAX_LENGTH {
238        None
239    } else {
240        let q = a / b;
241        let r = a % b;
242        if !r.is_zero() {
243            None
244        } else {
245            let pow = common_power(&q, b)?;
246            Some(pow + BigRational::one())
247        }
248    }
249}
250
251fn have_common_power(a: &BigRational, b: &BigRational) -> bool {
252    if !a.is_positive() || !b.is_positive() {
253        return false;
254    }
255
256    let na = a.numer();
257    let da = a.denom();
258    let nb = b.numer();
259    let db = b.denom();
260
261    if da == &BigInt::one() {
262        if db == &BigInt::one() {
263            return common_power(&na, &nb).is_some();
264        } else if nb == &BigInt::one() {
265            return common_power(&na, &db).is_some();
266        }
267    } else if na == &BigInt::one() {
268        if db == &BigInt::one() {
269            return common_power(&da, &nb).is_some();
270        } else if nb == &BigInt::one() {
271            return common_power(&da, &db).is_some();
272        }
273    }
274
275    let num_a_num_b = common_power(&na, &nb);
276    let num_a_den_b = common_power(&na, &db);
277    if let Some(power) = num_a_num_b {
278        if let Some(power2) = common_power(&da, &db) {
279            if power == power2 {
280                return true;
281            }
282        }
283    }
284    if let Some(power) = num_a_den_b {
285        if let Some(power2) = common_power(&da, &nb) {
286            if power == power2 {
287                return true;
288            }
289        }
290    }
291    false
292}
293
294/// Approximates the number of bits to the left of the binary point.
295/// Negative value indicates leading zeros to the right of the binary
296fn estimate_whole_bits(n: &BigRational) -> Option<i32> {
297    if n.is_zero() {
298        None
299    } else {
300        Some(n.numer().bits() as i32 - n.denom().bits() as i32)
301    }
302}
303
304/// Returns a string representation of the rational number, truncated to n
305/// decimal places. Rounded to zero.
306pub fn rational_truncated(r: &BigRational, n: u32) -> String {
307    let mut digits = (r.numer().abs() * BigInt::from(10).pow(n) / r.denom().abs()).to_string();
308    if digits.len() <= n as usize {
309        digits = format!("{:0>1$}", digits, n as usize + 1);
310    }
311    let len = digits.len();
312    let sign = if r.is_negative() { "-" } else { "" };
313    let dot = if n == 0 { "" } else { "." };
314    format!(
315        "{}{}{}{}",
316        sign,
317        &digits[..len - n as usize],
318        dot,
319        &digits[len - n as usize..]
320    )
321}
322
323/// Returns the number of digits to the right of the decimal point required to
324/// represent the rational exactly, if possible.
325pub fn rational_digits_required(n: &BigRational) -> Option<u32> {
326    if n.is_integer() {
327        return Some(0);
328    }
329    if n.denom().bits() > 10000 {
330        return None;
331    }
332    let mut denom = n.denom().clone();
333
334    let mut two_powers = 0;
335    while !denom.bit(0) {
336        denom >>= 1;
337        two_powers += 1;
338    }
339
340    let mut five_powers = 0;
341    while &denom % 5 == BigInt::zero() {
342        denom /= 5;
343        five_powers += 1;
344    }
345
346    if !denom.abs().is_one() {
347        return None;
348    }
349
350    Some(two_powers.max(five_powers))
351}
352
353/// Constructs a BigRational from a string representation.
354///
355/// This is needed because if we use a floating point number, we lose precision.
356///
357/// Can be in scientific notation, e.g., "1.23e-4" or standard notation,
358/// e.g., "0.1234".
359pub fn parse_rational(s: &str) -> Option<BigRational> {
360    let s = s.trim();
361    if s.is_empty() {
362        return None;
363    }
364
365    if let Some((base, exp)) = s.split_once('e').or_else(|| s.split_once('E')) {
366        let base = parse_rational(base)?;
367        let exp = i32::from_str(exp).ok()?;
368        let factor = BigRational::from(BigInt::from(10).pow(exp.abs() as u32));
369        return if exp.is_positive() {
370            Some(base * factor)
371        } else {
372            Some(base / factor)
373        };
374    }
375
376    if let Some((whole, frac)) = s.split_once('.') {
377        let frac_len = frac.len();
378        let whole = BigInt::from_str(whole).ok()?;
379        let frac = BigInt::from_str(frac).ok()?;
380        let scale = BigInt::from(10).pow(frac_len as u32);
381        return Some(BigRational::new(whole * &scale + frac, scale));
382    }
383
384    BigInt::from_str(s).ok().map(|n| BigRational::from(n))
385}
386
387pub fn reduce_trig_arg(n: &BigRational) -> BigRational {
388    if n >= &BigRational::new((-1).into(), 2.into()) && n < &BigRational::new(3.into(), 2.into()) {
389        return n.clone();
390    }
391    let floor = (n + BigRational::one()).floor().to_integer();
392    let offset = floor & !BigInt::one();
393    n - BigRational::from_integer(offset)
394}
unified_real/lib.rs

unified_real/
lib.rs
