A different `getGcd` implementation

[Rudxain]

The getGcd function defined at src/utils.ts is dealing with BigInts, so I think it might be a good idea to use Stein's algorithm only if both ints are too large. I'm aware the function is already optimized enough, and adding more code would make minification harder, but this could increase performance.

I already have a hybrid implementation here, but it's too complicated because it also uses Lehmer's algo, which I think is overkill for both libraries, so it shouldn't be included. As some source comments point out, computing the GCD of 2 Mersennes can be done calculating the bitlength, which we can get "for free" if we use an iterative isMersenne helper function, instead of the short expression approach my library uses.

If this is approved, I'll make a pull request with the changes. I'll start writing a small draft to post it here

[Rudxain]

This is the draft, I tried to follow the coding style:

export function bitLength(n: bigint): bigint {
  // this expects only positive input args
  // so it's undefined for n < 1, just like the logarithm
  let len = 1n
  while (n > 1n) {
    n >>= 1n
    len++
  }
  return len
}

// count binary trailing zeros
export function ctz(n: bigint, tuple?: boolean): bigint | bigint[] {
  let c = 0n
  while ( (n & 1n) === 0n ) {
    n >>= 1n
    c++
  }
  // `tuple` returns both the count and the input arg with all TZs removed
  return tuple ? [n, c] : c
}

// use Euclidean & Stein's algorithm , if a>=b, getGcd(a,b) = getGcd(b,a mod b)
export function getGcd(a: bigint, b: bigint): bigint {
  // params default value is positive, >=1n
  if (a === b) return a
  let i: bigint, j: bigint, k: bigint
  ;[a, i] = ctz(a, true) // I get 1 TS error here for some reason: "no iterator"
  ;[b, j] = ctz(b, true)
  k = i < j ? i : j // min

  if (a < b) [a, b] = [b, a]
  if (a - b === 1n) return 1n // IDK if this is necessary

  const LIMIT = 1n << 16n // this is arbitrary, we need benchmarks to find the "sweet spot"
  if (bitLength(a) < LIMIT || bitLength(b) < LIMIT) {
    do // if any or both are small, do Euclid
      [a, b] = [b, a % b]
    while (b)
    return a << k
  } else {
    // otherwise, use Stein's
    while (true)
    {
      if (a > b) [a, b] = [b, a]
      b -= a
      if (!b) return a << k
      b = ctz(b, true)[0] // remove all TZs
    }
  }
}

[Yaffle]

GCD should be avoided if possible as it is quite slow,
I am looking how you are using it, and seems, you could simply test if a * 10**e (where e is a large enough) divides b .

[Rudxain]

GCD should be avoided if possible as it is quite slow

What do you mean? Are you referring to the Euclidean algorithm or Euclid's algorithm? (According to Wikipedia, they should not be confused, because Euclid's uses repeated subtraction, while Euclidean uses modulo). I know computing GCDs with any algorithm is usually slower than division. I just wanted to contribute to the function implementation, in case it's needed somewhere (IDK where).

You could simply test if a * 10**e (where e is a large enough) divides b

Do you mean Lehmer's algo? I didn't know that, thanks for the info! That algorithm has so many variants, and it's so complicated that I don't know what I'm doing lol

[Yaffle]

Ah, sorry, I thought, you are the author,
for this particular library the only usage is at https://github.com/yukinotech/JSBD/blob/main/src/jsbd.ts#L100 and it can be avoided may be.
Btw, I have another library with GCD - https://github.com/Yaffle/bigint-gcd - the README has some performance numbers.

[Rudxain]

Oh, it's ok I understand. And I visited your repo and it's interesting and useful! Although I didn't understand some stuff, I'm still learning about the math behind different GCD algos