# lib/fixed64.bend -- Q8.56 signed fixed point on two U32 limbs.
#
# An Fx is a 64-bit two's-complement integer (hi:lo) read as value / 2^56:
# range about +-128, resolution 2^-56 (1.4e-17), against F32's 24-bit
# mantissa. It is built only from U32 operations the GPU already has:
# 32x32 products are assembled from 16-bit halves (Fx.mul_hi), 64-bit
# carries come from `sum < addend`, and the escape test of a Mandelbrot
# loop is one comparison on the high limb. demos/mandelbrot/deep uses it
# to zoom past the point where F32 cannot even represent the centre.
#
# Fx.mul truncates toward zero. Fx.div10 and Fx.shr are for non-negative
# values (the decimal parser and pixel geometry). The specification of
# every operation against the exact integers is stated as open claims in
# lib/claims/fixed64.bend, and tested against Python in lib/test_fixed64.py.

import Base

type Fx is Data:
  Fx{+hi: U32, +lo: U32}

def Fx.b2u(b: Bool) -> U32:
  match b:
    case True{}:
      1
    case False{}:
      0

def Fx.zero() -> Fx:
  Fx{0, 0}

# k as an Fx, for 0 <= k < 128
def Fx.from_int(k: U32) -> Fx:
  Fx{(k << 24n : U32), 0}

def Fx.hi(f: Fx) -> U32:
  match f:
    case Fx{hi, lo}:
      hi

def Fx.lo(f: Fx) -> U32:
  match f:
    case Fx{hi, lo}:
      lo

# two's complement on two limbs: ~x + 1, the +1 carrying into hi iff lo was 0
def Fx.neg_hi(+hi: U32, +lo: U32) -> U32:
  (U32.not(hi) + Fx.b2u(U32.is_eq(lo, 0)) : U32)

def Fx.neg_lo(lo: U32) -> U32:
  (U32.not(lo) + 1 : U32)

def Fx.neg(f: Fx) -> Fx:
  match f:
    case Fx{hi, lo}:
      Fx{Fx.neg_hi(hi, lo), Fx.neg_lo(lo)}

def Fx.add(a: Fx, b: Fx) -> Fx:
  match a b:
    case Fx{ah, al} Fx{bh, bl}:
      +lo = (al + bl : U32)
      Fx{(ah + bh + Fx.b2u(lo < al) : U32), lo}

def Fx.sub(a: Fx, b: Fx) -> Fx:
  Fx.add(a, Fx.neg(b))

def Fx.is_neg(f: Fx) -> Bool:
  match f:
    case Fx{hi, lo}:
      U32.is_eq((hi >> 31n : U32), 1)

# high 32 bits of the 64-bit product x * y (the low 32 are x * y itself)
def Fx.mul_hi(+x: U32, +y: U32) -> U32:
  +xh = (x >> 16n : U32)
  +xl = U32.and(x, 65535)
  +yh = (y >> 16n : U32)
  +yl = U32.and(y, 65535)
  +ll = (xl * yl : U32)
  +lh = (xl * yh : U32)
  +hl = (xh * yl : U32)
  +mid = ((ll >> 16n) + U32.and(lh, 65535) + U32.and(hl, 65535) : U32)
  (xh * yh + (lh >> 16n) + (hl >> 16n) + (mid >> 16n) : U32)

# unsigned (ah:al) * (bh:bl), a 128-bit product in limbs p0..p3, of which
# bits 56..119 are the Q8.56 result. p0 is never needed: nothing carries
# out of it.
def Fx.umul(+ah: U32, +al: U32, +bh: U32, +bl: U32) -> Fx:
  +h0 = Fx.mul_hi(al, bl)
  +s1 = (h0 + al * bh : U32)
  +c1 = (Fx.b2u(s1 < h0) : U32)
  +p1 = (s1 + ah * bl : U32)
  +c2 = (Fx.b2u(p1 < s1) : U32)
  +m1 = Fx.mul_hi(al, bh)
  +t1 = (m1 + Fx.mul_hi(ah, bl) : U32)
  +d1 = (Fx.b2u(t1 < m1) : U32)
  +t2 = (t1 + ah * bh : U32)
  +d2 = (Fx.b2u(t2 < t1) : U32)
  +p2 = (t2 + c1 + c2 : U32)
  +d3 = (Fx.b2u(p2 < t2) : U32)
  +p3 = (Fx.mul_hi(ah, bh) + d1 + d2 + d3 : U32)
  Fx{(U32.or(p2 >> 24n, p3 << 8n) : U32), (U32.or(p1 >> 24n, p2 << 8n) : U32)}

def Fx.mul_signs(+ah: U32, +al: U32, +bh: U32, +bl: U32, na: Bool, nb: Bool) -> Fx:
  match na nb:
    case False{} False{}:
      Fx.umul(ah, al, bh, bl)
    case True{} False{}:
      Fx.neg(Fx.umul(Fx.neg_hi(ah, al), Fx.neg_lo(al), bh, bl))
    case False{} True{}:
      Fx.neg(Fx.umul(ah, al, Fx.neg_hi(bh, bl), Fx.neg_lo(bl)))
    case True{} True{}:
      Fx.umul(Fx.neg_hi(ah, al), Fx.neg_lo(al), Fx.neg_hi(bh, bl), Fx.neg_lo(bl))

def Fx.mul(a: Fx, b: Fx) -> Fx:
  match a b:
    case Fx{ah, al} Fx{bh, bl}:
      Fx.mul_signs(ah, al, bh, bl, U32.is_eq((ah >> 31n : U32), 1), U32.is_eq((bh >> 31n : U32), 1))

def Fx.sq(+a: Fx) -> Fx:
  Fx.mul(a, a)

# f * d for a small non-negative integer d (any U32; the product must fit)
def Fx.mul_small(f: Fx, +d: U32) -> Fx:
  match f:
    case Fx{hi, lo}:
      Fx{(hi * d + Fx.mul_hi(lo, d) : U32), (lo * d : U32)}

# f / 2^n for non-negative f; m must be 32 - n
def Fx.shr(f: Fx, +n: Nat, +m: Nat) -> Fx:
  match f:
    case Fx{hi, lo}:
      Fx{(hi >> n : U32), (U32.or(lo >> n, hi << m) : U32)}

# f / 10 for non-negative f: long division in base 2^16
def Fx.div10(f: Fx) -> Fx:
  match f:
    case Fx{hi, lo}:
      +c3 = (hi >> 16n : U32)
      +c2 = U32.and(hi, 65535)
      +c1 = (lo >> 16n : U32)
      +c0 = U32.and(lo, 65535)
      +q3 = (c3 / 10 : U32)
      +x2 = ((c3 % 10) * 65536 + c2 : U32)
      +q2 = (x2 / 10 : U32)
      +x1 = ((x2 % 10) * 65536 + c1 : U32)
      +q1 = (x1 / 10 : U32)
      +x0 = ((x1 % 10) * 65536 + c0 : U32)
      +q0 = (x0 / 10 : U32)
      Fx{((q3 << 16n) + q2 : U32), ((q1 << 16n) + q0 : U32)}

# branchless select: m must be 0 or 1; answers a when m = 1, b when m = 0
def Fx.sel(+m: U32, a: Fx, b: Fx) -> Fx:
  match a b:
    case Fx{ah, al} Fx{bh, bl}:
      Fx{(m * ah + (1 - m) * bh : U32), (m * al + (1 - m) * bl : U32)}

# is a non-negative f greater than 4? (one comparison on the high limb)
def Fx.gt4(f: Fx) -> Bool:
  match f:
    case Fx{hi, lo}:
      U32.is_gt(hi, 67108864)

# to F32, for comparisons: hi / 2^24 + lo / 2^56 on the magnitude
def Fx.to_f32.mag(+hi: U32, +lo: U32) -> F32:
  F32.add(F32.mul(U32.to_f32(hi), 0.000000059604645), F32.mul(U32.to_f32(lo), 0.000000000000000013877788))

def Fx.to_f32.go(+hi: U32, +lo: U32, neg: Bool) -> F32:
  match neg:
    case False{}:
      Fx.to_f32.mag(hi, lo)
    case True{}:
      F32.neg(Fx.to_f32.mag(Fx.neg_hi(hi, lo), Fx.neg_lo(lo)))

def Fx.to_f32(f: Fx) -> F32:
  match f:
    case Fx{hi, lo}:
      Fx.to_f32.go(hi, lo, U32.is_eq((hi >> 31n : U32), 1))

# ---------------------------------------------------------------------
# decimal parsing: "-0.743643887037151" -> Fx. Digits beyond the
# resolution are truncated. The parser is a fold over the string with a
# small state record, so that the branching helpers need not recurse.
# ---------------------------------------------------------------------

type Fx.Dec is Data:
  Dec{acc: Fx, scale: Fx, frac: Bool, neg: Bool}

def Fx.dec.digit(acc: Fx, scale: Fx, frac: Bool, neg: Bool, +d: U32) -> Fx.Dec:
  match frac:
    case False{}:
      Dec{Fx.add(Fx.mul_small(acc, 10), Fx.from_int(d)), scale, False{}, neg}
    case True{}:
      +sc = scale
      Dec{Fx.add(acc, Fx.mul_small(sc, d)), Fx.div10(sc), True{}, neg}

def Fx.dec.step2(acc: Fx, scale: Fx, frac: Bool, neg: Bool, +d: U32, minus: Bool, dot: Bool) -> Fx.Dec:
  match minus dot:
    case True{} True{}:
      Dec{acc, scale, frac, neg}
    case True{} False{}:
      Dec{acc, scale, frac, True{}}
    case False{} True{}:
      Dec{acc, Fx.div10(Fx.from_int(1)), True{}, neg}
    case False{} False{}:
      Fx.dec.digit(acc, scale, frac, neg, d)

def Fx.dec.step(st: Fx.Dec, +c: U32) -> Fx.Dec:
  match st:
    case Dec{acc, scale, frac, neg}:
      Fx.dec.step2(acc, scale, frac, neg, (c - 48 : U32), U32.is_eq(c, 45), U32.is_eq(c, 46))

def Fx.dec.fin(st: Fx.Dec) -> Fx:
  match st:
    case Dec{acc, scale, frac, neg}:
      match neg:
        case False{}:
          acc
        case True{}:
          Fx.neg(acc)

def Fx.dec.go(s: String, st: Fx.Dec) -> Fx:
  match s:
    case SNil{}:
      Fx.dec.fin(st)
    case SCon{Chr{c}, t}:
      Fx.dec.go(t, Fx.dec.step(st, c))

def Fx.from_dec(s: String) -> Fx:
  Fx.dec.go(s, Dec{Fx.zero(), Fx.zero(), False{}, False{}})