Symbolica 2.0: programmable symbols for Python and Rust

# Symbolica 2.0: programmable symbols – Symbolica | Modern Computer Algebra Source: [https://symbolica.io/posts/symbolica_2_0_release/](https://symbolica.io/posts/symbolica_2_0_release/) ## Introduction Symbolica is a high\-performance symbolic computation framework for Python and Rust\. You can use it to manipulate symbolic expressions and turn them into fast numerical kernels for computing Jacobians, numerical optimization, integration, and much more\. Today marks the 2\.0 release of Symbolica, with many exciting new features and improvements\. The theme of this release is programmable symbols: more of Symbolica’s behavior can now be customized by the user\. This makes it possible to define mathematical objects that simplify, differentiate, expand, print, and evaluate like built\-ins\. Since 1\.0, Symbolica has accumulated improvements in several directions: - a simpler Rust API with more operator overloading and builder\-style APIs - a symbol registration system, with namespaces, aliases, tags, user data, and custom hooks - a redesigned evaluator interface that supports double\-float arithmetic and JIT compilation - richer output for notebooks and documents, including HTML output, graph and polynomial display, Typst output, colorized printing, and more structural multiline formatting - new built\-in mathematical functions, including gamma, polylogarithms, Bessel functions, Riemann zeta, and related series/evaluation hooks See the[migration guide](https://symbolica.io/docs/migration.html)to learn more about the changes and how to update your code\. ## Better output Symbolica gained an automatic line\-wrapping output mode, with colored brackets, similar to how code is styled\. This makes it easier to read large, nested expressions\. In notebooks, such as Jupyter or[Marimo](https://marimo.io/), the default output is a colorful HTML\-mode and you can easily switch to LaTeX mode\. Typst output is available now as well\.  ## Improved Rust API One of the most visible changes is that ordinary Rust programs need far fewer imports and fewer long type paths\. The new prelude collects the common traits, macros, domains, and evaluator types that most users need\. Rust ergonomics have also been improved with additional overloads, automatic type conversions, builder patterns, and a`call`method on symbols: - Symbolica 2\.0 - Symbolica 1\.0 ``` use symbolica::prelude::*; fn main() { let (x, y, f) = symbol!("x", "y", "f"); let e: Atom = 2 + (x + 1).pow(-2) + f.call((1 + y, y)) / 3; let s = e.series(x, 0, 1).unwrap(); println!("{e}"); // -> 2+1/3*f(1+y,y)+1/(1+x)^2 println!("{s}"); // -> 3+1/3*f(1+y,y)-2*x+𝒪(x^2) } ``` ``` use symbolica::{ atom::{Atom, AtomCore}, function, symbol, }; fn main() { let (x, y, f) = symbol!("x", "y", "f"); let e = Atom::num(2) + (Atom::var(x) + 1).npow(-2) + function!(f, Atom::num(1) + y, y) / 3; let s = e.series(x, Atom::num(0), 1.into(), true).unwrap(); println!("{e}"); // -> 2+1/3*f(1+y,y)+1/(1+x)^2 println!("{s}"); // -> 3+1/3*f(1+y,y)-2*x+𝒪(x^2) } ``` Structures used for settings and functions with many arguments now use a builder pattern\. For example, the construction of a high\-performance numerical evaluator now looks like this: - Symbolica 2\.0 - Symbolica 1\.0 ``` use symbolica::prelude::*; fn main() -> Result<(), EvaluationError> { let mut evaluator = parse!("x^2 + 2*x + 1 + f(x)") .evaluator(&[parse!("x")]) .add_function(symbol!("f"), vec![symbol!("y")], parse!("cos(y + 1)"))? .horner_iterations(2) .build()? .map_coeff(&|c| c.re.to_f64()); println!("{}", evaluator.evaluate_single(&[3.0])); Ok(()) } ``` ``` use symbolica::{atom::AtomCore, evaluate::{FunctionMap, OptimizationSettings}, parse, symbol}; fn main() { let mut fn_map = FunctionMap::new(); fn_map .add_function( symbol!("f"), "f".to_string(), vec![symbol!("y")], parse!("cos(y + 1)"), ) .unwrap(); let params = vec![parse!("x")]; let optimization_settings = OptimizationSettings { horner_iterations: 2, ..OptimizationSettings::default() }; let mut evaluator = parse!("x^2 + 2*x + 1 + f(x)") .evaluator(&fn_map, &params, optimization_settings) .unwrap() .map_coeff(&|c| c.re.to_f64()); println!("{}", evaluator.evaluate_single(&[3.0])); } ``` With the builder pattern, settings arguments that used to belong to different substructures \(`add\_function`belongs to`FunctionMap`,`horner\_iterations`belongs to`OptimizationSettings`\) can now be set together in a more fluent way\. Also note that fallible operations now return a dedicated error type\. ## Programmable symbols In Symbolica 1\.0, a symbol could already carry important algebraic attributes: - Python - Rust ``` from symbolica import * x, y = S("x", "y") dot = S("dot", is_symmetric=True, is_linear=True) print(dot(3*x + y, x + 2*y)) # -> 3*dot(x,x)+7*dot(x,y)+2*dot(y,y) ``` ``` use symbolica::prelude::*; fn main() { let dot = symbol!("dot"; Symmetric, Linear); println!("{}", parse!("dot(3*x+y, x+2*y)")); // -> 3*dot(x,x)+7*dot(x,y)+2*dot(y,y) } ``` In 2\.0, symbols can now install more hooks that run at specific points in the algebraic lifecycle: HookUse caseNormalizationRewrite a symbol or function when it is normalized\.PrintingOverride Symbolica, LaTeX, or custom\-mode output\.DerivativeDefine a custom derivative rule\.SeriesTeach Symbolica how a function behaves near a singular point\.EvaluationRegister numerical implementations for evaluators\.For example, let us define a[gamma function](https://en.wikipedia.org/wiki/Gamma_function)\. Since gamma has a pole at 0, we cannot Taylor expand around 0\. Using the series hook we can tell Symbolica how to regularize the function near that point by applying the identity\\\(\\Gamma\(a \+ 1\) = a \\Gamma\(a\)\\\): - Python - Rust ``` from typing import Sequence from symbolica import * x = S("x") def regularize_gamma(args: Sequence[Series]) -> tuple[Expression, Expression] | None: if args[0].get_coefficient(0) == 0: a = args[0].to_expression() return (1 / a, S("gamma")(a + 1)) gamma = S( "gamma", derivative=lambda t, _: t * E("digamma")(t[0]), series=regularize_gamma, ) print(gamma(x).series(x, 0, 0)) # gamma(1)*x^-1+gamma(1)*digamma(1)+𝒪(x^1) ``` ``` use symbolica::prelude::*; fn main() { let _gamma = symbol!( "gamma", der = |f, index, out| { if index == 0 { let arg = f.as_fun_view().unwrap().get(0); **out = symbol!("gamma").call(arg) * symbol!("digamma").call(arg); } }, series = |args| { if args[0].coefficient(0.into()) == Atom::Zero { let a = args[0].to_atom(); Some((1 / &a, symbol!("gamma").call(a + 1))) } else { None } } ); let x = symbol!("x"); println!("{}", parse!("gamma(x)").series(x, 0, 0).unwrap()); } ``` For this example, the built\-in gamma function simplifies further to`1/x\-γ`\. ## Evaluators: from expression trees to compiled kernels Expression evaluation has seen the largest amount of engineering since 1\.0\. The high\-level flow is still the same: Symbolica rewrites an expression into a small instruction program, optimizes it, then evaluates it many times for different numeric inputs\. See[this blog post](https://symbolica.io/posts/expression_evaluation/index.html)for the details of how this works\. ### Evaluators for symbols Let us define a custom symbol`cosh`and teach Symbolica how to evaluate it for different numeric domains\. We can do this by registering an evaluation hook: - Python - Rust ``` import cmath import math from symbolica import * x = S("x") cosh = S( "cosh", eval={ "float": lambda args: math.cosh(args[0]), "complex": lambda args: cmath.cosh(args[0]), "cpp": "template<typename T> T python_cosh2(T a) { return std::cosh(a); }", }, ) expr = cosh(1/2) + x expr.to_float() # -> 1.1276259652063807 + x ``` ``` use symbolica::prelude::*; fn main() { let cosh = symbol!( "cosh", eval = EvaluationInfo::new() .register(|args: &[f64]| args[0].cosh()) .register(|args: &[Complex<f64>]| args[0].cosh()) ); let x = symbol!("x"); let expr = cosh.call(Atom::num(1) / 2) + x; println!("{}", expr.to_float(16)); // -> 1.127625965206381+x } ``` Since it is known how to evaluate`cosh`, Symbolica looks up a suitable evaluator in the call to`to\_float`\. ### JIT compilation Besides generating custom ASM, C\+\+, and CUDA code, Symbolica can now compile evaluators just in time\. This uses the[symjit](https://github.com/siravan/symjit)crate by Shahriar Iravanian, with whom we worked closely to make the integration feel native to Symbolica\. The JIT path supports custom evaluator hooks and is now the default evaluation backend for Python\. In practice, it is competitive with the custom ASM backend while keeping compilation times under control\. ### Double\-float arithmetic Symbolica 2\.0 adds a double\-float evaluation path\. Double\-float arithmetic stores a number as the unevaluated sum of two`f64`s, giving roughly 106 bits of precision \(31 decimal digits compared to 16 for ordinary doubles\) while staying more than 3x faster than arbitrary\-precision`Float`arithmetic\. In Python,`evaluate\_with\_prec\(\.\.\., 32\)`takes this path automatically\. ## Special functions Another major development made possible by the new hook system is that Symbolica has grown a larger mathematical vocabulary\. Symbolica now supports polygamma functions, polylogarithms, Bessel functions, Riemann zeta, geometric functions and their inverses, with evaluation hooks and Laurent/Puiseux series behavior around poles\. Some special values normalize immediately: ``` from symbolica import * print(E("gamma(1/2)")) # 𝜋^(1/2) print(E("polylog(-3, z)")) # z*(1+4*z+z^2)/(1-z)^4 ``` All special functions can be evaluated: ``` g = E("gamma(5/6)") print(g) # gamma(5/6) print(g.to_float()) # 1.128787029908126 ``` Symbolica can use special functions inside optimized evaluators: ``` from symbolica import * import numpy as np x = S('x') e = E('pi + zeta(3) + gamma(x)') ev = e.evaluator([x]) print(e) # gamma(x)+zeta(3)+𝜋 print(ev.evaluate(np.array([[2.0]]))) # [[5.34364956]] ``` When exporting code, constants such as\\\(\\pi\\\)and\\\(\\zeta\(3\)\\\)are replaced by their numerical value, so the generated code does not need to depend on a special function library\. ## Faster manipulation Not all important changes are visible in the public API\. There have also been substantial improvements to the machinery underneath expression manipulation: - Pattern matching can prove structural mismatches and terminate early - Term sorting can now sort based on a byte slice comparison \(30% faster normalization\) - Horner schemes can be applied on linear expression structures \(\>2x memory reduction\) - Common pair elimination uses much less memory - Improvements to polynomial GCD computations, including a new modular algorithm These combined changes have increased performance by factors between 2 and 10,000x on real\-world use cases\. For these improvements, user feedback was invaluable\. If you have a slow computation, please share it with us so we can make Symbolica faster\. ## Technical deep dive: type\-erased evaluation callbacks This section is for the technically curious\. It explains how Symbolica’s evaluation hooks can support multiple numerical domains without sacrificing ergonomics or performance\. The evaluation\-hook API has an interesting internal problem: a single symbol may have several numerical implementations, each with a different Rust type that implements a specific trait with a member function of the form`fn\(&\[T\]\) \-\> T`for some numeric domain`T`\. For example, a function might have distinct fast implementations for: ``` |args: &[f64]| -> f64 |args: &[Complex<f64>]| -> Complex<f64> |args: &[Float]| -> Float |args: &[Complex<Float>]| -> Complex<Float> ``` Those function types cannot all be stored directly in one homogeneous field[1](https://symbolica.io/posts/symbolica_2_0_release/#fn1)\. Symbolica solves this by type\-erasing the callbacks at registration time and recovering the concrete type when an evaluator for a numeric domain is built\. 1An`enum`of all possible variants cannot work, because the user can define their own type`T`and register a callback for it\. The core shape is: ``` pub trait ExternalFunction<T>: Fn(&[T]) -> T + Send + Sync + DynClone {} pub type EvalFn<T> = Box<dyn ExternalFunction<T>>; pub struct EvaluationInfo { eval_fns: HashMap<TypeId, Box<dyn Any + Send + Sync>>, } ``` When a user writes: ``` use symbolica::prelude::*; let eval = EvaluationInfo::new() .register(|args: &[f64]| 2.0 * args[0]) .register(|args: &[Complex<f64>]| args[0] + args[0]); ``` each callback is stored under`TypeId::of::<T\>\(\)`, where`T`is the numeric domain carried by the callback argument slice\. The actual callback is boxed as`Any`, which removes the concrete type from the outer container\. ``` pub fn register<T: 'static>(mut self, f: impl ExternalFunction<T> + 'static) -> Self { let f: EvalFn<T> = Box::new(f); self.eval_fns.insert( TypeId::of::<T>(), Box::new(f) as Box<dyn Any + Send + Sync>, ); self } ``` Later, when an evaluator is specialized to a domain, Symbolica asks for the implementation of that exact type: ``` impl EvaluationInfo { pub fn get_evaluator<T: 'static>(&self) -> Option<EvalFn<T>> { let e = self.eval_fns.get(&TypeId::of::<T>())?; let f = e.downcast_ref::<EvalFn<T>>()?; Some(f.clone()) } } ``` The important design point is that the type erasure is not in the evaluator hot loop\. It is in the registration and specialization path\. Once the evaluator has been built for`f64`,`Complex<f64\>`,`Float`, or another domain, it stores concrete callbacks in the instruction stream or hands concrete calls to the JIT backend\. ### Choosing an evaluation domain The other half of the evaluator story is how Symbolica chooses the numeric type in the first place\. Internally this is expressed through the`EvaluationDomain`trait\. In simplified form, the trait looks like this: ``` pub trait EvaluationDomain: Sized + 'static { fn get_evaluator(info: &EvaluationInfo) -> Option<Box<dyn ExternalFunction<Self>>> { info.get_evaluator() } } ``` Thus, for a type`T`, the default implementation fetches the callback registered for`T`in the`EvaluationInfo`, if it exists\. However, sometimes a domain can provide a useful fallback\. For example, the SIMD domain can vectorize a scalar`f64`callback when no native SIMD callback was registered, and double\-float evaluation can fall back to an arbitrary\-precision implementation before converting the result back to double\-float\. ``` impl EvaluationDomain for f64 {} impl EvaluationDomain for wide::f64x4 { fn get_evaluator( info: &EvaluationInfo, ) -> Option<Box<dyn ExternalFunction<Self>>> { if let Some(f) = info.get_evaluator::<wide::f64x4>() { return Some(f); } // create a vectorized version of the scalar function if it exists if let Some(f) = f64::get_evaluator(info) { Some(Box::new(move |args| { let mut buffer = vec![0.; args.len()]; let mut res = [0.; 4]; for i in 0..4 { for (b, v) in buffer.iter_mut().zip(args) { *b = v.as_array()[i]; } res[i] = f(&buffer); } res.into() })) } else { None } } } ``` ## AI use AI is the elephant in the room of software development, so it is worth stating how it has fit into the development of Symbolica 2\.0\. Symbolica was developed mostly without AI assistance for more than two years, and it involved a lot of planning, optimization work, and design trade\-offs\. Initially, I was skeptical about the usefulness of AI for a high\-performance project like this\. However, if a tool can do what I can do faster, it is foolish to ignore it\. So I set up isolated experiments: porting small projects from Mathematica to Symbolica, doing website and documentation work, and trying more serious implementation and research tasks\. One experiment was to let Codex 5\.5 implement a GCD algorithm from a research paper, using already written code snippets as context\. The result was interesting: it identified the key points from the paper and implemented most of the algorithm correctly, though sometimes using nonsensical helper functions\. When it encountered bugs, however, it sometimes drew the wrong conclusions and started patching increasingly specific cases instead of identifying the root cause, which was that the test input did not satisfy an assumption and triggered an infinite loop[2](https://symbolica.io/posts/symbolica_2_0_release/#fn2)\. In the end, writing critical code and debugging hard cases is still a manual process\. 2When it couldn’t figure it out, it added a fallback routine to another GCD algorithm instead of solving the actual problem\! The largest gain from using AI has been in large refactors and peripheral tasks, such as the website redesign and checking whether examples on the website still work\. That has been a huge time saver\. ## Upcoming features An upcoming feature is making the GMP backend for arbitrary\-precision arithmetic optional \(using[malachite](https://www.malachite.rs/)and[astro\-float](https://github.com/stencillogic/astro-float)instead\), which leads to a WASM\-able Symbolica build \(and easier compilation on Windows\) at the expense of performance\. Stay tuned for more details\! Symbolica is**free**for hobbyists and a single\-core instance is**free**for non\-commercial use\. If you want new features, more documentation, or support for your workflow, your organization can purchase a license and directly support the project\. ## Closing Symbolica 2\.0 makes the symbolic system more programmable: users can teach symbols how to normalize, print, differentiate, expand around singularities, and evaluate numerically\. Let me know what you think in the comments\!