symplex 0.23.0

Exact symbolic mathematics for Rust: calculus, summation, solving, linear algebra, transforms, compile-time dimensional analysis, and Rust/C code generation
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

//! Arbitrary-precision complex numbers as `(re, im)` pairs of
//! [`BigFloat`]s, with the field operations every numeric algorithm needs
//! (`evalf`, the polynomial root finders, algebraic-number verification).
//!
//! This is the arbitrary-precision analogue of [`num_complex::Complex64`],
//! which the public API uses for `f64` results (`eval_complex64`,
//! `nroots`).  It stays a tuple because `num_complex::Complex<BigFloat>`
//! would need `BigFloat: Num`, and every operation here also threads an
//! explicit precision and rounding mode.
//!
//! Only the operations that need no transcendental constants live here;
//! `exp`, `ln`, the trigonometric family and `sqrt` stay in
//! `transforms::evalf`, next to the constant cache they draw on.

use astro_float::{BigFloat, RoundingMode};

/// A complex number at some working precision: `(real part, imaginary part)`.
///
/// Crate-internal, arbitrary-precision counterpart of
/// [`num_complex::Complex64`].
pub(crate) type Complex = (BigFloat, BigFloat);

/// `0 + 0i` at `prec` bits.
pub(crate) fn c_zero(prec: usize) -> Complex {
    (BigFloat::new(prec), BigFloat::new(prec))
}

/// `1 + 0i` at `prec` bits.
pub(crate) fn c_one(prec: usize) -> Complex {
    (BigFloat::from_i32(1, prec), BigFloat::new(prec))
}

/// `0 + 1i` at `prec` bits.
pub(crate) fn c_i(prec: usize) -> Complex {
    (BigFloat::new(prec), BigFloat::from_i32(1, prec))
}

/// `r + 0i`.
pub(crate) fn c_from_real(r: BigFloat, prec: usize) -> Complex {
    (r, BigFloat::new(prec))
}

pub(crate) fn c_add(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    (a.0.add(&b.0, prec, rm), a.1.add(&b.1, prec, rm))
}

pub(crate) fn c_sub(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    (a.0.sub(&b.0, prec, rm), a.1.sub(&b.1, prec, rm))
}

pub(crate) fn c_mul(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    // (a+bi)(c+di) = (ac-bd) + (ad+bc)i
    let ac = a.0.mul(&b.0, prec, rm);
    let bd = a.1.mul(&b.1, prec, rm);
    let ad = a.0.mul(&b.1, prec, rm);
    let bc = a.1.mul(&b.0, prec, rm);
    (ac.sub(&bd, prec, rm), ad.add(&bc, prec, rm))
}

pub(crate) fn c_div(a: &Complex, b: &Complex, prec: usize, rm: RoundingMode) -> Complex {
    // (a+bi)/(c+di) = ((ac+bd) + (bc-ad)i) / (c²+d²)
    let ac = a.0.mul(&b.0, prec, rm);
    let bd = a.1.mul(&b.1, prec, rm);
    let bc = a.1.mul(&b.0, prec, rm);
    let ad = a.0.mul(&b.1, prec, rm);
    let denom =
        b.0.mul(&b.0, prec, rm)
            .add(&b.1.mul(&b.1, prec, rm), prec, rm);
    let re = ac.add(&bd, prec, rm).div(&denom, prec, rm);
    let im = bc.sub(&ad, prec, rm).div(&denom, prec, rm);
    (re, im)
}

pub(crate) fn c_neg(a: &Complex) -> Complex {
    (a.0.neg(), a.1.neg())
}

/// `|z| = √(re² + im²)`.
pub(crate) fn c_abs(a: &Complex, prec: usize, rm: RoundingMode) -> BigFloat {
    let re2 = a.0.mul(&a.0, prec, rm);
    let im2 = a.1.mul(&a.1, prec, rm);
    let sum = re2.add(&im2, prec, rm);
    sum.sqrt(prec, rm)
}

/// `zⁿ` by binary powering.
pub(crate) fn c_powi(base: &Complex, n: usize, prec: usize, rm: RoundingMode) -> Complex {
    if n == 0 {
        return c_one(prec);
    }
    let mut result = c_one(prec);
    let mut b = base.clone();
    let mut exp = n;
    while exp > 0 {
        if exp & 1 == 1 {
            result = c_mul(&result, &b, prec, rm);
        }
        b = c_mul(&b, &b, prec, rm);
        exp >>= 1;
    }
    result
}