fibodevy/BigInts
GitHub: fibodevy/BigInts
为 Pascal 语言提供无依赖的任意精度整数和小数运算库，性能接近 GMP 并内置丰富的数论与数学函数。
Stars: 2 | Forks: 0
# BigInts
为 Pascal 提供的任意精度整数和小数，封装在单个自包含单元 `BigInts` ([bigints.pas](bigints.pas)) 中。包含三种值类型，支持普通数字拥有的所有运算符，除了可用内存外没有大小限制，除 RTL 外无任何依赖。
| 类型 | 描述 |
|---|---|
| `BigInt` | 有符号数；位运算符使用具有无限符号扩展的二进制补码语义，类似于 Python 整数 |
| `UBigInt` | 无符号数；任何会导致值降到零以下的操作都会引发 `ERangeError` |
| `BigDecimal` | 十进制浮点数：一个 `BigInt` 尾数乘以 10 的幂；精确的 `+ - *`，除法可选择精度 |
## 功能
- 完整的运算符覆盖：`+ - * div mod / ** shl shr and or xor not`，所有比较操作，`inc`/`dec`，复合赋值（`+=`, `*=` 等），一元 `+`/`-`
- 与任意一侧的普通整数进行混合运算，支持从 `Int64`/`QWord`/string 隐式转换，双向显式转换（包括 `Double`；整数转换永远不会通过浮点数进行四舍五入）
- 任意大小的字面量，支持 `_` 分隔符以及 `$ 0x % 0b & 0o` 前缀；支持 2 到 36 进制之间的解析和格式化
- 乘法：Schoolbook、Karatsuba 和 Toom-3，由可调阈值决定选择策略，并配有专门的平方计算路径
- 除法：Knuth 算法 D，外加针对大数的分治法进制转换
- 模运算：Montgomery `modPow`（加上常数时间复杂度的 `modPowSec`）、`modInverse`、`modSqrt` (Tonelli-Shanks)、`sqrtModN`、`crt`、`discreteLog`
- 质数：Miller-Rabin `isProbablePrime`（在 3.3e24 以下为确定性判定）、Baillie-PSW `isPrime`、`nextPrime`/`prevPrime`、`randomPrime`/`randomSafePrime`/`randomStrongPrime`、精确的 `primePi`/`primeCount`
- 因式分解：试除法加上 Pollard-Brent rho 算法，质因数及其指数被分组为 `(p, e)` 元组；由此衍生出乘法函数 `eulerPhi`、`carmichaelLambda`、`moebius`、`sigma`、`tau`、`divisors`、`radical`
- 数论与组合数学：Lehmer `gcd`、`gcdExt`、`jacobi`、`kronecker`、`continuedFraction`，以及 `factorial`、`fibonacci`、`lucas`、`binomial`、`multinomial`、`catalan`、`bell`、`stirling1`/`stirling2`、`bernoulli`、`partitions`、`subfactorial`、`primorial`
- 随机数：可插拔的生成器（xoshiro256**、PCG64、splitmix64、`System.Random`、OS 熵），确定性种子生成，均匀分布的 `randomBelow`/`randomRange`
- 互操作性：支持大端序和小端序的字节序列化、`hashCode`、带分组的数字输出
- 小数：精确的十进制运算（`0.1 + 0.2 = 0.3`），任意精度的除法和开方，六种舍入模式，双向的最短和精确浮点数转换，以及完整的分析工具箱 - `pi`、`exp`、`ln`、`log2`/`log10`/`logBase`、分数幂、三角函数、双曲函数、`gamma`、`erf`、`atan2`、`hypot`、任意精度的 `agm`（见下文的 BigDecimal 章节）
- 速度：在 x64 平台上，核心操作的速度实测为 GMP 的 1.3 到 4.6 倍（基准测试见下文）；内层循环使用汇编语言编写，并在 `USEASM` 定义下提供纯 Pascal 回退方案
## 快速入门
```
program quickstart;
{$mode unleashed}
uses BigInts;
begin
var a: BigInt := '123456789012345678901234567890';
var b: BigInt := '-0xDEAD_BEEF';
writeln($'{a * b}');
writeln((BigInt(2) ** 4096).digitCount); // 1234
var (q, r) := a.divMod(b);
writeln($'{q} rem {r}');
var p := UBigInt.randomPrime(256);
writeln(p.isProbablePrime); // TRUE
for var (f, e) in UBigInt(720).factorize do
write($'{f}^{e} '); // 2^4 3^2 5^1
writeln;
{$ifdef WINDOWS}readln;{$endif}
end.
```
## 方法
所有方法均采用 camelCase 命名，并可通过代码补全功能发现。除非另有说明，方法同时存在于两种类型上。
### 转换输入与输出
| 方法 | 备注 |
|---|---|
| `parse(s)`, `parse(s, base)` | 静态方法；自动检测 `$ 0x % 0b & 0o` 前缀，允许使用 `_` 分隔符和正负号 |
| `tryParse(s, out v)`, `tryParse(s, base, out v)` | 静态方法；遇到错误输入时不抛出异常 |
| `toString`, `toString(base)` | 支持 2..36 进制；负数在任何进制下均格式化为正负号加上绝对值 |
| `toHex`, `toBin`, `toOct` | 用于 16、2、8 进制的简写 |
| `toStringGrouped(sep = '_', groupSize = 3)` | `1_234_567` 样式的输出 |
| `toInt8`, `toUInt8`, `toInt16`, `toUInt16`, `toInt32`, `toUInt32`, `toInt64`, `toUInt64`, `toDouble` | 如果值超出范围，则引发 `ERangeError` |
| `toInt128`, `toUInt128` | 适用于编译器提供原生 128 位类型的目标平台 |
| `fitsInInt8`, `fitsInUInt8`, ... `fitsInInt64`, `fitsInUInt64`（以及在可用情况下的 `fitsInInt128`/`fitsInUInt128`） | 针对每种原生位宽的对应检查 |
| `toUBigInt` / `toBigInt` | 跨越有符号与无符号的边界；如果为负值，则引发 `ERangeError` |
| `toDecimal` | 扩展为 `BigDecimal`（精确转换，绝不进行四舍五入） |
| `toBytesLE`, `toBytesBE`, `fromBytesLE`, `fromBytesBE` | `UBigInt`：原始绝对值；`BigInt`：带符号位的最小二进制补码，类似于 Java 的 `toByteArray` |
### 谓词与符号
| 方法 | 备注 |
|---|---|
| `isZero`, `isOne`, `isEven`, `isOdd`, `isPowerOfTwo` | |
| `sign` | 返回 -1、0 或 1 |
| `isNegative`, `isPositive` | 仅限 `BigInt` |
| `abs`, `magnitude`, `negate` | 仅限 `BigInt`；`magnitude` 返回 `UBigInt` 类型的绝对值 |
### 位操作
| 方法 | 备注 |
|---|---|
| `bitLength`, `popCount`, `lowestSetBit` | |
| `testBit(i)`, `setBit(i)`, `clearBit(i)`, `flipBit(i)`, `bits[i]` | 在 `BigInt` 上，这些操作基于无限的二进制补码扩展 |
| `complement(width)` | 仅限 `UBigInt`：对低 `width` 位进行按位非操作；由于无符号值不存在无限的补码，因此 `UBigInt` 没有 `not` 运算符 |
### 比较与除法
| 方法 | 备注 |
|---|---|
| `compare`, `equals`, `min`, `max` | 外加全套比较运算符 |
| `divMod(d)` | 执行一次除法，返回 `(q, r)` 元组 |
| `floorDiv(d)`, `floorMod(d)` | 仅限 `BigInt`；类似于 Python，向负无穷方向舍入 |
| `ceilDiv(d)` | 向正无穷方向舍入 |
| `swap(other)`, `hashCode`, `digitCount` | |
### 数学运算
| 方法 | 备注 |
|---|---|
| `sqr`, `sqrt`, `nthRoot(n)`, `nthRootRem(n)` | 平方和整数（下取整）开方；`nthRootRem` 还会返回余数 |
| `isKthPower(k)` | 判断该值是否为完全的 `k` 次方 |
| `pow(e)`, `**` | 常规求幂 |
| `modPow(e, m)` | 对奇数 `m` 使用窗口指数的 Montgomery 算法；在 `BigInt` 上，模数必须为正数，结果落在 `0..m-1` 区间内，负指数通过模逆元进行计算 |
| `modPowSec(e, m)` | 结果与 `modPow` 相同，但运算序列不会根据指数的位产生分支（具备抗侧信道攻击能力，适用于保密指数） |
| `modInverse(m)` | 当逆元不存在时引发 `EBigIntError` |
| `gcd`, `lcm` | Lehmer gcd 算法 |
| `isProbablePrime(rounds = 24)` | Miller-Rabin 算法；在 3.3e24 以下使用确定性见证，大于该值使用随机轮次 |
| `isPrime` | Baillie-PSW（先进行确定性小范围测试，然后进行强基数 2 的 Miller-Rabin 测试以及强 Lucas 测试）；目前无已知的反例 |
| `nextPrime` | 大于自身值的下一个质数 |
### 常量与生成器（类函数）
| 方法 | 备注 |
|---|---|
| `zero`, `one`, `two`, `ten`, `minusOne` | `minusOne` 仅限 `BigInt` |
| `pow2(n)` | |
| `random(bits)` | 在 `2^bits` 以下均匀分布；生成器是可插拔的，详见扩展功能章节 |
| `factorial(n)` | 二进制拆分 |
| `fibonacci(n)` | 快速倍增 |
## 扩展功能
构建在核心算术之上的可选数学层。
### 随机数
| 方法 | 备注 |
|---|---|
| `randomBelow(bound)` | 在 `0..bound-1` 均匀分布，采用拒绝采样 |
| `randomRange(lo, hi)` | 在 `lo..hi` 均匀分布，包含两端；在 `BigInt` 上支持负数边界 |
| `randomPrime(bits, rounds = 24)` | 精确的位长度：最高位置 1、奇数、经过 Miller-Rabin 测试 |
| `randomSafePrime(bits)` | 安全质数 `p`，即 `(p-1)/2` 也是质数 |
| `randomStrongPrime(bits)` | Gordon 算法：`p-1` 和 `p+1` 都包含大的质因数 |
`random` 及其相关方法的后端由 `BigIntRngAlgo` 变量选择：
| 生成器 | 备注 |
|---|---|
| `rngXoshiro256ss` | 默认选项；xoshiro256** |
| `rngPcg64` | 带有参考乘数和流的 PCG XSL-RR 128/64 |
| `rngSplitMix64` | 小巧且快速；内部也用于扩展种子 |
| `rngSystem` | 基于历史 `RandSeed` 的 `System.Random` 流 |
| `rngOS` | 每次调用都获取新的 OS 熵（RtlGenRandom, /dev/urandom）；用于生成密钥时请选择此项 |
生成器状态为每个线程独立（`threadvar`）：线程中的第一次随机抽取会从 OS 熵进行自我播种，因此未设种子的值每次运行都会不同，且线程间拥有独立的流——没有共享状态，也没有锁。`BigIntRandomSeed(seed)` 使调用线程变得可复现（它同时会设置 `RandSeed`，因此 `rngSystem` 模式也会随之改变）；`BigIntRandomize` 会根据 OS 熵重新播种。如果需要跨线程的复现性，请分别对每个线程进行播种。`rngSystem` 模式保留了普通的 `System.Random` 契约（由 `RandSeed` 驱动，延迟自动播种机制对其不产生影响）。
### 数论
| 方法 | 备注 |
|---|---|
| `gcdExt(other)` | 仅限 `BigInt`；扩展欧几里得算法，返回 `(g, x, y)` 元组，满足 `a*x + b*y = g` |
| `jacobi(n)` | 雅可比符号，适用于奇数正整数 `n`，返回 -1、0 或 1 |
| `kronecker(n)` | 克罗内克符号，是雅可比符号对任意整数的完整扩展（处理因数 2 和负数参数） |
| `modSqrt(p)` | 模质数平方根（Tonelli-Shanks 算法）；如果是非二次剩余则引发 `EBigIntError` |
`sqrtModN(n)` | 模合数 `n` 的所有平方根（分解、提升、CRT）；需要满足 `gcd(self, n) = 1`，如果是非二次剩余则返回空数组 |
| `discreteLog(target, m)` | 小步大步算法：寻找满足 `self^x = target (mod m)` 的最小 `x`，若无则返回 -1；返回 `Int64`，适用于小规模实例 |
| `crt(remainders, moduli)` | `BigInt` 类函数；用于两两互质的正模数的中国剩余定理 |
| `isPerfectSquare` | 快速 mod-16 过滤器，然后进行精确的根检查 |
| `sqrtRem` | 返回 `(root, rem)` 元组，满足 `self = root^2 + rem` |
| `prevPrime` | 小于自身值的最大质数；如果 `self <= 2` 则引发异常 |
| `factorize` | 返回按质数升序排列的 `(p, e)` 元组数组；10^4 以下使用试除法，以上使用 Pollard-Brent rho 算法；`BigInt.factorize` 对绝对值进行因式分解 |
`factorize` 的运行时间随第二大质因数的平方根成正比增长，因此两个随机大质数的乘积将需要极长的时间才能分解——这是因式分解的本质，而不是 bug。
以下方法直接基于因式分解结果运行（因此其开销等同于 `factorize`）：
| 方法 | 备注 |
|---|---|
| `eulerPhi` | 欧拉函数，即不超过 `n` 且与 `n` 互质的整数个数 |
| `carmichaelLambda` | 群的指数：对于所有互质的 `a`，满足 `a^k = 1 (mod n)` 的最小 `k` |
| `moebius` | 莫比乌斯函数，返回 -1、0 或 1 |
| `sigma(k = 1)` | 约数的 `k` 次方之和；`sigma(0)` 即为 `tau` |
| `tau` | 约数的个数 |
| `radical` | 不同质因数的乘积 |
| `divisors` | 升序排列的所有约数 |
| `isSquarefree`, `isPerfect`, `isCarmichael` | 对应的谓词判定 |
质数计数和有理近似（`BigInt`/`UBigInt` 类函数）：
| 方法 | 备注 |
|---|---|
| `primePi(n)` | 使用分段筛法得出 `<= n` 的质数的精确数量（`QWord`，实际适用于约 1e10） |
| `primeCount(lo, hi)` | 在 `lo..hi` 范围内质数的精确数量 |
| `continuedFraction(num, den)` | `num/den` 连续分数的系数 |
| `fromContinuedFraction(cf)` | 将系数计算并还原为最简分数 `(num, den)`；对 `cf` 进行切片可获得渐近分数 |
### 组合数学
全部为类函数。
| 方法 | 备注 |
|---|---|
| `lucas(n)` | 斐波那契的伴随序列，采用一次快速倍增运算 |
| `binomial(n, k)` | 乘法形式，每个中间除法都是精确的 |
| `multinomial(ks)` | `(sum ks)! / prod(ks[i]!)`，以二项式之积表示 |
| `catalan(n)` | `binomial(2n, n) div (n + 1)` |
| `primorial(n)` | 所有不超过 `n` 的质数的乘积，使用奇数筛法加上平衡乘法 |
| `risingFactorial(x, n)`, `fallingFactorial(x, n)` | Pochhammer 符号；`BigInt` 接受负数底数 |
| `subfactorial(n)` | 错排数 `!n` |
| `bell(n)` | 贝尔数，通过贝尔三角形计算 |
| `stirling1(n, k)` | 第一类有符号斯特林数（`BigInt`） |
| `stirling2(n, k)` | 第二类斯特林数 |
| `partitions(n)` | 整数分拆数 `p(n)`，使用欧拉五边形定理递推 |
| `bernoulli(n)` | `BigInt` 类函数；伯努利数，返回精确的最简 `(num, den)` 分数 |
### 罗马数字与单词
`BigInt` 和 `UBigInt` 还支持将自身格式化为人类可读的形式：
| 方法 | 备注 |
|---|---|
| `toRoman` | 将 `1..3999` 之间的数值转换为罗马数字 |
| `toWords` | 英语短进位制单词（`one million two hundred thirty-four thousand ...`），最高达 `10^66` |
## BigDecimal
基于相同整数核心构建的任意精度十进制浮点数：一个值由一个 `BigInt` 尾数乘以 10 的幂构成，并保持规范形式（没有尾随零）。`0.1` 就是精确的 `0.1`，财务计算永远不会发生漂移，尾数也能完全利用整数引擎的速度。
- `+ - *` 总是精确的，`div`/`mod`（整数商和精确余数）也一样。
- 默认情况下，`/` 会四舍五入保留 18 位小数；`divide(b, precision)` 可以自定义。商总是保留完整的整数部分以及至少 `precision` 位有效数字，并包含一位会被 `toString` 舍去的隐藏保护位，因此 `(1/3) * 3` 会打印为 `1`。完全精确的商仍然是精确的：`1 / 8` 结果为 `0.125`。
- 比较运算是数值层面的（`0.5 = 5E-1`），且能看到存储的保护位，所以即使 `(1/3) * 3` 打印出来是 `1`，`(1/3) * 3 < 1` 仍然成立。
- 与整数、字符串、`BigInt` 和 `UBigInt` 的混合运算可隐式转换；浮点数只能通过显式转换或 `from*` 构造器进行转换，因此不会有未知的二进制舍入误差混入。
```
program decimals;
{$mode unleashed}
uses BigInts;
begin
var price: BigDecimal := '19.99';
writeln($'{price * 3}'); // 59.97
writeln($'{BigDecimal(1) / 3}'); // 0.333333333333333333
writeln($'{BigDecimal(1) / 3 * 3}'); // 1
writeln($'{BigDecimal(2).sqrt(30)}'); // 1.41421356237309504880168872421
writeln($'{BigDecimal.fromDouble(0.1)}'); // 0.1
writeln($'{BigDecimal.fromDoubleExact(0.1)}'); // 0.1000000000000000055511151231257827021181583404541015625
var pi: BigDecimal := '3.14159265';
writeln($'{pi.rounded(-2)} {pi.rounded(-2, bdrCeil)} {pi.trunc}'); // 3.14 3.15 3
writeln($'{BigDecimal('123456.789').toScientific}'); // 1.23456789E5
{$ifdef WINDOWS}readln;{$endif}
end.
```
解析层运行在相同的缩放整数核心之上：每个函数都接受一个 `precision` 参数（小数位数，默认为 18），并像除法一样通过隐藏的保护位对其显示的最后一位进行舍入。`pi` 来源于 Chudnovsky 二进制拆分并进行了缓存，巨大的三角函数参数会使用匹配精度的 pi 进行规约。
```
program analytic;
{$mode unleashed}
uses BigInts;
begin
writeln($'{BigDecimal.pi(50)}'); // 3.14159265358979323846264338327950288419716939937511
writeln($'{BigDecimal(2).ln(40)}'); // 0.6931471805599453094172321214581765680755
writeln($'{BigDecimal(2) ** BigDecimal('0.5')}'); // 1.414213562373095049
writeln($'{BigDecimal(1).sin(40)}'); // 0.8414709848078965066525023216302989996226
writeln($'{BigDecimal('1E6').logBase(BigDecimal(10))}'); // 6
writeln($'{BigDecimal('19.99').quantize(BigDecimal('0.05'))}'); // 20
var (num, den) := BigDecimal('0.375').toFraction;
writeln($'{num}/{den}'); // 3/8
writeln(BigDecimal('0.000123').toEngineering); // 123E-6
{$ifdef WINDOWS}readln;{$endif}
end.
```
| 方法 | 备注 |
|---|---|
| `parse(s)`, `tryParse(s, out v)`, `:=` from string | `[sign]digits[.digits][E[sign]digits]`，允许使用 `_` 分隔符 |
| `calc(s, precision = 18)`, `tryCalc(s, out v, precision = 18)` | 计算整个表达式字符串，详情见下方的计算器部分 |
| `toString`, `toScientific`, `toEngineering` | 普通格式 `-123.45` / 标准化格式 `-1.2345E2` / 指数为 3 的倍数的工程格式，如 `123E-6` |
| `toInt8`..`toInt64`, `toUInt8`..`toUInt64`, `toBigInt`, `toUBigInt`, `fitsInInt8`..`fitsInUInt64` | 精确转换：如果不属于整数或超出范围则引发 `ERangeError`（如果是负数，`toUInt*`/`toUBigInt` 也会报错） |
| `trunc`, `floor`, `ceil`, `round` | 转换为 `BigInt`：分别向零、向负无穷、向正无穷舍入，以及向最近的偶数舍入（如同 Pascal 的 `round`） |
| `frac` | `trunc` 丢弃的小数部分，满足 `self = trunc + frac` |
| `toFraction` | 转换为精确的分数形式 `(num, den)` 元组：`0.375` 变为 `(3, 8)` |
| `rounded(toDigit = 0, mode = bdrRound)` | 在任何小数位进行舍入：`0` = 整数，`-2` = 分，`3` = 千位；模式包括 `bdrTrunc bdrCeil bdrFloor bdrRound bdrHalfUp bdrHalfEven` |
| `quantize(step, mode = bdrRound)` | 舍入到任意步长的最近倍数，如 `0.05` |
| `divide(b, precision = 18)`, `divMod(d)` | 按指定精度进行除法运算 / 返回带有精确余数的整数商 |
| `fromDouble`, `fromSingle`, explicit float casts | 能够读取并还原为相同浮点数的最短小数表示：`0.1` 得到 `0.1` |
| `fromDoubleExact`, `fromSingleExact` | 精确的二进制值表示：`0.1` 会给出全部 55 位数字 |
| `toDouble`, `toSingle` | 正确地舍入到最近的浮点数，逢半向上；向上溢出得到无穷大，向下溢出得到零 |
| `toExtended`, `fromExtended`, `fromExtendedExact` | 适用于支持 80 位浮点数类型的目标平台 |
| `sqrt(precision = 18)`, `nthRoot(n, precision = 18)` | 保留 `precision` 位小数，带有与除法相同的隐藏保护位 |
| `pow(e)` | 当 `e >= 0` 的整数时完全精确；负数指数则按默认精度作除法处理 |
| `pow(y, precision = 18)`, `**` | 通过 `exp(y * ln x)` 计算分数指数幂 |
| `exp`, `ln`, `log2`, `log10`, `logBase(b)` | 均接受 `(precision = 18)` 参数；`log10` 对 10 的幂是精确的，`log2` 对 2 的幂是精确的 |
| `sin`, `cos`, `tan`, `arcsin`, `arccos`, `arctan` | 弧度制；大参数会在匹配精度下对 pi/2 进行模运算规约 |
| `sinh`, `cosh`, `tanh` | 在相同的指数核心上计算双曲函数 |
| `gamma`, `lnGamma` | Gamma 函数（通过反射公式覆盖负数）及其对数（正参数） |
| `factorial` | 实数阶乘 `x! = gamma(x+1)`，对于小的非负整数是精确的 |
| `erf`, `erfc` | 误差函数及其补函数；`erfc` 在参数较大时使用连续分数计算 |
| `continuedFraction(maxTerms = 0)` | 该值的（有限、精确的）连分数；非常适合用于寻找最佳有理近似 |
| `roundToSignificant(digits, mode = bdrRound)` | 舍入到指定的有效数字位数，而不是小数位数 |
| `pi(precision)`, `e(precision)` | 类函数；pi 的值在多次调用间会被缓存 |
| `atan2(y, x, precision)`, `hypot(x, y, precision)`, `agm(a, b, precision)` | 类函数：感知象限的反正切、欧几里得长度、算术几何平均值 |
| `gcd`, `lcm` | 基于十进制网格的运算：`gcd(0.25, 0.15) = 0.05` |
| `precision`, `mostSignificantExponent`, `getDigit(i)` | 有效数字位数、最高位数字的指数、位于 `10^i` 上的数字 |
| `shift10(n)`, `shifted10(n)` | 在不改变尾数的情况下乘以 10 的幂 |
| `isZero`, `isOne`, `isIntegral`, `isEven`, `isOdd`, `isNegative`, `isPositive`, `sign`, `abs`, `negate` | 谓词和符号辅助方法；带小数部分的值既不是偶数也不是奇数 |
| `compare`, `equals`, `approxEquals(other, eps)`, `min`, `max`, `hashCode`, `swap` | 外加完整的运算符和比较集合 |
| `zero`, `one`, `two`, `ten` | 类常量 |
### 计算器
`calc` 会根据您传入的精度对整个字符串表达式进行求值，并返回一个 `BigDecimal`。它是一个纯粹、无状态的求值器：没有变量，没有赋值，相同的字符串总是产生相同的结果。`tryCalc` 在遇到错误表达式时返回 `false` 而不是抛出异常；而 `calc` 在遇到语法错误时会抛出带有字符位置信息的 `EConvertError`，并让数学错误（`EDivByZero`、`ERangeError`、EBigIntError`）直接透传抛出。
```
writeln(BigDecimal.calc('2^100 / 3').toScientific); // 4.2...E29
writeln(BigDecimal.calc('sin(pi/6)', 30)); // 0.5
writeln(BigDecimal.calc('(1 + sqrt(5)) / 2', 50)); // the golden ratio, 50 digits
```
运算符（从最低优先级到最高优先级）：
| 运算符 | 含义 | 备注 |
|---|---|---|
| `+` `-` | 加、减 | 左结合 |
| `*` `/` `div` `mod` `%` | 乘法、实数除法、整数除法、余数 | `/` 对于 `10/4` 会给出 `2.5`；`div` 给出 `2`；`%` 等同于 `mod` |
| `-` `+` (一元) | 符号 | 优先级低于幂运算，因此 `-2^2 = -4` |
| `^` `**` | 幂运算 | 右结合，`2^3^2 = 512`；`2^-3 = 0.125` |
| `!` (后缀) | 阶乘 | 作用于非负整数 |
函数（名称不区分大小写）：
| 组别 | 函数 |
|---|---|
| 根与幂 | `sqrt(x)` `cbrt(x)` `root(x, n)` `pow(x, y)` `sqr(x)` |
| 指数与对数 | `exp(x)` `ln(x)` `log(x)`=log10, `log(x, b)`=以 b 为底, `log2(x)` `log10(x)` `logb(x, b)` |
| 三角函数 | `sin cos tan (x)`, `asin acos atan (x)` (也可写作 `arcsin`/`arccos`/`arctan`), `sinh cosh tanh (x)` |
| 舍入函数 | `floor(x)` `ceil(x)` `round(x)` `trunc(x)` |
| 特殊函数 | `gamma(x)` `lngamma(x)` `erf(x)` `erfc(x)` `factorial(x)` |
| 双参数函数 | `min max gcd lcm atan2 hypot agm (a, b)` |
| 常量 | `pi` `e` `tau` `phi` |
每个函数都会在工作精度下映射到同名的方法，因此 `calc('sin(1)', 40)` 等同于 `BigDecimal(1).sin(40)`。目前仅公开了接收数字并返回单个数字的函数；整数数论相关的方法和返回元组的方法仍保留在 API 层面。
## 值得了解的语义
- `div`/`mod` 像 Pascal 一样向零截断；`floorDiv`/`floorMod` 像 Python 一样向负无穷舍入；`ceilDiv` 向上舍入。
- `/` 是整数除法，与 `div` 相同（针对整数类型遵循 C 语言家族的惯例）。
- 负数 `BigInt` 的 `shr` 是算术右移（向负无穷舍入）；`shl` 会保留符号位。
- 负数 `BigInt` 的位运算使用带有无限符号扩展的二进制补码；`not x = -x-1`。
- 在任何进制下，负数的格式化形式均为正负号加上绝对值：`-255` 在十六进制下输出为 `-FF`。
- 复制这些值开销很低：小值（不超过 256 位）直接内联在值本身中，较大的值则共享一个引用计数的内存块，但进行修改的方法会优先取消共享，因此绝不会有变量在背后被意外改变。
- `0 ** 0 = 1`，除以零会引发 `EDivByZero`，无法容纳的转换会引发 `ERangeError`，解析错误会引发 `EConvertError`，定义域错误（负指数、无逆元、非二次剩余）会引发 `EBigIntError`。
## 性能
在 x86_64 平台上使用带有汇编内层循环的 64 位肢（用于乘法/带进位加法的行原语，并在 CPU 支持 ADX 时在运行时选用 mulx/adcx/adox 的 `addmul_1` 实现）；在其他平台上则使用可移植的 32 位 Pascal 肢。汇编代码隐藏在单元顶部的 `USEASM` 定义之下——将其注释掉即可获得完全可移植的纯 Pascal 版本（在 x64 上核心操作大约慢 4 到 8 倍）。实现了 Knuth 算法 D 除法、超过可调阈值（`BigIntKaratsubaThreshold`、`BigIntToom3Threshold`）的 Karatsuba 和 Toom-3 乘法与平方算法、带有窗口指数的 Montgomery modPow、分治法进制转换以及 Lehmer gcd。最大 256 位（`BIGINT_INLINE_LIMBS` 肢）的数值直接内联存储在变量中，无需分配堆内存；较大的结果则直接在热点路径中构建大小精确的堆缓冲区。在 x64 桌面平台上：计算 `factorial(50000)` 约需 ~12 ms，计算 `fibonacci(1000000)` 约需 ~10 ms。
### 对比 GMP 的基准测试
在一台 x64 桌面电脑上与 GMP 6.3.0（随 Git for Windows 附带的 64 位肢 `libgmp-10.dll`）进行对比测试，双方都启用了 `-O3`，结果为单次操作耗时。GMP 方面重用了其 mpz 目标变量，这符合 GMP 代码的常规编写方式；而 BigInts 方面为每次操作分配全新的值，这正是值语义所必须付出的代价。
| 操作 | BigInts | GMP | 比率 |
|---|---|---|---|
| 128位加法 | 0.015 us | 0.005 us | 3.0x |
| 1024位加法 | 0.037 us | 0.008 us | 4.6x |
| 16384位加法 | 0.148 us | 0.066 us | 2.2x |
| 262144位加法 | 1.78 us | 1.20 us | 1.5x |
| 128位乘法 | 0.019 us | 0.006 us | 3.2x |
| 1024位乘法 | 0.181 us | 0.13 us | 1.4x |
| 8192位乘法 | 6.23 us | 4.20 us | 1.5x |
| 65536位乘法 | 180 us | 82.1 us | 2.2x |
| 262144位乘法 | 1560 us | 548 us | 2.8x |
| 65536x1024位乘法 | 6.71 us | 8.64 us | 0.8x |
| 8192位平方 | 5.32 us | 2.59 us | 2.1x |
| 65536位平方 | 178 us | 54.9 us | 3.2x |
| 2048/1024位 divmod | 0.48 us | 0.283 us | 1.7x |
| 8192/4096位 divmod | 3.50 us | 2.63 us | 1.3x |
| 131072/65536位 divmod | 709 us | 208 us | 3.4x |
| 4096位 toString | 8.93 us | 4.25 us | 2.1x |
| 65536位 toString | 423 us | 229 us | 1.8x |
| 4096位 parse | 5.89 us | 3.83 us | 1.5x |
| 65536位 parse | 230 us | 132 us | 1.7x |
| 512位 modPow | 79.6 us | 37.9 us | 2.1x |
| 1024位 modPow | 431 us | 260 us | 1.7x |
| 2048位 modPow | 2710 us | 1960 us | 1.4x |
| 1024位 gcd | 10.6 us | 2.48 us | 4.3x |
| 16384位 gcd | 413 us | 114 us | 3.6x |
批量算术运算的性能达到了 GMP 的 1.3 至 4.6 倍。剩余的性能差距主要源于 GMP 采用的手写优化汇编、其在处理超大操作数时使用的更高阶 Toom 算法和 FFT，以及本单元尚未实现的次平方阶 gcd 和除法算法。最高 256 位的小数值直接内联存储而无需分配内存，因此一肢或两肢长度的加法或乘法运算耗时约为 15-20 纳秒，即便在这种尺寸下，其速度依然保持在 GMP 的 3 倍以内，否则对于这种每次操作都新建值的库来说，其耗时将主要耗费在内存分配上。
## 示例
以下每个代码块都是一个完整的程序：将其复制到 `.lpr` 文件中，并将 `bigints.pas` 放在同一路径下，即可直接编译并运行。
### 字面量与格式化
```
program literals;
{$mode unleashed}
uses BigInts;
begin
var a: UBigInt := '123_456_789_000_000_000_000_000';
var b: BigInt := '-0xDEAD_BEEF';
var c: UBigInt := '%1010_1010';
writeln(a.toStringGrouped); // 123_456_789_000_000_000_000_000
writeln(b.toString); // -3735928559
writeln(c.toString(36)); // 4Q
writeln(UBigInt.parse('zz', 36).toString); // 1295
{$ifdef WINDOWS}readln;{$endif}
end.
```
### 各种形式的除法
```
program division;
{$mode unleashed}
uses BigInts;
begin
var (q, r) := BigInt(-7).divMod(BigInt(2));
writeln($'{q} {r}'); // -3 -1 (truncated, like Pascal div/mod)
writeln(BigInt(-7).floorDiv(2).toString); // -4 (like Python)
writeln(BigInt(-7).floorMod(2).toString); // 1
writeln(UBigInt(7).ceilDiv(UBigInt(2)).toString); // 4
{$ifdef WINDOWS}readln;{$endif}
end.
```
### 二进制补码位运算
```
program bitwise;
{$mode unleashed}
uses BigInts;
begin
writeln((BigInt(-1) and BigInt($FF)).toString); // 255: -1 is an infinite run of ones
writeln((not BigInt(0)).toString); // -1
writeln((BigInt(-5) shr 1).toString); // -3: arithmetic shift
writeln(BigInt(-255).toHex); // -FF: sign plus magnitude in every base
{$ifdef WINDOWS}readln;{$endif}
end.
```
### 质数与一个简易 RSA 示例
```
program rsa;
{$mode unleashed}
uses BigInts;
begin
BigIntRandomize;
var p := UBigInt.randomPrime(512);
var q := UBigInt.randomPrime(512);
var n := p * q;
var e: UBigInt := 65537;
var d := e.modInverse((p - 1) * (q - 1));
var msg: UBigInt := '0x48656C6C6F21'; // "Hello!"
var cipher := msg.modPow(e, n);
writeln(cipher.modPow(d, n) = msg); // TRUE
{$ifdef WINDOWS}readln;{$endif}
end.
```
### 随机数套件
```
program random_suite;
{$mode unleashed}
uses BigInts;
begin
BigIntRngAlgo := rngPcg64;
BigIntRandomSeed(42); // reproducible from here on
writeln(UBigInt.random(128).toHex); // F6A4492CA8314B92F0D3403191F1E9AF
writeln(UBigInt.randomBelow(UBigInt.ten ** 20).toString);
writeln(BigInt.randomRange(-50, 50).toString);
{$ifdef WINDOWS}readln;{$endif}
end.
```
### 因式分解
```
program factor;
{$mode unleashed}
uses BigInts;
begin
var n: UBigInt := '123456789012345678';
for var (p, e) in n.factorize do
write($'{p}^{e} '); // 2^1 3^3 21491747^1 106377431^1
writeln;
{$ifdef WINDOWS}readln;{$endif}
end.
```
### 中国剩余定理
```
program crt_demo;
{$mode unleashed}
uses BigInts;
begin
// x = 2 (mod 3), x = 3 (mod 5), x = 2 (mod 7)
var x := BigInt.crt([BigInt(2), BigInt(3), BigInt(2)], [BigInt(3), BigInt(5), BigInt(7)]);
writeln(x.toString); // 23
{$ifdef WINDOWS}readln;{$endif}
end.
```
## 许可证
本项目采用 Mozilla Public License 2.0 (MPL-2.0) 进行授权。详情请参阅 LICENSE 文件。
`examples/` 目录下的示例程序可无限制免费使用；您可以自由地将其中代码复制到您自己的项目中，且无需承担任何 MPL 义务。
标签：Pascal, 任意精度数学库, 大整数计算, 密码学基础, 数论算法, 算法库