1pub 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}