Fix catastrophic cancellation in cot/coth/csc/csch near poles

[chatman-media]

Problem

cot, coth, csc and csch lose accuracy — and eventually return NaN — near their poles.

Each is computed as numerator / d where the denominator is a difference of two terms that both approach 1 as the argument approaches the pole:

fn	denominator	near the pole (a, b → 0)
cot	cos(2a) − cosh(2b)	(1 − …) − (1 + …)
coth	cosh(2a) − cos(2b)	(1 + …) − (1 − …)
csc	½cosh(2b) − ½cos(2a)	½(1 + …) − ½(1 − …)
csch	cos(2b) − cosh(2a)	(1 − …) − (1 + …)

Subtracting two nearly-equal numbers is catastrophic cancellation: the leading 1s cancel and only rounding noise survives. Once the terms round to exactly 1, the denominator becomes 0 and the result is NaN/Infinity.

const Complex = require("complex.js");

Complex(0, 1e-8).cot().toString();   // "-90071992.54740992i"   (true: -1e8 i, ~10% off)
Complex(0, 1e-10).cot().toString();  // "NaN + Infinityi"        (true: -1e10 i)
Complex(1e-8, 0).coth().toString();  // "90071992.54740992"      (true:  1e8)
Complex(1e-8, 0).csc().toString();   // "90071992.54740994"      (true: ~1e8)
Complex(1e-8, 0).csch().toString();  // "90071992.54740992"      (true: ~1e8)

These functions are real-valued on the real (resp. imaginary) axis, so the error is easy to confirm: coth(x) → 1/tanh(x), cot(x) → 1/tan(x), etc.

Fix

Rewrite each denominator with the algebraically-identical, cancellation-free identities

cosh(2x) − 1 = 2·sinh(x)²
1 − cos(2x)  = 2·sin(x)²

e.g. for cot: cos(2a) − cosh(2b) = −(2·sin(a)² + 2·sinh(b)²). No subtraction of near-equal terms, so full precision is retained right up to the pole. The numerators are unchanged.

This is the same class of numerical-stability fix as #54 (which addressed tanh overflow); the two are disjoint — #54 touches the +cosh family for large inputs, this touches the −cosh family for small inputs.

Verification

• Added regression tests for the four functions near their poles (they fail on main, with errors up to 10% / NaN, and pass with the fix). Reference values computed with mpmath at 25 digits.
• Fuzzed all four against mpmath over a grid plus random points (incl. near-pole): worst-case relative error drops from ~1e-1 (and NaN) to ~5e-16 (machine precision).
• Full existing suite passes. The two pinned coth/cot test values were updated — the new outputs match mpmath exactly (the old ones were off by 1 ulp), mirroring the test update in Improve numerical stability of tanh #54.
• dist/* rebuilt with crude-build.