commit 680a0698c60fe87b21c0d65c20751c2f4b044aed
parent 8e0da908a4b40fc350164b26de2f7ce3c8707868
Author: Jared Tobin <jared@jtobin.io>
Date: Fri, 24 Jan 2025 11:21:07 +0400
lib: pure quotrem_knuth
Diffstat:
2 files changed, 72 insertions(+), 1 deletion(-)
diff --git a/lib/Data/Word/Extended.hs b/lib/Data/Word/Extended.hs
@@ -97,6 +97,11 @@ data Word640 = Word640
{-# UNPACK #-} !Word64
deriving (Eq, Show, Generic)
+data Word1152 = Word1152 -- yikes
+ {-# UNPACK #-} !Word576
+ {-# UNPACK #-} !Word576
+ deriving (Eq, Show, Generic)
+
-- conversion -----------------------------------------------------------------
to_integer :: Word256 -> Integer
@@ -503,6 +508,44 @@ quotrem_by1_gen u ulen d =
!nacc = set576 acc j q_j
in loop (pred j) nacc r
+quotrem_knuth_gen
+ :: Word576
+ -> Int
+ -> Word256
+ -> Int
+ -> Word1152
+quotrem_knuth_gen u ulen d dlen = loop (ulen - dlen - 1) zero576 u where
+ !d_hi = sel256 d (dlen - 1)
+ !d_lo = sel256 d (dlen - 2)
+ !rec = recip_2by1 d_hi
+ loop j !qacc !uacc
+ | j < 0 = Word1152 qacc uacc
+ | otherwise =
+ let !u_2 = sel576 uacc (j + dlen)
+ !u_1 = sel576 uacc (j + dlen - 1)
+ !u_0 = sel576 uacc (j + dlen - 2)
+ !qhat
+ | u_2 >= d_hi = 0xffff_ffff_ffff_ffff
+ | otherwise =
+ let !(P qh rh) = quotrem_2by1 u_2 u_1 d_hi rec
+ !(P ph pl) = mul_c qh d_lo
+ in if ph > rh || (ph == rh && pl > u_0)
+ then qh - 1
+ else qh
+
+ !(Word640 u0 borrow) = sub_mul u j d dlen qhat
+ !u1 = set576 u0 (j + dlen) (u_2 - borrow)
+ in if u_2 < borrow
+ then
+ let !qh = qhat - 1
+ !(Word640 u2 r) = add_big u1 j d dlen
+ !u3 = set576 u2 (j + dlen) r
+ !q = set576 qacc j qh
+ in loop (pred j) q u3
+ else
+ let !q = set576 qacc j qhat
+ in loop (pred j) q u1
+
-- XX primarray
quotrem_knuth
:: PrimMonad m
@@ -543,6 +586,7 @@ quotrem_knuth quo u d = do
loop (pred j)
loop (lu - ld - 1)
+
-- XX needs work
-- quotrem_256
-- :: Word256
diff --git a/test/Main.hs b/test/Main.hs
@@ -232,6 +232,32 @@ quotrem_by1_gen_case1 = do
H.assertEqual "remainder matches" pec_rem r
H.assertEqual "quotient matches" pec_quo q
+quotrem_knuth_gen_case0 :: H.Assertion
+quotrem_knuth_gen_case0 = do
+ let !u = Word576
+ 2162362899639802732
+ 8848548347662387477
+ 13702897166684377657
+ 16799544643779908154
+ 1
+ 0 0 0 0
+ !d = Word256
+ 16950798510782491100
+ 2612788699139816405
+ 5146719872810836952
+ 14966148379609982000
+ !(Word1152 q nu) = quotrem_knuth_gen u 5 d 4
+ !pec_q = Word576 2 0 0 0 0 0 0 0 0
+ !pec_u = Word576
+ 5154254025493923764
+ 3622970949382754665
+ 3409457421062703753
+ 5313991958269495770
+ 0
+ 0 0 0 0
+ H.assertEqual "divisor matches" pec_u nu
+ H.assertEqual "quotient matches" pec_q q
+
quotrem_knuth_case0 :: H.Assertion
quotrem_knuth_case0 = do
let (q, u) = runST $ do
@@ -258,7 +284,7 @@ quotrem_knuth_case0 = do
5154254025493923764
, 3622970949382754665
, 3409457421062703753
- , 5313991958269495770 -- was off by 1 :|
+ , 5313991958269495770
, 0
]
H.assertEqual "divisor matches" pec_u u
@@ -388,6 +414,7 @@ main = defaultMain $
, H.testCase "quotrem_by1 matches case1" quotrem_by1_case1
, H.testCase "quotrem_by1_gen matches case0" quotrem_by1_gen_case0
, H.testCase "quotrem_by1_gen matches case1" quotrem_by1_gen_case1
+ , H.testCase "quotrem_knuth_gen matches case0" quotrem_knuth_gen_case0
, H.testCase "quotrem_knuth matches case0" quotrem_knuth_case0
, H.testCase "quotrem matches case0" quotrem_case0
, H.testCase "quotrem matches case1" quotrem_case1
