Symbolica 2.0: 面向Python和Rust的可编程符号

# Symbolica 2.0：可编程符号 – Symbolica | 现代计算机代数 来源：https://symbolica.io/posts/symbolica_2_0_release/ ## 引言 Symbolica 是一个用于 Python 和 Rust 的高性能符号计算框架。你可以用它来操作符号表达式，并将其转化为用于计算雅可比矩阵、数值优化、积分等场景的快速数值内核。今天，Symbolica 2.0 版本正式发布，带来了许多令人兴奋的新特性和改进。本次发布的主题是**可编程符号**：用户可以自定义 Symbolica 更多的行为。这使得定义像内置对象一样能够简化、求导、展开、打印和求值的数学对象成为可能。 自 1.0 版本以来，Symbolica 在多个方向上积累了改进： - 更简单的 Rust API，提供更多运算符重载和构建器风格的 API - 符号注册系统，支持命名空间、别名、标签、用户数据和自定义钩子 - 重新设计的求值器接口，支持双浮点算术和即时编译（JIT） - 笔记本和文档的丰富输出，包括 HTML 输出、图形和多项式展示、Typst 输出、彩色打印以及更具结构性的多行格式化 - 新增内置数学函数，包括 gamma、多对数、贝塞尔函数、黎曼 zeta 以及相关的级数/求值钩子 请参阅**迁移指南**（https://symbolica.io/docs/migration.html）了解具体变更以及如何更新代码。 ## 更好的输出 Symbolica 增加了自动换行输出模式，并带有彩色括号，类似于代码的样式。这使得阅读大型嵌套表达式更加容易。在笔记本中，例如 Jupyter 或 Marimo（https://marimo.io/），默认输出为彩色 HTML 模式，你可以轻松切换为 LaTeX 模式。现在也支持 Typst 输出。 ## 改进的 Rust API 最明显的变化之一是普通 Rust 程序所需的导入和冗长的类型路径大大减少。新的预导出版本收集了大多数用户需要的常见 trait、宏、域和求值器类型。Rust 的易用性也得到了改善，包括额外的重载、自动类型转换、构建器模式以及符号上的 `call` 方法： - Symbolica 2.0 - Symbolica 1.0 ```rust 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+O(x^2) } ``` ```rust 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+O(x^2) } ``` 现在，用于设置和接受多个参数的函数采用构建器模式。例如，构造一个高性能数值求值器现在的代码如下： - Symbolica 2.0 - Symbolica 1.0 ```rust 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(()) } ``` ```rust 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, ¶ms, optimization_settings) .unwrap() .map_coeff(&|c| c.re.to_f64()); println!("{}", evaluator.evaluate_single(&[3.0])); } ``` 通过构建器模式，原本属于不同子结构（`add_function` 属于 `FunctionMap`，`horner_iterations` 属于 `OptimizationSettings`）的设置参数现在可以更流畅地设置在一起。还要注意，可失败的操作现在返回一个专用的错误类型。 ## 可编程符号 在 Symbolica 1.0 中，符号已经可以携带重要的代数属性： - Python - Rust ```python 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) ``` ```rust 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) } ``` 在 2.0 版本中，符号现在可以安装更多钩子，这些钩子在代数生命周期的特定点执行： | 钩子 | 用途 | |------------|----------------------------------------------------------------------| | 规范化 | 在符号或函数被规范化时重写。 | | 打印 | 覆盖 Symbolica、LaTeX 或自定义模式的输出。 | | 导数 | 定义自定义导数规则。 | | 级数 | 教会 Symbolica 函数在奇点附近如何表现。 | | 求值 | 为求值器注册数值实现。 | 例如，让我们定义一个 gamma 函数（https://en.wikipedia.org/wiki/Gamma_function）。由于 gamma 在 0 处有一个极点，我们无法在 0 附近进行泰勒展开。使用级数钩子，我们可以告诉 Symbolica 如何通过恒等式 \( \Gamma(a+1) = a \Gamma(a) \) 在该点附近正则化函数： - Python - Rust ```python 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)+O(x^1) ``` ```rust 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()); } ``` 对于这个例子，内置的 gamma 函数会进一步化简为 `1/x - γ`。 ## 求值器：从表达式树到编译内核 自 1.0 版本以来，表达式求值经历了最大的工程投入。高级流程仍然相同：Symbolica 将表达式重写为一个小的指令程序，对其进行优化，然后针对不同的数值输入多次求值。请参阅这篇博客文章（https://symbolica.io/posts/expression_evaluation/index.html）了解其工作原理的详细信息。 ### 符号的求值器 让我们定义一个自定义符号 `cosh`，并教会 Symbolica 如何针对不同的数值域对其进行求值。我们可以通过注册一个求值钩子来实现： - Python - Rust ```python 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 ``` ```rust use symbolica::prelude::*; fn main() { let cosh = symbol!( "cosh", eval = EvaluationInfo::new() .register(|args: &[f64]| args[0].cosh()) .register(|args: &[Complex]| args[0].cosh()) ); let x = symbol!("x"); let expr = cosh.call(Atom::num(1) / 2) + x; println!("{}", expr.to_float(16)); // -> 1.127625965206381+x } ``` 由于已经知道如何求值 `cosh`，Symbolica 会在 `to_float` 调用中查找合适的求值器。 ### JIT 编译 除了生成自定义的 ASM、C++ 和 CUDA 代码外，Symbolica 现在还可以即时编译求值器。这使用了 symjit（https://github.com/siravan/symjit）crate，作者是 Shahriar Iravanian，我们与其密切合作，使集成在 Symbolica 中感觉非常自然。JIT 路径支持自定义求值器钩子，并且现在是 Python 的默认求值后端。在实践中，它在保持编译时间可控的同时，与自定义 ASM 后端具有竞争力。 ### 双浮点算术 Symbolica 2.0 增加了双浮点求值路径。双浮点算术将一个数存储为两个 `f64` 的未求和，提供大约 106 位精度（31 位十进制数字，而普通双精度为 16 位），同时比任意精度的 `Float` 算术快 3 倍以上。在 Python 中，`evaluate_with_prec(..., 32)` 会自动选择这条路径。 ## 特殊函数 新的钩子系统还使得 Symbolica 拥有更丰富的数学词汇表成为可能。Symbolica 现在支持 polygamma 函数、多对数函数、贝塞尔函数、黎曼 zeta、几何函数及其反函数，并带有求值钩子和在极点附近的 Laurent/Puiseux 级数行为。一些特殊值会立即化简： ```python from symbolica import * print(E("gamma(1/2)")) # π^(1/2) print(E("polylog(-3, z)")) # z*(1+4*z+z^2)/(1-z)^4 ``` 所有特殊函数都可以求值： ```python g = E("gamma(5/6)") print(g) # gamma(5/6) print(g.to_float()) # 1.128787029908126 ``` Symbolica 可以在优化的求值器中使用特殊函数： ```python 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]] ``` 在导出代码时，诸如 \( \pi \) 和 \( \zeta(3) \) 之类的常量会被其数值替换，因此生成的代码不需要依赖特殊函数库。 ## 更快的操作 并非所有重要的变化都体现在公共 API 中。表达式操作底层机制也进行了实质性改进： - 模式匹配可以证明结构性不匹配并提前终止 - 项排序现在可以基于字节切片比较（规范化速度提升 30%） - 霍纳方案可以应用于线性表达式结构（内存减少 >2 倍） - 公共对的消除使用更少的内存 - 改进的多项式 GCD 计算，包括一个新的模块化算法 这些综合变化在实际用例中使性能提升了 2 到 10,000 倍。对于这些改进，用户的反馈至关重要。如果你有计算慢的情况，请与我们分享，以便我们让 Symbolica 更快。 ## 技术深入探讨：类型擦除的求值回调 本节面向技术好奇的读者。它解释了 Symbolica 的求值钩子如何在不牺牲易用性或性能的情况下支持多个数值域。 求值钩子 API 内部有一个有趣的问题：一个符号可能有多个数值实现，每个实现对应一个不同的 Rust 类型，该类型实现了一个特定的 trait，其成员函数形式为 `fn(&[T]) -> T`，其中 `T` 是某个数值域。例如，一个函数可能有针对以下类型的独立快速实现： ```rust |args: &[f64]| -> f64 |args: &[Complex]| -> Complex |args: &[Float]| -> Float |args: &[Complex]| -> Complex ``` 这些函数类型不能直接存储在一个同质的字段中¹。Symbolica 通过在注册时擦除回调的类型，并在为某个数值域构建求值器时恢复具体类型来解决这个问题。 ¹ 所有可能变体的枚举是不可能的，因为用户可能定义自己的类型 `T` 并为其注册回调。 核心形状如下： ```rust 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>>, } ``` 当用户编写： ```rust use symbolica::prelude::*; let eval = EvaluationInfo::new() .register(|args: &[f64]| 2.0 * args[0]) .register(|args: &[Complex]| args[0] + args[0]); ``` 每个回调都存储在 `TypeId::of::<T>()` 下，其中 `T` 是回调参数切片携带的数值域。实际的回调被装箱为 `Any`，从而从外部容器中移除具体类型。 ```rust 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 } ``` 稍后，当求值器特化为某个域时，Symbolica 会请求该确切类型的实现： ```rust 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()) } } ``` 重要的设计点在于，类型擦除不在求值器的热循环中。它在注册和特化路径中。一旦求值器为 `f64`、`Complex`、`Float` 或其他域构建完成，它就会在指令流中存储具体的回调，或者将具体的调用交给 JIT 后端。 ### 选择求值域 求值器故事的另一个方面是 Symbolica 如何首先选择数值类型。内部通过 `EvaluationDomain` trait 来表达。简化形式如下： ```rust pub trait EvaluationDomain: Sized + 'static { fn get_evaluator<T: 'static>(info: &EvaluationInfo) -> Option<EvalFn<T>> { info.get_evaluator::<T>() } } ``` 因此，对于类型 `T`，默认实现会获取在 `EvaluationInfo` 中为 `T` 注册的回调（如果存在）。然而，有时一个域可以提供有用的回退。例如，SIMD 域可以在没有注册原生 SIMD 回调时向量化标量 `f64` 回调；双浮点求值可以回退到任意精度的实现，然后再将结果转换回双浮点。 ```rust impl EvaluationDomain for f64 {} impl EvaluationDomain for wide::f64x4 { fn get_evaluator<T: 'static>( info: &EvaluationInfo, ) -> Option<EvalFn<T>> { if let Some(f) = info.get_evaluator::<wide::f64x4>() { return Some(f); } // 如果存在标量函数，则创建一个向量化版本 if let Some(f) = f64::get_evaluator::<f64>(info) { Some(Box::new(move |args: &[f64x4]| { let mut buffer = vec![0.; args.len() * 4]; 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 的使用 AI 是软件开发中房间里的大象，因此值得说明它在 Symbolica 2.0 开发中的角色。Symbolica 的开发在两年多的时间里基本上没有借助 AI，这涉及大量的规划、优化工作和设计权衡。起初，我对于 AI 在这样一个高性能项目中的实用性持怀疑态度。然而，如果一个工具可以比我更快地完成我能做的事情，忽视它是不明智的。因此，我设置了一些独立的实验：将小型项目从 Mathematica 移植到 Symbolica，进行网站和文档工作，并尝试更严肃的实现和研究任务。 一个实验是让 Codex 5.5 根据一篇研究论文实现一个 GCD 算法，使用已经写好的代码片段作为上下文。结果很有趣：它识别了论文中的关键点，并正确实现了大部分算法，尽管有时会使用无意义的辅助函数。然而，当遇到错误时，它有时会得出错误的结论，并开始修补越来越特定的情况，而不是识别根本原因——即测试输入不满足某个假设，导致无限循环²。最终，编写关键代码和调试困难问题仍然是一个手动过程。 ² 当它无法解决时，它添加了