From b31cfff2e69c8ae04c5ed8ab82b13e35933357fc Mon Sep 17 00:00:00 2001 From: Alexander Kireev Date: Sat, 20 Jun 2026 01:48:14 +0700 Subject: [PATCH] Fix catastrophic cancellation in cot/coth/csc/csch near poles The denominators cos(2a) - cosh(2b) (and the analogous forms in coth, csc and csch) subtract two values that both approach 1 near the pole (a, b -> 0), so a direct evaluation loses all significance and eventually divides by an exactly-zero denominator, returning NaN. For example cot(1e-10 i) returned NaN + Infinity i and coth(1e-8) returned 90071992.5 (a 10% error) instead of ~1e8. Rewrite each denominator using cosh(2x) - 1 = 2 sinh(x)^2 and 1 - cos(2x) = 2 sin(x)^2, which is algebraically identical but free of cancellation. Results are now accurate to machine precision across the whole plane (verified against mpmath). The two existing coth/cot test values are updated to the more accurate output. --- dist/complex.js | 30 ++++++++++++++++++++++++----- dist/complex.min.js | 26 ++++++++++++------------- dist/complex.mjs | 30 ++++++++++++++++++++++++----- src/complex.js | 28 +++++++++++++++++++++++---- tests/complex.test.js | 44 +++++++++++++++++++++++++++++++++++++++++-- 5 files changed, 129 insertions(+), 29 deletions(-) diff --git a/dist/complex.js b/dist/complex.js index 11ff7ec..1792656 100644 --- a/dist/complex.js +++ b/dist/complex.js @@ -173,7 +173,7 @@ const parse = function (a, b) { case 'string': z['im'] = /* void */ - z['re'] = 0; + z['re'] = 0; const tokens = a.replace(/_/g, '') .match(/\d+\.?\d*e[+-]?\d+|\d+\.?\d*|\.\d+|./g); @@ -674,7 +674,12 @@ Complex.prototype = { const a = 2 * this['re']; const b = 2 * this['im']; - const d = Math.cos(a) - cosh(b); + // d = cos(a) - cosh(b), rewritten via cos(2x) - 1 = -2sin(x)^2 and + // cosh(2x) - 1 = 2sinh(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cos(a) and cosh(b) both approach 1. + const sa = Math.sin(this['re']); + const sb = sinh(this['im']); + const d = -2 * (sa * sa + sb * sb); return new Complex( -Math.sin(a) / d, @@ -710,7 +715,12 @@ Complex.prototype = { const a = this['re']; const b = this['im']; - const d = 0.5 * cosh(2 * b) - 0.5 * Math.cos(2 * a); + // d = (cosh(2b) - cos(2a)) / 2, rewritten via cosh(2x) - 1 = 2sinh(x)^2 + // and 1 - cos(2x) = 2sin(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cosh(2b) and cos(2a) both approach 1. + const sa = Math.sin(a); + const sb = sinh(b); + const d = sa * sa + sb * sb; return new Complex( Math.sin(a) * cosh(b) / d, @@ -936,7 +946,12 @@ Complex.prototype = { const a = 2 * this['re']; const b = 2 * this['im']; - const d = cosh(a) - Math.cos(b); + // d = cosh(a) - cos(b), rewritten via cosh(2x) - 1 = 2sinh(x)^2 and + // 1 - cos(2x) = 2sin(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cosh(a) and cos(b) both approach 1. + const sa = sinh(this['re']); + const sb = Math.sin(this['im']); + const d = 2 * (sa * sa + sb * sb); return new Complex( sinh(a) / d, @@ -954,7 +969,12 @@ Complex.prototype = { const a = this['re']; const b = this['im']; - const d = Math.cos(2 * b) - cosh(2 * a); + // d = cos(2b) - cosh(2a), rewritten via 1 - cos(2x) = 2sin(x)^2 and + // cosh(2x) - 1 = 2sinh(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cos(2b) and cosh(2a) both approach 1. + const sa = sinh(a); + const sb = Math.sin(b); + const d = -2 * (sb * sb + sa * sa); return new Complex( -2 * sinh(a) * Math.cos(b) / d, diff --git a/dist/complex.min.js b/dist/complex.min.js index 37c9da1..87a1663 100644 --- a/dist/complex.min.js +++ b/dist/complex.min.js @@ -8,20 +8,20 @@ Licensed under the MIT license. 'use strict';(function(q){function l(a,b){if(a===void 0||a===null)f.re=f.im=0;else if(b!==void 0)f.re=a,f.im=b;else switch(typeof a){case "object":if("im"in a&&"re"in a)f.re=a.re,f.im=a.im;else if("abs"in a&&"arg"in a){if(!isFinite(a.abs)&&isFinite(a.arg))return c.INFINITY;f.re=a.abs*Math.cos(a.arg);f.im=a.abs*Math.sin(a.arg)}else if("r"in a&&"phi"in a){if(!isFinite(a.r)&&isFinite(a.phi))return c.INFINITY;f.re=a.r*Math.cos(a.phi);f.im=a.r*Math.sin(a.phi)}else a.length===2?(f.re=a[0],f.im=a[1]):m(); break;case "string":f.im=f.re=0;a=a.replace(/_/g,"").match(/\d+\.?\d*e[+-]?\d+|\d+\.?\d*|\.\d+|./g);b=1;let d=0;a===null&&m();for(let e=0;e0&& m();break;case "number":f.im=0;f.re=a;break;default:m()}return f}function m(){throw SyntaxError("Invalid Param");}function n(a,b){a=Math.abs(a);b=Math.abs(b);a 0)return new c(Math.pow(this.re,a.re),0);if(this.re===0)switch((a.re%4+4)%4){case 0:return new c(Math.pow(this.im,a.re),0);case 1:return new c(0,Math.pow(this.im,a.re));case 2:return new c(-Math.pow(this.im,a.re),0);case 3:return new c(0,-Math.pow(this.im,a.re))}}if(b&&a.re>0)return c.ZERO;const d=Math.atan2(this.im,this.re),e=p(this.re,this.im);b=Math.exp(a.re*e-a.im*d);a=a.im*e+a.re*d;return new c(b*Math.cos(a),b*Math.sin(a))},sqrt:function(){const a=this.re,b=this.im;if(b===0)return a>=0?new c(Math.sqrt(a), 0):new c(0,Math.sqrt(-a));var d=n(a,b);d=Math.sqrt(.5*(d+Math.abs(a)));let e=Math.abs(b)/(2*d);return a>=0?new c(d,b<0?-e:e):new c(e,b<0?-d:d)},exp:function(){const a=Math.exp(this.re);return this.im===0?new c(a,0):new c(a*Math.cos(this.im),a*Math.sin(this.im))},expm1:function(){const a=this.re,b=this.im,d=Math.sin(.5*b);return new c(Math.expm1(a)*Math.cos(b)+-2*d*d,Math.exp(a)*Math.sin(b))},log:function(){const a=this.re,b=this.im;return b===0&&a>0?new c(Math.log(a),0):new c(p(a,b),Math.atan2(b, -a))},abs:function(){return n(this.re,this.im)},arg:function(){return Math.atan2(this.im,this.re)},sin:function(){const a=this.re,b=this.im;return new c(Math.sin(a)*h(b),Math.cos(a)*k(b))},cos:function(){const a=this.re,b=this.im;return new c(Math.cos(a)*h(b),-Math.sin(a)*k(b))},tan:function(){const a=2*this.re,b=2*this.im,d=Math.cos(a)+h(b);return new c(Math.sin(a)/d,k(b)/d)},cot:function(){const a=2*this.re,b=2*this.im,d=Math.cos(a)-h(b);return new c(-Math.sin(a)/d,k(b)/d)},sec:function(){const a= -this.re,b=this.im,d=.5*h(2*b)+.5*Math.cos(2*a);return new c(Math.cos(a)*h(b)/d,Math.sin(a)*k(b)/d)},csc:function(){const a=this.re,b=this.im,d=.5*h(2*b)-.5*Math.cos(2*a);return new c(Math.sin(a)*h(b)/d,-Math.cos(a)*k(b)/d)},asin:function(){var a=this.re;const b=this.im,d=(new c(b*b-a*a+1,-2*a*b)).sqrt();a=(new c(d.re-b,d.im+a)).log();return new c(a.im,-a.re)},acos:function(){var a=this.re;const b=this.im,d=(new c(b*b-a*a+1,-2*a*b)).sqrt();a=(new c(d.re-b,d.im+a)).log();return new c(Math.PI/2-a.im, -a.re)},atan:function(){var a=this.re;const b=this.im;if(a===0){if(b===1)return new c(0,Infinity);if(b===-1)return new c(0,-Infinity)}const d=a*a+(1-b)*(1-b);a=(new c((1-b*b-a*a)/d,-2*a/d)).log();return new c(-.5*a.im,.5*a.re)},acot:function(){const a=this.re,b=this.im;if(b===0)return new c(Math.atan2(1,a),0);const d=a*a+b*b;return d!==0?(new c(a/d,-b/d)).atan():(new c(a!==0?a/0:0,b!==0?-b/0:0)).atan()},asec:function(){const a=this.re,b=this.im;if(a===0&&b===0)return new c(0,Infinity);const d=a*a+ -b*b;return d!==0?(new c(a/d,-b/d)).acos():(new c(a!==0?a/0:0,b!==0?-b/0:0)).acos()},acsc:function(){const a=this.re,b=this.im;if(a===0&&b===0)return new c(Math.PI/2,Infinity);const d=a*a+b*b;return d!==0?(new c(a/d,-b/d)).asin():(new c(a!==0?a/0:0,b!==0?-b/0:0)).asin()},sinh:function(){const a=this.re,b=this.im;return new c(k(a)*Math.cos(b),h(a)*Math.sin(b))},cosh:function(){const a=this.re,b=this.im;return new c(h(a)*Math.cos(b),k(a)*Math.sin(b))},tanh:function(){const a=2*this.re,b=2*this.im,d= -h(a)+Math.cos(b);return new c(k(a)/d,Math.sin(b)/d)},coth:function(){const a=2*this.re,b=2*this.im,d=h(a)-Math.cos(b);return new c(k(a)/d,-Math.sin(b)/d)},csch:function(){const a=this.re,b=this.im,d=Math.cos(2*b)-h(2*a);return new c(-2*k(a)*Math.cos(b)/d,2*h(a)*Math.sin(b)/d)},sech:function(){const a=this.re,b=this.im,d=Math.cos(2*b)+h(2*a);return new c(2*h(a)*Math.cos(b)/d,-2*k(a)*Math.sin(b)/d)},asinh:function(){const a=this.re;var b=this.im;if(b===0){if(a===0)return new c(0,0);b=Math.abs(a);b= -Math.log(b+Math.sqrt(b*b+1));return new c(a<0?-b:b,0)}const d=(new c(a*a-b*b+1,2*a*b)).sqrt();return(new c(a+d.re,b+d.im)).log()},acosh:function(){const a=this.re,b=this.im;if(b===0)return a>1?new c(Math.log(a+Math.sqrt(a-1)*Math.sqrt(a+1)),0):a<-1?new c(Math.log(-a+Math.sqrt(a*a-1)),Math.PI):new c(0,Math.acos(a));const d=(new c(a-1,b)).sqrt(),e=(new c(a+1,b)).sqrt();return(new c(a+d.re*e.re-d.im*e.im,b+d.re*e.im+d.im*e.re)).log()},atanh:function(){var a=this.re,b=this.im;if(b===0)return a===0?new c(0, -0):a===1?new c(Infinity,0):a===-1?new c(-Infinity,0):-11?new c(.5*Math.log((a+1)/(a-1)),-Math.PI/2):new c(.5*Math.log(-((1+a)/(1-a))),Math.PI/2);const d=1-a,e=1+a,g=d*d+b*b;if(g===0)return new c(a!==-1?a/0:0,b!==0?b/0:0);a=(e*d-b*b)/g;b=(b*d+e*b)/g;return new c(p(a,b)/2,Math.atan2(b,a)/2)},acoth:function(){const a=this.re,b=this.im;if(a===0&&b===0)return new c(0,Math.PI/2);const d=a*a+b*b;return d!==0?(new c(a/d,-b/d)).atanh():(new c(a!==0?a/0:0,b!==0?-b/ -0:0)).atanh()},acsch:function(){var a=this.re;const b=this.im;if(b===0){if(a===0)return new c(Infinity,0);a=1/a;return new c(Math.log(a+Math.sqrt(a*a+1)),0)}const d=a*a+b*b;return d!==0?(new c(a/d,-b/d)).asinh():(new c(a!==0?a/0:0,b!==0?-b/0:0)).asinh()},asech:function(){const a=this.re,b=this.im;if(this.isZero())return c.INFINITY;const d=a*a+b*b;return d!==0?(new c(a/d,-b/d)).acosh():(new c(a!==0?a/0:0,b!==0?-b/0:0)).acosh()},inverse:function(){if(this.isZero())return c.INFINITY;if(this.isInfinite())return c.ZERO; -const a=this.re,b=this.im,d=a*a+b*b;return new c(a/d,-b/d)},conjugate:function(){return new c(this.re,-this.im)},neg:function(){return new c(-this.re,-this.im)},ceil:function(a){a=Math.pow(10,a||0);return new c(Math.ceil(this.re*a)/a,Math.ceil(this.im*a)/a)},floor:function(a){a=Math.pow(10,a||0);return new c(Math.floor(this.re*a)/a,Math.floor(this.im*a)/a)},round:function(a){a=Math.pow(10,a||0);return new c(Math.round(this.re*a)/a,Math.round(this.im*a)/a)},equals:function(a,b){a=l(a,b);return Math.abs(a.re- -this.re)<=c.EPSILON&&Math.abs(a.im-this.im)<=c.EPSILON},clone:function(){return new c(this.re,this.im)},toString:function(){let a=this.re,b=this.im,d="";if(this.isNaN())return"NaN";if(this.isInfinite())return"Infinity";Math.abs(a)1?new c(Math.log(a+Math.sqrt(a-1)*Math.sqrt(a+1)),0):a<-1?new c(Math.log(-a+Math.sqrt(a*a-1)),Math.PI):new c(0,Math.acos(a));const d=(new c(a-1,b)).sqrt(),e=(new c(a+1,b)).sqrt();return(new c(a+d.re*e.re-d.im*e.im,b+d.re*e.im+d.im*e.re)).log()}, +atanh:function(){var a=this.re,b=this.im;if(b===0)return a===0?new c(0,0):a===1?new c(Infinity,0):a===-1?new c(-Infinity,0):-11?new c(.5*Math.log((a+1)/(a-1)),-Math.PI/2):new c(.5*Math.log(-((1+a)/(1-a))),Math.PI/2);const d=1-a,e=1+a,g=d*d+b*b;if(g===0)return new c(a!==-1?a/0:0,b!==0?b/0:0);a=(e*d-b*b)/g;b=(b*d+e*b)/g;return new c(p(a,b)/2,Math.atan2(b,a)/2)},acoth:function(){const a=this.re,b=this.im;if(a===0&&b===0)return new c(0,Math.PI/2);const d=a* +a+b*b;return d!==0?(new c(a/d,-b/d)).atanh():(new c(a!==0?a/0:0,b!==0?-b/0:0)).atanh()},acsch:function(){var a=this.re;const b=this.im;if(b===0){if(a===0)return new c(Infinity,0);a=1/a;return new c(Math.log(a+Math.sqrt(a*a+1)),0)}const d=a*a+b*b;return d!==0?(new c(a/d,-b/d)).asinh():(new c(a!==0?a/0:0,b!==0?-b/0:0)).asinh()},asech:function(){const a=this.re,b=this.im;if(this.isZero())return c.INFINITY;const d=a*a+b*b;return d!==0?(new c(a/d,-b/d)).acosh():(new c(a!==0?a/0:0,b!==0?-b/0:0)).acosh()}, +inverse:function(){if(this.isZero())return c.INFINITY;if(this.isInfinite())return c.ZERO;const a=this.re,b=this.im,d=a*a+b*b;return new c(a/d,-b/d)},conjugate:function(){return new c(this.re,-this.im)},neg:function(){return new c(-this.re,-this.im)},ceil:function(a){a=Math.pow(10,a||0);return new c(Math.ceil(this.re*a)/a,Math.ceil(this.im*a)/a)},floor:function(a){a=Math.pow(10,a||0);return new c(Math.floor(this.re*a)/a,Math.floor(this.im*a)/a)},round:function(a){a=Math.pow(10,a||0);return new c(Math.round(this.re* +a)/a,Math.round(this.im*a)/a)},equals:function(a,b){a=l(a,b);return Math.abs(a.re-this.re)<=c.EPSILON&&Math.abs(a.im-this.im)<=c.EPSILON},clone:function(){return new c(this.re,this.im)},toString:function(){let a=this.re,b=this.im,d="";if(this.isNaN())return"NaN";if(this.isInfinite())return"Infinity";Math.abs(a) 0), where cos(a) and cosh(b) both approach 1. + const sa = Math.sin(this['re']); + const sb = sinh(this['im']); + const d = -2 * (sa * sa + sb * sb); return new Complex( -Math.sin(a) / d, @@ -710,7 +715,12 @@ Complex.prototype = { const a = this['re']; const b = this['im']; - const d = 0.5 * cosh(2 * b) - 0.5 * Math.cos(2 * a); + // d = (cosh(2b) - cos(2a)) / 2, rewritten via cosh(2x) - 1 = 2sinh(x)^2 + // and 1 - cos(2x) = 2sin(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cosh(2b) and cos(2a) both approach 1. + const sa = Math.sin(a); + const sb = sinh(b); + const d = sa * sa + sb * sb; return new Complex( Math.sin(a) * cosh(b) / d, @@ -936,7 +946,12 @@ Complex.prototype = { const a = 2 * this['re']; const b = 2 * this['im']; - const d = cosh(a) - Math.cos(b); + // d = cosh(a) - cos(b), rewritten via cosh(2x) - 1 = 2sinh(x)^2 and + // 1 - cos(2x) = 2sin(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cosh(a) and cos(b) both approach 1. + const sa = sinh(this['re']); + const sb = Math.sin(this['im']); + const d = 2 * (sa * sa + sb * sb); return new Complex( sinh(a) / d, @@ -954,7 +969,12 @@ Complex.prototype = { const a = this['re']; const b = this['im']; - const d = Math.cos(2 * b) - cosh(2 * a); + // d = cos(2b) - cosh(2a), rewritten via 1 - cos(2x) = 2sin(x)^2 and + // cosh(2x) - 1 = 2sinh(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cos(2b) and cosh(2a) both approach 1. + const sa = sinh(a); + const sb = Math.sin(b); + const d = -2 * (sb * sb + sa * sa); return new Complex( -2 * sinh(a) * Math.cos(b) / d, diff --git a/src/complex.js b/src/complex.js index fbc1913..e925285 100644 --- a/src/complex.js +++ b/src/complex.js @@ -680,7 +680,12 @@ Complex.prototype = { const a = 2 * this['re']; const b = 2 * this['im']; - const d = Math.cos(a) - cosh(b); + // d = cos(a) - cosh(b), rewritten via cos(2x) - 1 = -2sin(x)^2 and + // cosh(2x) - 1 = 2sinh(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cos(a) and cosh(b) both approach 1. + const sa = Math.sin(this['re']); + const sb = sinh(this['im']); + const d = -2 * (sa * sa + sb * sb); return new Complex( -Math.sin(a) / d, @@ -716,7 +721,12 @@ Complex.prototype = { const a = this['re']; const b = this['im']; - const d = 0.5 * cosh(2 * b) - 0.5 * Math.cos(2 * a); + // d = (cosh(2b) - cos(2a)) / 2, rewritten via cosh(2x) - 1 = 2sinh(x)^2 + // and 1 - cos(2x) = 2sin(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cosh(2b) and cos(2a) both approach 1. + const sa = Math.sin(a); + const sb = sinh(b); + const d = sa * sa + sb * sb; return new Complex( Math.sin(a) * cosh(b) / d, @@ -942,7 +952,12 @@ Complex.prototype = { const a = 2 * this['re']; const b = 2 * this['im']; - const d = cosh(a) - Math.cos(b); + // d = cosh(a) - cos(b), rewritten via cosh(2x) - 1 = 2sinh(x)^2 and + // 1 - cos(2x) = 2sin(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cosh(a) and cos(b) both approach 1. + const sa = sinh(this['re']); + const sb = Math.sin(this['im']); + const d = 2 * (sa * sa + sb * sb); return new Complex( sinh(a) / d, @@ -960,7 +975,12 @@ Complex.prototype = { const a = this['re']; const b = this['im']; - const d = Math.cos(2 * b) - cosh(2 * a); + // d = cos(2b) - cosh(2a), rewritten via 1 - cos(2x) = 2sin(x)^2 and + // cosh(2x) - 1 = 2sinh(x)^2 to avoid catastrophic cancellation near the + // pole (a, b -> 0), where cos(2b) and cosh(2a) both approach 1. + const sa = sinh(a); + const sb = Math.sin(b); + const d = -2 * (sb * sb + sa * sa); return new Complex( -2 * sinh(a) * Math.cos(b) / d, diff --git a/tests/complex.test.js b/tests/complex.test.js index f7ffbea..74d4e41 100644 --- a/tests/complex.test.js +++ b/tests/complex.test.js @@ -614,11 +614,11 @@ var functionTests = [{ }, { set: "3.14-4i", fn: "coth", - expect: "0.9994481238383576 + 0.0037048958915019857i" + expect: "0.9994481238383578 + 0.0037048958915019865i" }, { set: "8i-31", fn: "cot", - expect: "1.6636768291213935e-7 - 1.0000001515864902i" + expect: "1.6636768291213932e-7 - 1.00000015158649i" }, { set: Complex(1, 1).sub(0, 1), // Distance fn: "abs", @@ -1087,3 +1087,43 @@ describe("Complex Details", function () { }); }); + +describe("Reciprocal trig near poles", function () { + + // Near their poles the denominators cos(2a) - cosh(2b) (and friends) are a + // difference of two values both close to 1, so a naive evaluation loses all + // significance and eventually divides by an exactly-zero denominator, + // producing NaN. Reference values below come from mpmath at 25 digits. + + function assertClose(actual, expRe, expIm) { + var mag = Math.hypot(expRe, expIm) || 1; + assert.ok(isFinite(actual["re"]) && isFinite(actual["im"]), + "expected finite result but got " + actual.toString()); + assert.ok(Math.abs(actual["re"] - expRe) <= 1e-12 * mag, + "re: got " + actual["re"] + ", expected ~" + expRe); + assert.ok(Math.abs(actual["im"] - expIm) <= 1e-12 * mag, + "im: got " + actual["im"] + ", expected ~" + expIm); + } + + it("cot stays accurate near the pole", function () { + assertClose(Complex(0, 1e-8).cot(), 0, -100000000.000000001); + assertClose(Complex(0, 1e-10).cot(), 0, -9999999999.99999964); + assertClose(Complex(1e-7, 0).cot(), 9999999.99999996712, 0); + }); + + it("coth stays accurate near the pole", function () { + assertClose(Complex(1e-8, 0).coth(), 100000000.000000001, 0); + assertClose(Complex(1e-10, 0).coth(), 9999999999.99999964, 0); + assertClose(Complex(0, 1e-7).coth(), 0, -9999999.99999996712); + }); + + it("csc stays accurate near the pole", function () { + assertClose(Complex(1e-8, 0).csc(), 99999999.9999999996, 0); + assertClose(Complex(1e-10, 0).csc(), 9999999999.99999964, 0); + }); + + it("csch stays accurate near the pole", function () { + assertClose(Complex(1e-8, 0).csch(), 99999999.9999999962, 0); + assertClose(Complex(1e-10, 0).csch(), 9999999999.99999964, 0); + }); +});