moonlight-triangulation-0.1.0.0: src-core/Moonlight/Triangulation/Internal/Dyadic.hs
{-# 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