josdejong / mathjs

An extensive math library for JavaScript and Node.js
https://mathjs.org
Apache License 2.0
14.41k stars 1.24k forks source link

Overflow encountered with tanh and coth functions #2967

Open richinex opened 1 year ago

richinex commented 1 year ago

Hi, I work with impedance spectroscopy and I sometimes have to deal with large complex numbers. I have experience NaN results. I think the mathjs tanh implementation has an overflow issue as I do not experience this problem with Python. Take for instance the problem below:

const s = math.multiply(math.i, 2 * Math.PI * 15000);

    console.log("s:", s);
    const Rd = math.complex(1, 0);
    console.log("Rd:", Rd);
    const Cd = math.complex(10, 0);
    console.log("Cd:", Cd);
    const RdtimesSCd = math.multiply(Rd, math.multiply(s, Cd));
    // const RdOverSCd = Rd/(math.multiply(s, Cd));
    console.log("Rd * (s * Cd):", RdtimesSCd);
    const sqrtRdtimesSCd = math.sqrt(RdtimesSCd);
    console.log("sqrt(Rd * (s * Cd)):", sqrtRdtimesSCd);
    const tanhSqrtRdtimesSCd = math.tanh(sqrtRdtimesSCd);
    console.log("tanh(sqrt(Rd * (s * Cd))):", tanhSqrtRdtimesSCd);
    console.log("test using custom Tanh is: ", (custom_tanh(math.sqrt(math.multiply(Rd, math.multiply(s, Cd))))));

// Results
s: Complex {re: 0, im: 94247.7796076938}
Simulator.js:421 Rd: Complex {re: 1, im: 0}
Simulator.js:423 Cd: Complex {re: 10, im: 0}
Simulator.js:426 Rd * (s * Cd): Complex {re: 0, im: 942477.7960769379}
Simulator.js:428 sqrt(Rd * (s * Cd)): Complex {re: 686.4684246478267, im: 686.4684246478267}
Simulator.js:430 tanh(sqrt(Rd * (s * Cd))): Complex {re: NaN, im: -0}  //mathjs gives NaNs
Simulator.js:431 test using custom Tanh is:  Complex {re: 1, im: -0} // correct result

So I thought I'd share my implementation below. Maybe it helps. It is a translation of numba's implementation to JavaScript https://github.com/numba/numba/blob/c04edc02deb40b84c3b21ed2a9835353181a5916/numba/targets/cmathimpl.py#L343-L367:

import * as math from 'mathjs';

export function custom_tanh(z) {
    if (math.isComplex(z)) {
      const x = z.re;
      const y = z.im;

      if (Number.isFinite(x)) {
        const tx = math.tanh(x);
        const ty = Math.tan(y);
        const cx = 1 / math.cosh(x);
        const txty = tx * ty;

        const denom = 1 + txty * txty;
        const newRealPart = tx * (1 + ty * ty) / denom;
        const newImaginaryPart = ((ty / denom) * cx) * cx;

        return math.complex(newRealPart, newImaginaryPart);
      } else {
        const real = Math.sign(x);
        const imag = Number.isFinite(y) ? Math.sign(Math.sin(2 * y)) : 0;

        return math.complex(real, imag);
      }
    } else {
      return math.tanh(z);
    }
  }

 export function custom_coth(z) {
    if (math.isComplex(z)) {
      // Avoid division by zero
      if (z.re === 0 && z.im === 0) {
        return math.complex(Infinity, Infinity);
      }

      const tanhZ = custom_tanh(z);
      const one = math.complex(1, 0);
      return math.divide(one, tanhZ);
    } else {
      // Avoid division by zero
      if (z === 0) {
        return Infinity;
      }

      return 1 / math.tanh(z);
    }
  }
josdejong commented 1 year ago

Thanks for sharing. There is indeed an overflow occurring in the calculation of tanh for example. Here a minimal demo:

math.tanh(math.complex(1, 1))       // { re: 1.0839233273386948, im: 0.27175258531951174 }
math.tanh(math.complex(100, 100))   // { re: 1, im: -2.4171061928460553e-87 }
math.tanh(math.complex(1000, 1000)) // { re: NaN, im: -0 }

The tanh and coth functions for Complex numbers are currently straight from the Complex.js library, the calculation contains cosh which goes to Infinity for large values:

https://github.com/infusion/Complex.js/blob/9bcf22ac333b1c2225164857f221cddb10f979dd/complex.js#L942-L958

Before solving this in mathjs, it will be good to raise this issue in the https://github.com/infusion/Complex.js project, ideally it would be solved there. @richinex can you open an issue there?

richinex commented 1 year ago

Yea sure. I have raised an issue there.