zenith-float 1.0.0

Arbitrary-precision software floating-point arithmetic on integer limbs.
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
# zenith-float library inventory

This is a complete inventory of what the **`zenith-float` crate** exposes as a software math library. Arithmetic uses integer limbs only. Hardware IEEE binary interchange formats are not part of this crate; convert to those formats in the caller if needed.

Related docs: [getting started](GETTING_STARTED.md), [help](HELP.md), [error bounds](README.md), [`expr!` rounding](EXPR.md).

License: MIT OR Apache-2.0.

---

## 1. Crate layout

| Crate | Role |
| --- | --- |
| **`zenith-float`** | Public package. Re-exports the numeric kernel and the macros `expr!`, `cexpr!`, `exact!`, `fbig!`. |
| **`zenith-float-num`** | Numeric kernel. Applications should depend on `zenith-float`, not this crate. |
| **`zenith-float-macro`** | Procedural macros. Not a direct application dependency. |
| **`zenith-float-compare`** | Workspace comparison benches against other float crates. Not part of the public math API. |

`no_std` is supported when a memory allocator is available (`default-features = false` to drop `std`).

---

## 2. Cargo features (`zenith-float`)

| Feature | Default | What it enables |
| --- | --- | --- |
| `std` | yes | `Display` / radix format traits, `FromStr`, `std::error::Error` for `Error`, `SharedConsts`, serde when `serde` is on. |
| `random` | no | `random_uniform` / `random_gaussian` / `random_exponential` / `random_fill`; existing `random_normal(p, exp_from, exp_to)` mantissa draw; `seeded_random`, `reseed_random`. |
| `serde` | no | `Serialize` / `Deserialize` for `ExactNum` / `ExactComplex` / `ExactRational` / `ExactInt` / arrays / `Ball`. Decimal strings carry `@p=`. Implies `std`. |
| `hdf5` | no | `to_hdf5` / `from_hdf5` on IEEE and `ExactNum` arrays via crates.io `hdf5-rust`. Implies `std`. No `libhdf5`. |
| `mpfr-tests` | no | Optional MPFR / GMP comparison tests in the kernel crate (Linux x86_64, `rug`). Not a runtime math engine. |

---

## 3. Numeric model

A finite `ExactNum` is a **sign**, a **binary mantissa** stored as an array of words, and an **`Exponent` (`i32`)**. Precision is a bit count, always rounded **up to the word size**.

| Item | Value |
| --- | --- |
| `Word` | `u64` except on 32-bit targets, where it is `u32` |
| `WORD_BIT_SIZE` | 64 or 32 accordingly |
| `WORD_MAX` | `Word::MAX` |
| `WORD_BASE` | `WORD_MAX + 1` as a double-width integer |
| `WORD_SIGNIFICANT_BIT` | high bit of a word |
| `INLINE_WORDS` | `2` — small mantissas stored inline |
| `EXPONENT_BIT_SIZE` | 32 |
| `EXPONENT_MAX` | `i32::MAX` on non-32-bit pointers; `i32::MAX / 4` on 32-bit |
| `EXPONENT_MIN` | `i32::MIN` on non-32-bit pointers; `i32::MIN / 4` on 32-bit |
| Default operator precision | 128 bits (`ToEven`) for `+ − × ÷` on `ExactNum` (internal default; not a public constant) |
| `MAX_PREC_RETRY` | `256` — extra-precision budget for correct-rounding retries (`try_set_precision`, `ziv_round`, transcendentals) |

Special values: **`+Inf`**, **`-Inf`**, **`NaN`** (optional associated `Error`). Finite values may be **subnormal** at `EXPONENT_MIN`. Results may be marked **inexact** when bits were rounded or an argument was inexact.

Public sentinels: `NAN`, `INF_POS`, `INF_NEG`.

---

## 4. Rounding

`RoundingMode` (copy, eq, debug):

| Variant | Meaning |
| --- | --- |
| `None` | Skip the final round; extra bits may remain. |
| `Up` | Round half toward +∞ |
| `Down` | Round half toward −∞ |
| `ToZero` | Round half toward zero |
| `FromZero` | Round half away from zero |
| `ToEven` | Round half to even |
| `ToOdd` | Round half to odd |

Methods that take a mode other than `None` round to the requested precision. `expr!` raises working precision for cancellation; it does not itself guarantee correct rounding. See [EXPR.md](EXPR.md).

---

## 5. Errors

`Error`:

| Variant | Typical public mapping |
| --- | --- |
| `ExponentOverflow(Sign)` | `±Inf` |
| `DivisionByZero` | `NaN` |
| `InvalidArgument` | `NaN` |
| `PrecisionRetryExhausted` | `NaN` — Ziv/`MAX_PREC_RETRY` budget, not a domain error |
| `MemoryAllocation` | `NaN` |

`ExactNum::err()` returns the associated error on `NaN`. `Error` implements `Display`; with `std` it implements `std::error::Error`. `From<TryReserveError>` and `From<LayoutError>` → `MemoryAllocation`.

---

## 6. Sign and radix

**`Sign`:** `Pos`, `Neg`. Methods: `invert`, `is_positive`, `is_negative`, `to_int` (`i8`: +1 / −1). `Hash`, `Eq`, `Copy`.

**`Radix`:** bases **2 through 36**. Named constants `Bin` (2), `Oct` (8), `Dec` (10), `Hex` (16). Methods: `try_new(u8)`, `value`, `commensurable_shift` (power-of-two bases), `uses_underscore_exponent` (bases > 10 use `_e` so `e` can be a digit), `bits_per_digit`. `From<Radix> for u8`.

---

## 7. Public types and re-exports

From `zenith_float` / `zenith_float_num`:

- `ExactNum`, `ExactComplex`, `ExactRational`, `ExactInt`, `ExactNumPoly`
- `chebyshev_coeffs`, `chebyshev_eval`, `chebyshev_error_bound`, `clenshaw`
- `gauss_legendre`, `tanh_sinh`, `gauss_laguerre`, `gauss_hermite`
- `bisect`, `newton`, `brent`, `illinois`, `root_default_tol`
- `rk4`, `rk45_adaptive`, `euler`, `ode_min_step`
- `dct`, `idct`, `dst`, `idst`, `fft_real`, `ifft_real`
- `hann_window`, `hamming_window`, `blackman_window`, `kaiser_window`, `rectangular_window`
- `mod_pow`, `mod_inv`, `miller_rabin`, `pollard_rho`
- `sha256`, `sha512`, `hmac_sha256`, `constant_time_eq`
- `Ieee32`, `Ieee64`, `Ieee32Array`, `Ieee64Array`, `ExactNumArray`
- `Consts`, `ConstCache` (alias of `Consts`), `ConstCacheInfo`, `CachedFBig`, `SharedConsts` (`std` only)
- `Context` (module `zenith_float::ctx`), trait `Contextable`
- `Ball`, `ComplexBall`, `ziv_round`, `ziv_round_vec`
- `RadixFloat`
- `FromExt`
- `RoundingMode`, `Radix`, `Sign`, `Error`, `Exponent`, `Word`
- Word/exponent constants listed in §3
- `MAX_PREC_RETRY`, `INLINE_WORDS`
- `NAN`, `INF_POS`, `INF_NEG`
- `POLY_COMPANION_CLOSED_DEG`, `CHEBYSHEV_MAX_DEGREE`, `ORTHOPOLY_N_MAX`, `QUADRATURE_MAX_NODES`, `TANH_SINH_LEVELS_MAX`, `ROOT_MAX_ITER`, `ROOT_DEFAULT_TOL`, `ODE_MAX_STEPS`, `ODE_MIN_STEP`, `DSP_MAX_POINTS`, `IEEE_SIMD_LANE_WIDTH`, `POLLARD_RHO_ITER_MAX`, `BINARY_INLINE_LEN`, `BINARY_INLINE_MANT_BITS`, `BINARY_FORMAT_VERSION`, `BINARY_MAX_U32`, `BINARY_MAX_ELEMS`, `CSV_MAX_ROWS`, `CSV_MAX_COLS`, `PROPTEST_CASES`
- Feature `random`: `random_seed`, `reseed_random`, `seeded_random`, `DEFAULT_RANDOM_SEED`, `RandomDist`

Module `ctx` is public. `macro_util` is `#[doc(hidden)]` and exists for `expr!` / `cexpr!` expansion (`check_exponent_range`, `check_complex_exponent_range`, `complex_cancel_bits`, `compute_added_err`, `ErrAlgo`, `TrigFun`, …). Do not treat it as application API.

Not re-exported: internal `Mantissa`, `WordBuf`, series helpers, `DEFAULT_P`.

---

## 7b. Software IEEE widths and arrays

Hardware floating-point is not used. `Ieee32` / `Ieee64` store IEEE-754 binary32 / binary64 as `u32` / `u64` (`from_bits` / `to_bits`). Arithmetic is integer, to-nearest ties-to-even.

| Type | Notes |
| --- | --- |
| `Ieee32` / `Ieee64` | add/sub/mul/div/sqrt/`mul_add`, classify, `next_up`/`next_down`/`next_after`, `frexp`, `from_i32`, `to_exact` / `from_exact` |
| `Ieee32Array` / `Ieee64Array` | Row-major dense bits (1-D is shape `(1, n)`); elementwise `+ − × ÷` (matching shape), scalar broadcast, `sum`, `dot`, `sqrt`, `fma`, software `matmul`; add/sub/mul/div/sqrt/fma use integer SIMD (`u32`/`u64` lanes, SSE2/NEON) and are **bit-identical** to the scalar integer kernel; specials via widen-to-`ExactNum` (named methods; `cc` is an argument, not a global) |
| `ExactNumArray` | Row-major `ExactNum` at a stored default `p`. Named elementwise methods match the scalar (`int`/`floor`/…, roots, logs, circular/hyperbolic, §15 specials). Methods that take `p`/`rm`/`cc` on `ExactNum` take the same arguments here — including `Consts` when the `expr!` leaf needs a cache. Software `matmul`. `lu_decomp(p, rm)` → `(L, U, P)` or `None` if singular or a workspace reserve fails. `qr_decomp(p, rm)` → `(Q, R)` or `None` on failed reserve; rank-deficient → zero \(R_{kk}\). `svd_decomp(p, rm)` → `(U, Σ, V^T)` (thin, \(\sigma\) descending) or `None` if empty, non-finite, failed reserve, or not converged in `SVD_ITER_MAX` sweeps per value. `eigen_decomp(p, rm)` → `(Λ, V)` for real symmetric \(A\) (\(\lambda\) descending) or `None` if non-square, non-symmetric, empty, non-finite, or not converged in `EIGEN_ITER_MAX` sweeps per value. `fft`/`ifft(p, rm, cc)` — radix-2 Cooley–Tukey; `(1,n)` or `(n,1)` real, `(2,n)` complex; unnormalized forward; `ifft` divides by `n`; `n` a power of two ≤ `FFT_MAX_POINTS`. `random_fill(shape, dist, p, rm, cc)` (`random` feature / tests) samples `RandomDist::Uniform` / `Normal` / `Exponential`. |
| `ExactRational` | Exact `num/den` with integer-valued `ExactNum` parts, reduced to lowest terms (`den > 0`). `new` / `from_i64` / `from_ints`; `add`/`sub`/`mul`/`div`; `to_exact_num(p, rm)`; `is_integer`; `floor`/`ceil`/`round`; `partial_cmp` by cross-multiply; `parse_exact` / `format_exact`. Zero `den` is `NaN`. Not a float. `0.1` is `1/10`. |
| `ExactInt` | Signed limb integer (little-endian `Word`s). `from_i64`/`from_u64`/`from_i128`/`from_u128`; `add`/`sub`/`mul`; `div_rem` (truncated, zero divisor `None`); `gcd`; `pow`; `to_exact_num` / `from_exact_num`; `bit_length`. Not a truncated `ExactNum`. |
| `ExactNumPoly` | Dense univariate, coefficients lowest degree first. `from_coeffs` / `from_i64_coeffs`; `eval` (Horner/`polyval`); `add`/`sub`/`mul`; `div_rem` (zero divisor `None`); `gcd` (Euclidean, integer content, monic); `compose`; `derivative`; `integral` (constant 0); `companion_matrix`; `roots_real` through `POLY_COMPANION_CLOSED_DEG=2`. |
| Chebyshev | `chebyshev_coeffs(f, n, a, b, p, rm, cc)` interpolates at `n` Chebyshev–Gauss nodes (`n ≤ CHEBYSHEV_MAX_DEGREE`). `chebyshev_eval` maps `[a,b]` then `clenshaw`. `chebyshev_error_bound` is the ℓ¹ tail `Σ_{k≥1}|c_k|`. |
| Orthogonal polys | `ExactNum::{hermite_he,hermite_h,laguerre,gen_laguerre,chebyshev_t,chebyshev_u,gegenbauer}(n, …)` via three-term recurrences. `n ≤ ORTHOPOLY_N_MAX`. Standard Gegenbauer `C_n^{(λ)}` (`C_n^{(1/2)}=P_n`, `C_n^{(1)}=U_n`). |
| Quadrature | `gauss_legendre(f,a,b,n,…)`; `tanh_sinh(f,a,b,…)`; `gauss_laguerre(f,n,…)`; `gauss_hermite(f,n,…)`. `n ≤ QUADRATURE_MAX_NODES`. Tanh–sinh uses `h=2π/(p ln 2)` and up to `TANH_SINH_LEVELS_MAX` halvings. |
| Root finding | `bisect`/`newton`/`brent`/`illinois` on `ExactNum`. Opposite signs required for bracket methods (`None` otherwise). `ROOT_MAX_ITER=256`; `root_default_tol` is `2^{ROOT_DEFAULT_TOL}`. |
| ODE | `rk4`/`rk45_adaptive`/`euler` return `(t, y)` row `ExactNumArray`s. Dormand–Prince 5(4) with `atol`/`rtol`; `h_min=2^{ODE_MIN_STEP}`; cap `ODE_MAX_STEPS`. |
| DSP | Type-II `dct`/`dst` via a `2N` FFT; type-III inverses scaled so the round-trip is the identity. `fft_real`/`ifft_real` wrap the radix-2 DFT. Symmetric Hann / Hamming / Blackman / Kaiser / rectangular windows. Real length a power of two ≤ `DSP_MAX_POINTS` for transforms; windows allow any `n` in `1..=DSP_MAX_POINTS`. |
| Modular `ExactInt` | `mod_pow` / `mod_inv` / `miller_rabin` / Brent `pollard_rho`. Zero modulus and non-units are `None`. Cap `POLLARD_RHO_ITER_MAX`. |
| Hash | FIPS 180-4 `sha256`/`sha512`; HMAC-SHA-256; `constant_time_eq` always scans both slices. |
| Binary I/O | `ExactNum` / `ExactNumArray`: `to_bytes` / `from_bytes` / `write_bytes`. Inline 16-byte big-endian record when the mantissa is ≤ `BINARY_INLINE_MANT_BITS` (`to_inline_bytes` / `write_inline_bytes` / `InlineBinaryBuffer`). Wider values use a heap record of `u32` limbs. Invalid input is `Err`, not a panic. |
| CSV | `Ieee64Array` cells are unsigned binary64 bit patterns (empty / `nan` → `NAN`). `ExactNumArray` cells are `Display@p=` (`std`). `to_csv` / `from_csv` take a path; `to_csv_string` / `from_csv_str` are in-memory. Caps: `CSV_MAX_ROWS`, `CSV_MAX_COLS`. |
| HDF5 | Feature `hdf5`: `Ieee64Array` / `Ieee32Array` / `ExactNumArray` `to_hdf5` / `from_hdf5` (`path`, `dataset`). Native binary64 / binary32 bit lanes; `ExactNum` cells are opaque (`u32` BE length + `to_bytes`, padded). `Ieee64Array::append_hdf5` grows the first axis. `write_hdf5` / `read_hdf5` take an open `Hdf5File`. Re-exports `Hdf5File`, `HDF5DType`, `HDF5Error`. Wrong dataset name is `Err`. CSV is unchanged. |

These arrays are not NumPy-fast. `matmul` is a sequential triple loop (IEEE or `ExactNum` mul-then-add per term). No BLAS. Integer SIMD is not an FPU; a SIMD bit pattern that differs from the scalar kernel is a bug. `expr!` stays scalar.

---

## 8. `ExactNum` — constructors and classification

| Method / item | Signature / notes |
| --- | --- |
| `new` | `new(p: usize)` — zero at precision `p` |
| `nan` | `nan(err: Option<Error>)` |
| `from_word` | `from_word(d: Word, p)` |
| `from_i8`, `from_i16`, `from_i32`, `from_i64`, `from_i128` | signed integers + precision |
| `from_u8`, `from_u16`, `from_u32`, `from_u64`, `from_u128` | unsigned integers + precision |
| `From` | `i8`/`u8`/`i16`/`u16`/`i32`/`u32`/`i64`/`u64`/`i128`/`u128` at default 128-bit precision |
| `Default` | `new` at default 128-bit precision |
| `from_raw_parts` | mantissa slice, used length, sign, exponent, inexact flag |
| `from_words` | mantissa slice, sign, exponent |
| `max_value` / `min_value` | largest / most-negative finite at precision `p` |
| `min_positive` / `min_positive_normal` | smallest positive (incl. subnormal) / smallest normal |
| `random_normal` | feature `random`: random finite with exponent in `[exp_from, exp_to]` (not a Gaussian) |
| `random_uniform` / `random_gaussian` / `random_exponential` | feature `random` / tests: `[a,b]`; Box–Muller \(N(\mu,\sigma^2)\); Exp(\(\lambda\)) |
| `is_inf_pos` / `is_inf_neg` / `is_inf` | infinities |
| `is_nan` | |
| `is_int` | integer-valued finite |
| `is_zero` / `is_positive` / `is_negative` / `is_subnormal` | |
| `err` | `Option<Error>` on NaN |
| `classify` | `core::num::FpCategory` (`Zero`, `Subnormal`, `Normal`, `Infinite`, `Nan`) |
| `inexact` / `set_inexact` | exactness flag |
| `Clone`, `Debug` | |
| `PartialEq`, `Eq` | equality via `cmp`; NaN is not equal |
| `PartialOrd` | `None` if either is NaN |
| `Neg` | `ExactNum` and `&ExactNum` |

No `Hash`, no total `Ord`, no `AddAssign` / `MulAssign`.

---

## 9. `ExactNum` — arithmetic (explicit precision)

Unless noted, results round with `(p, rm)`.

| Method | Meaning |
| --- | --- |
| `add` / `sub` / `mul` / `div` | four operations |
| `add_full_prec` / `sub_full_prec` / `mul_full_prec` | no precision reduction (exact product/sum of finite values when it fits) |
| `fma` / `mul_add` | fused `a*b+c` with one final round; `mul_add` is an alias of `fma` |
| `two_sum` / `two_product` | `(hi, lo)`; `p`/`rm` round **hi** only; finite `hi + lo` equals the exact limb sum/product |
| `fused_sum` / `fused_dot` | extra-precision accumulation, one final round |
| `polyval` | Horner; `coeffs[0]` is the constant term (lowest degree first) |
| `rem` | remainder (`self` rem `d2`); no `p`/`rm` |
| `reciprocal` | `1/self` |
| `neg` | copy with inverted sign (also `Neg` / `inv_sign` in place) |
| `pow` | `self^n` (arbitrary `ExactNum` exponent); needs `Consts` |
| `powi` | `self^n` for `usize` `n` |
| `powsi` | `self^n` for `isize` `n` |

Operator traits **`Add` `Sub` `Mul` `Div`** for all combinations of `ExactNum` and `&ExactNum`, always at **128 bits, `ToEven`**. Use the named methods for other precision/rounding.

There is **no** `%` operator trait; remainder is `rem` or `expr!` `%`.

---

## 10. `ExactNum` — comparison, clamp, sign

| Method | Meaning |
| --- | --- |
| `cmp` | `Option<SignedWord>` — negative / zero / positive; `None` if unordered (NaN). `SignedWord` is `i128` on 64-bit targets and `i64` on 32-bit; it is the return type of `cmp`/`abs_cmp` even though it is not re-exported at the crate root. |
| `abs_cmp` | compare absolute values |
| `min` / `max` | |
| `clamp` | |
| `signum` | −1, 0, or +1 as `ExactNum` |
| `abs` | absolute value |
| `copysign` | magnitude of `self`, sign of `sign`; sign applied before precision reduction |
| `next_after` | next representable toward `toward` at precision `p` |

---

## 11. `ExactNum` — exponent, mantissa, precision

| Method | Meaning |
| --- | --- |
| `exponent` | `Option<Exponent>` |
| `precision` | `Option<usize>` mantissa bit length |
| `sign` | `Option<Sign>` |
| `as_raw_parts` | `Option<(&[Word], usize, Sign, Exponent, bool)>` |
| `mantissa_digits` | `Option<&[Word]>` |
| `mantissa_max_bit_len` | |
| `is_inline` | mantissa fits in `INLINE_WORDS` |
| `set_exponent` | |
| `set_sign` | |
| `set_precision` | round to `p`; `Result<(), Error>` |
| `try_set_precision` | Ziv-style: round if uniquely determined given working precision `s`; `bool` |
| `frexp` | `(significand in [0.5, 1), exponent)` as `(ExactNum, Exponent)` |
| `ldexp` / `scalb` | `self × 2^n` (`n: Exponent`) |
| `logb` | `floor(log2(|self|))` as a float |
| `ilogb` | same as integer `Option<Exponent>` |

`frexp` and `ilogb` are **methods only** (not `expr!` leaves).

---

## 12. `ExactNum` — integer / rounding of the value

| Method | Meaning |
| --- | --- |
| `int` | integer part |
| `fract` | fractional part |
| `ceil` / `floor` | |
| `round` | round with `n` binary fractional bits and mode `rm` |

---

## 13. `ExactNum` — roots, logs, exponentials

Need `Consts` where marked `cc`.

| Method | `expr!` leaf? |
| --- | --- |
| `sqrt(p, rm)` | `sqrt(x)` |
| `cbrt(p, rm)` | `cbrt(x)` |
| `nth_root(n, p, rm)` | `root(x, n)` — `n=0` → NaN; `n=2`/`3` delegate to sqrt/cbrt; even `n` of negative → NaN |
| `ln` / `log2` / `log10` | yes |
| `log(self, base, p, rm, cc)` | `log(x, b)` |
| `log1p` | `log1p(x)` |
| `exp` / `exp2` / `exp10` | yes |
| `expm1` | yes |
| `pow` | `pow(b, x)` |

---

## 14. `ExactNum` — circular and hyperbolic

All take `(p, rm, cc)` except `hypot` (no cache).

| Method | `expr!` |
| --- | --- |
| `sin` / `cos` / `tan` | yes |
| `sin_cos` → `(sin, cos)` | **no** (shared argument reduction) |
| `asin` / `acos` / `atan` | yes |
| `atan2(y=self, x, …)` | `atan2(y, x)` |
| `hypot` | `hypot(x, y)` |
| `rem_pi` | reduce modulo `2π` into `(-2π, 2π)` |
| `sinh` / `cosh` / `tanh` | yes |
| `sinh_cosh` → `(sinh, cosh)` | **no** |
| `asinh` / `acosh` / `atanh` | yes |

---

## 15. `ExactNum` — special functions

Each method’s rustdoc has a **`# Precision`** section (algorithm, region, named thresholds, ULP / Ziv, MPFR oracle).

| Method | Notes | `expr!` |
| --- | --- | --- |
| `erf` / `erfc` | `erfc = 1 − erf` | yes |
| `gamma` | poles at non-positive integers → NaN; +Inf at 0 | yes |
| `ln_gamma` | log-gamma for positive `self` | yes |
| `digamma` | \(\psi(z)\); reflection for \(z<0\); poles at non-positive integers | yes |
| `gammainc(s, x)` | lower \(\gamma(s,x)\); \(s>0\), \(x\ge 0\) | `gammainc(s, x)` |
| `gammainc_upper(s, x)` | upper \(\Gamma(s,x)\) | `gammainc_upper(s, x)` |
| `ei` | Cauchy PV for \(x<0\); \(x=0\) is a pole | yes |
| `si` | odd; all real | yes |
| `ci` | `self > 0`; otherwise NaN | yes |
| `li` | `self > 0`, `self ≠ 1`; `Ei(ln self)` | yes |
| `fresnel_s` / `fresnel_c` | odd; series or auxiliary \(f,g\) | yes |
| `ai` / `bi` / `ai_prime` / `bi_prime` | Airy; series for \(\lvert x\rvert<\texttt{AIRY\_SERIES\_THRESHOLD}\); asymptotic otherwise | `ai`, `bi` |
| `bessel_j(n, p, rm, cc)` | `J_n(self)`, integer order; Miller for large \(n\) | `bessel_j(x, n)` |
| `bessel_j_nu` / `bessel_y` / `bessel_i` / `bessel_k` | Real order; \(K\): \(x>0\) | yes |
| `elliptic_k` / `elliptic_e_complete` / `elliptic_f` / `elliptic_e` / `elliptic_pi_complete` / `elliptic_pi` | Carlson; \(m=k^2\), \(x=\sin\varphi\); \(K(1)=+\infty\); \(K(m>1)=m^{-1/2}K(1/m)\) | `elliptic_k`, `elliptic_e`, `elliptic_f`, `elliptic_e_inc`, `elliptic_pi`, `elliptic_pi_inc` |
| `legendre_p` / `assoc_legendre_p` | Integer \(n\); \(p+O(n)\) recurrence; Condon–Shortley | `legendre_p(x, n)`, `legendre_p_assoc(x, n, m)` |
| `hypergeom_2f1` | Series / Gauss / Pfaff; real \(z\le -1\) when defined; non-real \(z>1\) → NaN | `hypergeom_2f1(a,b,c,z)` |
| `betainc` | Regularized \(I_x(a,b)\) | `betainc(a,b,x)` |
| `normal_pdf` / `normal_cdf` | \(\sigma>0\); CDF via `erf` | `normal_pdf(x,μ,σ)` |
| `gamma_pdf` / `beta_pdf` | Scale \(\beta\); \(B\) via \(\Gamma\) | yes |
| `poisson_pmf` / `binomial_pmf` | Non-negative integer \(k\) | yes |
| `chi_squared_cdf` / `student_t_pdf` | Regularized \(\gamma(k/2,x/2)/\Gamma(k/2)\); \(t\) via \(\Gamma\) | yes |

---

## 16. `ExactNum` — parse, format, radix conversion

| API | Notes |
| --- | --- |
| `parse(s, rdx, p, rm, cc)` | returns `ExactNum` (not `Result`); Inf / NaN / `err()` on failure; scientific `1.5e3` |
| `parse_exact(s)` | dyadic decimals only (`0.5`, `0.125`); `0.1` is `None` — use `ExactRational::parse_exact` |
| `format_exact(rdx)` | min digits for a terminating expansion; `None` if the radix cannot represent `self` exactly |
| `format(rdx, rm, cc)` | `Result<String, Error>`; Inf / −Inf / NaN / Err strings |
| `with_radix(self, radix)` | wrap as `RadixFloat` |
| `convert_from_radix(sign, digits, e, rdx, p, rm, cc)` | digit bytes, scientific exponent in that radix |
| `convert_to_radix(rdx, rm, cc)` | `(Sign, Vec<u8>, Exponent)` |
| `FromStr` (`std`) | decimal parse, `ToEven`, unbounded precision then value |
| `Display` (`std`) | decimal |
| `LowerExp` / `UpperExp` (`std`) | decimal scientific (`e` / `E`) |
| `Binary` / `Octal` (`std`) | |
| `UpperHex` / `LowerHex` (`std`) | hex; specials stay `Inf` / `NaN` |

For radices where `e` is a digit, the exponent is written with `_e`.

---

## 17. Macros

Public macros (crate root): `expr!`, `cexpr!`, `exact!`, `fbig!`. Import them like any other item, e.g. `use zenith_float::{expr, cexpr, exact};`, or call them as `zenith_float::exact!("1.25")`.

### `expr!(expression, context)`

**Context** (`Contextable`): `Context`, or tuples

- `(precision, RoundingMode, &mut Consts)` — exponent limits `EXPONENT_MIN`/`MAX`
- `(precision, RoundingMode, &mut Consts, emin, emax)` — `emin` clamped to `[EXPONENT_MIN, 0]`, `emax` to `[0, EXPONENT_MAX]`

**Syntax allowed:** paths (variables), integer / float / string literals, unary `-`, binary `+ − * / %`, calls, parentheses. Inputs are treated as exact (`set_inexact(false)`). Working precision is raised; final result uses the context rounding mode. Overflow of `emin`/`emax` → `±Inf`.

**Binary operators:** `+`, `-`, `*`, `/`, `%` (maps to `ExactNum::rem`).

**Function leaves (complete list):**

`recip`, `sqrt`, `cbrt`, `root`, `ln`, `log2`, `log10`, `log`, `log1p`, `exp`, `exp2`, `exp10`, `expm1`, `pow`, `rem_pi`, `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `atan2`, `hypot`, `fma`, `mul_add`, `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh`, `erf`, `erfc`, `gamma`, `ln_gamma`, `digamma`, `gammainc`, `gammainc_upper`, `ei`, `si`, `ci`, `li`, `fresnel_s`, `fresnel_c`, `ai`, `bi`, `bessel_j`, `bessel_j_nu`, `bessel_y`, `bessel_i`, `bessel_k`, `elliptic_k`, `elliptic_e`, `elliptic_e_inc`, `elliptic_f`, `elliptic_pi`, `elliptic_pi_inc`, `legendre_p`, `legendre_p_assoc`, `hypergeom_2f1`, `betainc`, `normal_pdf`, `normal_cdf`, `gamma_pdf`, `beta_pdf`, `poisson_pmf`, `binomial_pmf`, `chi_squared_cdf`, `student_t_pdf`, `ldexp`, `scalb`, `logb`.

**Named constants in the expression:** `pi`, `e`, `ln_2`, `ln_10`, `sqrt2`, `phi`, `euler_gamma`.

### `cexpr!(expression, context)`

Same **context** as `expr!` (`Context` or the two tuples above). Same extra-precision loop: working precision starts at `p + WORD_BIT_SIZE`, internals use `RoundingMode::None`, the root is rounded once with `(p, rm)`, then each part is clamped with `check_complex_exponent_range`.

**`I`:** expands to `ExactComplex::i(p_wrk)` (`0 + 1i`). Lowercase `i` remains a variable. Do not name a variable `e` (`e` is Euler’s number, lifted as `e + 0i`).

**Cancellation:** add/sub use **two** `errs[]` slots, one for the real part and one for the imaginary part (`complex_cancel_bits` returns `(Option<usize>, Option<usize>)`). A subtraction can lose bits on both parts at different rates; a shared max would underestimate working precision. `p_wrk` is `p_rnd + sum(errs)`.

**Operators:** `+`, `-`, `*`, `/`. No `%`.

**Function leaves (complete list):**

`recip`, `sqrt`, `cbrt`, `root`, `ln`, `log2`, `log10`, `log`, `log1p`, `exp`, `exp2`, `exp10`, `expm1`, `pow`, `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `hypot`, `fma`, `mul_add`, `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh`, `abs`, `arg`, `conj`, `ldexp`, `scalb`, `logb`, `erf`, `erfc`, `gamma`, `ln_gamma`, `digamma`, `ei`, `si`, `ci`, `li`, `fresnel_s`, `fresnel_c`, `ai`, `bi`, `bessel_j_nu`, `bessel_y`, `bessel_i`, `bessel_k`, `elliptic_k`, `elliptic_e`, `elliptic_e_inc`, `elliptic_f`, `elliptic_pi`, `elliptic_pi_inc`, `hypergeom_2f1`.

**Named constants:** `I`, `pi`, `e`, `ln_2`, `ln_10`, `sqrt2`, `phi`, `euler_gamma` (reals except `I` are `x + 0i`).

**Reals lifted as `x + 0i`:** `abs`, `arg`, `logb`, literals, and the real constants. `ldexp` / `scalb` scale **both** parts by `2^n`. `hypot(z, w)` is principal `sqrt(z² + w²)`. `fma` / `mul_add` are extra-precision `a*b+c`, then one round. `root(z, n)` is `nth_root` (integer `n`).

**Branch cuts (principal values, observable):**

| Family | Cut / range |
| --- | --- |
| `arg` | `atan2(im, re)`, values in (−π, π] |
| `ln`, `log2`, `log10`, `log`, `log1p` | cut on (−∞, 0]; `ln(−1) = iπ` |
| `sqrt`, `cbrt`, `root` | cut on (−∞, 0]; `Re(sqrt) ≥ 0`; `sqrt(−1) = +i` |
| `pow` | `exp(w · ln(z))`, so the `ln` cut on the **base** |
| `asin` / `acos` / `atan` / `asinh` / `acosh` / `atanh` | principal branches of the usual identities |
| `ei` / `si` / `ci` | `Ei` cut on (−∞, 0]; `Si`/`Ci` via `Ei(±iz)`; `Ci(0)` → NaN |
| `li` | `Ei(ln z)`; cut on (−∞, 1]; pole at 1 → NaN |
| `fresnel_s` / `fresnel_c` | via `erf`; entire |
| `bessel_j_nu` / `bessel_i` | entire for integer ν; cut on (−∞, 0] otherwise; \(z=0\) non-integer ν → NaN |
| `bessel_y` / `bessel_k` | cut on (−∞, 0]; \(z=0\) → NaN |
| `elliptic_k` | cut on \([1,+\infty)\) in \(m\); \(K(1)=+\infty\); \(K(2+\varepsilon i)\) and \(K(2-\varepsilon i)\) are conjugates |
| `hypergeom_2f1` | cut on \([1,+\infty)\) in \(z\); \({}_2F_1(1,1;2;2+\varepsilon i)\) and \({}_2F_1(1,1;2;2-\varepsilon i)\) differ |
| incomplete \(F,E,\Pi\) | Carlson cuts when \(1-x^2\), \(1-mx^2\), or \(1-nx^2\) lies on \((-\infty,0]\) |

**Trig / hyperbolic:** `sin(x+iy) = sin(x)cosh(y) + i cos(x)sinh(y)` (and the matching identities). The **complex** argument is never passed to `rem_pi`. Only a real component uses real `sin_cos` / `sinh_cosh` (those may reduce that real).

**Explicitly not in `cexpr!`:**

| Leaf | Why |
| --- | --- |
| `%` / `rem_pi` | remainder and π-reduction are real; complex trig uses the identities above |
| `atan2` | no standard two-complex analogue; use `arg` for `atan2(im, re)` |

Literals and real constants are lifted as `x + 0i`. Variables may be `ExactComplex` or anything `FromExt` can wrap as a real.

**Not in `expr!`:** `frexp`, `ilogb`, `sin_cos`, `sinh_cosh`, `int`/`fract`/`ceil`/`floor`/`round`, `min`/`max`/`clamp`, `cmp`, `copysign`, `next_after`, `powi`/`powsi`, `abs`/`signum`, parse/format, raw parts, `nth_root` under the name `nth_root` (use `root`). Complex values use `cexpr!`, not `expr!`.

### `exact!("…")` and `fbig!("…")`

Compile-time parse of a **string literal** into an exact `ExactNum` (same expansion). Decimal literal syntax as implemented by the macro parser.

---

## 18. Constants cache

**`Consts::new() -> Result<Self, Error>`** — allocate progressive caches.

| Method | Constant |
| --- | --- |
| `pi(p, rm)` | π |
| `e(p, rm)` | e |
| `ln_2(p, rm)` | ln 2 |
| `ln_10(p, rm)` | ln 10 |
| `sqrt2(p, rm)` | √2 |
| `phi(p, rm)` | φ = (1+√5)/2 |
| `euler_gamma(p, rm)` | Euler–Mascheroni γ |
| `cache_info()` | `ConstCacheInfo` { `pi`, `e`, `ln2`, `ln10`, `sqrt2`, `phi`, `euler` } bit lengths |

**`CachedFBig`:** wrap a precomputed `ExactNum`; `new`, `cached_bit_len`, `inner`, `round(p, rm)`.

**`SharedConsts` (`std`):** mutex around `Consts`. `new()`, `with(|cc: &mut Consts| …)` (recovers poisoned mutex).

---

## 19. `Context` and `Contextable`

**`zenith_float::ctx::Context`:** holds precision, rounding mode, `Consts`, `emin`, `emax`.

| Method | |
| --- | --- |
| `new(p, rm, cc, emin, emax)` | |
| `to_raw_parts` | `(p, rm, Consts, emin, emax)` |
| `with_rounding_mode` | run a closure with a temporary rounding mode, then restore |
| `set_precision` / `set_rounding_mode` / `set_consts` / `set_emin` / `set_emax` | |
| `precision` / `rounding_mode` / `consts` / `emin` / `emax` | |
| `const_pi` / `const_e` / `const_ln2` / `const_ln10` / `const_sqrt2` / `const_phi` / `const_euler_gamma` | at context `p`/`rm` |
| `clone` | `Result<Self, Error>` (deep-clone cache) |

**`Contextable`:** `precision`, `rounding_mode`, `consts`, the seven `const_*` methods, `emin`, `emax`. Implemented for `Context` and the two tuples in §17.

---

## 20. `ExactComplex`

Cartesian `re + i·im` as two `ExactNum`s.

| API | |
| --- | --- |
| `new(re, im)` | |
| `re` / `im` | references |
| `zero(p)` / `one(p)` / `i(p)` | |
| `is_nan` | either part NaN |
| `conj` | |
| `abs(p, rm)` | modulus (`hypot`) |
| `arg(p, rm, cc)` | argument (`atan2`) |
| `add` / `sub` / `mul` / `div` | with `(p, rm)` |
| `exp` / `ln` / `sin` / `cos` / `tan` | with `(p, rm, cc)` |
| `sinh` / `cosh` / `tanh` | |
| `sqrt` / `pow` | principal branch |
| `asin` / `acos` / `atan` / `asinh` / `acosh` / `atanh` | principal branches |
| `log2` / `log10` / `log` / `log1p` | principal; `log2`/`log10`/`log` via `ln` ratios |
| `exp2` / `exp10` / `expm1` | via `exp` |
| `ldexp` / `scalb` | scale both parts by `2^n` |
| `cbrt` / `nth_root` | principal; `exp(ln / n)` |
| `hypot` | principal `sqrt(z² + w²)` |
| `fma` / `mul_add` | extra-precision `a*b+c` |
| `logb` | `logb(\|z\|)` as a real |
| `Add` `Sub` `Mul` `Div` | 128-bit `ToEven` like reals |
| `erf` / `erfc` | Faddeeva \(w(z)\); entire; NaN in → NaN out |
| `gamma` / `ln_gamma` / `digamma` | Stirling + reflection; \(\ln\Gamma\) cut on \((-\infty,0]\); poles → NaN |
| `ei` / `si` / `ci` / `li` | series or asymptotic `Ei`; `Si`/`Ci` via `Ei(±iz)`; `li=Ei(ln z)` |
| `fresnel_s` / `fresnel_c` | via `erf`; entire |
| `ai` / `bi` | series for small \(\lvert z\rvert\); asymptotic or \(\omega\)-connection for large \(\lvert z\rvert\) |
| `bessel_j_nu` / `bessel_y` / `bessel_i` / `bessel_k` | series or Hankel; \(I_ν=i^{-ν}J_ν(iz)\); \(K_ν=(\pi/2)i^{ν+1}H_ν^{(1)}(iz)\) |
| `elliptic_k` / `elliptic_e_complete` / `elliptic_f` / `elliptic_e` / `elliptic_pi_complete` / `elliptic_pi` | Carlson in \(\mathbb{C}\); \(m=k^2\), \(x=\sin\varphi\); \(K(1)=+\infty\); \(K\) cut on \([1,+\infty)\) |
| `hypergeom_2f1` | Series / Euler / Pfaff / Kummer; cut on \([1,+\infty)\) in \(z\); non-positive integer \(c\) → NaN |
| Other specials (`_2F1`, …) | not yet; backlog (software limbs, principal branches) |

No `expr!` for complexes — use `cexpr!`. Serde: struct `{ "re", "im" }` of decimal strings.

---

## 21. `Ball` and `ziv_round`

**`Ball`:** enclosure `mid ± rad` (`rad` stored non-negative).

| Method | |
| --- | --- |
| `new(mid, rad)` | |
| `mid` / `rad` | |
| `add` / `mul` | interval arithmetic with an extra rounding ulp in the radius |
| `exp` / `sin` / `cos` / `ln` / `sqrt` / `erf` / `bessel_j0` / `bessel_j1` | Lipschitz (or \(\lvert e^m\rvert(e^r-1)\)) plus `BALL_TRANSCENDENTAL_ERROR_TERMS` ulps |
| `contains(x, p)` | `x` in `[mid−rad, mid+rad]`; NaN/Inf never contained |

**`ComplexBall`:** disk `mid` (an `ExactComplex`) with radius `rad`.

| Method | |
| --- | --- |
| `new(mid, rad)` / `mid` / `rad` | |
| `add` / `mul` | disk arithmetic plus `BALL_TRANSCENDENTAL_ERROR_TERMS` slack |
| `exp` / `ln` / `sin` / `cos` | Lipschitz on the disk; `ln` requires the disk to exclude \(0\) |
| `contains(z, p)` | \(\lvert z-\mathrm{mid}\rvert\le\mathrm{rad}\) |

**`ziv_round(p, rm, compute)`:** call `compute(working_p)` and `try_set_precision` until the rounding is unique or `MAX_PREC_RETRY` is exhausted (then NaN / `PrecisionRetryExhausted`).

**`ziv_round_vec(p, rm, inputs, compute)`:** lift each input to `p_wrk`, then the same loop. `compute(p_wrk, &[ExactNum])`. Golds: `hypot(3,4)=5`; `atan2(1,1)=π/4` at 256 bits.

---

## 22. `RadixFloat`

Binary kernel plus a **parse/format radix** tag.

| Method | |
| --- | --- |
| `new(value, radix)` | validates 2..=36 |
| `with_radix` | no extra check (for named radices) |
| `value` / `radix` / `into_inner` | |
| `parse(s, radix, p, rm, cc)` | |
| `format(rm, cc)` | using the stored radix |
| `Deref` / `DerefMut` | to `ExactNum` |
| `ExactNum::with_radix` | constructor |

---

## 23. `FromExt`

`FromExt<T>::from_ext(v, p, rm, cc)`:

- For types with `From<T> for ExactNum`: convert then `set_precision`.
- For `&str`: `parse` in decimal.

Used by `expr!` to lift variables and literals.

---

## 24. Binary interchange

Network byte order. Inline record is 16 bytes: flag, version, `n_sig`, inexact, `i32` exponent, two `u32` limbs (least-significant first, matching kernel word order). Specials use flags `0x02`–`0x0A`. Mantissas wider than 64 bits use flag `0x10`/`0x11` plus `u32` limbs. Arrays use flag `0x20` and a shape/`p` header.

---

## 25. Serde (`serde` feature)

- **`ExactNum`:** `"<Display>@p=<bits>"`. JSON integers use `DEFAULT_P`. Rehydration parses at the stored bit count.
- **`ExactComplex`:** `re` / `im` encoded as above.
- **`ExactRational`:** `num` / `den` encoded strings.
- **`ExactInt`:** signed decimal.
- **`ExactNumArray`:** `shape`, `data`, `p`, `rm`. Shape product must equal `data.len()`.
- **`Ieee32Array` / `Ieee64Array`:** integer IEEE bit patterns (not hardware JSON floats).
- **`Ball`:** `mid` / `rad` encoded strings.

---

## 26. Random (`random` feature)

| Item | |
| --- | --- |
| `DEFAULT_RANDOM_SEED` | `0x5EED_CAFE_BADC_0D00` |
| `random_seed()` | current seed |
| `reseed_random(seed)` | identical sequences after the same seed |
| `seeded_random::<T>()` | draw from the crate RNG |
| `ExactNum::random_normal` | random normalized mantissa in an exponent range (testing) |
| `ExactNum::random_uniform` / `random_gaussian` / `random_exponential` | `[a,b]`; Box–Muller; inverse-CDF exponential |
| `RandomDist` / `ExactNumArray::random_fill` | `Uniform` / `Normal` / `Exponential` |
| `ZENITH_TEST_SEED` | env override for tests / replay |

Unit tests default to the deterministic seed. Without reseeding, non-test `random` feature uses OS entropy.

---

## 27. Implementation notes (not extra public functions)

The kernel uses integer add/mul (schoolbook, Karatsuba/Toom, FFT at large sizes), Newton division, series and argument reduction for elementary/special functions, and a progressive constant cache. None of those algorithms are separate public types.

MPFR/`rug` appear only in **tests** (`mpfr-tests`), not as the evaluation engine.

---

## 28. What is not in this crate

These are load-bearing product choices, not a backlog:

- Hardware binary interchange types or converters (including a feature-gated module). Callers or a separate crate pack bits from `frexp` / `ilogb`.
- Using `expr!` for complex values (that macro is real-valued). Complex expressions are `cexpr!`.
- Symbolic CAS, formula rewriting, or host-application IR
- Remainder as a Rust `%` operator on `ExactNum`
- `Hash` / total order including NaN
- Assigning operators (`+=` …)

Callers choose precision `p` and rounding; this crate supplies the mechanism, not per-domain “formula class” policies.