# lib/fx4.bend -- Q8.120 signed fixed point on four flat U32 limbs.
#
# GENERATED by lib/gen_fx4.py -- edit the generator, not this file.
#
# An Fx4 is a 128-bit two's-complement integer (a:b:c:d, a most
# significant) read as value / 2^120: range about +-128, resolution
# 2^-120 (7.5e-37), against lib/fixed64.bend's 2^-56.  That is a factor
# of 2^64 more precision, and it moves the Mandelbrot zoom wall from
# about 1e14 to about 1e33.
#
# The point of the FLAT record is the GPU.  lib/bigfix.bend does the same
# arithmetic on a list of limbs, and is correct, but it allocates a cons
# cell per limb per operation, so it runs on one core and the Metal
# watchdog kills it.  Fx4 has no allocation and no recursion in its
# arithmetic: every operation is straight-line U32, and the 16 partial
# products of the multiply are unrolled.
#
# Fx4.mul truncates toward zero.  Fx4.div10 and Fx4.shr are for
# non-negative values (the decimal parser and pixel geometry).
# Tested against Python big integers by lib/test_fx4.py.

import Base

type Fx4 is Data:
  Fx4{+a: U32, +b: U32, +c: U32, +d: U32}

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

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

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

# high 32 bits of the 64-bit product x * y (the low 32 are x * y itself)
def Fx4.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)

# two's complement on four limbs: ~x + 1.  The +1 carries out of limb i
# exactly when every limb below it, and limb i itself, was zero.
def Fx4.neg(x: Fx4) -> Fx4:
  match x:
    case Fx4{xa, xb, xc, xd}:
      +kd = Fx4.b2u(U32.is_eq(xd, 0))
      +kc = (kd * Fx4.b2u(U32.is_eq(xc, 0)) : U32)
      +kb = (kc * Fx4.b2u(U32.is_eq(xb, 0)) : U32)
      Fx4{(U32.not(xa) + kb : U32), (U32.not(xb) + kc : U32), (U32.not(xc) + kd : U32), (U32.not(xd) + 1 : U32)}

# ripple carry, LSB first.  A carry is `sum < addend`; adding the
# incoming carry can itself overflow, but never both at once (the two
# would need a + b >= 2^33 - 1).
def Fx4.add(x: Fx4, y: Fx4) -> Fx4:
  match x y:
    case Fx4{xa, xb, xc, xd} Fx4{ya, yb, yc, yd}:
      +rd = (xd + yd : U32)
      +k0 = (Fx4.b2u(rd < xd) : U32)
      +t1 = (xc + yc : U32)
      +rc = (t1 + k0 : U32)
      +k1 = (Fx4.b2u(t1 < xc) + Fx4.b2u(rc < t1) : U32)
      +t2 = (xb + yb : U32)
      +rb = (t2 + k1 : U32)
      +k2 = (Fx4.b2u(t2 < xb) + Fx4.b2u(rb < t2) : U32)
      Fx4{(xa + ya + k2 : U32), rb, rc, rd}

def Fx4.sub(x: Fx4, y: Fx4) -> Fx4:
  Fx4.add(x, Fx4.neg(y))

def Fx4.is_neg(x: Fx4) -> Bool:
  match x:
    case Fx4{xa, xb, xc, xd}:
      U32.is_eq((xa >> 31n : U32), 1)

# unsigned (a:b:c:d) * (a:b:c:d): the 256-bit product accumulated column
# by column (comba), r0 the column word and r1, r2 its carries.  Only
# columns 3..7 are kept -- the Q8.120 result is bits 120..247 -- but the
# lower columns must still be accumulated for their carries.
def Fx4.umul(x: Fx4, y: Fx4) -> Fx4:
  match x y:
    case Fx4{xa, xb, xc, xd} Fx4{ya, yb, yc, yd}:
      +lo1 = (xd * yd : U32)
      +hi2 = Fx4.mul_hi(xd, yd)
      +s3 = (0 + lo1 : U32)
      +y4 = (Fx4.b2u(s3 < lo1) : U32)
      +s5 = (0 + y4 : U32)
      +y6 = (Fx4.b2u(s5 < y4) : U32)
      +s7 = (s5 + hi2 : U32)
      +y8 = (Fx4.b2u(s7 < hi2) : U32)
      +s9 = (0 + y6 + y8 : U32)
      +lo10 = (xd * yc : U32)
      +hi11 = Fx4.mul_hi(xd, yc)
      +s12 = (s7 + lo10 : U32)
      +y13 = (Fx4.b2u(s12 < lo10) : U32)
      +s14 = (s9 + y13 : U32)
      +y15 = (Fx4.b2u(s14 < y13) : U32)
      +s16 = (s14 + hi11 : U32)
      +y17 = (Fx4.b2u(s16 < hi11) : U32)
      +s18 = (0 + y15 + y17 : U32)
      +lo19 = (xc * yd : U32)
      +hi20 = Fx4.mul_hi(xc, yd)
      +s21 = (s12 + lo19 : U32)
      +y22 = (Fx4.b2u(s21 < lo19) : U32)
      +s23 = (s16 + y22 : U32)
      +y24 = (Fx4.b2u(s23 < y22) : U32)
      +s25 = (s23 + hi20 : U32)
      +y26 = (Fx4.b2u(s25 < hi20) : U32)
      +s27 = (s18 + y24 + y26 : U32)
      +lo28 = (xd * yb : U32)
      +hi29 = Fx4.mul_hi(xd, yb)
      +s30 = (s25 + lo28 : U32)
      +y31 = (Fx4.b2u(s30 < lo28) : U32)
      +s32 = (s27 + y31 : U32)
      +y33 = (Fx4.b2u(s32 < y31) : U32)
      +s34 = (s32 + hi29 : U32)
      +y35 = (Fx4.b2u(s34 < hi29) : U32)
      +s36 = (0 + y33 + y35 : U32)
      +lo37 = (xc * yc : U32)
      +hi38 = Fx4.mul_hi(xc, yc)
      +s39 = (s30 + lo37 : U32)
      +y40 = (Fx4.b2u(s39 < lo37) : U32)
      +s41 = (s34 + y40 : U32)
      +y42 = (Fx4.b2u(s41 < y40) : U32)
      +s43 = (s41 + hi38 : U32)
      +y44 = (Fx4.b2u(s43 < hi38) : U32)
      +s45 = (s36 + y42 + y44 : U32)
      +lo46 = (xb * yd : U32)
      +hi47 = Fx4.mul_hi(xb, yd)
      +s48 = (s39 + lo46 : U32)
      +y49 = (Fx4.b2u(s48 < lo46) : U32)
      +s50 = (s43 + y49 : U32)
      +y51 = (Fx4.b2u(s50 < y49) : U32)
      +s52 = (s50 + hi47 : U32)
      +y53 = (Fx4.b2u(s52 < hi47) : U32)
      +s54 = (s45 + y51 + y53 : U32)
      +lo55 = (xd * ya : U32)
      +hi56 = Fx4.mul_hi(xd, ya)
      +s57 = (s52 + lo55 : U32)
      +y58 = (Fx4.b2u(s57 < lo55) : U32)
      +s59 = (s54 + y58 : U32)
      +y60 = (Fx4.b2u(s59 < y58) : U32)
      +s61 = (s59 + hi56 : U32)
      +y62 = (Fx4.b2u(s61 < hi56) : U32)
      +s63 = (0 + y60 + y62 : U32)
      +lo64 = (xc * yb : U32)
      +hi65 = Fx4.mul_hi(xc, yb)
      +s66 = (s57 + lo64 : U32)
      +y67 = (Fx4.b2u(s66 < lo64) : U32)
      +s68 = (s61 + y67 : U32)
      +y69 = (Fx4.b2u(s68 < y67) : U32)
      +s70 = (s68 + hi65 : U32)
      +y71 = (Fx4.b2u(s70 < hi65) : U32)
      +s72 = (s63 + y69 + y71 : U32)
      +lo73 = (xb * yc : U32)
      +hi74 = Fx4.mul_hi(xb, yc)
      +s75 = (s66 + lo73 : U32)
      +y76 = (Fx4.b2u(s75 < lo73) : U32)
      +s77 = (s70 + y76 : U32)
      +y78 = (Fx4.b2u(s77 < y76) : U32)
      +s79 = (s77 + hi74 : U32)
      +y80 = (Fx4.b2u(s79 < hi74) : U32)
      +s81 = (s72 + y78 + y80 : U32)
      +lo82 = (xa * yd : U32)
      +hi83 = Fx4.mul_hi(xa, yd)
      +s84 = (s75 + lo82 : U32)
      +y85 = (Fx4.b2u(s84 < lo82) : U32)
      +s86 = (s79 + y85 : U32)
      +y87 = (Fx4.b2u(s86 < y85) : U32)
      +s88 = (s86 + hi83 : U32)
      +y89 = (Fx4.b2u(s88 < hi83) : U32)
      +s90 = (s81 + y87 + y89 : U32)
      +lo91 = (xc * ya : U32)
      +hi92 = Fx4.mul_hi(xc, ya)
      +s93 = (s88 + lo91 : U32)
      +y94 = (Fx4.b2u(s93 < lo91) : U32)
      +s95 = (s90 + y94 : U32)
      +y96 = (Fx4.b2u(s95 < y94) : U32)
      +s97 = (s95 + hi92 : U32)
      +y98 = (Fx4.b2u(s97 < hi92) : U32)
      +s99 = (0 + y96 + y98 : U32)
      +lo100 = (xb * yb : U32)
      +hi101 = Fx4.mul_hi(xb, yb)
      +s102 = (s93 + lo100 : U32)
      +y103 = (Fx4.b2u(s102 < lo100) : U32)
      +s104 = (s97 + y103 : U32)
      +y105 = (Fx4.b2u(s104 < y103) : U32)
      +s106 = (s104 + hi101 : U32)
      +y107 = (Fx4.b2u(s106 < hi101) : U32)
      +s108 = (s99 + y105 + y107 : U32)
      +lo109 = (xa * yc : U32)
      +hi110 = Fx4.mul_hi(xa, yc)
      +s111 = (s102 + lo109 : U32)
      +y112 = (Fx4.b2u(s111 < lo109) : U32)
      +s113 = (s106 + y112 : U32)
      +y114 = (Fx4.b2u(s113 < y112) : U32)
      +s115 = (s113 + hi110 : U32)
      +y116 = (Fx4.b2u(s115 < hi110) : U32)
      +s117 = (s108 + y114 + y116 : U32)
      +lo118 = (xb * ya : U32)
      +hi119 = Fx4.mul_hi(xb, ya)
      +s120 = (s115 + lo118 : U32)
      +y121 = (Fx4.b2u(s120 < lo118) : U32)
      +s122 = (s117 + y121 : U32)
      +y123 = (Fx4.b2u(s122 < y121) : U32)
      +s124 = (s122 + hi119 : U32)
      +y125 = (Fx4.b2u(s124 < hi119) : U32)
      +s126 = (0 + y123 + y125 : U32)
      +lo127 = (xa * yb : U32)
      +hi128 = Fx4.mul_hi(xa, yb)
      +s129 = (s120 + lo127 : U32)
      +y130 = (Fx4.b2u(s129 < lo127) : U32)
      +s131 = (s124 + y130 : U32)
      +y132 = (Fx4.b2u(s131 < y130) : U32)
      +s133 = (s131 + hi128 : U32)
      +y134 = (Fx4.b2u(s133 < hi128) : U32)
      +s135 = (s126 + y132 + y134 : U32)
      +lo136 = (xa * ya : U32)
      +hi137 = Fx4.mul_hi(xa, ya)
      +s138 = (s133 + lo136 : U32)
      +y139 = (Fx4.b2u(s138 < lo136) : U32)
      +s140 = (s135 + y139 : U32)
      +y141 = (Fx4.b2u(s140 < y139) : U32)
      +s142 = (s140 + hi137 : U32)
      +y143 = (Fx4.b2u(s142 < hi137) : U32)
      +s144 = (0 + y141 + y143 : U32)
      Fx4{(U32.or(s138 >> 24n, s142 << 8n) : U32), (U32.or(s129 >> 24n, s138 << 8n) : U32), (U32.or(s111 >> 24n, s129 << 8n) : U32), (U32.or(s84 >> 24n, s111 << 8n) : U32)}

# signs handled on the magnitudes, so the truncation is toward zero
def Fx4.mul_signs(x: Fx4, y: Fx4, nx: Bool, ny: Bool) -> Fx4:
  match nx ny:
    case False{} False{}:
      Fx4.umul(x, y)
    case True{} False{}:
      Fx4.neg(Fx4.umul(Fx4.neg(x), y))
    case False{} True{}:
      Fx4.neg(Fx4.umul(x, Fx4.neg(y)))
    case True{} True{}:
      Fx4.umul(Fx4.neg(x), Fx4.neg(y))

def Fx4.sign(x: Fx4) -> Bool:
  match x:
    case Fx4{xa, xb, xc, xd}:
      U32.is_eq((xa >> 31n : U32), 1)

def Fx4.mul(+x: Fx4, +y: Fx4) -> Fx4:
  Fx4.mul_signs(x, y, Fx4.sign(x), Fx4.sign(y))

def Fx4.sq(+x: Fx4) -> Fx4:
  Fx4.mul(x, x)

# x * m for a small non-negative integer m (the product must fit)
def Fx4.mul_small(x: Fx4, +m: U32) -> Fx4:
  match x:
    case Fx4{xa, xb, xc, xd}:
      +rd = (xd * m : U32)
      +k0 = Fx4.mul_hi(xd, m)
      +t1 = (xc * m : U32)
      +rc = (t1 + k0 : U32)
      +k1 = (Fx4.mul_hi(xc, m) + Fx4.b2u(rc < t1) : U32)
      +t2 = (xb * m : U32)
      +rb = (t2 + k1 : U32)
      +k2 = (Fx4.mul_hi(xb, m) + Fx4.b2u(rb < t2) : U32)
      Fx4{(xa * m + k2 : U32), rb, rc, rd}

# x / 2^n for non-negative x; m must be 32 - n
def Fx4.shr(x: Fx4, +n: Nat, +m: Nat) -> Fx4:
  match x:
    case Fx4{xa, xb, xc, xd}:
      Fx4{(xa >> n : U32), (U32.or(xb >> n, xa << m) : U32), (U32.or(xc >> n, xb << m) : U32), (U32.or(xd >> n, xc << m) : U32)}

# x / 10 for non-negative x: long division in base 2^16, eight digits
def Fx4.div10(x: Fx4) -> Fx4:
  match x:
    case Fx4{xa, xb, xc, xd}:
      +h0 = (xa >> 16n : U32)
      +h1 = U32.and(xa, 65535)
      +h2 = (xb >> 16n : U32)
      +h3 = U32.and(xb, 65535)
      +h4 = (xc >> 16n : U32)
      +h5 = U32.and(xc, 65535)
      +h6 = (xd >> 16n : U32)
      +h7 = U32.and(xd, 65535)
      +q0 = (h0 / 10 : U32)
      +r0 = (h0 % 10 : U32)
      +x1 = (r0 * 65536 + h1 : U32)
      +q1 = (x1 / 10 : U32)
      +r1 = (x1 % 10 : U32)
      +x2 = (r1 * 65536 + h2 : U32)
      +q2 = (x2 / 10 : U32)
      +r2 = (x2 % 10 : U32)
      +x3 = (r2 * 65536 + h3 : U32)
      +q3 = (x3 / 10 : U32)
      +r3 = (x3 % 10 : U32)
      +x4 = (r3 * 65536 + h4 : U32)
      +q4 = (x4 / 10 : U32)
      +r4 = (x4 % 10 : U32)
      +x5 = (r4 * 65536 + h5 : U32)
      +q5 = (x5 / 10 : U32)
      +r5 = (x5 % 10 : U32)
      +x6 = (r5 * 65536 + h6 : U32)
      +q6 = (x6 / 10 : U32)
      +r6 = (x6 % 10 : U32)
      +x7 = (r6 * 65536 + h7 : U32)
      +q7 = (x7 / 10 : U32)
      +r7 = (x7 % 10 : U32)
      Fx4{((q0 << 16n) + q1 : U32), ((q2 << 16n) + q3 : U32), ((q4 << 16n) + q5 : U32), ((q6 << 16n) + q7 : U32)}

# branchless select: m must be 0 or 1; answers x when m = 1, y when m = 0
def Fx4.sel(+m: U32, x: Fx4, y: Fx4) -> Fx4:
  match x y:
    case Fx4{xa, xb, xc, xd} Fx4{ya, yb, yc, yd}:
      +i = (1 - m : U32)
      Fx4{(m * xa + i * ya : U32), (m * xb + i * yb : U32), (m * xc + i * yc : U32), (m * xd + i * yd : U32)}

# is a non-negative x greater than 4? (one comparison on the top limb:
# the integer part is the top 8 bits, exactly as in Q8.56)
def Fx4.gt4(x: Fx4) -> Bool:
  match x:
    case Fx4{xa, xb, xc, xd}:
      U32.is_gt(xa, 67108864)

# to F32, for comparisons: only the top two limbs can matter
def Fx4.to_f32.mag(+xa: U32, +xb: U32) -> F32:
  F32.add(F32.mul(U32.to_f32(xa), 0.000000059604645), F32.mul(U32.to_f32(xb), 0.000000000000000013877788))

def Fx4.to_f32.hi(x: Fx4) -> F32:
  match x:
    case Fx4{xa, xb, xc, xd}:
      Fx4.to_f32.mag(xa, xb)

def Fx4.to_f32.go(x: Fx4, neg: Bool) -> F32:
  match neg:
    case False{}:
      Fx4.to_f32.hi(x)
    case True{}:
      F32.neg(Fx4.to_f32.hi(Fx4.neg(x)))

def Fx4.to_f32(+x: Fx4) -> F32:
  Fx4.to_f32.go(x, Fx4.sign(x))

# decimal parsing: "-0.7436438870371587047521915061147" -> Fx4, to as
# many digits as the type holds (about 36).  A fold over the string with
# a small state record, exactly as lib/fixed64.bend does it.
type Fx4.Dec is Data:
  Dec4{acc: Fx4, scale: Fx4, frac: Bool, neg: Bool}

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

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

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

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

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

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