Skip to content

feat(float): implement arbitrary-precision trigonometric functions and constants - #60

Merged
cmpute merged 3 commits into
cmpute:masterfrom
CokieMiner:feat/trigonometry
Jun 8, 2026
Merged

feat(float): implement arbitrary-precision trigonometric functions and constants#60
cmpute merged 3 commits into
cmpute:masterfrom
CokieMiner:feat/trigonometry

Conversation

@CokieMiner

Copy link
Copy Markdown
Contributor

Overview

This PR introduces a comprehensive suite of trigonometric primitives and constants to the dashu-float crate. The implementation follows the library's architectural roadmap by consolidating advanced functions into a new math module and adopting a non-panicking FpResult pattern for robust error handling.

Key Features

  • Complete Trigonometric Suite: Implements sin, cos, tan, asin, acos, atan, atan2, and an optimized joint sin_cos evaluation.
  • *High-Performance $\pi: Implements the Chudnovsky algorithm with Binary Splitting, providing state-of-the-art performance for high-precision constant generation.
  • Non-Panicking API: All functions return FpResult<B> (wrapping Normal, NaN, Infinite, etc.), ensuring safety for domain errors and singularities without library crashes.
  • Mathematically Rigorous:
    • IEEE 754 Compliance: Full support for signed infinities and edge cases in atan2.
    • Numerical Stability: Magnitude-aware range reduction and dynamic guard-digit calculation ( + \lceil\log(x)\rceil + 50$) to prevent catastrophic cancellation.
    • Base-Agnostic Precision: Dynamic bit-mapping using EstimatedLog2 ensures correctness across all bases (binary, decimal, etc.).

Internal Improvements

  • Robust Rounding: Patched a bug in split_at_point_internal where intermediate results with high digit counts could trigger assertion failures for values $|x| &lt; 1.0$.
  • API Ergonomics: Exposed functions via both the Context API and inherent methods on FBig, with clear documentation on panic conditions and return types.

Verification Results

  • Fuzz Testing: Validated against the MPFR library with 1,500+ randomized iterations covering extreme precision ranges (up to 1,000 digits) and exponents (up to $\pm 500$).
  • Unit Tests: Comprehensive coverage in tests/trig.rs for all primitives and IEEE 754 infinity logic.
  • Standards: Verified clean by cargo clippy and maintained no_std compatibility by using pure integer and core-compliant arithmetic.

Performance Note

The Chudnovsky implementation leverages dashu's fast multiplication for (M(n)\log^2 n)$ complexity. A // TODO has been added for future static caching of $\pi$ to further optimize repeated calls at common precisions.

@cmpute

cmpute commented Apr 23, 2026

Copy link
Copy Markdown
Owner

Thanks for the PR! I will try to find a time to review and merge 👍

@CokieMiner

Copy link
Copy Markdown
Contributor Author

I tried to follow more or less what was there already if you don't agree with any design decision just ask for the changes

@CokieMiner

CokieMiner commented Apr 23, 2026

Copy link
Copy Markdown
Contributor Author

Just a note I forgot to leave there are still some panick paths that I left that in principle should never happen but if you want to discuss how to handle them I will patche to crate design preferences

@CokieMiner
CokieMiner force-pushed the feat/trigonometry branch 2 times, most recently from fb8c49a to adce916 Compare April 23, 2026 19:38
@CokieMiner

Copy link
Copy Markdown
Contributor Author

Update: Enhanced Precision Stability & Performance Optimizations

I've just pushed a series of updates to address some edge-case precision issues and optimize memory usage:

  • Precision Guarding: Refined the guard digit calculation in compute_work_context to scale with input magnitude ($50 + \text{mag}/10$). This ensures $100%$ accuracy during argument reduction even for massive exponents (e.g., $10^{2000}$), which previously suffered from catastrophic cancellation.
  • Tangent Accuracy: Fixed a 1-ulp precision leak in tan() by computed the intermediate sin/cos division within the high-precision work context before the final rounding.
  • Zero-Clone Optimization: Refactored the range reduction logic to be "zero-clone" by reordering operations to consume owned values only after their final use as references.
  • Memory Efficiency: Optimized the Chudnovsky algorithm base case in consts.rs to use reference-based multiplication, significantly reducing allocations during $\pi$ constant generation.
  • Aggressive Fuzzing: Hardened the test suite with a 2000-iteration random fuzzing pass against rug, covering extreme exponent ranges and Pythagorean identity checks.

The implementation is now significantly leaner and more robust against extreme scale inputs.

@cmpute
cmpute force-pushed the feat/trigonometry branch from adce916 to fe969f4 Compare May 31, 2026 14:32
@CokieMiner

CokieMiner commented May 31, 2026

Copy link
Copy Markdown
Contributor Author

Any attention needed here from my part?

@cmpute

cmpute commented May 31, 2026

Copy link
Copy Markdown
Owner

I still need sometime to review the whole PR, I do have some feedback from Deepseek, I do agree with its suggestions:

Issues and Suggestions

  1. rug added as a dev-dependency — potential concern

float/Cargo.toml:64 adds rug = "1.24" which depends on GMP/MPIR via gmp-mp-fr-sys. This is a native C dependency that:

  • Requires system libraries (gmp, mpfr, mpc) or builds from source
  • May fail to compile on some CI environments or platforms
  • Adds significant build time

Consider gating behind a feature flag (e.g., dev-utils) so it's only compiled when explicitly needed.

  1. Float arithmetic in pi() bit calculation

consts.rs:33: (self.precision as f64 * log2_b_ub as f64) as usize + 1 uses f64 to compute ceil(precision * log2(B)). This could be off by 1 for large
precision values where f64 loses integer precision (above 2^53). For typical use this is fine, but since this is an arbitrary-precision library, consider
integer ceiling division instead.

  1. tan missing assert_limited_precision check

trig.rs:273-289 — tan delegates to sin_cos which calls compute_work_context using self.precision. If self.precision == 0, compute_work_context produces
saturating_add(0 + ...) which could be wrong. Other functions like sin, cos, atan call assert_limited_precision but tan does not. Should add it for
consistency.

  1. Duplicated range reduction code

sin, cos, and sin_cos each contain ~20 nearly identical lines for:

  1. compute_work_context → round input
  2. Compute π, half_pi
  3. Divide by half_pi, round, compute remainder
  4. Extract k, k % 4

This should be extracted into a helper like reduce_to_quadrant(x) -> (FBig, i8) to avoid the duplication and risk of divergence.

  1. atan2 infinite case — deeply nested branching

trig.rs:782-823 — The infinity handling in atan2 has 4 levels of nesting. Consider a lookup table or early-return pattern:
let (y_sign, x_sign) = (y.sign(), x.sign());
let (sy, sx) = (y_sign == Sign::Positive, x_sign == Sign::Positive);
let res = match (y.is_infinite(), x.is_infinite(), sy, sx) { ... };

  1. FBig::from_repr usage in test_pythagorean_identity_fuzz

trig_random.rs:330-331 — Uses FBig::from_repr(s_r.value(), dashu_ctx). If from_repr is not a public API, this test may break. Verify this is intentionally
available.

  1. sin_internal convergence threshold

trig.rs:524: let threshold = sum.sub_ulp() is computed from the initial sum (x), not from the running sum. For x near 0, this threshold is very small — good.
But for x close to π/4 where sin(x) ≈ 0.7, the threshold from sub_ulp() on the initial value may be too tight or too loose depending on the precision. The
Taylor series should converge rapidly regardless, so this is a minor concern.

  1. Missing blank line

trig.rs:700 — Missing blank line before /// Calculate the arctangent doc comment after the closing brace of acos.

  1. compute_work_context and assert_limited_precision ordering

sin and cos call assert_limited_precision(self.precision) after the infinite check but before compute_work_context. If self.precision == 0, the function
panics. But atan calls assert_limited_precision after the infinite check AND the zero check. This is fine since infinite inputs return early, but the
inconsistency in ordering is worth noting.

  1. No #[cfg(test)] guard on test utility

trig_random.rs:1089 — mod helper_macros; is included in the test file but the random_dbig function is not #[cfg(test)] gated (it's in a test file so this is
automatically test-only, but worth confirming).

@CokieMiner

Copy link
Copy Markdown
Contributor Author

1 - So maybe gatting the fuzz tests behind some kind of feature, I never used features for dev dependencies and I don't think it is a good idea to create features users are never gona use do you recommend trying to add similar fuzz and property based testing similar to the Pythagorean identity one? there are some I can try to add with some simple expression manipulation on trig identities but i think fuzzing against a well established lib is always better as this can hide some symmetric bugs on the algorithms.

2 - ok fair I had to keep in attention that in another project I'm doing with rug and agree, didn't think it was a problem at the time.

3 - ok fair

4 - seems not that much code but fair also

5 - that one I need to see the code

6 - most of these seem mostly irrelevant

I can try fixing them Tuesday as I need to study for eletromagnetism 1 and modern physics finals, any other finding let me know.

@CokieMiner

Copy link
Copy Markdown
Contributor Author

Done

@CokieMiner

Copy link
Copy Markdown
Contributor Author

Mb let those two slip in the middle of the other from other parts of the codebase

@cmpute

cmpute commented Jun 6, 2026

Copy link
Copy Markdown
Owner

Thanks for your help on editing. First, regarding the fuzz test, it seems that the fuzz test takes a lot of CI time. So my suggestion would be split the fuzz tests into a separate crate (just like why I have the benchmark crate in a separate folder as a stand alone crate. Maybe create a 'fuzz' folder and put all fuzzy related code in there.

Besides, I have run AI code review again and here's its feedback. Please review (you can directly feed them into your AI):

[
    {
      "file": "float/src/math/trig.rs",
      "line": 440,
      "summary": "atan2 calls self.atan() (original precision) instead of work_context.atan(), wasting the 50 guard digits it allocated — inconsistent with asin which correctly uses
  work_context.atan()",
      "failure_scenario": "atan2(3, -4) at precision 20: work_context has precision 70, y_f/x_f is computed at 70 digits, but self.atan() re-rounds to precision 20 before computing
  the arctangent, losing the extra precision"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 401,
      "summary": "atan2 is missing assert_limited_precision(self.precision), unlike sin/cos/tan/atan — with precision=0, infinity inputs succeed but finite inputs panic inside
  self.atan()",
      "failure_scenario": "Context::new(0).atan2(infinity, one) returns pi/4, but Context::new(0).atan2(one, one) panics inside atan at assert_limited_precision — inconsistent API
  behavior"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 262,
      "summary": "asin and acos are also missing assert_limited_precision(self.precision), same defect as atan2 — they silently accept precision=0 instead of panicking like their
  sibling functions",
      "failure_scenario": "Context::new(0).asin(x) silently computes with work_precision=50 and returns a result, while Context::new(0).sin(x) panics — inconsistent contract across
  trig functions"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 244,
      "summary": "tan calls compute_work_context(x) to inflate precision, then passes that work_context to sin_cos which calls compute_work_context again — roughly doubling the work
  precision and making the Chudnovsky pi computation much more expensive than needed",
      "failure_scenario": "tan(1e1000) at precision 20: tan's work_context has precision ~1150, sin_cos's internal compute_work_context inflates it again to ~2300, requiring pi
  computed to 2300 digits instead of 1150"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 312,
      "summary": "acos → asin → atan call chain stacks 3 layers of +50 guard digits (150 total), because asin calls the public atan() which re-validates and re-creates its own work
  context instead of calling atan_internal directly",
      "failure_scenario": "acos(0.5) at precision P: acos creates P+50, asin creates P+100, atan creates P+150 — deepest series computation runs at 150 extra digits instead of the
  intended 50-100"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 412,
      "summary": "atan2 unconditionally computes pi and half_pi via the expensive Chudnovsky algorithm before checking the infinity match arms — for (finite_y, +inf_x) the result is
  ZERO and pi is entirely wasted",
      "failure_scenario": "atan2(1, infinity) computes pi to 50+ digits only to return ZERO — the Chudnovsky binary splitting is O(n log²n) and dominates the function cost"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 375,
      "summary": "atan_internal Euler series has only linear convergence (ratio → 1/2 when x≈1 after range reduction), requiring O(P) iterations for P digits — much slower than
  sin/cos's O(√P) quadratic convergence",
      "failure_scenario": "atan(1) at precision 10000 in base 10: ~10000 iterations of multi-precision arithmetic, compared to sin/cos which need only ~100 iterations for similar
  precision"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 19,
      "summary": "x.exponent + x.digits() as isize can overflow isize for extreme inputs — saturating_add would be safer",
      "failure_scenario": "If a Repr somehow has exponent near isize::MAX and significand with 20+ digits, the addition wraps to a negative value and x_mag becomes 0, producing an
  undersized work context"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 428,
      "summary": "atan2(finite_y, +inf_x) returns unsigned ZERO regardless of y's sign, violating IEEE 754 which requires signed zero — though Repr cannot represent negative zero,
  this is a spec deviation worth documenting",
      "failure_scenario": "atan2(-1, +inf) == atan2(+1, +inf) == 0, but IEEE 754 requires atan2(-1, +inf) == -0"
    },
    {
      "file": "float/src/math/trig.rs",
      "line": 64,
      "summary": "sin/cos check x.significand.is_zero() instead of x.is_zero() for the zero fast-path — currently safe because infinity is handled earlier, but fragile if code is
  reordered",
      "failure_scenario": "If someone reorders the infinity check after the zero check, infinities (which have significand==0) would incorrectly return ZERO instead of NaN"
    }
  ]

@CokieMiner

Copy link
Copy Markdown
Contributor Author

Bro I'm not a openclaw 😭😭, at least format it in bullet points so it is readable I read every response/thinking tokens my ai outputs so I can make sure they don't deviate from my vision, and can correct them in the middle of thinking and implementation if they start to drift.

Also like some of those are ok some don't make sense if the algorithm as linear convergence and I did research and didn't find a better one how the fuck does that AI want me to improve it write a full paper and proof with a new algorithm that does it, get a better LLM.

@cmpute

cmpute commented Jun 6, 2026

Copy link
Copy Markdown
Owner

@CokieMiner Sorry about that.. Well I can fix those issues later myself, but can you separate the fuzz test to a separate crate? That will make CI and local testing easier and cleaner.

@CokieMiner

Copy link
Copy Markdown
Contributor Author

already did all just doing a final ceck on a detail

@CokieMiner
CokieMiner force-pushed the feat/trigonometry branch 2 times, most recently from 2ddf0e2 to 3b7bc33 Compare June 6, 2026 15:58
@CokieMiner

Copy link
Copy Markdown
Contributor Author

everything fixed and clean history

@CokieMiner

Copy link
Copy Markdown
Contributor Author

Ok one of them is my bad, the other is do to fuzz workspace having rug as a dep that doesn't build in arm you will need to exclude fuzz compilation for arm targets

@CokieMiner
CokieMiner force-pushed the feat/trigonometry branch 2 times, most recently from f391b22 to 31351c2 Compare June 6, 2026 16:54
@CokieMiner

Copy link
Copy Markdown
Contributor Author

Ok noticed main add a lot of updates so rebased on top

@cmpute cmpute left a comment

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for your efforts, but please read these suggestions to make the coding style of the crate more consistent

Comment thread fuzz/Cargo.toml Outdated
Comment thread fuzz/Cargo.toml
Comment thread Cargo.toml Outdated
Comment thread float/src/math/consts.rs
Comment thread float/src/math/consts.rs Outdated
Comment thread float/src/math/trig.rs Outdated
Comment thread float/src/math/trig.rs Outdated
Comment thread float/src/round_ops.rs Outdated
Comment thread float/tests/trig.rs
Comment thread float/Cargo.toml
@cmpute

cmpute commented Jun 7, 2026

Copy link
Copy Markdown
Owner

The CI thing was broken because fuzz crate is included in the build. Exclude it from the workspace will fix.

@CokieMiner
CokieMiner force-pushed the feat/trigonometry branch from 61878da to 6e00ea5 Compare June 7, 2026 05:47
@CokieMiner

Copy link
Copy Markdown
Contributor Author

Any confusion left just say

@CokieMiner

Copy link
Copy Markdown
Contributor Author

Also while thinking of a cache strategy, I came up with specific design that could be useful, don't know what other libs use but:

After calculating pi for the first time, we cache it. Every time pi is requested again, we compare the requested precision with the stored one. If the requested precision is smaller, we simply round the cached value and return it. If it's higher, we could potentially use the cached value as a starting point to compute the missing digits, and then update the cache with this new, higher-precision value.

For example for pi we: Store the final FBig and the raw Chudnovsky integers (P, Q, T, terms) alongside it.

Since binary splitting is associative, if a higher precision is later requested, we only need to compute chudnovsky_bs(cached_terms, n_new) and combine with the cached (P, Q, T) using the standard recurrence, avoiding a full recomputation from scratch. For lower precisions, we just round the cached FBig down, which is nearly free.

This gives three O(1)-or-better cases:

  • precision < cache → round cached FBig
  • precision == cache → return directly
  • precision > cache → extend from cached_terms, update cache

A problem you might have with this design is that the cache will probably have to be constant specific.

I can try to implement it but can only do it Tuesday.

@cmpute

cmpute commented Jun 7, 2026

Copy link
Copy Markdown
Owner

@CokieMiner Yes I did think about adding a caching strategy, and I did note that in my TODOs. However, it's a little bit complex in Rust, we need to consider threading safety. The best way I can think is to associate the cache with the context.

I will suggest leave it for the next PR (and for the next big version of dashu)

Comment thread float/src/math/trig.rs
);
};

let k_mod_4_big = k.rem_euclid(IBig::from(4));

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

actually you can directly call k.rem_euclid(4u8), iirc which directly produces an u8

@CokieMiner CokieMiner Jun 7, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just tried that even commited without testing, but IBig doesn't implement the RemEuclid < u8 > trait. I have to use IBig::from(4) and downcast it as a workaround.

Comment thread float/src/round_ops.rs Outdated
@CokieMiner
CokieMiner force-pushed the feat/trigonometry branch from 2daeb1d to ecdf1de Compare June 7, 2026 15:43
@cmpute

cmpute commented Jun 7, 2026

Copy link
Copy Markdown
Owner

Well there is currently no easy way to create const UBig from u64 outside dashu-int crate, I think for the bigger constant we can stick to UBig::from()

@CokieMiner

Copy link
Copy Markdown
Contributor Author

I think anything other than a global cache defeats the purpose of the caching, but happy to defer to your judgment on the design.

@CokieMiner

Copy link
Copy Markdown
Contributor Author

I just checked and UBig::from does the optimization you asked for inside.

@CokieMiner
CokieMiner force-pushed the feat/trigonometry branch from ecdf1de to dab788f Compare June 7, 2026 16:15
@CokieMiner

Copy link
Copy Markdown
Contributor Author

Any other problem you find just leave a comment. Think I responded to all questions and all details you pointed out.

@cmpute
cmpute merged commit 59bb3c5 into cmpute:master Jun 8, 2026
13 checks passed
@cmpute

cmpute commented Jun 8, 2026

Copy link
Copy Markdown
Owner

Thanks for your efforts on this PR! I think we can merge this now and make future improvements step by step

@cmpute

cmpute commented Jun 8, 2026

Copy link
Copy Markdown
Owner

Regarding constants caching, we can discuss in #15, BTW

@CokieMiner
CokieMiner deleted the feat/trigonometry branch June 8, 2026 22:27
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants