{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE UnboxedTuples #-}

module Moonlight.Triangulation.Internal.Dyadic
  ( exactOrientDet
  , exactOrientSignDouble
  , exactInCircleDet
  , exactBarycentricDeterminants
  , exactDiametralDot
  , integerRatioToDouble
  ) where

import Data.Bits (shiftL, shiftR)
#if WORD_SIZE_IN_BITS == 64
import GHC.Exts
  ( Double (D#)
  , Double#
  , Int#
  , Word#
  , and#
  , castDoubleToWord64#
  , eqWord#
  , gtWord#
  , int2Word#
  , isTrue#
  , or#
  , plusWord#
  , plusWord2#
  , subWordC#
  , timesWord2#
  , uncheckedShiftL#
  , uncheckedShiftRL#
  , word2Int#
  , word64ToWord#
  , (*#)
  , (+#)
  , (-#)
  , (<=#)
  , (==#)
  , (>#)
  , (>=#)
  )
#endif

-- Every finite binary64 value is a dyadic rational. Aligning all mantissas to
-- one exponent gives exact integer predicates without constructing Rational
-- expression trees.
type Decoded = (Integer, Int)

commonExponent :: [Decoded] -> Int
commonExponent = foldl' step 0
 where
  step :: Int -> Decoded -> Int
  step !current (!mantissa, !power)
    | mantissa == 0 = current
    | otherwise = min current power
{-# INLINE commonExponent #-}

alignDecoded :: Int -> Decoded -> Integer
alignDecoded !power (!mantissa, !sourcePower)
  | mantissa == 0 = 0
  | otherwise = mantissa `shiftL` (sourcePower - power)
{-# INLINE alignDecoded #-}

aligned6
  :: Double -> Double -> Double -> Double -> Double -> Double
  -> (Integer, Integer, Integer, Integer, Integer, Integer)
aligned6 a b c d e f =
  let !da = decodeFloat a
      !db = decodeFloat b
      !dc = decodeFloat c
      !dd = decodeFloat d
      !de = decodeFloat e
      !df = decodeFloat f
      !power = commonExponent [da, db, dc, dd, de, df]
   in ( alignDecoded power da
      , alignDecoded power db
      , alignDecoded power dc
      , alignDecoded power dd
      , alignDecoded power de
      , alignDecoded power df
      )

aligned8
  :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
  -> (Integer, Integer, Integer, Integer, Integer, Integer, Integer, Integer)
aligned8 a b c d e f g h =
  let !da = decodeFloat a
      !db = decodeFloat b
      !dc = decodeFloat c
      !dd = decodeFloat d
      !de = decodeFloat e
      !df = decodeFloat f
      !dg = decodeFloat g
      !dh = decodeFloat h
      !power = commonExponent [da, db, dc, dd, de, df, dg, dh]
   in ( alignDecoded power da
      , alignDecoded power db
      , alignDecoded power dc
      , alignDecoded power dd
      , alignDecoded power de
      , alignDecoded power df
      , alignDecoded power dg
      , alignDecoded power dh
      )

exactOrientDet
  :: Double -> Double -> Double -> Double -> Double -> Double -> Integer
exactOrientDet ax ay bx by cx cy =
  let (!iax, !iay, !ibx, !iby, !icx, !icy) = aligned6 ax ay bx by cx cy
      !acx = iax - icx
      !acy = iay - icy
      !bcx = ibx - icx
      !bcy = iby - icy
   in acx * bcy - acy * bcx

exactInCircleDet
  :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
  -> Integer
exactInCircleDet ax ay bx by cx cy dx dy =
  let (!iax, !iay, !ibx, !iby, !icx, !icy, !idx, !idy) =
        aligned8 ax ay bx by cx cy dx dy
      !adx = iax - idx
      !ady = iay - idy
      !bdx = ibx - idx
      !bdy = iby - idy
      !cdx = icx - idx
      !cdy = icy - idy
      !abdet = adx * bdy - bdx * ady
      !bcdet = bdx * cdy - cdx * bdy
      !cadet = cdx * ady - adx * cdy
      !alift = adx * adx + ady * ady
      !blift = bdx * bdx + bdy * bdy
      !clift = cdx * cdx + cdy * cdy
   in alift * bcdet + blift * cadet + clift * abdet

exactDiametralDot
  :: Double -> Double -> Double -> Double -> Double -> Double -> Integer
exactDiametralDot ax ay bx by px py =
  let (!iax, !iay, !ibx, !iby, !ipx, !ipy) = aligned6 ax ay bx by px py
      !pax = ipx - iax
      !pay = ipy - iay
      !pbx = ipx - ibx
      !pby = ipy - iby
   in pax * pbx + pay * pby

exactBarycentricDeterminants
  :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
  -> (Integer, Integer, Integer, Integer)
exactBarycentricDeterminants ax ay bx by cx cy qx qy =
  let (!iax, !iay, !ibx, !iby, !icx, !icy, !iqx, !iqy) =
        aligned8 ax ay bx by cx cy qx qy
      determinant :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer -> Integer
      determinant px py rx ry sx sy =
        let !psx = px - sx
            !psy = py - sy
            !rsx = rx - sx
            !rsy = ry - sy
         in psx * rsy - psy * rsx
      !denominator = determinant iax iay ibx iby icx icy
      !weightA = determinant iqx iqy ibx iby icx icy
      !weightB = determinant iax iay iqx iqy icx icy
      !weightC = determinant iax iay ibx iby iqx iqy
   in (denominator, weightA, weightB, weightC)

integerRatioToDouble :: Integer -> Integer -> Double
integerRatioToDouble numerator denominator
  | denominator == 0 = 0 / 0
  | numerator == 0 = 0
  | otherwise =
      let !precision = floatDigits (0 :: Double)
          !numeratorMagnitude = abs numerator
          !denominatorMagnitude = abs denominator
          !numeratorBits = integerBitLength numeratorMagnitude
          !denominatorBits = integerBitLength denominatorMagnitude
          !numeratorShift = max 0 (numeratorBits - precision)
          !denominatorShift = max 0 (denominatorBits - precision)
          !scaledNumerator = fromInteger (numeratorMagnitude `shiftR` numeratorShift)
          !scaledDenominator = fromInteger (denominatorMagnitude `shiftR` denominatorShift)
          !magnitude = scaleFloat (numeratorShift - denominatorShift) (scaledNumerator / scaledDenominator)
          !sameSign = (numerator < 0) == (denominator < 0)
       in if sameSign then magnitude else negate magnitude

integerBitLength :: Integer -> Int
integerBitLength = go 0
 where
  go !bits value
    | value <= 0xffffffff = bits + wordBitLength value
    | otherwise = go (bits + 32) (value `shiftR` 32)

  wordBitLength = count 0
  count :: Int -> Integer -> Int
  count !bits 0 = bits
  count !bits value = count (bits + 1) (value `shiftR` 1)

-- ---------------------------------------------------------------------------
-- Fixed-precision exact orient sign for Double.
--
-- The generic dyadic path answers every exact query with arbitrary-precision
-- Integers: six decodes, one alignment, and two multiplies, each allocating.
-- Practical inputs have an exponent spread small enough that the determinant's
-- exact sign is decided by 128-bit differences and 256-bit products in machine
-- words, without a single heap object. 'exactOrientSignDouble' takes that path
-- and falls back to 'exactOrientDet' the moment an operand is non-finite or an
-- alignment shift would outgrow the fixed width. The two agree by
-- construction: both compute the sign of the same integer determinant.
--
-- The fixed-width worker reads a Double as one machine word and aligns
-- mantissas across a 128-bit pair, so it is only meaningful where a machine
-- word is 64 bits wide. On a narrower target the same sign is taken from the
-- arbitrary-precision determinant directly, which is the branch this path
-- already falls back to whenever an alignment shift would outgrow the width.

#if WORD_SIZE_IN_BITS == 64

exactOrientSignDouble
  :: Double -> Double -> Double -> Double -> Double -> Double -> Ordering
exactOrientSignDouble ax ay bx by cx cy =
  case ax of
    D# axw ->
      case ay of
        D# ayw ->
          case bx of
            D# bxw ->
              case by of
                D# byw ->
                  case cx of
                    D# cxw ->
                      case cy of
                        D# cyw ->
                          case orientSignWorker axw ayw bxw byw cxw cyw of
                            2# -> compare (exactOrientDet ax ay bx by cx cy) 0
                            0# -> EQ
                            sign ->
                              case sign ># 0# of
                                1# -> GT
                                _ -> LT
{-# NOINLINE exactOrientSignDouble #-}

-- Decode a Double into sign bit (0/1), mantissa, and power-of-two exponent
-- with value = (-1)^sign * mantissa * 2^exponent. Zero decodes to a zero
-- mantissa; subnormals decode without a hidden bit. The fourth component is 1
-- when the value is finite and 0 when it is not.
decodeExact :: Double# -> (# Int#, Word#, Int#, Int# #)
decodeExact d =
  case word64ToWord# (castDoubleToWord64# d) of
    bits ->
      let neg = word2Int# (uncheckedShiftRL# bits 63#)
          exponentField = word2Int# (and# (uncheckedShiftRL# bits 52#) 2047##)
          mantissaField = and# bits 4503599627370495##
       in case exponentField of
            0# -> (# neg, mantissaField, -1074#, 1# #)
            2047# -> (# neg, mantissaField, 0#, 0# #)
            raw -> (# neg, or# mantissaField 4503599627370496##, raw -# 1075#, 1# #)

-- A mantissa of at most 53 bits shifted left by at most 73 bits: the pair
-- (high, low) of a value below 2^126.
shiftMantissa :: Word# -> Int# -> (# Word#, Word# #)
shiftMantissa mantissa k =
  case k >=# 64# of
    1# -> (# uncheckedShiftL# mantissa (k -# 64#), 0## #)
    _ ->
      case k ==# 0# of
        1# -> (# 0##, mantissa #)
        _ ->
          (#
            uncheckedShiftRL# mantissa (64# -# k),
            uncheckedShiftL# mantissa k
          #)

-- The exact sign and 128-bit magnitude of sa*ma*2^ea - sc*mc*2^ec, aligned
-- to the caller-supplied floor exponent, as (sign, high, low, status) with
-- sign in {-1, 0, 1} and status 1 when an alignment shift outgrows the fixed
-- width. The floor never exceeds the exponent of a nonzero operand, so every
-- shift is non-negative; one shared floor is what makes the four differences
-- of one determinant comparable after multiplication.
differenceExact
  :: Int# -> Word# -> Int# -> Int# -> Word# -> Int# -> Int# -> (# Int#, Word#, Word#, Int# #)
differenceExact nega ma ea negc mc ec emin =
  case ma of
    0## ->
      case mc of
        0## -> (# 0#, 0##, 0##, 0# #)
        _ -> aligned (negateSign (positiveSign negc)) mc (ec -# emin)
    _ ->
      case mc of
        0## -> aligned (positiveSign nega) ma (ea -# emin)
        _ ->
          case ea -# emin of
            da ->
              case da ># 73# of
                1# -> (# 0#, 0##, 0##, 1# #)
                _ ->
                  case ec -# emin of
                    dc ->
                      case dc ># 73# of
                        1# -> (# 0#, 0##, 0##, 1# #)
                        _ ->
                          case shiftMantissa ma da of
                            (# ahi, alo #) ->
                              case shiftMantissa mc dc of
                                (# chi, clo #) ->
                                  case nega ==# negc of
                                    1# ->
                                      -- Same operand signs: subtract magnitudes.
                                      case compareWord2 ahi alo chi clo of
                                        0# -> (# 0#, 0##, 0##, 0# #)
                                        1# ->
                                          case subtractWord2 ahi alo chi clo of
                                            (# hi, lo #) -> (# positiveSign nega, hi, lo, 0# #)
                                        _ ->
                                          case subtractWord2 chi clo ahi alo of
                                            (# hi, lo #) -> (# negateSign (positiveSign nega), hi, lo, 0# #)
                                    _ ->
                                      -- Opposite operand signs: add magnitudes.
                                      case addWord2 ahi alo chi clo of
                                        (# hi, lo #) -> (# positiveSign nega, hi, lo, 0# #)
  where
    aligned sign mantissa k =
      case k ># 73# of
        1# -> (# 0#, 0##, 0##, 1# #)
        _ ->
          case shiftMantissa mantissa k of
            (# hi, lo #) -> (# sign, hi, lo, 0# #)
    positiveSign neg = case neg of
      1# -> -1#
      _ -> 1#
    negateSign sign = case sign of
      1# -> -1#
      _ -> 1#

-- Lexicographic comparison of 128-bit magnitudes: 1, 0, or -1.
compareWord2 :: Word# -> Word# -> Word# -> Word# -> Int#
compareWord2 ahi alo chi clo =
  case eqWord# ahi chi of
    1# ->
      case eqWord# alo clo of
        1# -> 0#
        _ ->
          case gtWord# alo clo of
            1# -> 1#
            _ -> -1#
    _ ->
      case gtWord# ahi chi of
        1# -> 1#
        _ -> -1#

-- 128-bit difference of magnitudes, first operand at least the second.
subtractWord2 :: Word# -> Word# -> Word# -> Word# -> (# Word#, Word# #)
subtractWord2 ahi alo chi clo =
  case subWordC# alo clo of
    (# low, borrow #) ->
      case subWordC# ahi chi of
        (# high0, _ #) ->
          case subWordC# high0 (int2Word# borrow) of
            (# high, _ #) -> (# high, low #)

-- 128-bit sum of magnitudes each below 2^126: the total stays below 2^127 and
-- the final carry is empty by the shift bound.
addWord2 :: Word# -> Word# -> Word# -> Word# -> (# Word#, Word# #)
addWord2 ahi alo chi clo =
  case plusWord2# alo clo of
    (# carry0, low #) ->
      case plusWord2# ahi chi of
        (# _, high0 #) ->
          case plusWord2# high0 carry0 of
            (# _, high #) -> (# high, low #)

-- 128-bit by 128-bit exact product, (r3, r2, r1, r0), most significant first.
-- Each factor stays below 2^127, so the product stays below 2^254 and the top
-- accumulation cannot overflow.
multiplyWord2 :: Word# -> Word# -> Word# -> Word# -> (# Word#, Word#, Word#, Word# #)
multiplyWord2 ahi alo bhi blo =
  case timesWord2# alo blo of
    (# h00, l00 #) ->
      case timesWord2# alo bhi of
        (# h01, l01 #) ->
          case timesWord2# ahi blo of
            (# h10, l10 #) ->
              case timesWord2# ahi bhi of
                (# h11, l11 #) ->
                  case plusWord2# h00 l01 of
                    (# carryA, sumA #) ->
                      case plusWord2# sumA l10 of
                        (# carryB, r1 #) ->
                          case plusWord# carryA carryB of
                            carry2 ->
                              case plusWord2# h01 h10 of
                                (# carryC, sumC #) ->
                                  case plusWord2# sumC l11 of
                                    (# carryD, sumD #) ->
                                      case plusWord2# sumD carry2 of
                                        (# carryE, r2 #) ->
                                          case plusWord# (plusWord# carryC carryD) carryE of
                                            carry3 ->
                                              case plusWord# h11 carry3 of
                                                r3 -> (# r3, r2, r1, l00 #)

-- Lexicographic comparison of 256-bit magnitudes: 1, 0, or -1.
compareWord4
  :: Word# -> Word# -> Word# -> Word# -> Word# -> Word# -> Word# -> Word# -> Int#
compareWord4 a3 a2 a1 a0 b3 b2 b1 b0 =
  case compareWord2 a3 a2 b3 b2 of
    0# -> compareWord2 a1 a0 b1 b0
    answer -> answer

-- The exponent a mantissa contributes to the alignment floor: a zero
-- mantissa is exact at any floor and votes for the impossibly high sentinel.
floorExp :: Word# -> Int# -> Int#
floorExp mantissa power =
  case mantissa of
    0## -> 2000000#
    _ -> power

minExp :: Int# -> Int# -> Int#
minExp a b = if isTrue# (a <=# b) then a else b

orientSignWorker :: Double# -> Double# -> Double# -> Double# -> Double# -> Double# -> Int#
orientSignWorker ax ay bx by cx cy =
  case decodeExact ax of
    (# negax, max_, eax, okax #) ->
      case decodeExact ay of
        (# negay, may, eay, okay #) ->
          case decodeExact bx of
            (# negbx, mbx, ebx, okbx #) ->
              case decodeExact by of
                (# negby, mby, eby, okby #) ->
                  case decodeExact cx of
                    (# negcx, mcx, ecx, okcx #) ->
                      case decodeExact cy of
                        (# negcy, mcy, ecy, okcy #) ->
                          case okax +# okay +# okbx +# okby +# okcx +# okcy of
                            6# ->
                              -- The alignment floor is the least exponent
                              -- among nonzero mantissas; a zero mantissa is
                              -- exact at any floor and must not drag it down.
                              case floorExp max_ eax `minExp` floorExp may eay `minExp` floorExp mbx ebx `minExp` floorExp mby eby `minExp` floorExp mcx ecx `minExp` floorExp mcy ecy of
                                emin ->
                                  case differenceExact negax max_ eax negcx mcx ecx emin of
                                    (# s1, d1hi, d1lo, f1 #) ->
                                      case differenceExact negby mby eby negcy mcy ecy emin of
                                        (# s2, d2hi, d2lo, f2 #) ->
                                          case differenceExact negay may eay negcy mcy ecy emin of
                                            (# s3, d3hi, d3lo, f3 #) ->
                                              case differenceExact negbx mbx ebx negcx mcx ecx emin of
                                                (# s4, d4hi, d4lo, f4 #) ->
                                                  case f1 +# f2 +# f3 +# f4 of
                                                    0# ->
                                                      combineSigns
                                                        (s1 *# s2) d1hi d1lo d2hi d2lo
                                                        (s3 *# s4) d3hi d3lo d4hi d4lo
                                                    _ -> 2#
                            _ -> 2#
  where
    -- det = leftSign * leftProduct - rightSign * rightProduct
    combineSigns leftSign d1hi d1lo d2hi d2lo rightSign d3hi d3lo d4hi d4lo =
      case leftSign of
        0# ->
          case rightSign of
            0# -> 0#
            _ -> 0# -# rightSign
        _ ->
          case rightSign of
            0# -> leftSign
            _ ->
              case leftSign ==# rightSign of
                1# ->
                  case multiplyWord2 d1hi d1lo d2hi d2lo of
                    (# p3, p2, p1, p0 #) ->
                      case multiplyWord2 d3hi d3lo d4hi d4lo of
                        (# q3, q2, q1, q0 #) ->
                          case compareWord4 p3 p2 p1 p0 q3 q2 q1 q0 of
                            0# -> 0#
                            1# -> leftSign
                            _ -> 0# -# leftSign
                _ -> leftSign

#else

-- The narrow-word answer to the same question. 'exactOrientDet' is the
-- determinant the fixed-width worker exists to avoid allocating, not a
-- different quantity, so the two branches agree by construction.
exactOrientSignDouble
  :: Double -> Double -> Double -> Double -> Double -> Double -> Ordering
exactOrientSignDouble ax ay bx by cx cy =
  compare (exactOrientDet ax ay bx by cx cy) 0
{-# NOINLINE exactOrientSignDouble #-}

#endif