Skip to main content

simsym/
lib.rs

1//! simsym — a simple symbolic computation library.
2//!
3//! ## Cargo features
4//!
5//! - **`simplify`** (default): algebraic simplification via [`Expr::simplify`]
6//! - **`diff`** (default): symbolic differentiation, [`Expr::gradient`], [`Expr::hessian`]
7//! - **`integrate`** (default, implies `diff`): symbolic integration and [`Expr::integrate_definite`]
8//!
9//! Build a minimal core (expression AST, rationals, evaluation only):
10//!
11//! ```bash
12//! cargo build --no-default-features
13//! ```
14//!
15//! ```rust
16//! use simsym::prelude::*;
17//! let x = symbol("x");
18//! let f = x.pow(2) + rational(2, 1) * x;
19//! let val = f.eval(&[(x, rational(1, 2))]).unwrap();
20//! assert_eq!(val, rational(5, 4));
21//! ```
22
23pub mod constant;
24pub mod display;
25pub mod eval;
26pub mod expr;
27mod num_convert;
28#[cfg(feature = "bigint")]
29mod num_convert_big;
30mod ops_ext;
31pub mod poly;
32pub mod rational;
33pub mod symbol;
34
35#[cfg(feature = "simplify")]
36pub mod simplify;
37#[cfg(not(feature = "simplify"))]
38pub mod simplify_nop;
39#[cfg(not(feature = "simplify"))]
40pub use simplify_nop as simplify;
41
42#[cfg(any(feature = "diff", feature = "integrate"))]
43pub mod calculus;
44
45#[cfg(feature = "bigint")]
46pub mod rational_big;
47
48pub use eval::EvalError;
49pub use expr::Expr;
50pub use constant::Constant;
51pub use num_convert::RationalConvertError;
52pub use rational::{rational, rational_from_i32, Rational};
53#[cfg(feature = "bigint")]
54pub use rational_big::{big_rational, BigInt, BigRational};
55#[cfg(feature = "bigint")]
56pub use constant::big_const;
57pub use symbol::{symbol, Symbol};
58
59#[cfg(any(feature = "diff", feature = "integrate"))]
60pub use calculus::{integrate_numeric, DefiniteIntegralError, NumericOptions};
61
62#[cfg(feature = "integrate")]
63pub use calculus::{integrate, integrate_definite, IntegrateError};
64
65#[cfg(feature = "diff")]
66pub use calculus::{diff, diff_without_simplify, gradient, hessian};
67
68pub use simsym_macros::expr;
69
70macro_rules! unary_fn {
71    ($($name:ident),* $(,)?) => {
72        $(pub fn $name(e: impl Into<Expr>) -> Expr {
73            expr::$name(e.into())
74        })*
75    };
76}
77
78unary_fn!(
79    sin, cos, tan, cot, sec, csc, asin, acos, atan, acot, asec, acsc, sinh, cosh, tanh, coth,
80    sech, csch, asinh, acosh, atanh, acoth, asech, acsch, exp, ln,
81);
82
83pub mod prelude {
84    pub use crate::{
85        acos, acosh, acot, acoth, acsc, acsch, asec, asech, asin, asinh, atan, atanh, cos, cosh,
86        cot, coth, csc, csch, exp, expr, ln, rational, rational_from_i32, sec, sech, sin, sinh,
87        symbol, tan, tanh, EvalError, Expr, Rational, Symbol,
88    };
89    #[cfg(any(feature = "diff", feature = "integrate"))]
90    pub use crate::{DefiniteIntegralError, NumericOptions, integrate_numeric};
91    #[cfg(feature = "integrate")]
92    pub use crate::{integrate, integrate_definite, IntegrateError};
93    #[cfg(feature = "diff")]
94    pub use crate::{diff_without_simplify, gradient, hessian};
95}
96
97#[cfg(test)]
98mod tests {
99    use super::prelude::*;
100
101    #[test]
102    fn eval_exact() {
103        let x = symbol("x");
104        let f = x.pow(2) + rational(2, 1) * x;
105        let v = f.eval(&[(x, rational(1, 2))]).unwrap();
106        assert_eq!(v, rational(5, 4));
107    }
108
109    #[test]
110    #[cfg(feature = "bigint")]
111    fn bigint_i128_constant() {
112        use crate::rational_big::BigRational;
113        let e: Expr = i128::MAX.into();
114        assert_eq!(e.to_string(), format!("{}", i128::MAX));
115        let br = BigRational::from(i128::MAX);
116        assert_eq!(br.numer().to_string(), i128::MAX.to_string());
117        assert!(Rational::try_from_big(&br).is_none());
118        assert!((e.eval_f64(&[]).unwrap() - i128::MAX as f64).abs() < 1.0);
119    }
120
121    #[test]
122    fn integer_and_float_conversions() {
123        assert_eq!(Rational::try_from(42u8).unwrap(), rational(42, 1));
124        assert_eq!(Rational::try_from(-3i8).unwrap(), rational(-3, 1));
125        assert!(Rational::try_from(i128::MAX).is_err());
126        let r = Rational::try_from(0.25f32).unwrap();
127        assert!((r.to_f32() - 0.25f32).abs() < 1e-6);
128        let x = symbol("x");
129        let f = 2u16 * x.to_expr() + Expr::from(1i8);
130        assert_eq!(f.eval(&[(x, rational(3, 1))]).unwrap(), rational(7, 1));
131    }
132
133    #[test]
134    fn eval_f32_works() {
135        let x = symbol("x");
136        let f = sin(x) + cos(x);
137        let v64 = f.clone().eval_f64(&[(x, 1.0)]).unwrap();
138        let v32 = f.eval_f32(&[(x, 1.0f32)]).unwrap();
139        assert!((v64 - f64::from(v32)).abs() < 1e-5);
140    }
141
142    #[test]
143    fn extended_trig_eval_f64() {
144        let x = symbol("x");
145        let env = [(x, 0.5)];
146        assert!((cot(x).eval_f64(&env).unwrap() - 1.0 / 0.5f64.tan()).abs() < 1e-10);
147        assert!((sinh(x).eval_f64(&env).unwrap() - 0.5f64.sinh()).abs() < 1e-10);
148        assert_eq!(cot(x).to_string(), "cot(x)");
149        assert_eq!(sinh(x).to_string(), "sinh(x)");
150    }
151
152    #[test]
153    #[cfg(all(feature = "diff", feature = "simplify"))]
154    fn diff_power() {
155        let x = symbol("x");
156        let f = x.pow(2);
157        let df = f.diff(x).simplify();
158        let expected = rational(2, 1) * x.to_expr();
159        assert_eq!(df.simplify().to_string(), expected.simplify().to_string());
160    }
161
162    #[test]
163    #[cfg(all(feature = "integrate", feature = "diff"))]
164    fn diff_without_simplify_matches_numeric() {
165        let x = symbol("x");
166        let f = sin(x).pow(3) * cos(x).pow(2);
167        let fx = f.clone().integrate(x).unwrap();
168        let df = fx.diff_without_simplify(x);
169        let v = f.eval_f64(&[(x, 0.5)]).unwrap();
170        let rv = df.eval_f64(&[(x, 0.5)]).unwrap();
171        assert!((v - rv).abs() < 1e-5);
172    }
173
174    #[test]
175    #[cfg(all(feature = "integrate", feature = "simplify"))]
176    fn integrate_polynomial() {
177        let x = symbol("x");
178        let f = x.pow(2);
179        let antideriv = f.integrate(x).unwrap().simplify();
180        let expected = x.pow(3) / rational(3, 1);
181        assert_eq!(antideriv.to_string(), expected.simplify().to_string());
182    }
183
184    #[test]
185    #[cfg(feature = "integrate")]
186    fn integrate_exp_manual() {
187        let x = symbol("x");
188        let f = exp(x) * (x.pow(3) + rational(2, 1) * x);
189        assert!(f.integrate(x).is_ok());
190    }
191
192    #[test]
193    #[cfg(all(feature = "diff", feature = "simplify"))]
194    fn diff_exp_times_polynomial() {
195        let x = symbol("x");
196        let f = exp(x) * (x.pow(3) + rational(2, 1) * x);
197        let df = f.clone().diff(x).simplify();
198        let x0 = rational(1, 2);
199        let f_val = f.clone().eval_f64(&[(x, x0.to_f64())]).unwrap();
200        let df_val = df.eval_f64(&[(x, x0.to_f64())]).unwrap();
201        let h = 1e-6;
202        let numeric = (f.clone().eval_f64(&[(x, x0.to_f64() + h)]).unwrap()
203            - f.eval_f64(&[(x, x0.to_f64() - h)]).unwrap())
204            / (2.0 * h);
205        assert!((df_val - numeric).abs() < 1e-4, "f'={df_val} numeric={numeric}");
206        assert!((df_val - f_val).abs() > 0.0);
207    }
208
209    #[test]
210    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
211    fn integrate_exp_times_polynomial() {
212        let x = symbol("x");
213        let f = exp(x) * (x.pow(3) + rational(2, 1) * x);
214        let antideriv = f.clone().integrate(x).unwrap().simplify();
215        let recovered = antideriv.diff(x).simplify();
216        let x0 = rational(1, 3);
217        let env = [(x, x0.to_f64())];
218        let f_val = f.eval_f64(&env).unwrap();
219        let recovered_val = recovered.eval_f64(&env).unwrap();
220        assert!(
221            (f_val - recovered_val).abs() < 1e-6,
222            "d/dx(F) should match f: {recovered_val} vs {f_val}"
223        );
224    }
225
226    #[test]
227    #[cfg(feature = "simplify")]
228    fn polynomial_normal_form_orders_terms() {
229        let x = symbol("x");
230        let messy = x.pow(3) - rational(3, 1) * x.pow(2) + rational(6, 1) * x + (-rational(6, 1))
231            + rational(2, 1) * x
232            + (-rational(2, 1));
233        let neat =
234            x.pow(3) - rational(3, 1) * x.pow(2) + rational(8, 1) * x + (-rational(8, 1));
235        assert_eq!(messy.simplify().to_string(), neat.simplify().to_string());
236    }
237
238    #[test]
239    #[cfg(all(feature = "integrate", feature = "simplify"))]
240    fn integrate_cancels_coefficients_in_quotient() {
241        let x = symbol("x");
242        let f = x.pow(3) + rational(2, 1) * x;
243        let antideriv = f.integrate(x).unwrap().simplify();
244        assert_eq!(
245            antideriv.to_string(),
246            (x.pow(4) / rational(4, 1) + x.pow(2)).simplify().to_string()
247        );
248    }
249
250    #[test]
251    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
252    fn sin_diff_and_integrate() {
253        let x = symbol("x");
254        let f = sin(x);
255        let df = f.clone().diff(x).simplify();
256        assert_eq!(df.to_string(), cos(x).simplify().to_string());
257        let antideriv = f.integrate(x).unwrap().simplify();
258        assert_eq!(antideriv.to_string(), (-cos(x)).simplify().to_string());
259    }
260
261    #[test]
262    #[cfg(all(feature = "diff", feature = "simplify"))]
263    fn gradient_two_vars() {
264        let x = symbol("x");
265        let y = symbol("y");
266        let f = x * y;
267        let g = f.gradient(&[x, y]);
268        assert_eq!(g[0].clone().simplify().to_string(), y.to_string());
269        assert_eq!(g[1].clone().simplify().to_string(), x.to_string());
270    }
271
272    #[test]
273    #[cfg(all(feature = "diff", feature = "simplify"))]
274    fn hessian_bilinear() {
275        let x = symbol("x");
276        let y = symbol("y");
277        let f = x * y;
278        let h = f.hessian(&[x, y]);
279        assert_eq!(h[0][0].clone().simplify().to_string(), "0");
280        assert_eq!(h[0][1].clone().simplify().to_string(), "1");
281        assert_eq!(h[1][0].clone().simplify().to_string(), "1");
282        assert_eq!(h[1][1].clone().simplify().to_string(), "0");
283    }
284
285    #[test]
286    #[cfg(all(feature = "integrate", feature = "simplify"))]
287    fn integrate_tan_and_ln() {
288        let x = symbol("x");
289        let tan_int = crate::expr::tan(x.to_expr()).integrate(x).unwrap().simplify();
290        assert_eq!(
291            tan_int.to_string(),
292            (-crate::expr::ln(crate::expr::cos(x.to_expr()))).simplify().to_string()
293        );
294        let ln_int = crate::expr::ln(x.to_expr()).integrate(x).unwrap().simplify();
295        let x_expr = x.to_expr();
296        let expected_ln = x_expr.clone() * crate::expr::ln(x_expr.clone()) - x_expr;
297        assert_eq!(ln_int.to_string(), expected_ln.simplify().to_string());
298    }
299
300    #[test]
301    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
302    fn integrate_sin_squared() {
303        let x = symbol("x");
304        let f = sin(x).pow(2);
305        let f_int = f.clone().integrate(x).unwrap().simplify();
306        let recovered = f_int.diff(x).simplify();
307        let v = f.eval_f64(&[(x, 0.3)]).unwrap();
308        let rv = recovered.eval_f64(&[(x, 0.3)]).unwrap();
309        assert!((v - rv).abs() < 1e-6);
310    }
311
312    #[test]
313    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
314    fn integrate_sin_cos_product() {
315        let x = symbol("x");
316        let f = sin(x) * cos(x);
317        let f_int = f.clone().integrate(x).unwrap().simplify();
318        let recovered = f_int.diff(x).simplify();
319        let v = f.eval_f64(&[(x, 0.4)]).unwrap();
320        let rv = recovered.eval_f64(&[(x, 0.4)]).unwrap();
321        assert!((v - rv).abs() < 1e-5);
322    }
323
324    #[test]
325    #[cfg(all(feature = "integrate", feature = "simplify"))]
326    fn integrate_partial_fractions() {
327        use crate::expr::const_;
328        let x = symbol("x");
329        let f = const_(Rational::one())
330            / (x.to_expr() - const_(rational(1, 1)))
331            / (x.to_expr() - const_(rational(2, 1)));
332        let antideriv = f.integrate(x).unwrap().simplify();
333        assert!(
334            antideriv.to_string().contains("ln"),
335            "expected ln terms, got {antideriv}"
336        );
337    }
338
339    #[test]
340    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
341    fn integrate_one_over_x_squared_plus_one() {
342        use crate::expr::const_;
343        let x = symbol("x");
344        let f = const_(Rational::one()) / (x.to_expr().pow(2) + rational(1, 1));
345        let antideriv = f.clone().integrate(x).unwrap().simplify();
346        assert!(
347            antideriv.to_string().contains("atan"),
348            "expected atan, got {antideriv}"
349        );
350        let recovered = antideriv.diff(x).simplify();
351        let v = f.eval_f64(&[(x, 0.5)]).unwrap();
352        let rv = recovered.eval_f64(&[(x, 0.5)]).unwrap();
353        assert!((v - rv).abs() < 1e-6, "{rv} vs {v}");
354    }
355
356    #[test]
357    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
358    fn integrate_rational_polynomial_quotient() {
359        use crate::expr::const_;
360        let x = symbol("x");
361        let f = (x.to_expr().pow(2) + rational(1, 1))
362            / (x.to_expr() - const_(rational(1, 1)));
363        let antideriv = f.clone().integrate(x).unwrap().simplify();
364        let recovered = antideriv.diff(x).simplify();
365        let v = f.eval_f64(&[(x, 2.0)]).unwrap();
366        let rv = recovered.eval_f64(&[(x, 2.0)]).unwrap();
367        assert!((v - rv).abs() < 1e-5, "{rv} vs {v}");
368    }
369
370    #[test]
371    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
372    fn integrate_sin_fourth_power() {
373        let x = symbol("x");
374        let f = sin(x).pow(4);
375        let antideriv = f.clone().integrate(x).unwrap().simplify();
376        let recovered = antideriv.diff(x).simplify();
377        let v = f.eval_f64(&[(x, 0.25)]).unwrap();
378        let rv = recovered.eval_f64(&[(x, 0.25)]).unwrap();
379        assert!((v - rv).abs() < 1e-5, "{rv} vs {v}");
380    }
381
382    #[test]
383    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
384    fn integrate_atan_x() {
385        let x = symbol("x");
386        let f = crate::expr::atan(x.to_expr());
387        let antideriv = f.clone().integrate(x).unwrap().simplify();
388        assert!(antideriv.to_string().contains("atan"));
389        let recovered = antideriv.diff(x).simplify();
390        let v = f.eval_f64(&[(x, 1.0)]).unwrap();
391        let rv = recovered.eval_f64(&[(x, 1.0)]).unwrap();
392        assert!((v - rv).abs() < 1e-5);
393    }
394
395    #[test]
396    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
397    fn integrate_exp_sin() {
398        let x = symbol("x");
399        let f = exp(x) * sin(x);
400        let antideriv = f.clone().integrate(x).unwrap().simplify();
401        let recovered = antideriv.diff(x).simplify();
402        let v = f.eval_f64(&[(x, 0.7)]).unwrap();
403        let rv = recovered.eval_f64(&[(x, 0.7)]).unwrap();
404        assert!((v - rv).abs() < 1e-5, "{rv} vs {v}");
405    }
406
407    #[test]
408    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
409    fn integrate_x_squared_ln() {
410        let x = symbol("x");
411        let f = x.to_expr().pow(2) * crate::expr::ln(x.to_expr());
412        let antideriv = f.clone().integrate(x).unwrap().simplify();
413        let recovered = antideriv.diff(x).simplify();
414        let v = f.eval_f64(&[(x, 2.0)]).unwrap();
415        let rv = recovered.eval_f64(&[(x, 2.0)]).unwrap();
416        assert!((v - rv).abs() < 1e-5);
417    }
418
419    #[test]
420    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
421    fn integrate_tan_cubed() {
422        let x = symbol("x");
423        let f = crate::expr::tan(x.to_expr()).pow(3);
424        let antideriv = f.clone().integrate(x).unwrap().simplify();
425        let recovered = antideriv.diff(x).simplify();
426        let v = f.eval_f64(&[(x, 0.3)]).unwrap();
427        let rv = recovered.eval_f64(&[(x, 0.3)]).unwrap();
428        assert!((v - rv).abs() < 1e-5);
429    }
430
431    #[test]
432    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
433    fn integrate_one_over_x_squared_minus_one() {
434        use crate::expr::const_;
435        let x = symbol("x");
436        let f = const_(Rational::one()) / (x.to_expr().pow(2) - const_(rational(1, 1)));
437        let antideriv = f.clone().integrate(x).unwrap().simplify();
438        assert!(antideriv.to_string().contains("ln"));
439        let recovered = antideriv.diff(x).simplify();
440        let v = f.eval_f64(&[(x, 2.0)]).unwrap();
441        let rv = recovered.eval_f64(&[(x, 2.0)]).unwrap();
442        assert!((v - rv).abs() < 1e-5);
443    }
444
445    #[test]
446    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
447    fn integrate_exp_sin_with_phase() {
448        let x = symbol("x");
449        let f = exp(rational(2, 1) * x.to_expr() + rational(1, 1))
450            * sin(rational(3, 1) * x.to_expr() + rational(1, 4));
451        let antideriv = f.clone().integrate(x).unwrap().simplify();
452        let recovered = antideriv.diff(x).simplify();
453        let v = f.eval_f64(&[(x, 0.2)]).unwrap();
454        let rv = recovered.eval_f64(&[(x, 0.2)]).unwrap();
455        assert!((v - rv).abs() < 1e-5, "{rv} vs {v}");
456    }
457
458    #[test]
459    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
460    fn integrate_sin_squared_cos_squared() {
461        let x = symbol("x");
462        let f = sin(x).pow(2) * cos(x).pow(2);
463        let antideriv = f.clone().integrate(x).unwrap().simplify();
464        let recovered = antideriv.diff(x).simplify();
465        let v = f.eval_f64(&[(x, 0.4)]).unwrap();
466        let rv = recovered.eval_f64(&[(x, 0.4)]).unwrap();
467        assert!((v - rv).abs() < 1e-5);
468    }
469
470    #[test]
471    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
472    fn integrate_sin_cubed_cos_squared() {
473        let x = symbol("x");
474        let f = sin(x).pow(3) * cos(x).pow(2);
475        let antideriv = f.clone().integrate(x).unwrap().simplify();
476        let recovered = antideriv.diff(x).simplify();
477        let v = f.eval_f64(&[(x, 0.35)]).unwrap();
478        let rv = recovered.eval_f64(&[(x, 0.35)]).unwrap();
479        assert!((v - rv).abs() < 1e-5);
480    }
481
482    #[test]
483    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
484    fn integrate_sec_fourth() {
485        let x = symbol("x");
486        let f = cos(x).pow(-4);
487        let antideriv = f.clone().integrate(x).unwrap().simplify();
488        let recovered = antideriv.diff(x).simplify();
489        let v = f.eval_f64(&[(x, 0.25)]).unwrap();
490        let rv = recovered.eval_f64(&[(x, 0.25)]).unwrap();
491        assert!((v - rv).abs() < 1e-5, "{rv} vs {v}");
492    }
493
494    #[test]
495    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
496    fn integrate_sin_cos_with_phase() {
497        let x = symbol("x");
498        let f = sin(x.to_expr() + rational(1, 4)) * cos(rational(2, 1) * x.to_expr());
499        let antideriv = f.clone().integrate(x).unwrap().simplify();
500        let recovered = antideriv.diff(x).simplify();
501        let v = f.eval_f64(&[(x, 0.5)]).unwrap();
502        let rv = recovered.eval_f64(&[(x, 0.5)]).unwrap();
503        assert!((v - rv).abs() < 1e-5);
504    }
505
506    #[test]
507    #[cfg(all(feature = "integrate", feature = "diff", feature = "simplify"))]
508    fn integrate_x_exp_by_parts() {
509        let x = symbol("x");
510        let f = x.to_expr() * exp(x);
511        let f_int = f.clone().integrate(x).unwrap().simplify();
512        let recovered = f_int.diff(x).simplify();
513        let v = f.eval_f64(&[(x, 1.0)]).unwrap();
514        let rv = recovered.eval_f64(&[(x, 1.0)]).unwrap();
515        assert!((v - rv).abs() < 1e-6);
516    }
517
518    #[test]
519    #[cfg(feature = "integrate")]
520    fn numeric_integral_sin() {
521        let x = symbol("x");
522        let f = sin(x);
523        let pi = std::f64::consts::PI;
524        let val = integrate_numeric(&f, x, 0.0, pi, &[], NumericOptions::default()).unwrap();
525        assert!((val - 2.0).abs() < 1e-6);
526    }
527}