# Eigenmath-Plus

**Repository Path**: rabbitc999/eigenmath-plus

## Basic Information

- **Project Name**: Eigenmath-Plus
- **Description**: Eigenmath-Plus是一个轻量级、可嵌入的符号计算系统，支持代数运算、微积分、线性代数、常微分方程求解、数值计算及多种工程数学变换。它采用 C 语言编写，提供简洁的交互式命令行界面和脚本执行能力，适用于教学、科研以及作为复杂系统的计算后端。
- **Primary Language**: Unknown
- **License**: MIT
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2026-08-30
- **Last Updated**: 2026-08-31

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# Eigenmath-Plus

**Eigenmath-Plus** 是一个轻量级、可嵌入的符号计算系统（纯 C 实现），在经典 **Eigenmath** 的基础上扩充了代数、微积分、线性代数、常微分方程、数值计算与工程变换能力。它既能作为桌面/命令行计算工具，也能嵌入到 MCU（单片机）计算器、仪器仪表和教学设备中运行：不依赖外部库、数字运算不占用堆内存、代码与内存占用可控，非常适合资源受限的嵌入式环境。

---

## 特性概览

- **符号代数**：多项式展开、因式分解（实域/复域/多元立方和差、$x^n \pm 1$ 分圆多项式分解、符号失败时自动 `nroots` 数值回退）、多元多项式 Gröbner 基（Buchberger 算法）、方程求解、表达式化简、数论与组合（gcd/lcm、素数判定、斐波那契、二项式系数、欧拉函数等）、模算术与有限域（模幂、模逆、中国剩余定理、多项式取模）、精确复数运算（欧拉公式、复对数、复三角/双曲函数）。
- **微积分**：符号求导（包括隐函数）、不定积分（换元、分部积分、有理函数）、定积分（含反常积分、变限积分、含参变量积分及莱布尼茨求导）、多元微积分（梯度/雅可比/海森矩阵/散度/旋度）、数值积分。
- **极限与级数**：极限计算（洛必达法则支持）、泰勒/洛朗展开、有限及无穷级数求和、带符号上下限的闭式求和、超几何函数 ${}_2F_1$。
- **线性代数**：矩阵行化简（RREF）、秩、零空间、行列式、逆矩阵、LU/QR 分解、特征分解（`eigen`）、奇异值分解（`svd`）、若尔当分解（`jordan`）、矩阵范数与条件数（`norm`/`cond`）、伪逆、矩阵指数/对数/平方根。
- **常微分方程**：符号求解（可分离变量、一阶线性、二阶常系数齐次/部分非齐次）及数值求解（RK4）。
- **数值工具**：割线法求根、多项式数值求根（`nroots`，含复根/重根）、隐函数数据生成、函数采样（一维/二维网格）、描述统计（均值/方差/标准差）与正态分布（PDF/CDF/分位数）。
- **变换**：拉普拉斯变换、傅里叶变换、Z 变换、信号函数（sinc/阶跃/冲激/卷积/角度单位）。
- **假设系统**：支持对变量添加正性假设，简化绝对值等表达式。
- **数值稳健性**：整数内部统一为 64 位无符号，四则/乘方/阶乘溢出自动转双精度浮点；数值求逆带条件数检查，奇异/病态返回 `nil` 并警告。
- **脚本支持**：通过 `.in` 文件批量执行命令，方便自动化计算和回归测试。

---

## MCU 用途与设计取舍

本项目在面向单片机计算器场景时做了针对性设计，目标是在 Cortex-M 等 MCU 上以有限的 RAM / Flash 提供“够用”的科学计算能力。

### 适用场景

- **手持 / 桌面计算器固件**：交互输入表达式，输出分数、矩阵、微积分结果；
- **仪器仪表与教学设备**：作为嵌入式计算内核，通过按键 / 串口 / 显示屏调用；
- **离线批处理**：在设备上执行预置脚本（查表、单位换算、工程计算等）。

### 面向 MCU 的设计

- **数字运算不依赖堆内存**：整数与有理数直接以 64 位无符号整数内联存储，不再为每个数字分配动态内存，避免小堆碎片化与分配失败；
- **资源占用确定**：删除任意精度大整数后，Flash 与 RAM 占用显著下降，行为可预测，适合裸机或 RTOS 环境；
- **超范围自动转浮点**：四则运算、乘方、阶乘等结果超过 $0 \sim 2^{64}-1$ 时自动转为双精度浮点继续计算；数论函数（`gcd`、`powermod`、`isprime` 等）仍要求精确整数，不会用浮点给出错误结果；
- **纯 C、单线程、无外部依赖**：只依赖标准 C 库与数学库，容易用 `arm-none-eabi-gcc` 等交叉编译器移植。

### MCU 交叉编译示例

```sh
arm-none-eabi-gcc -mcpu=cortex-m4 -mthumb -Os \
  -ffunction-sections -fdata-sections \
  main.c eigenmath_*.c eigenmath_data.c -lm -o eigenmath.elf
```

嵌入时通过 C API `EigenMath_Calculate(in, out, outSize)` 以字符串方式执行脚本并捕获输出（返回 0 成功 / 1 出错，错误信息写入 `out`），也可使用 `EigenMath_Run(buf)` / `EigenMath_RunBuf(buf)` 直接执行。需要进一步裁剪 Flash 时，可按需精简 `eigenmath_data.c` 中的函数表项（例如积分表、LaTeX 映射表）。

---

## 构建与安装

### 依赖项
- C 编译器（如 `gcc` 或 `clang`）
- 标准数学库（`libm`）
- `make`（可选，用于便捷构建）

### 编译步骤

```sh
git clone https://gitee.com/rabbitc999/eigenmath-plus.git   # 替换为实际仓库地址
cd eigenmath
make          # 编译生成 ./eigenmath
```

若需手动编译：
```sh
cc -Wall -Wextra -O2 main.c eigenmath_*.c eigenmath_data.c -lm -o eigenmath
```

---

## 使用方式

### 配置与内存池

需要按目标平台调整的参数集中在 `eigenmath_config.h` 中：

- 引擎参数：`STACKSIZE`（求值栈容量）、`BLOCKSIZE`/`MAXBLOCKS`（原子节点内存块）、`BUCKETSIZE`（符号表桶容量）、`STRBUFLEN`、`MAXDIM`、`EIGENMATH_MAX_ASSUMPTIONS`；
- 内存池参数：`EIGENMATH_MEM_POOL_SIZE`（池总大小）、`EIGENMATH_MEM_BLOCK_SIZE`（块大小）、`EIGENMATH_MEM_TABLE_TYPE`（管理表类型）。

引擎的全部动态内存（原子节点、字符串、输出缓冲等）均由 `eigenmath_mem_pool.c` 提供的固定容量内存池分配，不依赖系统 `malloc`/`free`，便于移植到 MCU 裸机环境。修改配置后重新编译即可；池过小会导致计算中途 `Stop: out of memory`，请按实际负载调整（默认 8 MB 可覆盖全量回归与万能公式积分用例；MCU 可按板载内存调小）。

### 三角函数万能公式与缓冲区清理

- `weierstrass(expr, x)`：把表达式中的 `sin(x)/cos(x)/tan(x)/sec(x)/csc(x)/cot(x)` 按万能公式（半角代换 `t = tan(x/2)`）化为 `t` 的有理式：
  `sin(x) = 2t/(1+t²)`、`cos(x) = (1-t²)/(1+t²)`、`tan(x) = 2t/(1-t²)` 等；
- 积分器已内置该代换：对 `sin/cos/tan/sec/csc/cot(x)` 的有理函数被积式，先作 `t = tan(x/2)` 代换化为有理函数积分，再代回，例如
  `integral(1/(1+sin(x)),x)` → `-tan(1/4 pi - 1/2 x)`，`integral(1/(sin(x)+cos(x)),x)` → `log(...)`；
- `EigenMath_BufferReset()`：清理输出缓冲区 `outbuf` 与词法缓冲区 `token_buf`，把内存归还内存池；每次 `EigenMath_Run` 开始时自动调用，嵌入式反复调用运行时可避免缓冲区累积。
- `EigenMath_Calculate(in, out, outSize)`：字符串输入/字符串输出接口，内部自行建立错误跳转环境并捕获输出，不暴露 `jmp_buf`；
- `EigenMath_EnvSetup()` / `EigenMath_EnvCleanup()`：主运行环境（setjmp/longjmp）的封装函数，供需要批量错误恢复的主机程序使用。
- `EigenMath_SetUserVariable(name, expr)` / `EigenMath_GetUserVariable(name, out, outSize)`：用户变量的设置（解析并求值表达式后绑定）与获取（格式化为字符串）。

### 大整数、大浮点数与有效位数

- `USE_BIG_INT`：`eigenmath_config.h` 中的可配置宏（默认 0）。置 1 后，超出 64 位的整数自动使用任意精度大整数，乘法内核采用 NTT（数论变换，双素数 CRT 精确重建）加速；`eigenmath_bigint.c/.h` 提供独立的 `BigInt_T`/`BigFloat_T` 库接口（字符串转换、加减乘除、幂、比较等）；
- 大整数示例：`2^1000`、`factorial(50)`、超大整数字面量及其加减乘、比较均精确；
- 大有理数（BIGRATIONAL）：置 1 后，普通 64 位有理数（原子 `RATIONAL`，`u.q.a/u.q.b` 为 `uint64_t`）在四则运算/乘方结果超出 64 位时自动升级为大有理数（原子 `BIGRATIONAL`，分子分母均为 `BigInt_T *`），例如 `1/(2^100)`、`1/999999999999999999999999999999`、`(2^100+1)/(2^50+1)`、`2^-100`、`(3/2)^(-100)` 均精确表示并按 `num/den` 输出；负指数直接生成大有理数，不再退化为浮点；`numerator`/`denominator`/`floor`/`ceiling`/`round`/`numden` 等函数与精确比较（`<`、`>`、`==`）均支持大有理数；
- `digital`：浮点输出有效位数变量（默认 6）。`digital=30` 后浮点结果按 30 位有效数字输出；`eigenmath_bigint.h` 的 `BigFloat_T` 支持任意精度十进制浮点（尾数 × 10^指数，加减乘除按 `digits` 舍入）。

### 交互模式

直接在终端运行 `./eigenmath`，进入交互式提示符（`?`）。输入表达式或命令即可立即计算。

```text
? sin(pi/2)
> 1
? A = [[1,2],[3,4]]
> [[1,2],[3,4]]
? det(A)
> -2
```

退出交互模式：输入 `exit` 或按 `Ctrl+D`。

### 脚本模式

将命令序列写入 `.in` 文件，每行一条命令，通过重定向执行：

```sh
./eigenmath myscript.in
```

脚本执行后所有输出会打印到标准输出，常用于自动化测试或批量数据处理。

### 输出控制

- 结果以纯文本形式输出，矩阵、向量、分数等使用线性表示。
- 设置 `latex=1` 后结果以 LaTeX 代码输出，`latex=0` 恢复线性输出。

---

## 命令速查表

以下所有命令遵循 `command(arg1, arg2, ...)` 的 Lisp 风格语法，参数在括号内用逗号分隔。

### 1. 代数与多项式

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `solve(expr, var)` | `expr`: 方程或表达式（隐式 $=0$）；`var`: 符号变量 | 返回根列表（实数/复数），或解区间矩阵（不等式）。 |
| `factor(expr, var)` | `expr`: 多项式；`var`: 变量（可选，多元时自动识别） | 实数域因式分解。自动识别 $x^3 \pm y^3$ 与 $x^n \pm 1$（分圆多项式）。 |
| `cfactor(expr, var)` | 同上 | 复数域因式分解（引入虚数单位 $i$）。 |
| `collect(expr, var)` | `expr`: 表达式；`var`: 变量 | 按 `var` 的幂次合并同类项（如 $x + 2x \to 3x$）。 |
| `combine(expr)` | `expr`: 表达式 | 合并对数和指数（$\log a + \log b \to \log(ab)$）。 |
| `gcd(a, b, ...)` | 两个或多个整数 | 最大公约数。 |
| `lcm(a, b)` | 两个整数 | 最小公倍数。 |
| `isprime(n)` | `n`: 整数 | 返回 `1`（素数）或 `0`。 |
| `nextprime(n)` | `n`: 整数 | 返回大于 $n$ 的下一个素数。 |
| `choose(n, k)` / `binomial(n, k)` | $n, k$: 整数或符号 | 二项式系数 $C(n, k)$（$n$ 为符号时支持展开）。 |
| `perm(n, k)` | $n, k$: 整数或符号 | 排列数 $P(n, k)$。 |
| `fibonacci(n)` | $n$: 非负整数 | 第 $n$ 个斐波那契数。 |
| `totient(n)` | $n$: 正整数 | 欧拉 $\varphi$ 函数（小于 $n$ 且与 $n$ 互质的个数）。 |
| `divisors(n)` | $n$: 正整数 | 返回所有正因子的向量。 |
| `bernoulli(n)` | $n$: 非负整数 | 伯努利数（$B_1 = -1/2$）。 |
| `harmonic(n)` | $n$: 非负整数 | 调和数 $H_n$（$n$ 为符号时保留原式）。 |
| `powermod(a, e, m)` | 整数 $a, e, m$ | 模幂 $a^e \bmod m$。 |
| `invmod(a, m)` | 整数 $a, m$（互质） | 模逆元 $a^{-1} \bmod m$。 |
| `crt(a1, m1, a2, m2)` | 整数对 | 中国剩余定理（解同余式组）。 |
| `polymod(P, M)` | $P, M$: 多项式表达式 | 长除法余式 $P \bmod M$（$\mathbb{Q}[x]$ 上的除法）。 |
| `groebner(polys, vars)` | `polys`: 多项式向量；`vars`: 变量向量 | 返回约化 Gröbner 基（lex 序）。 |
| `trigexpand(expr)` | `expr`: 三角表达式 | 展开角度和差（$\sin(x+y)$）。 |
| `trigreduce(expr)` | `expr`: 三角表达式 | 降幂（$\sin^2 x \to \frac{1-\cos 2x}{2}$）。 |
| `trigsimplify(expr)` | `expr`: 三角表达式 | 应用 $\sin^2+\cos^2=1$ 等恒等式化简。 |

### 2. 微积分与极限

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `d(expr, var)` 或 `derivative(expr, var)` | `expr`: 函数；`var`: 变量 | 一阶导数。支持嵌套求高阶导（如 `d(d(sin(x),x),x)`）。 |
| `implicitdiff(F, x, y)` 或 `idiff(F, x, y)` | $F$: 隐函数方程；$x$: 自变量；$y$: 因变量 | 返回 $\frac{dy}{dx} = -\frac{F_x}{F_y}$。 |
| `grad(F, X)` | $F$: 标量函数；$X$: 变量向量 $[x,y]$ | 梯度向量。 |
| `jacobian(F, X)` | $F$: 函数向量；$X$: 变量向量 | 雅可比矩阵。 |
| `hessian(F, X)` | $F$: 标量函数；$X$: 变量向量 | 海森矩阵。 |
| `div(F, X)` | $F$: 向量场；$X$: 变量向量 | 散度。 |
| `curl(F, X)` | $F$: 3D 向量场；$X$: $[x,y,z]$ | 旋度。 |
| `integral(expr, var)` | `expr`: 被积函数；`var`: 积分变量 | 不定积分（结果不含常数 $C$）。 |
| `defint(expr, var, a, b)` | `expr`: 被积函数；`var`: 变量；$a,b$: 上下限 | 定积分（支持符号上下限和无穷 $\infty$）。 |
| `nintegrate(expr, var, a, b)` | `expr`: 数值函数；`var`: 变量；$a,b$: 浮点数 | 自适应 Simpson 数值积分。 |
| `limit(expr, var, point)` | `expr`: 函数；`var`: 变量；`point`: 趋近点（可为 $\infty$） | 极限值（自动洛必达）。 |
| `taylor(expr, var, order)` | `expr`: 函数；`var`: 变量；`order`: 展开阶数 | 在 $x=0$ 处的泰勒级数（有限项）。 |
| `laurent(expr, var, point, order)` | `expr`: 函数；`var`: 变量；`point`: 展开点；`order`: 阶数 | 洛朗级数展开。 |
| `sum(expr, var, a, b)` | `expr`: 通项；`var`: 求和变量；$a,b$: 起止（可为 $\infty$） | 闭式求和（幂和、等比）或保留原式。 |

### 3. 线性代数

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `rref(A)` | $A$: 矩阵 | 行最简阶梯形。 |
| `rank(A)` | $A$: 矩阵 | 秩。 |
| `nullspace(A)` | $A$: 矩阵 | 零空间基（每列一个向量）。 |
| `det(A)` | $A$: 方阵 | 行列式。 |
| `inv(A)` | $A$: 方阵 | 逆矩阵；奇异或病态（条件数 $>10^6$）时返回 `nil` 并警告。 |
| `lu(A)` | $A$: 方阵 | 返回 $[L, U]$（部分选主元）。 |
| `qr(A)` | $A$: 矩阵 | 返回 $[Q, R]$（Gram-Schmidt）。 |
| `eigenvals(A)` | $A$: 数值方阵 | 特征值向量（实部降序）。 |
| `eigenvecs(A)` | $A$: 数值方阵 | 特征向量矩阵（每列一个向量）。 |
| `eigen(A)` | $A$: 数值方阵 | 返回 $[\text{特征值向量}, \text{特征向量矩阵}]$。 |
| `jordan(A)` | $A$: 数值方阵（$\le 5\times 5$） | 若尔当分解，返回 $[P, J]$ 满足 $A = P J P^{-1}$。 |
| `svd(A)` | $A$: 数值矩阵 | 返回 $[U, S, V]$（经济分解，降序）。 |
| `linsolve(A, b)` | $A$: 系数矩阵；$b$: 向量或矩阵 | 解 $Ax=b$（数值 SVD 或符号 RREF）。 |
| `cross(u, v)` | $u, v$: 3维向量 | 叉积。 |
| `dot(u, v)` 或 `inner(u, v)` | $u, v$: 向量或矩阵 | 点积 / 矩阵乘法。 |
| `diag(a, b, ...)` 或 `diag([a,b,...])` | 标量列表或向量 | 对角矩阵。 |
| `charpoly(A)` 或 `charpoly(A, x)` | $A$: 方阵；$x$: 可选变量名 | 特征多项式 $\det(A - xI)$。 |
| `pseudoinverse(A)` | $A$: 矩阵 | Moore-Penrose 伪逆（满秩情形）。 |
| `norm(v)` / `norm(v, p)` | $v$: 向量或矩阵；$p$: $1, 2, \infty$ | 范数：向量默认欧几里得，矩阵默认 Frobenius；$1$/$\infty$ 为行列和范数，$2$ 为谱范数（最大奇异值）。 |
| `cond(A)` / `cond(A, p)` | $A$: 数值方阵；$p$: $1, 2$（默认）, $\infty$ | 条件数：默认谱条件数，$1$/$\infty$ 为相应范数条件数；奇异矩阵返回 $\infty$。 |
| `kron(A, B, ...)` 或 `kronecker(A, B, ...)` | 矩阵或向量 | 克罗内克积（链式支持）。 |
| `rowswap(A, i, j)` | $A$: 矩阵；$i, j$: 行号（1-based） | 交换两行，返回新矩阵。 |
| `rowadd(A, i, j, k)` | $A$: 矩阵；$i, j$: 行号；$k$: 系数 | $R_i \leftarrow R_i + k R_j$，返回新矩阵。 |
| `rowmul(A, i, k)` | $A$: 矩阵；$i$: 行号；$k$: 系数 | $R_i \leftarrow k R_i$，返回新矩阵。 |
| `colswap`, `coladd`, `colmul` | 参数同上（列操作） | 同理，操作列。 |

### 4. 常微分方程（ODE）

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `dsolve(ode, x)` | `ode`: ODE 方程（使用 `d(y(x),x)`）；`x`: 自变量 | 通解（含 $C_1, C_2$ 常数）。支持初值条件：`dsolve(ode, x, y(a)=b, eval(d(y(x),x),x,a)=b)`。 |
| `dsolve_system(d(Y,x) = A*Y, x)` | $Y$: 向量 $[y_1,y_2]$；$A$: 2×2 矩阵 | 线性系统通解（含实/复/重根）。 |
| `odesolve(f, x, x0, x1, y0, n)` | $f$: 函数或向量函数；$x$: 变量；$x0,x1$: 区间端点；$y0$: 初值；$n$: 步数 | RK4 数值解，返回矩阵（每行 $x, y_1, \ldots$）。 |
| `odesolve_adaptive(f, x, x0, x1, y0, tol)` | $f$: 函数；`tol`: 误差容忍度 | 自适应 RK45，返回自适应点矩阵。 |
| `slopefield(f, x, a, b, nx, y, c, d, ny)` | $f$: $\frac{dy}{dx}$；$x,y$: 变量；$a,b,c,d$: 边界；$nx,ny$: 网格数 | 返回方向场向量 $(x, y, dx, dy)$。 |
| `expm(A)` 或 `expm(A, x)` | $A$: 方阵；$x$: 可选标量 | 矩阵指数 $e^A$ 或 $e^{Ax}$（2×2 符号，高阶数值）。 |
| `logm(A)` | $A$: 数值方阵 | 矩阵对数（主分支）。 |
| `sqrtm(A)` | $A$: 数值方阵 | 矩阵平方根（主分支）。 |

### 5. 数值方法与数据生成

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `findroot(f, var, guess)` | $f$: 函数；`var`: 变量；`guess`: 浮点数初值 | 割线法求根。 |
| `nroots(poly, var)` | `poly`: 多项式；`var`: 变量 | 返回全部数值复根向量（含重根）。 |
| `fsolve(F, X, G)` | $F$: 函数向量；$X$: 变量向量；$G$: 初值向量 | 多维牛顿法求解 $F(X)=0$。 |
| `implist(F, x, a, b, n, y_guess)` | $F$: 隐函数；$x$: 自变量；$a,b$: 区间；$n$: 点数；`y_guess`: 初值 | 生成隐函数离散点列表。 |
| `plotdata(f, x, a, b, n)` | $f$: 函数；$x$: 变量；$a,b$: 区间；$n$: 步数 | 生成 $(x, f(x))$ 点矩阵（$n+1$ 个点）。 |
| `surfdata(F, x, a, b, nx, y, c, d, ny)` | $F$: 二元函数；$x,y$: 变量；$a,b,c,d$: 边界；$nx,ny$: 网格 | 生成 2D 网格采样矩阵 $(ny+1) \times (nx+1)$。 |
| `range(a, b, step)` | $a,b$: 数值；`step`: 步长（默认 $\pm 1$） | 等差数列向量。 |
| `seq(expr, var, start, end)` | `expr`: 表达式；`var`: 整数变量；`start`,`end`: 整数 | 序列向量。 |
| `diffseq(v)` | $v$: 向量 | 前向差分 $[v_2-v_1, v_3-v_2, \ldots]$。 |
| `cumsum(v)` | $v$: 向量 | 累积和。 |
| `recurrence(init, f, n)` | `init`: 初始向量；$f$: 二元函数；$n$: 项数 | 递推生成序列。 |

### 6. 工程变换与信号

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `laplace(f, t, s)` | $f$: 时域函数；$t$: 时域变量；$s$: 复频域变量 | 拉普拉斯变换。 |
| `ilaplace(F, s, t)` | $F$: 频域函数；$s$: 频域变量；$t$: 时域变量 | 逆拉普拉斯变换。 |
| `fourier(f, t, w)` | $f$: 时域函数；$t$: 变量；$w$: 频域变量 | 傅里叶变换。 |
| `ifourier(F, w, t)` | $F$: 频域函数；$w$: 频域变量；$t$: 时域变量 | 逆傅里叶变换。 |
| `ztransform(f, n, z)` | $f$: 序列；$n$: 离散变量；$z$: 复变量 | Z 变换。 |
| `iztransform(F, z, n)` | $F$: Z 域函数；$z$: 变量；$n$: 离散变量 | 逆 Z 变换。 |
| `sinc(x)` | $x$: 数值或符号 | $\frac{\sin x}{x}$（$x=0$ 时为 1）。 |
| `heaviside(x)` | $x$: 数值或符号 | 阶跃函数（$x=0$ 时为 $1/2$）。 |
| `dirac(x)` | $x$: 数值或符号 | 狄拉克 $\delta$ 函数（$x \ne 0$ 时为 0）。 |
| `convolve(u, v)` | $u, v$: 一维向量 | 离散卷积。 |
| `deg` / `rad` | 无参数（常量） | $\mathrm{deg} = \pi/180$；$\mathrm{rad} = 1$（角度单位）。 |

### 7. 复数与复分析

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `cplotdata(f, z, a, b, c, d, nx, ny)` | $f$: 复函数；$z$: 复变量；$a,b,c,d$: 实虚部边界；$nx,ny$: 网格 | 返回 $(x, y, \Re f, \Im f)$ 数据矩阵。 |
| `residue(f, z, z0)` | $f$: 复函数；$z$: 变量；$z_0$: 极点位置 | 留数（本质奇点保留原式）。 |
| `poles(f, z)` | $f$: 复函数；$z$: 变量 | 极点列表（解分母 $=0$）。 |
| `zeros(f, z)` | $f$: 复函数；$z$: 变量 | 零点列表（解分子 $=0$）。 |
| `rect(expr)` | `expr`: 复数表达式 | 转换为 $a + b i$ 矩形形式。 |
| `polar(expr)` | `expr`: 复数表达式 | 转换为 $r e^{i\theta}$ 极坐标形式。 |
| `real(expr)` | `expr`: 复数表达式 | 提取实部。 |
| `imag(expr)` | `expr`: 复数表达式 | 提取虚部。 |
| `conj(expr)` | `expr`: 复数表达式 | 共轭。 |
| `arg(expr)` | `expr`: 复数表达式 | 辐角（主值，$(-\pi, \pi]$）。 |
| `mag(expr)` | `expr`: 复数表达式 | 模长。 |

### 8. 特殊函数与统计

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `hypergeom(a, b, c, x)` | $a,b,c,x$: 参数 | 超几何函数 ${}_2F_1$。 |
| `zeta(s)` | $s$: 数值或整数 | 黎曼 $\zeta$ 函数（整数精确值）。 |
| `besselj(n, x)` | $n$: 整数阶；$x$: 数值 | 第一类贝塞尔函数。 |
| `bessely(n, x)` | $n$: 整数阶；$x$: 数值 | 第二类贝塞尔函数。 |
| `ellipk(k)` | $k$: 数值 | 第一类完全椭圆积分。 |
| `ellipe(k)` | $k$: 数值 | 第二类完全椭圆积分。 |
| `laguerre(x, n, m)` | $x$: 变量；$n$: 阶数；$m$: 参数 | 广义拉盖尔多项式。 |
| `legendre(f, n, m, x)` | $f$: 求值点；$n,m$: 阶数；$x$: 变量 | 连带勒让德函数。 |
| `hermite(x, n)` | $x$: 变量；$n$: 阶数 | 厄米多项式。 |
| `mean(v)` | $v$: 数值向量 | 算术均值。 |
| `variance(v)` | $v$: 数值向量 | 总体方差。 |
| `stddev(v)` | $v$: 数值向量 | 总体标准差。 |
| `normal(x)` / `normal(x, mu, sigma)` | $x$: 数值；$\mu$: 均值（默认 0）；$\sigma$: 标准差（默认 1） | 正态分布 PDF。 |
| `normalcdf(x)` / `normalcdf(x, mu, sigma)` | 同上 | 正态分布 CDF。 |
| `normalinv(p)` / `normalinv(p, mu, sigma)` | $p$: 概率 $(0,1)$ | 正态分布分位数（数值求根）。 |

### 9. 假设、格式与系统控制

| 命令格式 | 参数说明 | 返回值 / 说明 |
| :--- | :--- | :--- |
| `assume(cond)` | `cond`: $x > 0$, $x < 0$, `assume(n, integer)`, `assume(y, real)`, `assume(z, nonzero)` | 添加假设（影响 $\lvert x\rvert$, $\sin(n\pi)$ 等）。 |
| `clear()` | 无参数 | 清除所有假设和用户定义变量（恢复内置常量）。 |
| `eval(expr)` | `expr`: 表达式 | 强制再次求值（用于 `quote` 内部）。 |
| `quote(expr)` | `expr`: 表达式 | 返回表达式本身而不求值。 |
| `float(expr)` | `expr`: 表达式 | 将所有有理数和 $\pi$ 转为浮点数。 |
| `simplify(expr)` | `expr`: 表达式 | 尝试化简（含指数/三角/有理化）。 |
| `expform(expr)` | `expr`: 表达式 | 将三角函数转为指数形式。 |
| `logform(expr)` | `expr`: 表达式 | 将反三角函数转为对数形式。 |
| `infixform(expr)` | `expr`: 表达式 | 返回中缀格式字符串（用于调试）。 |
| `prefixform(expr)` | `expr`: 表达式 | 返回前缀格式（Lisp 风格）字符串。 |
| `latex` | 变量 | 控制结果输出格式：`latex=0` 线性输出，`latex=1` 输出 LaTeX。 |
| `print(expr1, expr2, ...)` | 多个表达式 | 计算并打印每个表达式（不压栈）。 |
| `status()` | 无参数 | 打印内存使用、GC 次数等统计信息。 |
| `stop()` | 无参数 | 触发运行时错误（测试用）。 |
| `exit()` | 无参数 | 退出进程。 |
| `for(var, start, end, body...)` | `var`: 变量；`start,end`: 整数；`body`: 语句块 | 循环执行（可嵌套）。 |
| `do(body...)` | `body`: 语句序列 | 顺序执行并返回最后一个结果。 |
| `loop(body...)` | `body`: 语句序列 | 无限循环（用 `break` 退出）。 |
| `break()` | 无参数 | 跳出 `loop` 或 `for` 循环。 |
| `test(c1, e1, c2, e2, ..., default)` | 条件-表达式对 | 类 C 的 `?:` 条件链。 |
| `testeq(a, b)` | $a, b$: 表达式 | 判断 $a = b$（返回 1/0）。 |
| `testge`, `testgt`, `testle`, `testlt` | $a, b$: 数值 | 数值比较（返回 1/0）。 |
| `nil` | 常量 | 空值 / 无结果。 |
| `inf` | 常量 | 正无穷大。 |
| `pi` | 常量 | 圆周率。 |
| `i` | 常量 | 虚数单位 $\sqrt{-1}$。 |
| `e` | 常量 | 自然对数的底；`2e` 表示 $2e$，而 `1e-3` 表示科学计数法 $0.001$。 |

---

## 命令参考

以下按功能分类列出所有常用命令，并附详细说明和示例。

### 基础语法

- **赋值**：`变量 = 表达式`  
  例如：`a = 3*x + 1`
- **等式**：`表达式1 = 表达式2`（用于方程求解）
- **向量**：`[e1, e2, ..., en]`  
  例如：`[1, 2, 3]`
- **矩阵**：`[[行1],[行2],...]`  
  例如：`[[1,2],[3,4]]`
- **范围字面量**：`[a..b]` 生成等差数列（步长 $\pm 1$），如 `[1..5]` $\to [1,2,3,4,5]$
- **切片索引**：`A[i:j]`（含端点）、`A[i:]`、`A[:j]`，多维 `A[i:j,k:l]`，整数索引混合使用自动降维
- **块拼接**：`vstack(A,B,...)` 沿行拼接、`hstack(A,B,...)` 沿列拼接，`[A,B]` 嵌套即沿第一维拼接

#### 序列工具：`range` / `seq` / `diffseq` / `cumsum` / `recurrence`
- `range(a, b, step)`：等差数列（含端点），`step` 默认 $\pm 1$；
- `seq(expr, var, start, end)`：表达式在变量整数取值下的序列；
- `diffseq(v)`：一阶前向差分 $v_{i+1}-v_i$；
- `cumsum(v)`：累积和 $[v_1, v_1+v_2, \ldots]$；
- `recurrence(init, f, n)`：递推生成，`init` 为初始向量，$f$ 为二元用户函数（如 `f(a,b)=a+b`），生成 $n$ 项。

```text
? range(1, 10, 3)
> [1,4,7,10]
? seq(x^2, x, 1, 4)
> [1,4,9,16]
? diffseq([1,4,9,16])
> [3,5,7]
? cumsum([1,2,3,4])
> [1,3,6,10]
? f(a,b)=a+b
? recurrence([0,1], f, 6)
> [0,1,1,2,3,5]
```
- **常量**：`pi`（圆周率）、`i`（虚数单位）、`inf`（无穷大）
- **注释**：`# 单行注释`（仅脚本文件内有效）
- **整数范围**：精确整数运算限于 $0 \sim 2^{64}-1$；四则运算、乘方、阶乘等结果超出该范围时自动转为浮点数，数论函数（`gcd`、`mod`、`powermod`、`isprime` 等）仍要求精确整数。

---

### 代数运算

#### `solve(equ, var)`
解方程（多项式、三角、指数、超越方程）与不等式，返回根列表或区间。

- **多项式**：返回所有根（实数或复数，最高 8 次，符号系数二次用求根公式）；
- **因式分解**：乘积/公因子自动分解后逐因子求解（如 $\sin x (x-1)=0$）；
- **三角方程**：$\sin$ / $\cos$ / $\tan$ 的线性参数方程返回主值解；
- **指数方程**：$e^u = c$ 给出 $u = \log c$；
- **超越方程**：无法解析时数值求根（单一根）；
- **不等式**（$<$、$\le$、$>$、$\ge$）：返回解区间列表，无穷用 $\infty$，无解返回 `nil`。

```text
? solve(2*x + 4 = 0, x)
> -2
? solve(x^2 - 5*x + 6, x)
> [2,3]
? solve(x^2 + 1, x)
> [-i,i]
? solve(x = cos(x), x)        # 超越方程：数值根
> [0.739085]
? solve(sin(x) = 1/2, x)      # 三角方程主值解
> [1/6*pi, 5/6*pi]
? solve(exp(x) = 2, x)
> [log(2)]
? solve(x^2 < 4, x)           # 不等式区间解
> [-2, 2]
? solve(x^2 > 4, x)
> [[-inf, -2], [2, inf]]
```

#### `factor(expr, var)` 与 `cfactor(expr, var)`
- `factor`：实数域因式分解。
- `cfactor`：复数域因式分解（含复数根）。
- 额外支持多元立方和差公式 $x^3 \pm y^3$ 的自动识别（无需指定变量）。
- 对 $x^n \pm 1$（$n \le 32$）自动给出有理数域不可约分解：$x^n - 1 = \prod_{d|n} \Phi_d(x)$、$x^n + 1 = \prod_{d|2n, d\nmid n} \Phi_d(x)$，因子用分圆多项式表示。
- 符号分解失败时（次数超过 8 或系数含浮点）自动用 `nroots` 数值求根回退：`factor` 给出实数域分解（共轭复根配成二次因子），`cfactor` 给出全部线性复因子。

```text
? factor(x^2 - 2, x)
> (x - sqrt(2))*(x + sqrt(2))
? factor(x^2 + 1, x)
> x^2 + 1
? cfactor(x^2 + 1, x)
> (x - i)*(x + i)
? factor(x^3 + y^3)
> (x + y) (x^2 - x y + y^2)
? factor(x^3 - y^3)
> (x - y) (x^2 + x y + y^2)
? factor(x^9 - 1)
> (x - 1) (x^2 + x + 1) (x^6 + x^3 + 1)
? factor(x^8 - 1)
> (x - 1) (x + 1) (x^2 + 1) (x^4 + 1)
```

#### `collect(expr, var)`
按指定变量的幂次合并同类项。

```text
? collect(x + x^2 + 2*x + 3*x^2, x)
> 4*x^2 + 3*x
? collect((x+1)^2 + x, x)
> x^2 + 3*x + 1
```

#### `combine(expr)`
合并对数和指数项（将 $\log a + \log b$ 合并为 $\log(ab)$，$e^a e^b$ 合并为 $e^{a+b}$）。

```text
? combine(log(x) + log(y))
> log(x*y)
? combine(exp(x)*z*exp(y))
> z*exp(x+y)
```

---

### 三角恒等变换

#### `trigexpand(expr)`
将三角函数的角度和差展开（如 $\sin(x+y) \to \sin x \cos y + \cos x \sin y$）。

```text
? trigexpand(sin(x+y))
> cos(x)*sin(y) + cos(y)*sin(x)
? trigexpand(cos(x-y))
> cos(x)*cos(y) + sin(x)*sin(y)
```

#### `trigreduce(expr)`
将三角函数的幂次降阶（如 $\sin^2 x \to \frac{1-\cos 2x}{2}$）。

```text
? trigreduce(sin(x)^2)
> 1/2 - 1/2*cos(2*x)
? trigreduce(cos(x)^2)
> 1/2 + 1/2*cos(2*x)
```

#### `trigsimplify(expr)`
尝试用基本三角恒等式化简表达式（如 $\sin^2+\cos^2=1$）。

```text
? trigsimplify(sin(x)^2 + cos(x)^2)
> 1
? trigsimplify(sin(x)^4 + 2*sin(x)^2*cos(x)^2 + cos(x)^4)
> 1
```

#### 余割 / 正割 / 余切及其反函数
`csc`、`sec`、`cot` 与 `acsc`、`asec`、`acot` 补齐三角家族，支持精确特殊值、数值求值、求导与积分。

```text
? csc(pi/6)
> 2
? sec(pi/3)
> 2
? cot(pi/4)
> 1
? acsc(2)
> 1/6*pi
? d(sec(x), x)
> sec(x)*tan(x)
? integral(sec(x), x)
> log(sec(x) + tan(x))
```

---

### 数论与组合

#### 最大公约数与最小公倍数
`gcd(a, b, ...)` 与 `lcm(a, b)` 支持 64 位无符号整数范围内的整数。

```text
? gcd(12, 18)
> 6
? lcm(4, 6)
> 12
```

#### 素数
`isprime(n)` 判断素数，`nextprime(n)` 求下一个素数。

```text
? isprime(7)
> 1
? isprime(9)
> 0
? nextprime(13)
> 17
```

#### 组合数
`choose(n,k)` 或 `binomial(n,k)` 二项式系数、`perm(n,k)` 排列数，$n$ 可为符号。

```text
? choose(10, 3)
> 120
? choose(x, 2)
> 1/2*x^2 - 1/2*x
? binomial(5, 2)
> 10
? perm(5, 2)
> 20
```

#### 正交多项式：`laguerre` / `legendre` / `hermite`
`laguerre(x,n,m)` 广义拉盖尔多项式 $L_n^{(m)}(x)$、`legendre(f,n,m,x)` 连带勒让德函数（$f$ 为求值点，$x$ 为变量）、`hermite(x,n)` 厄米多项式 $H_n(x)$。阶数 $n$ 为非负整数时精确展开，符号阶数保留原式。

```text
? laguerre(x,2,0)
> 1/2*x^2 - 2*x + 1
? legendre(0,2,0,x)
> -1/2
? hermite(x,3)
> 8*x^3 - 12*x
? hermite(1,3)          # 数值自变量按代入求值
> -4
```

#### 其他常用数论函数
`fibonacci(n)`、`totient(n)`（欧拉函数）、`divisors(n)`（正因子列表）、`bernoulli(n)`（伯努利数，$B_1=-1/2$）、`harmonic(n)`（调和数）。

```text
? fibonacci(20)
> 6765
? totient(12)
> 4
? divisors(12)
> [1,2,3,4,6,12]
? bernoulli(4)
> -1/30
? harmonic(4)
> 25/12
```

#### 模算术与有限域
- `powermod(a, e, m)`：模幂 $a^e \bmod m$（二进制快速幂）。
- `invmod(a, m)`：模逆元 $a^{-1} \bmod m$（扩展欧几里得，要求 $\gcd(a,m)=1$）。
- `crt(a1, m1, a2, m2)`：中国剩余定理，解 $x \equiv a_1 \pmod{m_1}$、$x \equiv a_2 \pmod{m_2}$ 的最小非负解。
- `polymod(P, M)`：一元多项式长除法取模，即 $P \bmod M$（$\mathbb{Q}[x]$ 上的除法余式，也可用于有限域多项式运算）。

```text
? powermod(2, 10, 1000)
> 24
? powermod(3, 100, 7)
> 4
? invmod(3, 7)
> 5
? crt(2, 3, 3, 5)
> 8
? polymod(x^3 + 1, x^2 + 1)
> -x + 1
? polymod(x^2 + 1, x + 1)
> 2
```

#### Gröbner 基：`groebner(polys, vars)`
对多项式组（中括号向量）按字典序（lex）计算约化 Gröbner 基（Buchberger 算法），可用于求解多元多项式方程组、理想成员判定等。

```text
? groebner([x^2 - y, x^3 - z], [x, y, z])
> [x^2 - y, x y - z, x z - y^2, y^3 - z^2]
? groebner([x^2 + y^2 - 1, x - y], [x, y])
> [x - y, y^2 - 1/2]
```

---

### 微积分

#### `d(expr, var)` 或 `derivative(expr, var)`
求导，支持高阶导数（通过嵌套）。

```text
? d(sin(x), x)
> cos(x)
? d(x^3, x)
> 3*x^2
? d(d(sin(x),x),x)   # 二阶导
> -sin(x)
```

#### `implicitdiff(F, x, y)` 或 `idiff(F, x, y)`
对隐函数 $F(x,y)=0$ 求 $\frac{dy}{dx}$，使用公式 $\frac{dy}{dx} = -\frac{F_x}{F_y}$。

```text
? implicitdiff(x^2 + y^2 - 1, x, y)
> -x/y
? idiff(y^2 - x, x, y)
> 1/(2*y)
```

#### 多元微积分：`grad` / `jacobian` / `hessian` / `div` / `curl`
变量以向量形式给出，如 $(x,y,z)$。

```text
? grad(x^2 + y^2, [x,y])
> [2*x, 2*y]
? jacobian([x^2 + y, x*y], [x,y])
> [[2*x,1],[y,x]]
? hessian(x^2*y, [x,y])
> [[2*y,2*x],[2*x,0]]
? div([x^2, y^2], [x,y])
> 2*x + 2*y
? curl([0,0,x], [x,y,z])
> [0,-1,0]
```

#### `integral(expr, var)`
不定积分，内置换元、分部积分和有理函数分解能力。

```text
? integral(2*x*cos(x^2) + sin(x^2), x)
> sin(x^2) - 1/2*cos(x^2)   # 注意结果可能需简化
? integral(x^3*sin(x), x)
> (-x^3 + 6*x)*cos(x) + (3*x^2 - 6)*sin(x)
? integral((x^3+1)/(x^2+1), x)
> 1/2*x^2 - ln(x^2+1)/2 + atan(x)
? integral((3*x+5)/(x^2+3*x+2), x)
> 2*ln(x+1) + ln(x+2)
```

#### `defint(expr, var, a, b)`
定积分、反常积分，以及变限/含参变量积分。上下限可以是符号表达式；当原函数无法解析求出时，结果保留为未求值的 `defint(expr, var, a, b)` 表达式。

```text
? defint(1/x^2, x, 1, inf)
> 1
? defint(1/sqrt(x), x, 0, 1)
> 2
? defint(exp(-x), x, 0, inf)
> 1
? defint(t, t, 1, x)                # 变上限积分
> x^2/2 - 1/2
? defint(sin(a*t), t, 0, x)         # 含参数 a
> 1/a - cos(a*x)/a
? defint(exp(t^3), t, x, x^2)       # 无法解析求解，保留原式
> defint(exp(t^3), t, x, x^2)
```

对 `defint` 求导使用完整的莱布尼茨公式：

$$
\frac{d}{dx}\int_{A(x)}^{B(x)} f(x,t)\,dt = \int_{A(x)}^{B(x)} \frac{\partial f}{\partial x}\,dt + f(x,B(x))B'(x) - f(x,A(x))A'(x)
$$

```text
? d(defint(exp(t^3), t, x, x^2), x)
> 2*x*exp(x^6) - exp(x^3)
? d(defint(exp(a*t^3), t, 0, 1), a)   # 含参积分对参数求导
> defint(t^3*exp(a*t^3), t, 0, 1)
```

能解析求解的变限/含参积分仍会先求出原函数再代入上下限，结果与莱布尼茨公式一致。

#### `nintegrate(expr, var, a, b)`
数值积分（自适应 Simpson 或 Gauss 型，用于无法求初等原函数的场合）。

```text
? nintegrate(exp(-x^2), x, 0, 1)
> 0.746824
? nintegrate(sin(x), x, 0, pi)
> 2.0
```

---

### 极限、级数和展开

#### `limit(expr, var, point)`
计算极限。优先采用泰勒展开主项，必要时自动应用洛必达法则；无穷极限会进行倒数代换。

```text
? limit(sin(x)/x, x, 0)
> 1
? limit((1-cos(x))/x^2, x, 0)
> 1/2
? limit(log(x)/x, x, inf)
> 0
? limit((1+1/x)^x, x, inf)
> exp(1)
? limit(x^(1/x), x, inf)
> 1
? limit((x^2+3*x)/(2*x^2-1), x, inf)
> 1/2
```

#### `taylor(expr, var, order)`
在 $x=0$ 处进行泰勒展开到指定阶数（内置常用函数的展开表）。

```text
? taylor(exp(x), x, 6)
> 1 + x + 1/2*x^2 + 1/6*x^3 + 1/24*x^4 + 1/120*x^5 + 1/720*x^6
? taylor(sin(x), x, 9)
> x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9
? taylor(log(1+x), x, 6)
> x - 1/2*x^2 + 1/3*x^3 - 1/4*x^4 + 1/5*x^5 - 1/6*x^6
```

#### `laurent(expr, var, point, order)`
在指定点处展开洛朗级数，`order` 为展开阶数。

```text
? laurent(1/(x^2*(1-x)), x, 0, 2)
> x^(-2) + x^(-1) + 1 + x + O(x^2)
? laurent((x+1)/(x-1)^2, x, 1, 2)
> 2/(x-1)^2 + 1/(x-1) + 0 + O([x-1]^2)
```

#### `sum(expr, var, a, b)`
求有限和、带符号上下限的闭式求和，或部分无穷级数（若 $b=\infty$ 则计算收敛级数和）。

符号求和会尝试按以下规则求闭式 $S(n)=\sum_{x=1}^{n} f(x)$，结果为 $S(b)-S(a-1)$：

- 与求和变量无关的项：$f\cdot (b-a+1)$；
- 加法线性拆分：$\sum(f+g)=\sum f+\sum g$；
- 乘法：提取不含求和变量的因子；
- 幂和：$x^p$（$p$ 为非负整数）用 Faulhaber 公式（伯努利数）求任意次幂；
- 等比：$r^{a x+b}$ 用 $r^{a+b}\frac{r^{a n}-1}{r^a-1}$；
- 均不匹配时原样返回 `sum(expr, var, a, b)`。

```text
? sum(1, x, 1, k)
> k
? sum(x, x, 1, k)
> 1/2*k^2 + 1/2*k
? sum(x^2, x, 1, k)
> 1/3*k^3 + 1/2*k^2 + 1/6*k
? sum(x^2, x, 1, k+1)
> 1/3*k^3 + 3/2*k^2 + 13/6*k + 1
? sum(2*x + 1, x, 1, k)
> k^2 + 2*k
? sum(2^x, x, 0, k)
> 2*2^k - 1
? sum(1/x, x, 1, k)
> sum(1/x, x, 1, k)   # 无闭式，保留原式
? sum(1/x^2, x, 1, inf)
> pi^2/6
? sum(1/x, x, 1, inf)
> inf   # 发散
? sum(1/(x*(x+1)), x, 1, inf)
> 1
? sum((1/2)^x, x, 0, inf)
> 2
```

对未求值的 `sum` 表达式求导使用离散莱布尼茨公式：

$$
\frac{d}{dp}\sum_{x=A}^{B} f(x,p) = \sum_{x=A}^{B} \frac{\partial f}{\partial p}
$$

```text
? d(sum(1/x, x, 1, k), k)
> 1/k
? d(sum(1/x^2, x, 1, k), k)
> 1/k^2
```

注意：能求出闭式的和会先求值，此时求导得到的是闭式关于变量的导数，例如 `d(sum(x^2,x,1,k),k) = k^2 + k + 1/6`。求和签名优先采用 `sum(F, X, A, B)`；旧的 `sum(X, A, B, F)` 写法在第二个参数不是符号时仍然兼容。

---

### 线性代数

#### `rref(A)`
将矩阵 $A$ 化为行最简阶梯形（RREF）。

```text
? A = [[1,2,1],[2,4,2],[1,1,0]]
? rref(A)
> [[1,0,0],[0,1,0],[0,0,0]]
```

#### `rank(A)` 和 `nullspace(A)`
- `rank`：矩阵的秩。
- `nullspace`：零空间的一组基（返回矩阵，每列是一个基向量）。

```text
? rank([[1,2,3],[2,4,6]])
> 1
? nullspace([[1,2,3],[2,4,6]])
> [[-2,1,0],[-3,0,1]]
```

#### `det(A)` 和 `inv(A)`
- `det`：行列式（通过行消元计算）。
- `inv`：矩阵的逆。数值矩阵先做条件数检查（$\infty$-范数估计，阈值 $10^6$）：奇异或病态时返回 `nil` 并打印警告，不再产生数值爆炸的“大数”结果；精确矩阵奇异时同样返回 `nil` 并警告。

```text
? det([[1,2],[3,4]])
> -2
? inv([[1,2],[3,4]])
> [[-2,1],[3/2,-1/2]]
? inv([[1.0,2],[2,4]])   # 奇异：返回 nil + 警告
> nil
```

#### `lu(A)` 和 `qr(A)`
- `lu`：返回 $(L, U, P)$ 三元组（部分选主元）。
- `qr`：返回 $(Q, R)$ 正交分解（Gram-Schmidt 过程）。

```text
? lu([[2,1],[4,3]])
> [[[1,0],[2,1]], [[2,1],[0,1]], [[1,0],[0,1]]]
? qr([[2,0],[0,3]])
> [[[1,0],[0,1]], [[2,0],[0,3]]]
```

#### 特征分解：`eigenvals` / `eigenvecs` / `eigen`
一般实矩阵（非对称/非 Hermitian）的数值特征分解（Hessenberg + Wilkinson 移位 QR 算法）。输入需为数值矩阵，特征值按实部降序排列，特征向量列与之一一对应；复特征值（及复特征向量）用内置虚数单位 $i$ 表示。

- `eigenvals(A)`：特征值向量；
- `eigenvecs(A)`：特征向量矩阵（每列一个单位特征向量）；
- `eigen(A)`：返回 $[\text{特征值向量}, \text{特征向量矩阵}]$。

```text
? eigenvals([[1,2],[3,4]])
> [5.37228,-0.372281]
? eigenvals([[0,-1],[1,0]])
> [i,-i]
? eigenvecs([[1,2],[2,1]])
> [[0.707107,0.707107],[0.707107,-0.707107]]
? eigen([[2,0],[0,3]])
> [[3,2],[[0,1],[1,0]]]
```

#### 若尔当分解：`jordan(A)`
数值若尔当分解（支持 $\le 5\times 5$ 数值方阵），返回 $[P, J]$ 满足 $A = P J P^{-1}$。$J$ 为若尔当标准形：对角块 $\lambda$，若有广义特征向量则对应超对角元素为 1（特征值按实部降序分组，复特征值以 $i$ 表示）。奇异/病态输入返回 `nil` 并警告。

```text
? jordan([[2,1],[0,2]])
> [[[1,0],[0,1]],[[2,1],[0,2]]]
? jordan([[0,-1],[1,0]])
> [[[i,-i],[1,1]],[[i,0],[0,-i]]]
```

#### 奇异值分解：`svd(A)`
数值 SVD（one-sided Jacobi 算法），返回 $[U, S, V]$，满足 $A \approx U\,S\,V^{T}$（经济分解：$U$ 为 $m\times k$、$S$ 为 $k\times k$ 对角、$V$ 为 $n\times k$，$k=\min(m,n)$），奇异值降序排列。可用于稳健伪逆、低秩近似、降维等。

```text
? svd([[0,2],[3,0]])
> [[[0,1],[1,0]], [[3,0],[0,2]], [[1,0],[0,1]]]
? r = svd([[1,2],[3,4],[5,6]])
? norm(dot(dot(r[1], r[2]), transpose(r[3])) - [[1,2],[3,4],[5,6]]) < 0.000000001
> 1
```

#### 线性方程组：`linsolve(A, b)`
求解 $Ax=b$。右端 $b$ 可为向量（返回解向量）或矩阵（逐列求解，返回矩阵）。
- 矩阵含小数（double）时走数值分支：基于 SVD 的最小二乘/最小范数解，对奇异或病态矩阵同样稳健；
- 全整数/有理数/符号矩阵走符号分支：增广矩阵行化简（RREF）精确求解；无解报 `inconsistent`，解不唯一报 `singular`。

```text
? linsolve([[2,1],[1,3]],[5,6])
> [9/5,7/5]
? linsolve([[2.0,1],[1,3]],[5,6])
> [1.8,1.4]
? linsolve([[1.0,1],[1,1]],[2,3])    # 奇异：最小范数最小二乘
> [1.25,1.25]
? linsolve([[1.0,1],[1,2],[1,3]],[1,2,2])   # 超定：最小二乘
> [0.666667,0.5]
```

#### 行/列基本操作（初等变换）
高斯消元/教学演示用的原子命令，索引为 1-based（与 `A[1,2]` 一致），均返回新矩阵（不改动原矩阵），支持符号元素与符号系数。

- `rowswap(A, i, j)` / `colswap(A, i, j)`：交换第 $i$、$j$ 行/列；
- `rowadd(A, i, j, k)` / `coladd(A, i, j, k)`：第 $i$ 行/列 $\leftarrow$ 第 $i$ 行/列 $+ k \times$ 第 $j$ 行/列；
- `rowmul(A, i, k)` / `colmul(A, i, k)`：第 $i$ 行/列 $\times k$。

```text
? rowswap([[1,2,3],[4,5,6],[7,8,9]], 1, 3)
> [[7,8,9],[4,5,6],[1,2,3]]
? rowadd([[1,2],[3,4]], 2, 1, 2)
> [[1,2],[5,8]]
? rowmul([[1,2],[3,4]], 1, 3)
> [[3,6],[3,4]]
? coladd([[1,2],[3,4]], 2, 1, a)     # 符号系数
> [[1,a + 2],[3,3 a + 4]]
```

分步高斯消元示例：

```text
? A = [[1,2,1],[2,4,2],[1,1,0]]
? A = rowadd(A, 2, 1, -2)
? A = rowadd(A, 3, 1, -1)
? A = rowswap(A, 2, 3)
? A
> [[1,2,1],[0,-1,-1],[0,0,0]]
```

#### 克罗内克积：`kron(A, B)` / `kronecker(A, B)`
克罗内克积（张量积）：$A\,(m\times n)$ 与 $B\,(p\times q)$ 的积为 $(mp)\times(nq)$ 块矩阵 $[a_{ij}B]$。支持矩阵与向量（1D 按列向量展开）、符号元素，可链式传入多个因子（`kron(A, B, C, ...)`）。常用于构造大规模稀疏矩阵及矩阵方程（Sylvester / Lyapunov 方程）的向量化。

```text
? kron([[1,2],[3,4]], [[0,5],[6,7]])
> [[0,5,0,10],[6,7,12,14],[0,15,0,20],[18,21,24,28]]
? kron([1,2], [3,4])
> [3,4,6,8]
? kron([[1]], [[2]], [[3]])
> [6]
```

构造 Sylvester 算符 $A\otimes I + I\otimes A$ 的示例：

```text
? A = [[1,2],[3,4]]
? I2 = [[1,0],[0,1]]
? kron(A, I2) + kron(I2, A)
> [[2,2,2,0],[3,5,0,2],[3,0,5,2],[0,3,3,8]]
```

#### 其他线性代数工具
- `cross(u, v)`：三维向量叉积；
- `dot(u, v)` 或 `inner(u, v)`：显式向量点积 $\sum_i u_i v_i$（对矩阵等价矩阵乘法）；注意 `u*v` 对张量是逐元素（Hadamard）乘法，点积/矩阵乘法请用 `dot`；
- `diag(a, b, ...)` 或 `diag((a,b,...))`：构造对角矩阵；
- `charpoly(A)` 或 `charpoly(A, x)`：特征多项式 $\det(A - xI)$；
- `pseudoinverse(A)`：Moore-Penrose 伪逆（满秩情形）；
- `norm(A)` / `norm(A, p)`：矩阵范数。默认 Frobenius 范数 $\|A\|_F$；`norm(A, 1)` 为最大列绝对值和，`norm(A, inf)` 为最大行绝对值和，`norm(A, 2)` 为谱范数（最大奇异值）。向量用法不变：默认欧几里得，$1$ 为绝对值之和，$\infty$ 为最大绝对值。
- `cond(A)` / `cond(A, p)`：矩阵条件数。默认谱条件数 $\kappa_2(A)=\sigma_{\max}/\sigma_{\min}$（SVD 计算）；`cond(A, 1)` / `cond(A, inf)` 为 $\|A\|_p\|A^{-1}\|_p$。奇异矩阵返回 $\infty$。

```text
? cross([1,0,0], [0,1,0])
> [0,0,1]
? dot([1,2,3], [4,5,6])
> 32
? diag(1,2,3)
> [[1,0,0],[0,2,0],[0,0,3]]
? charpoly([[1,2],[3,4]])
> x^2 - 5*x - 2
? pseudoinverse([[1,0],[0,1],[0,0]])
> [[1,0,0],[0,1,0]]
? norm([3,4])
> 5
? norm([[1,2],[3,4]],1)
> 6
? norm([[3,0],[0,4]],2)
> 4
? cond([[2,0],[0,3]])
> 1.5
? cond([[1,2],[2,4]])
> inf
```

#### 初等函数补全
`cbrt`（立方根）、`hypot`（$\sqrt{x^2+y^2}$）、`round`（四舍五入）、`trunc`（向零取整）、`min`/`max`（极值）。

```text
? cbrt(-8)
> -2
? hypot(3,4)
> 5
? round(5/2)
> 3
? trunc(-7/2)
> -3
? max(3,7,2)
> 7
```

#### 张量分析：`tensorrank` / `christoffel` / `covdiff`
- `tensorrank(T)`：张量阶数（标量 0、向量 1、矩阵 2、三维张量 3）；
- `contract(A)` 或 `contract(A, i, j, ...)`：指标缩并（迹等）；
- `christoffel(g, x)`：由度规 $g$（对称矩阵）与坐标向量 $x$ 计算克里斯托费尔符号 $\Gamma^i_{jk}$（三维张量）；
- `covdiff(T, g, x)`：协变导数 $\nabla T$（标量→梯度向量，向量→矩阵，二阶逆变张量→三维张量），约定 $T$ 为逆变分量。

```text
? tensorrank([[1,2],[3,4]])
> 2
? contract([[1,2],[3,4]])
> 5
? christoffel([[1,0],[0,r^2]], [r,t])[1,2,2]    # 二维极坐标
> -r
? christoffel([[1,0],[0,r^2]], [r,t])[2,1,2]
> 1/r
? covdiff(r^2*cos(t), [[1,0],[0,r^2]], [r,t])   # 标量 → 梯度
> [2*r*cos(t), -r^2*sin(t)]
? covdiff([r,0], [[1,0],[0,r^2]], [r,t])[2,2]   # 向量协变导数
> 1
```

---

### 常微分方程

#### `dsolve(ode, var)`
符号求解一阶或二阶常微分方程，支持：
- 可分离变量、一阶线性、**伯努利方程**（$y'+P(x)y=Q(x)y^n$）、二阶常系数齐次及部分非齐次、**欧拉-柯西方程**（$x^2y''+axy'+by=0$，含实根/重根/复根）；
- **初值条件**：`dsolve(ode, x, y(a)=b, eval(d(y(x),x),x,a)=b)` 自动代入通解解出常数，得到特解。

```text
? dsolve(d(y(x),x) = x*y(x), x)
> y(x) = c*exp(1/2*x^2)
? dsolve(d(y(x),x) + y(x) = 0, x)
> y(x) = c*exp(-x)
? dsolve(d(y(x),x) + y(x) = x, x)
> y(x) = x - 1 + c*exp(-x)
? dsolve(d(y(x),x,2) - y(x) = 0, x)
> y(x) = c1*exp(x) + c2*exp(-x)
? dsolve(d(y(x),x) = x*y(x)^2, x)
> y(x) = -2/(x^2 + c)
? dsolve(d(y(x),x,2) + y(x) = sin(x), x)
> y(x) = c1*cos(x) + c2*sin(x) - 1/2*x*cos(x)
? dsolve(d(y(x),x) = y(x), x, y(0)=1)          # IVP 特解
> exp(x)
? dsolve(d(y(x),x,2) + y(x) = 0, x, y(0)=1, eval(d(y(x),x),x,0)=0)
> cos(x)
? dsolve(d(y(x),x) + y(x)/x = x*y(x)^2, x)     # 伯努利
> 1/(C1*x - x^2)
? dsolve(x^2*d(y(x),x,2) + x*d(y(x),x) + y(x) = 0, x)   # 欧拉-柯西
> C1*cos(log(x)) + C2*sin(log(x))
```

#### `odesolve(f, x, x0, x1, y0, n)`
数值求解 ODE，使用经典 RK4 方法。一阶标量返回 $n+1$ 个点的函数值列表；**传入函数向量**（如 `[y2, -y1]` 表示 $y''=-y$）则求解高阶/方程组，返回 $(n+1)\times(m+1)$ 矩阵（每行 $x, y_1, \dots, y_m$）。

```text
? odesolve(-y, x, 0, 1, 1, 10)   # 解 dy/dx = -y, y(0)=1, 步数10
> [[0,1], [0.1,0.9048], ...]
? odesolve([y2, -y1], x, 0, 1, [1,0], 5)   # y'' = -y
> ((0,1,0),(0.2,0.980067,-0.198667),...)
```

#### 自适应步长：`odesolve_adaptive`
`odesolve_adaptive(f, x, x0, x1, y0, tol)` 使用 **RK45（Dormand-Prince）** 自动调整步长保证局部误差不超过 `tol`，返回自适应插值点矩阵（每行 $x, y_1, \dots, y_m$），支持标量与向量系统，限制最大步数（10000）。

```text
? odesolve_adaptive(-y, x, 0, 1, 1, 0.000001)
> ((0,1),(0.1,0.904837),(0.3,0.740818),...,(1,0.367879))
```

#### 方向场：`slopefield`
`slopefield(f, x, a, b, nx, y, c, d, ny)` 在矩形网格上生成 $\frac{dy}{dx}=f(x,y)$ 的方向场数据，返回每点 $(x, y, dx, dy)$（长度 1 归一化），配合 `plotdata` 可绘制斜率场。

```text
? slopefield(y, x, 0, 1, 1, y, 0, 1, 1)
> ((0,0,1,0),(1,0,1,0),(0,1,0.707107,0.707107),(1,1,0.707107,0.707107))
```

#### 符号线性系统与矩阵指数
- `dsolve_system(d([y1,y2],x) = A*[y1,y2], x)`：常系数线性系统 $Y'=AY$（含常数非齐次 $AY+B$）的符号通解，支持实/复/重特征值；
- `expm(A)` / `expm(A, x)`：矩阵指数 $e^A$ / $e^{Ax}$，任意阶方阵。2×2 符号矩阵给出精确符号结果（含 $x$ 参数时保持 $e^{Ax}$ 形式）；数值矩阵（含小数或高阶）走特征分解数值路径，支持复特征值。
- `logm(A)` / `sqrtm(A)`：矩阵对数与矩阵平方根（主分支），任意阶数值方阵；谱含非正实部时给出复结果（如负特征值的对数/平方根含虚数单位 $i$）。

```text
? dsolve_system(d([y1,y2],x) = [[0,1],[-1,0]]*[y1,y2], x)
> [C1*cos(x) - C2*sin(x), -C1*sin(x) - C2*cos(x)]
? expm([[0,1],[-1,0]], x)
> [[cos(x),sin(x)],[-sin(x),cos(x)]]
? expm([[1,2],[3,4]])        # 任意阶数值矩阵
> [[51.969,74.7366],[112.105,164.074]]
? sqrtm([[4,0],[0,9]])
> [[2,0],[0,3]]
? logm([[2.0,1],[1,2]])      # 与 expm 互逆
> [[0.549306,0.549306],[0.549306,0.549306]]
```

---

### 数值求根和数据生成

#### `findroot(f, var, guess)`
利用割线法求方程 $f(var)=0$ 在 `guess` 附近的根。

```text
? findroot(x^2 - 2, x, 1)
> 1.41421356
? findroot(sin(x) - 0.5, x, 0.5)
> 0.52359878
```

#### `nroots(poly, var)`：多项式数值求根
返回多项式 $poly(var)$ 的全部数值复根向量（含重根，个数等于次数）。`factor` / `cfactor` 在符号分解失败时自动调用它做数值回退：实数域把共轭复根配对成实二次因子，复数域输出全部线性因子。

```text
? nroots(x^3 - 1, x)
> [1,-0.5 - 0.866025 i,-0.5 + 0.866025 i]
? factor(x^16 - 2*x^8 - 1)   # 符号无法分解 → 数值回退
> (x + 1.11647) (x - 1.11647) (x^2 + 1.2465) ...
```

#### `fsolve(F, X, G)`：非线性方程组（多维牛顿法）
求方程组 $F(X)=0$ 的数值解。$F$ 为函数向量（或单个函数）、$X$ 为变量向量、$G$ 为初始猜测向量；每一步用中心差分构造雅可比矩阵，行化简求解线性步长，$\lVert F\rVert < 10^{-12}$ 收敛。

```text
? fsolve([x^2 - 2, y^2 - 3], [x, y], [1, 1])
> [1.41421, 1.73205]
? fsolve([x^2 + y^2 - 1, x - y], [x, y], [0.7, 0.7])
> [0.707107, 0.707107]
? fsolve(x^2 - 2, x, 1)          # 单方程（标量 F/X/G）
> [1.41421]
```

#### `implist(F, x, a, b, n, y_guess)`
对隐函数 $F(x,y)=0$ 生成离散数据点。对 $x$ 在 $[a,b]$ 区间内均匀取 $n$ 个点，在每个点上用牛顿法求解 $y$，初值沿用上一点解得的 $y$（首个点用 `y_guess`）。返回点列表。

```text
? implist(x^2 + y^2 - 1, x, -1, 1, 20, 0.5)
> [[ -1, 0], [-0.8947, 0.4472], ...]
```

#### `plotdata(f, x, a, b, n)`
生成函数 $y=f(x)$ 在 $[a,b]$ 区间上的采样点矩阵（共 $n+1$ 个点），用于外部绘图。

```text
? plotdata(sin(x), x, 0, 2*pi, 30)
> [[0,0], [0.2094,0.2079], ...]
? plotdata(x^2, x, 0, 2, 3)
> [[0,0], [0.6667,0.4444], [1.333,1.7778], [2,4]]
```

#### `surfdata(F, x, a, b, nx, y, c, d, ny)`
生成二元函数 $F(x,y)$ 在矩形域 $[a,b]\times[c,d]$ 上的网格采样矩阵（大小为 $(ny+1)\times(nx+1)$），可用于三维绘图。

```text
? surfdata(sin(x)*cos(y), x, -3, 3, 20, y, -3, 3, 20)
> [[0,0.14,...], ...]
? surfdata(x+y, x, 0, 1, 1, y, 0, 1, 1)
> [[0,1],[1,2]]
```

---

### 工程数学变换

#### `laplace(f, t, s)`
计算函数 $f(t)$ 的拉普拉斯变换 $F(s)=\int_0^\infty f(t)e^{-st}\,dt$。

```text
? laplace(sin(t), t, s)
> 1/(s^2 + 1)
? laplace(cos(t), t, s)
> s/(s^2 + 1)
? laplace(t^2, t, s)
> 2/s^3
? laplace(exp(2*t), t, s)
> 1/(s - 2)
? laplace(sinh(t), t, s)
> 1/(s^2 - 1)
? laplace(cosh(t), t, s)
> s/(s^2 - 1)
```

#### `fourier(f, t, w)`
计算傅里叶变换 $F(w)=\int_{-\infty}^{\infty} f(t)e^{-i w t}\,dt$。

```text
? fourier(exp(-abs(t)), t, w)
> 2/(w^2 + 1)
```

#### `ztransform(f, n, z)`
计算 Z 变换 $F(z)=\sum_{n=0}^{\infty} f(n) z^{-n}$。

```text
? ztransform(1, n, z)
> z/(z - 1)
? ztransform((1/2)^n, n, z)
> z/(z - 1/2)
? ztransform(n, n, z)
> z/(z - 1)^2
```

#### 逆变换：`ilaplace` / `ifourier` / `iztransform`
与三个正向变换对应的逆变换，支持常见变换对（常数/冲激、幂、指数、三角、双曲、带衰减的正弦/余弦、平移与线性组合），无法匹配时保留原式。

```text
? ilaplace(1/(s^2 + 1), s, t)
> sin(t)
? ilaplace(s/(s^2 + 4), s, t)
> cos(2*t)
? ilaplace(1/((s-1)^2 + 4), s, t)
> 1/2*exp(t)*sin(2*t)
? ilaplace((s-1)/((s-1)^2 + 4), s, t)
> exp(t)*cos(2*t)
? ifourier(2/(w^2 + 1), w, t)
> exp(-abs(t))
? iztransform(z/(z - 1/2), z, n)
> 2^(-n)
```

#### 信号与工程函数：`sinc` / `heaviside` / `dirac` / `deg` / `rad` / `convolve`
- `sinc(x)`：$\operatorname{sinc}(x)=\frac{\sin x}{x}$，$\operatorname{sinc}(0)=1$；
- `heaviside(x)`：阶跃函数（$x<0$ 为 0，$x>0$ 为 1，$x=0$ 为 $1/2$）；
- `dirac(x)`：狄拉克 $\delta$ 函数（$x \ne 0$ 时为 0），求导 $\frac{d}{dx}\operatorname{heaviside}(x)=\operatorname{dirac}(x)$；
- `deg`、`rad`：角度单位常量（$\mathrm{deg}=\pi/180$），如 $\sin(90^\circ)=1$；
- `convolve(u, v)`：两个一维向量的离散卷积。

```text
? sinc(0)
> 1
? d(heaviside(x), x)
> dirac(x)
? sin(90*deg)
> 1
? convolve([1,2,3], [4,5,6])
> [4,13,28,27,18]
```

---

### 复数运算

复数以 $a + b i$ 形式精确表示，四则运算（$+$、$-$、$\times$、$/$）、整数幂与开方自动按复数规则计算；$\exp$、$\log$、$\sin$、$\cos$、$\tan$、$\sinh$、$\cosh$、$\tanh$ 及六个反函数支持精确复数参数（欧拉公式/对数公式展开），一般复数幂 $z^w$ 按主值 $\exp(w\log z)$ 求值（含负底、虚底，如 $(-2)^i$、$i^i$），正实数底另有 $a^{c+di}=a^c(\cos(d\ln a)+i\sin(d\ln a))$ 的简化路径。全程不依赖 `complex.h`（纯符号与数值实现）。$e$ 与 $i$ 均为内置常量，可被用户赋值覆盖，`clear()` 后恢复。

```text
? e^i
> cos(1) + i*sin(1)
? e^(i*pi)
> -1
? sin(i)
> i*sinh(1)
? cos(1+i)
> cos(1)*cosh(1) - i*sin(1)*sinh(1)
? log(1+i)
> 1/2*log(2) + 1/4*i*pi
? log(-1)
> i*pi
? 2^i
> cos(log(2)) + i*sin(log(2))
? i^i
> exp(-1/2*pi)
? arcsin(2)
> 1/2*pi - i*log(3^(1/2) + 2)
? arctan(2*i)
> 1/2*pi + 1/2*i*log(3)
? arccosh(1/2)
> i*arctan(1/2*3^(1/2), 1/2)
? (1+2*i)*(3+4*i)
> -5 + 10*i
```

#### 复数分析：`cplotdata` / `residue` / `poles` / `zeros`
- `cplotdata(f, z, a, b, c, d, nx, ny)`：在矩形区域 $x\in[a,b]$、$y\in[c,d]$ 上按 $nx\times ny$ 网格采样复函数 $f(x+yi)$，返回 $(nx+1)(ny+1)\times 4$ 数据矩阵，每行 $(x, y, \Re f, \Im f)$；
- `residue(f, z, z0)`：计算 $f$ 在 $z_0$ 的留数（自动判定极点阶数，本质奇点保留原式）；
- `poles(f, z)` / `zeros(f, z)`：返回有理函数 $f$ 的极点/零点列表（含复数根，无则返回 `nil`）。

```text
? cplotdata(z, z, 0, 1, 0, 1, 1, 1)
> [[0,0,0,0],[1,0,1,0],[0,1,0,1],[1,1,1,1]]
? residue(1/(z^2+1), z, i)
> -1/2*i
? residue(sin(z)/z^2, z, 0)
> 1
? residue(exp(1/z), z, 0)      # 本质奇点，保留原式
> residue(exp(1/z), z, 0)
? poles(1/(z^2+1), z)
> [-i, i]
? zeros((z^2+1)/(z-3), z)
> [-i, i]
```

---

### 特殊函数

#### `hypergeom(a, b, c, x)`
高斯超几何函数 ${}_2F_1(a,b;c;x)$，支持符号特殊值（$x=0$、$b=c$ 或 $a=c$ 时的 $(1-x)^{-a}$、$a=0$ 或 $b=0$）、负整数参数的有限项多项式精确展开、以及数值级数求和。

```text
? hypergeom(1,1,1,x)
> 1/(1 - x)
? hypergeom(-2,3,4,x)
> 3/5*x^2 - 3/2*x + 1
? d(hypergeom(a,b,c,x), x)
> a*b*hypergeom(a+1, b+1, c+1, x)/c
```

#### 高级特殊函数：`zeta` / `besselj` / `bessely` / `ellipk` / `ellipe`
- `zeta(s)`：黎曼 $\zeta$ 函数。整数参数给出精确值（$\zeta(2)=\pi^2/6$、$\zeta(0)=-1/2$、$\zeta(-1)=-1/12$、$\zeta(-2n)=0$），一般实参数用 Euler-Maclaurin 求和 + 反射公式数值计算；
- `besselj(n, x)` / `bessely(n, x)`：整数阶第一/二类贝塞尔函数（级数 + 递推，支持负阶）；
- `ellipk(k)` / `ellipe(k)`：第一/二类完全椭圆积分（AGM 算法），$k=0$ 给出 $\pi/2$。

```text
? zeta(2)
> 1/6*pi^2
? zeta(3)
> 1.20206
? besselj(0, 1)
> 0.765198
? bessely(1, 1)
> -0.781213
? ellipk(1/2)
> 1.68575
```

---

### 统计

#### 描述统计：`mean` / `variance` / `stddev`
对一维向量计算均值、总体方差（除以 $n$）与标准差。

```text
? mean([1,2,3,4])
> 5/2
? variance([1,2,3,4])
> 5/4
? stddev([1,2,3,4])
> 5^(1/2)/2
```

#### 正态分布：`normal` / `normalcdf` / `normalinv`
- `normal(x)` 或 `normal(x, mu, sigma)`：概率密度 $\varphi(x)=\frac{1}{\sigma\sqrt{2\pi}}e^{-(x-\mu)^2/(2\sigma^2)}$；
- `normalcdf(x)` 或 `normalcdf(x, mu, sigma)`：累积分布（用 $\operatorname{erf}$ 表示）；
- `normalinv(p)` 或 `normalinv(p, mu, sigma)`：分位数（数值求解）。

```text
? normal(0)
> 1/(2^(1/2)*pi^(1/2))
? normalcdf(1.0)
> 0.841345
? normalinv(0.9)
> 1.28155
```

---

### 假设系统

#### `assume(condition)` 和 `clear()`
- `assume`：对变量施加类型/符号假设，支持：
  - `assume(x > 0)` / `assume(x < 0)`：正/负假设，影响 $\lvert x\rvert$（$\to x$ 或 $-x$）；
  - `assume(n, integer)`：整数假设，影响 $\sin(n\pi)$（$\to 0$）、$\cos(n\pi)$（$\to (-1)^n$）、$\tan(n\pi)$（$\to 0$）、`floor`/`ceiling`/`round`（$\to$ 自身）；
  - `assume(y, real)`：实数假设，影响 $\operatorname{conj}(y)$（$\to y$）、$\operatorname{imag}(y)$（$\to 0$）；
  - `assume(z, nonzero)`：非零假设（记录备用）。
- `clear()`：清除所有假设。

```text
? abs(x)
> abs(x)
? assume(x > 0)
> nil
? abs(x)
> x
? assume(n, integer)
> nil
? sin(n*pi)
> 0
? cos(n*pi)
> [-1]^n
? floor(n)
> n
? clear()
> nil
? abs(x)
> abs(x)
```

---

## 绘图与可视化

Eigenmath 本身不直接绘制图形，但提供了丰富的数据生成命令（`plotdata`、`surfdata`、`implist`），您可以将这些数据导出（通过脚本或重定向）并配合 Gnuplot、Matplotlib 等工具绘图。

例如，生成 $\sin(x)$ 数据并保存为 CSV：
```sh
echo "plotdata(sin(x),x,0,2*pi,100)" > sin.in
./eigenmath sin.in > sin.csv
```

---

## 测试与验证

项目包含回归测试套件，位于 `tests/` 目录下，通过 shell 脚本 `run_*.sh` 执行。

```sh
make test        # 运行全部测试
```

您也可以单独运行某个测试脚本：
```sh
./tests/run_common_funcs.sh
```

测试脚本会读取 `*_cases.txt` 中的 `表达式 => 期望输出` 用例，比对实际输出并报告差异。

---

## 贡献与开发

欢迎贡献代码、报告问题或提出改进建议。请确保：

- 遵循现有代码风格（缩进、命名等）。
- 为新功能添加测试用例。
- 更新此文档以反映变更。

### 项目结构

```
├── main.c              # 交互/脚本入口
├── eigenmath_runtime.c # 核心运行时（内存/栈/求值驱动/符号表）
├── eigenmath_scan.c    # 词法扫描与解析
├── eigenmath_arith.c   # 基本算术运算
├── eigenmath_trig.c    # 三角函数
├── eigenmath_algebra.c # 代数化简/因式分解/假设
├── eigenmath_calculus.c# 微积分/微分方程/积分变换
├── eigenmath_matrix.c  # 矩阵与线性代数
├── eigenmath_tensor.c  # 张量运算
├── eigenmath_numeric.c # 数值计算
├── eigenmath_special.c # 特殊函数
├── eigenmath_output.c  # 输出格式化
├── eigenmath_solve.c   # 方程求解/多项式
├── eigenmath_sequence.c# 求和/序列
├── eigenmath_mem_pool.c# 固定容量内存池分配器
├── eigenmath_data.c    # 数学常量与内置函数表
├── eigenmath_core.h    # 核心定义（Atom_T / EigenMath_Global_T / 符号表索引）
├── eigenmath_config.h  # 可配置参数（栈/符号表/内存池大小等）
├── eigenmath_mem_pool.h# 内存池 API
├── eigenmath_*.h       # 各功能模块头（runtime/scan/arith/.../sequence）
├── eigenmath.h         # 万能头文件（包含全部核心与模块声明）
├── tests/              # 回归测试脚本及预期输出
├── Makefile            # 构建控制
└── README.md           # 本文档
```

### 开发调试

使用 `make debug`（若支持）可生成带有调试信息的可执行文件，便于 `gdb` 调试。

---

## 更新记录

- **矩阵若尔当分解**：新增 `jordan(A)`，返回 $[P, J]$ 满足 $A=PJP^{-1}$，支持数值方阵（$\le 5\times 5$），复用特征值 QR 算法并自动构造广义特征向量（固定数组实现，无动态内存）。
- **多项式数值求根与因式分解回退**：`nroots(poly, var)` 返回全部数值复根；`factor`/`cfactor` 在符号分解失败（次数超过 8 或系数含浮点）时自动回退为数值分解，实数域配成实二次因子、复数域输出全部线性因子。
- **分圆多项式分解**：`factor` 对 $x^n \pm 1$（$n \le 32$）直接给出有理数域不可约分解，因子用分圆多项式表示。
- **矩阵求逆条件数检查**：`inv` 对数值矩阵先做 $\infty$-范数条件数估计（阈值 $10^6$），奇异或病态时返回 `nil` 并打印警告，不再输出数值爆炸的“大数”；精确矩阵奇异时同样返回 `nil` 并警告。
- **矩阵范数与条件数**：`norm` 支持矩阵 1-范数、$\infty$-范数与谱范数（最大奇异值）；新增 `cond` 谱条件数及 $1$/$\infty$ 范数条件数，奇异矩阵返回 $\infty$。
- **整数模型调整**：`USE_BIG_INT=0`（默认）时，有理数分子/分母为 64 位无符号整数，四则运算、乘方、阶乘等结果超出 $0 \sim 2^{64}-1$ 时自动转为双精度浮点继续计算；`USE_BIG_INT=1` 时，超出 64 位的整数自动升级为任意精度大整数，分数分子/分母自动升级为大有理数（`BIGRATIONAL`），加减乘除、乘方、比较、取整等保持精确，乘法内核采用 NTT 加速。
- **科学计数法**：`1e-3` 表示科学计数法（$0.001$），`2e` 表示 $2e$（自然常数），大写 $E$ 同样可作为指数标记。
- **输出控制**：移除终端 Unicode 绘图（Box Drawing 二维排版），结果统一为纯文本线性输出；新增 `latex` 变量，`latex=1` 时结果以 LaTeX 代码输出。

---

## 作者

- **最爱小兔糖糖（Rabbit Candy）**：项目作者与维护者，负责本项目全部设计与实现。

---

## 许可证

本项目源代码采用 BSD 2-Clause License 授权（见各源文件头部版权声明），允许自由使用、修改和分发。
版权归 **最爱小兔糖糖（Rabbit Candy）** 所有，保留所有权利。

---