Use robust range reduction for trigonometric functions at large finite angles

[LunaStev]

wrap_angle_pi_f64 currently reduces an angle with floating-point remainder against the stored 2π constant:

value % MATH_TWO_PI_F64

This keeps the intermediate numerically bounded, but it does not preserve enough phase information for very large finite f64 inputs.

Once the magnitude becomes large enough, binary64 no longer contains enough low-order quotient information for a simple remainder against an approximated 2π constant. sin_f64, cos_f64, and tan_f64 then evaluate their polynomials at an incorrectly reduced angle.

This is separate from polynomial accuracy near the origin: the dominant error is the large-argument range-reduction step.

Representative inputs to verify against a trusted reference include:

• 1e10
• 1e16
• 1e20
• large negative values
• values near large integer multiples of π/2 and 2π
• finite values close to the upper range of f64

Code evidence:

• std/math/trig.wave — wrap_angle_pi_f64
• std/math/trig.wave — sin_f64, cos_f64, tan_f64

Acceptance:

• Introduce a range-reduction strategy that remains meaningful for large finite f64 values.
• Preserve the existing NaN behavior for non-finite input.
• Add accuracy tests against a trusted reference over small, medium, and large magnitudes.
• Include difficult values near quadrant boundaries.
• Verify both f64 and the f32 wrappers.
• Document the intended error tolerance or supported numerical contract so future changes can be regression-tested.

A Cody-Waite style fast path plus a higher-precision reduction path, or another established range-reduction strategy, would fit this issue; the implementation choice is intentionally left open.