Complex Analysis and Special Functions
New in 0.2. symplex is a CAS over ℂ: a symbol with no assumptions may be complex, and the library refuses to pretend otherwise.
re, im, conjugate, arg
Ex::{re, im, conjugate, arg} are constructed with as much evaluation as the structure and assumptions allow. For a symbol declared Real the parts are immediate; for an unassumed symbol you get an unevaluated re(z) node (0.1 silently assumed every symbol real — a wrong answer for conjugate(z)).
use symplex::prelude::*;
fn main() {
let ctx = Context::new();
let x = ctx.symbol_with("x", &[Assumption::Real]);
let y = ctx.symbol_with("y", &[Assumption::Real]);
let z = ctx.symbol("z"); // may be complex
let i = ctx.i_unit();
let w = &x + &i * &y;
println!("{}", w.re()); // x
println!("{}", w.im()); // y
println!("{}", w.conjugate()); // x - y*I
println!("{}", w.abs_squared()); // x^2 + y^2
println!("{}", w.arg()); // atan2(y, x)
println!("{}", z.re()); // re(z) — not assumed real
println!("{}", z.conjugate()); // conjugate(z)
println!("{}", (&z.powi(2) + 1).conjugate()); // conjugate(z)^2 + 1 (distributes)
println!("{}", z.exp().conjugate()); // exp(conjugate(z)) (commutes with exp)
println!("{}", (&i * &z).re()); // -im(z)
println!("{}", z.exp().re()); // cos(im(z))*exp(re(z))
println!("{:?} {:?}", (&x.powi(2) + 1).is_real_valued(), z.is_real_valued()); // Some(true) None
}Splitting into real and imaginary parts
as_real_imag() returns (re, im) as a pair; expand_complex() rewrites the expression as re + im·I; polar() returns a Polar { modulus, argument } struct (symplex::expr_complex::Polar).
use symplex::prelude::*;
fn main() {
let ctx = Context::new();
let x = ctx.symbol_with("x", &[Assumption::Real]);
let y = ctx.symbol_with("y", &[Assumption::Real]);
let w = &x + &ctx.i_unit() * &y;
let (re, im) = w.powi(2).as_real_imag();
println!("({re}) + ({im})i"); // (x^2 - y^2) + (2*x*y)i
let (re, im) = w.exp().as_real_imag();
println!("({re}) + ({im})i"); // (cos(y)*exp(x)) + (sin(y)*exp(x))i
let (re, im) = w.sin().as_real_imag();
println!("({re}) + ({im})i"); // (sin(x)*cosh(y)) + (cos(x)*sinh(y))i
let (re, im) = (1 / &w).as_real_imag();
println!("({re}) + ({im})i"); // (x/(x^2 + y^2)) + (-y/(x^2 + y^2))i
println!("{}", (&ctx.i_unit() * &x).cos().expand_complex()); // cosh(x)
}Concrete complex numbers
Gaussian rationals are exact. Use expand() (or expand_complex()) to multiply out and as_real_imag() to divide; abs_squared() avoids the square root.
use symplex::prelude::*;
fn main() {
let ctx = Context::new();
let i = ctx.i_unit();
let a = ctx.complex(&ctx.int(3), &ctx.int(4)); // 3 + 4i
let b = ctx.int(1) - &i * 2; // 1 − 2i
println!("{}", (&a * &b).expand()); // -2*I + 11
let (re, im) = (&a / &b).as_real_imag();
println!("{re} + {im}i"); // -1 + 2i
println!("{}", a.abs_squared()); // 25
println!("{}", (ctx.int(1) + &i).arg().eval()); // 1/4*pi
println!("{}", (ctx.int(1) + &i).powi(8).expand()); // 16
println!("{}", (&i * &ctx.pi()).exp().eval()); // -1
println!("{}", ctx.int(-1).ln().eval()); // pi*I
println!("{}", ctx.int(-4).sqrt()); // 2*I
println!("{}", i.exp().eval_complex64().unwrap()); // 0.5403…+0.8414…i (a `Complex64`)
}Current gaps, stated plainly: abs(3 + 4i) is not folded to 5 by eval/simplify (use abs_squared()), and ln(i) / i^i stay symbolic (evaluate with eval_decimal).
Complex infinity
1/0 evaluates to complex infinity zoo (Context::complex_infinity()), distinct from the signed real infinities oo and -oo. The parser accepts zoo.
New constants
| Constant | Context method | Notes |
|---|---|---|
| γ (Euler–Mascheroni) | euler_gamma() | rationality unknown → assumption system leaves it open |
| G (Catalan) | catalan() | positive, real, finite |
| φ (golden ratio) | golden_ratio() | algebraic; nsimplify recognises φ² − φ − 1 = 0, simplify does not yet |
| ∞̃ | complex_infinity() | zoo |
All evaluate to arbitrary precision with eval_decimal(digits).
Special functions
0.2 adds si, ci, ei, li, zeta, polygamma(n, x) and kronecker_delta(i, j), with exact special values, derivative rules and evalf; Digamma folds at positive integers and half-integers; Bessel I/K and orthogonal polynomials of any degree evaluate numerically; all Bessel functions and orthogonal polynomials have derivative rules.
use symplex::prelude::*;
fn main() {
let ctx = Context::new();
let x = ctx.symbol_with("x", &[Assumption::Real]);
println!("{}", ctx.int(1).digamma().eval()); // -EulerGamma
println!("{}", ctx.rational(1, 2).digamma().eval()); // -EulerGamma - 2*ln(2)
println!("{}", ctx.int(1).polygamma(&ctx.int(1)).eval()); // 1/6*pi^2
println!("{}", ctx.int(1).polygamma(&ctx.int(2)).eval()); // -2*zeta(3)
println!("{}", x.digamma().diff(&x)); // polygamma(1, x)
println!("{} {} {} {}", ctx.int(2).zeta().eval(), ctx.int(4).zeta().eval(),
ctx.int(0).zeta().eval(), ctx.int(-1).zeta().eval()); // 1/6*pi^2 1/90*pi^4 -1/2 -1/12
println!("{}", ctx.int(3).zeta().eval_decimal(30).unwrap()); // 1.20205690315959428539973816151
println!("{} {}", ctx.int(0).si().eval(), ctx.infinity().si().eval()); // 0 1/2*pi
println!("{} {}", x.si().diff(&x), x.li().diff(&x)); // sin(x)/x 1/ln(x)
println!("{}", ctx.int(1).ei().eval_decimal(25).unwrap()); // 1.89511781635593675546652
println!("{} {}", x.kronecker_delta(&x), ctx.int(1).kronecker_delta(&ctx.int(2))); // 1 0
println!("{}", ctx.int(2).bessel_k(&ctx.int(0)).eval_decimal(20).unwrap()); // 0.11389387274953343565
println!("{}", ctx.rational(1, 3).legendre(&ctx.int(10)).eval()); // 13597/59049
}The parser accepts all of these by name (re, im, conjugate/conj, arg, si, ci, ei, li, zeta, polygamma, zoo), and compile() / to_rust_fn / to_c_fn support Γ, lnΓ, ψ, erf/erfc, W, B, Bessel and orthogonal polynomials.
See cargo run --example complex_analysis for the full tour.
More special functions (0.9)
0.9 adds the remaining SymPy special functions that show up as integration
results and in physics. All of them are Apply nodes with SymPy’s names, so
they print, parse (ctx.parse("erfi(x)")) and serialise like besselj; the
method lives on the argument and takes the parameters first, as with
x.bessel_j(&nu).
| Function | Constructor | Exact values (eval) | Derivative |
|---|---|---|---|
| erfi, erf⁻¹, erfc⁻¹ | x.erfi(), x.erfinv(), x.erfcinv() | erfi(0) = 0, odd; erfinv(±1) = ±∞; erfcinv(1) = 0 | 2e^{x²}/√π; (√π/2) e^{erfinv²} |
| Eₙ, E₁ | x.expint(&n), x.e1() | Eₙ(0) = 1/(n−1), E₀(x) = e^{−x}/x, Eₙ(∞) = 0 | −Eₙ₋₁(x) |
| Shi, Chi | x.shi(), x.chi() | Shi(0) = 0, odd; Chi(0) = −∞ | sinh x/x, cosh x/x |
| Fresnel S, C | x.fresnels(), x.fresnelc() | S(0) = 0, S(±∞) = ±1/2, odd | sin(πx²/2), cos(πx²/2) |
| γ(s,x), Γ(s,x) | x.lowergamma(&s), x.uppergamma(&s) | closed forms for integer and half-integer s (Γ(1,x) = e^{−x}, Γ(0,x) = E₁(x), Γ(½,x) = √π erfc(√x), …) | ±x^{s−1}e^{−x} |
| Liₛ(z) | z.polylog(&s) | Liₛ(0) = 0, Liₛ(1) = ζ(s), Liₛ(−1) = −η(s), Li₁ = −ln(1−z), Li₀ = z/(1−z), Li₋ₙ rational, Li₂(½) | Liₛ₋₁(z)/z |
| η(s) | s.dirichlet_eta() | η(1) = ln 2, η(s) = (1−2^{1−s})ζ(s) when ζ(s) folds | formal |
| Ai, Bi, Ai′, Bi′ | x.airyai(), x.airybi(), x.airyaiprime(), x.airybiprime() | values at 0 in terms of Γ(⅓), Γ(⅔); limits at ±∞ | Ai′ = airyaiprime, Ai″ = x·Ai |
| K(m), E(m) | m.elliptic_k(), m.elliptic_e() | K(0) = E(0) = π/2, E(1) = 1, K(1) = z∞, K(½) = Γ(¼)²/(4√π) | (E − (1−m)K)/(2m(1−m)), (E − K)/(2m) |
| F(φ|m), Π(n|m) | phi.elliptic_f(&m), n.elliptic_pi(&m) | F(0|m) = 0, F(φ|0) = φ, F(π/2|m) = K(m); Π(0|m) = K(m), Π(n|0) = π/(2√(1−n)), Π(n|n) = E(n)/(1−n) | ∂φF = 1/√(1−m sin²φ) (∂ₘF formal); both partials of Π |
| Cₙ^(a), Pₙ^(a,b), Pₙ^m, Lₙ^(a) | x.gegenbauer(&n, &a), x.jacobi(&n, &a, &b), x.assoc_legendre(&n, &m), x.assoc_laguerre(&n, &a) | explicit polynomials for integer n (symbolic a, b allowed); Cₙ^(½) = Pₙ, Cₙ^(1) = Uₙ, Pₙ^(0,0) = Pₙ^0 = Pₙ, Lₙ^(0) = Lₙ | 2a Cₙ₋₁^(a+1), (n+a+b+1)/2 Pₙ₋₁^(a+1,b+1), (nxPₙ^m − (n+m)Pₙ₋₁^m)/(x²−1), −Lₙ₋₁^(a+1) |
assoc_legendre uses the Condon–Shortley phase like SymPy (P₁¹(x) = −√(1−x²));
the elliptic integrals take the parameter m = k².
Everything evaluates numerically to arbitrary precision (eval_decimal(40)
agrees with mpmath) for real arguments in the real domain: the error-function
family by series / asymptotics / Halley iteration, Eₙ, γ, Γ by series
and Legendre’s continued fraction, Liₛ by the direct series or the
expansion in ln z (any real z ≠ 1 for s ≤ 0, |z| ≤ 1 otherwise),
Airy by Maclaurin series with guard bits or the large-|x| asymptotic
expansions, and the elliptic integrals by Carlson’s symmetric forms
(n < 1, m < 1). Arguments outside these domains (complex Liₛ, K(m > 1),
erfinv(|y| ≥ 1), …) return Err(Unevaluable).
use symplex::prelude::*;
fn main() {
let ctx = Context::new();
let (x, z, a) = (ctx.symbol("x"), ctx.symbol("z"), ctx.symbol("a"));
println!("{}", x.erfi().diff(&x)); // 2*exp(x^2)/sqrt(pi)
println!("{}", ctx.rational(7, 10).erfi().eval_decimal(20).unwrap()); // 0.94028293383350747659
println!("{}", x.uppergamma(&ctx.int(3)).eval()); // x^2*exp(-x) + 2*x*exp(-x) + 2*exp(-x)
println!("{}", ctx.int(1).polylog(&ctx.int(2)).eval()); // 1/6*pi^2
println!("{}", z.polylog(&ctx.int(-1)).eval()); // z*(-z + 1)^(-2)
println!("{}", ctx.int(0).elliptic_k().eval()); // 1/2*pi
println!("{}", ctx.rational(1, 2).elliptic_k().eval_decimal(20).unwrap()); // 1.8540746773013719184
println!("{}", x.gegenbauer(&ctx.int(2), &a).eval()); // 2*a^2*x^2 + 2*a*x^2 - a
println!("{}", x.airyaiprime().diff(&x)); // x*airyai(x)
println!("{}", x.polylog(&ctx.int(2)).to_latex()); // \operatorname{Li}_{2}\left(x\right)
}to_lean returns NotImplemented for all of these (Mathlib has no standard
spelling), and compile / to_rust_fn / to_c_fn report the missing runtime
rather than generating code.