axionbuster · March 17, 2025 17:04 · Mar 17, 2025 · Mar 17, 2025 · Mar 17, 2025 · Mar 17, 2025
diff --git a/Hilbert3D.hs b/Hilbert3D.hs
@@ -119,11 +119,10 @@ w .>>. s = w `uncheckedShiftRLWord32#` word2Int# (word8ToWord# s)
 (<) :: Word32# -> Word32# -> Int#
 (<) = ltWord32#
 
-last, db, trp, octal, pred :: Word32# -> Word32#
+last, db, octal, pred :: Word32# -> Word32#
 last = (.&. 1#Word32)
 octal = (.&. 7#Word32)
 db = (.<<. 1#Word8)
-trp w = (w .<<. 1#Word8) + w
 pred w = w - 1#Word32
 
 w32to8 :: Word32# -> Word8#

diff --git a/Hilbert3D.hs b/Hilbert3D.hs
@@ -83,7 +83,7 @@ shil# m s__ = case magic .>>. (trp8 $ octal8 $ w32to8 s__) of
 
 bsiz :: Word32# -> Word32# -> Word32# -> Word8#
 bsiz x y z =
-  let max = case x < y of
+  let !max = case x < y of
         1# -> case y < z of
           1# -> z
           _ -> y

diff --git a/Hilbert3D.hs b/Hilbert3D.hs
@@ -26,7 +26,7 @@ shil (W8# m) (W32# s) = case shil# m s of
 -- by David Walker, U. Tennessee at Chattanooga, Aug 2023.
 
 hils# :: Word8# -> Word32# -> Word32# -> Word32# -> Word32#
-hils# m _ _ _ | isTrue# (m <=! 0#Word8) = 0#Word32
+hils# m _ _ _ | isTrue# (m ==! 0#Word8) = 0#Word32
 hils# m _ _ _ | isTrue# (m >=! 32#Word8) = 0#Word32
 hils# m x__ y__ z__ = case pred (1#Word32 .<<. m) of
   w -> case (# x__ .&. w, y__ .&. w, z__ .&. w #) of
@@ -97,8 +97,8 @@ bsiz x y z =
 (-!) :: Word8# -> Word8# -> Word8#
 (-!) = subWord8#
 
-(<=!), (>=!) :: Word8# -> Word8# -> Int#
-(<=!) = leWord8#
+(==!), (>=!) :: Word8# -> Word8# -> Int#
+(==!) = eqWord8#
 (>=!) = geWord8#
 
 trp8, octal8, rem3 :: Word8# -> Word8#

diff --git a/Hilbert3D.hs b/Hilbert3D.hs
@@ -26,7 +26,8 @@ shil (W8# m) (W32# s) = case shil# m s of
 -- by David Walker, U. Tennessee at Chattanooga, Aug 2023.
 
 hils# :: Word8# -> Word32# -> Word32# -> Word32# -> Word32#
-hils# m _ _ _ | isTrue# ((m <=! 0#Word8) || (m >=! 32#Word8)) = 0#Word32
+hils# m _ _ _ | isTrue# (m <=! 0#Word8) = 0#Word32
+hils# m _ _ _ | isTrue# (m >=! 32#Word8) = 0#Word32
 hils# m x__ y__ z__ = case pred (1#Word32 .<<. m) of
   w -> case (# x__ .&. w, y__ .&. w, z__ .&. w #) of
     (# x, y, z #) -> case bsiz x y z of
@@ -40,7 +41,7 @@ hils# m x__ y__ z__ = case pred (1#Word32 .<<. m) of
     go s i x y z = case last (x .>>. i)
       .|. (last (y .>>. i) .<<. 1#Word8)
       .|. (last (z .>>. i) .<<. 2#Word8) of
-      p -> case octal (magic .>>. w32to8 (trp p)) of
+      p -> case octal (magic .>>. trp8 (w32to8 p)) of
         o -> case (s .<<. 3#Word8) + o of
           s' -> case 1#Word32 .<<. i of
             w -> do
@@ -81,7 +82,6 @@ shil# m s__ = case magic .>>. (trp8 $ octal8 $ w32to8 s__) of
         _ -> k (pred $ w - y) (pred $ db w - z) x
 
 bsiz :: Word32# -> Word32# -> Word32# -> Word8#
-bsiz 0#Word32 0#Word32 0#Word32 = 1#Word8
 bsiz x y z =
   let max = case x < y of
         1# -> case y < z of
@@ -90,7 +90,9 @@ bsiz x y z =
         _ -> case x < z of
           1# -> z
           _ -> x
-   in 32#Word8 -! (wordToWord8# $ clz32# $ word32ToWord# max)
+   in case max of
+        0#Word32 -> 1#Word8
+        _ -> 32#Word8 -! (wordToWord8# $ clz32# $ word32ToWord# max)
 
 (-!) :: Word8# -> Word8# -> Word8#
 (-!) = subWord8#
@@ -127,9 +129,6 @@ pred w = w - 1#Word32
 w32to8 :: Word32# -> Word8#
 w32to8 w = wordToWord8# (word32ToWord# w)
 
-(||) :: Int# -> Int# -> Int#
-(||) = orI#
-
 infixl 9 +
 
 infixl 9 -
diff --git a/Hilbert3D.hs b/Hilbert3D.hs
@@ -0,0 +1,135 @@
+{-# LANGUAGE BlockArguments, MagicHash, UnboxedTuples #-}
+{-# LANGUAGE ExtendedLiterals, NoImplicitPrelude, ViewPatterns #-}
+
+-- | Convert between 3D Hilbert curve coordinates and distance.
+--
+-- The chosen curve goes from (0,0,0) to (1,0,0) to (1,0,1)
+-- for the first three points (listed in x,y,z order).
+module Data.Hilbert3D (hils, hils#, shil, shil#) where
+
+import GHC.Base (($))
+import GHC.Exts
+import GHC.Word
+
+-- | Given order (in range [1, 31]), x, y, and z, respectively,
+-- compute the distance of the point along the Hilbert curve.
+hils :: Word8 -> Word32 -> Word32 -> Word32 -> Word32
+hils (W8# m) (W32# x) (W32# y) (W32# z) = W32# (hils# m x y z)
+
+-- | Given order (in range [1, 31]), and distance along the
+-- Hilbert curve, return its coordinates in (x, y, z) order.
+shil :: Word8 -> Word32 -> (Word32, Word32, Word32)
+shil (W8# m) (W32# s) = case shil# m s of
+  (# x, y, z #) -> (W32# x, W32# y, W32# z)
+
+-- from Algorithms for Encoding and Decoding 3D Hilbert Orderings
+-- by David Walker, U. Tennessee at Chattanooga, Aug 2023.
+
+hils# :: Word8# -> Word32# -> Word32# -> Word32# -> Word32#
+hils# m _ _ _ | isTrue# ((m <=! 0#Word8) || (m >=! 32#Word8)) = 0#Word32
+hils# m x__ y__ z__ = case pred (1#Word32 .<<. m) of
+  w -> case (# x__ .&. w, y__ .&. w, z__ .&. w #) of
+    (# x, y, z #) -> case bsiz x y z of
+      i -> case i -! 1#Word8 of
+        j -> case rem3 (m -! i) of
+          1#Word8 -> go 0#Word32 j z x y
+          2#Word8 -> go 0#Word32 j y z x
+          _ -> go 0#Word32 j x y z
+  where
+    magic = 0o54236710#Word32 -- octal lookup table. zyx to octant.
+    go s i x y z = case last (x .>>. i)
+      .|. (last (y .>>. i) .<<. 1#Word8)
+      .|. (last (z .>>. i) .<<. 2#Word8) of
+      p -> case octal (magic .>>. w32to8 (trp p)) of
+        o -> case (s .<<. 3#Word8) + o of
+          s' -> case 1#Word32 .<<. i of
+            w -> do
+              let k = case i of
+                    0#Word8 -> \_ _ _ -> s'
+                    _ -> go s' (i -! 1#Word8)
+              case o of
+                0#Word32 -> k z x y
+                1#Word32 -> k y z (x - w)
+                2#Word32 -> k y (z - w) (x - w)
+                3#Word32 -> k (pred $ w - x) y (pred $ db w - z)
+                4#Word32 -> k (pred $ w - x) (y - w) (pred $ db w - z)
+                5#Word32 -> k (pred $ db w - y) (pred $ db w - z) (x - w)
+                6#Word32 -> k (pred $ db w - y) (pred $ w - z) (x - w)
+                _ -> k z (pred $ w - x) (pred $ db w - y)
+
+shil# :: Word8# -> Word32# -> (# Word32#, Word32#, Word32# #)
+shil# m s__ = case magic .>>. (trp8 $ octal8 $ w32to8 s__) of
+  o -> case s__ .>>. 3#Word8 of
+    s ->
+      go 2#Word32 s (last o) (last (o .>>. 1#Word8)) (last (o .>>. 2#Word8))
+  where
+    magic = 0o23764510#Word32 -- inverse table. octant to zyx.
+    go _ 0#Word32 x y z = case rem3 (m -! bsiz x y z) of
+      1#Word8 -> (# y, z, x #)
+      2#Word8 -> (# z, x, y #)
+      _ -> (# x, y, z #)
+    go w s x y z = do
+      let k = go (w .<<. 1#Word8) (s .>>. 3#Word8)
+      case octal s of
+        0#Word32 -> k y z x
+        1#Word32 -> k (z + w) x y
+        2#Word32 -> k (z + w) x (y + w)
+        3#Word32 -> k (pred $ w - x) y (pred $ db w - z)
+        4#Word32 -> k (pred $ w - x) (y + w) (pred $ db w - z)
+        5#Word32 -> k (z + w) (pred $ db w - x) (pred $ db w - y)
+        6#Word32 -> k (z + w) (pred $ db w - x) (pred $ w - y)
+        _ -> k (pred $ w - y) (pred $ db w - z) x
+
+bsiz :: Word32# -> Word32# -> Word32# -> Word8#
+bsiz 0#Word32 0#Word32 0#Word32 = 1#Word8
+bsiz x y z =
+  let max = case x < y of
+        1# -> case y < z of
+          1# -> z
+          _ -> y
+        _ -> case x < z of
+          1# -> z
+          _ -> x
+   in 32#Word8 -! (wordToWord8# $ clz32# $ word32ToWord# max)
+
+(-!) :: Word8# -> Word8# -> Word8#
+(-!) = subWord8#
+
+(<=!), (>=!) :: Word8# -> Word8# -> Int#
+(<=!) = leWord8#
+(>=!) = geWord8#
+
+trp8, octal8, rem3 :: Word8# -> Word8#
+trp8 w = (w `uncheckedShiftLWord8#` 1#) `plusWord8#` w
+octal8 = (`andWord8#` 7#Word8)
+rem3 = (`remWord8#` 3#Word8)
+
+(-), (+), (.&.), (.|.) :: Word32# -> Word32# -> Word32#
+(-) = subWord32#
+(+) = plusWord32#
+(.&.) = andWord32#
+(.|.) = orWord32#
+
+(.<<.), (.>>.) :: Word32# -> Word8# -> Word32#
+w .<<. s = w `uncheckedShiftLWord32#` word2Int# (word8ToWord# s)
+w .>>. s = w `uncheckedShiftRLWord32#` word2Int# (word8ToWord# s)
+
+(<) :: Word32# -> Word32# -> Int#
+(<) = ltWord32#
+
+last, db, trp, octal, pred :: Word32# -> Word32#
+last = (.&. 1#Word32)
+octal = (.&. 7#Word32)
+db = (.<<. 1#Word8)
+trp w = (w .<<. 1#Word8) + w
+pred w = w - 1#Word32
+
+w32to8 :: Word32# -> Word8#
+w32to8 w = wordToWord8# (word32ToWord# w)
+
+(||) :: Int# -> Int# -> Int#
+(||) = orI#
+
+infixl 9 +
+
+infixl 9 -