Javascript elliptical point multiplication algorit

Question

The following is an implementation of elliptical curve Point Multiplication, but it's not working as expected (using recent Chrome / Node with BigInt for illustration):

const bi0 = BigInt(0)
const bi1 = BigInt(1)
const bi2 = BigInt(2)
const bi3 = BigInt(3)

const absMod = (n, p) => n < bi0 ? (n % p) + p : n % p

export function pointAdd (xp, yp, xq, yq, p) {
  const lambda = (yq - yp) / (xq - xp)
  const x = absMod(lambda ** bi2 - xp - xq, p)
  const y = absMod(lambda * (xp - x) - yp, p)
  return { x, y }
}

export function pointDouble (xp, yp, a, p) {
  const numer = bi3 * xp ** bi2 + a
  const denom = (bi2 * yp) ** (p - bi2)
  const lambda = (numer * denom) % p
  const x = absMod(lambda ** bi2 - bi2 * xp, p)
  const y = absMod(lambda * (xp - x) - yp, p)
  return { x, y }
}

export function pointMultiply (d, xp, yp, a, p) {
  const add = (xp, yp, { x, y }) => pointAdd(xp, yp, x, y, p)
  const double = (x, y) => pointDouble(x, y, a, p)
  const recur = ({ x, y }, n) => {
    if (n === bi0) { return { x: bi0, y: bi0 } }
    if (n === bi1) { return { x, y } }
    if (n % bi2 === bi1) { return add(x, y, recur({ x, y }, n - bi1)) }
    return recur(double(x, y), n / bi2)
  }
  return recur({ x: xp, y: yp }, d)
}

Given a known curve with properties, the above succeeds for points P2 - P5, but starts failing at P6 onwards:

const p = BigInt('17')
const a = BigInt('2')
const p1 = { x: BigInt(5), y: BigInt(1) }

const d = BigInt(6)
const p6 = pointMultiply(d, p1.x, p1.y, a, p)
p6.x === BigInt(16) // incorrect value of 8 was returned
p6.y === BigInt(13) // incorrect value of 14 was returned

The known curve has values:

P   X    Y
—————————— 
1   5    1
2   6    3
3  10    6
4   3    1
5   9   16
6  16   13
7   0    6
8  13    7
9   7    6
10  7   11

I'm not sure if I misunderstand the algorithm or I've made an error in the implementation.

Answer 1

I don't know javascript very well, so the code confuses me. But ...

In function pointAdd, and everywhere else, "division" must be done mod p. You're evidently doing it correctly in pointDouble, but not in pointAdd. In pointAdd, use the same pattern: Instead of

const lambda = (yq - yp) / (xq - xp)

do

const numer = yq - yp
const denom = (xq - xp) ** (p - bi2)
const lambda = (numer * denom) % p

though I would think it would much clearer and less error prone to simply have a modular inverse function instead of computing X^p-2 everywhere you need it.