josdejong/mathjs

Overflow encountered with tanh and coth functions

Open

#2,967 建立於 2023年6月8日

在 GitHub 查看
 (2 留言) (0 反應) (0 負責人)JavaScript (1,298 fork)batch import
bugcategory:numericaldependencieshelp wanted

倉庫指標

Star
 (13,832 star)
PR 合併指標
 (30 天內沒有已合併 PR)

描述

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);
    }
  }

貢獻者指南