axionbuster · December 21, 2024 03:39 · Dec 21, 2024 · Dec 21, 2024 · Dec 21, 2024 · Dec 21, 2024
diff --git a/Lib.hs b/Lib.hs
@@ -41,6 +41,9 @@ march start direction = runST do
         | otherwise = min a b -- if both are NaN, then pick either
       minimum_ = foldr1 minnonan
       sig = floor . signum <$> direction
+      func 1 = floor
+      func (-1) = ceiling
+      func _ = floor -- direction is zero, so it doesn't matter
   cur <- new start
   com <- new $ pure 0 -- Kahan sum compensator
   fix \this -> do
@@ -51,8 +54,9 @@ march start direction = runST do
     let (!) = index
         t cur' i =
           -- solve for time to next intersection in dimension i
-          let s = fi (truncate (cur' ! i) + sig ! i) - cur' ! i
-           in s / direction ! i
+          let s = sig ! i
+              u = fi (func s (cur' ! i) + s) - cur' ! i
+           in u / direction ! i
         add c x y =
           -- Kahan's compensated sum (x += y)
           let y' = y - c

diff --git a/Lib.hs b/Lib.hs
@@ -40,9 +40,9 @@ march start direction = runST do
         | isNaN b = a
         | otherwise = min a b -- if both are NaN, then pick either
       minimum_ = foldr1 minnonan
+      sig = floor . signum <$> direction
   cur <- new start
   com <- new $ pure 0 -- Kahan sum compensator
-  let sig = floor . signum <$> direction
   fix \this -> do
     -- mechanism:
     -- using the parametric equation of the line segment

diff --git a/Lib.hs b/Lib.hs
@@ -42,7 +42,7 @@ march start direction = runST do
       minimum_ = foldr1 minnonan
   cur <- new start
   com <- new $ pure 0 -- Kahan sum compensator
-  let sig = signum <$> direction
+  let sig = floor . signum <$> direction
   fix \this -> do
     -- mechanism:
     -- using the parametric equation of the line segment
@@ -51,8 +51,8 @@ march start direction = runST do
     let (!) = index
         t cur' i =
           -- solve for time to next intersection in dimension i
-          let s = fi (floor $ cur' ! i) + sig ! i
-           in (s - cur' ! i) / direction ! i
+          let s = fi (truncate (cur' ! i) + sig ! i) - cur' ! i
+           in s / direction ! i
         add c x y =
           -- Kahan's compensated sum (x += y)
           let y' = y - c

diff --git a/Lib.hs b/Lib.hs
@@ -71,3 +71,4 @@ march start direction = runST do
     write cur newcur
     write com newcom
     ((tim, newcur) :) <$> this
+{-# INLINEABLE march #-}
diff --git a/Lib.hs b/Lib.hs
@@ -50,7 +50,7 @@ march start direction = runST do
     -- then use the 'time' to get the coordinates of the intersection
     let (!) = index
         t cur' i =
-          -- solve for time to next intersection
+          -- solve for time to next intersection in dimension i
           let s = fi (floor $ cur' ! i) + sig ! i
            in (s - cur' ! i) / direction ! i
         add c x y =

diff --git a/Lib.hs b/Lib.hs
@@ -39,18 +39,19 @@ march start direction = runST do
         | isNaN a = b
         | isNaN b = a
         | otherwise = min a b -- if both are NaN, then pick either
+      minimum_ = foldr1 minnonan
   cur <- new start
   com <- new $ pure 0 -- Kahan sum compensator
-  sig <- new $ signum <$> direction
+  let sig = signum <$> direction
   fix \this -> do
     -- mechanism:
     -- using the parametric equation of the line segment
     -- find the closest intersection with the grid -> get 'time' value
     -- then use the 'time' to get the coordinates of the intersection
     let (!) = index
-        t cur' sig' i =
+        t cur' i =
           -- solve for time to next intersection
-          let s = fi (floor $ cur' ! i) + sig' ! i
+          let s = fi (floor $ cur' ! i) + sig ! i
            in (s - cur' ! i) / direction ! i
         add c x y =
           -- Kahan's compensated sum (x += y)
@@ -60,10 +61,8 @@ march start direction = runST do
            in (s, c')
     com' <- read com
     cur' <- read cur
-    tim <- do
-      sig' <- read sig
-      pure $ foldr1 minnonan $ tabulate @f (t cur' sig')
-    let vadd v w = tabulate @f \i ->
+    let tim = minimum_ $ tabulate @f $ t cur'
+        vadd v w = tabulate @f \i ->
           -- elementwise error-compensated vector addition
           add (com' ! i) (v ! i) (w ! i)
         s = vadd cur' $ direction <&> (* tim)

diff --git a/Lib.hs b/Lib.hs
@@ -43,6 +43,10 @@ march start direction = runST do
   com <- new $ pure 0 -- Kahan sum compensator
   sig <- new $ signum <$> direction
   fix \this -> do
+    -- mechanism:
+    -- using the parametric equation of the line segment
+    -- find the closest intersection with the grid -> get 'time' value
+    -- then use the 'time' to get the coordinates of the intersection
     let (!) = index
         t cur' sig' i =
           -- solve for time to next intersection
@@ -55,15 +59,14 @@ march start direction = runST do
               c' = (s - x) - y'
            in (s, c')
     com' <- read com
-    let vadd v w = tabulate @f \i ->
-          -- elementwise compensated vector addition
-          add (com' ! i) (v ! i) (w ! i)
+    cur' <- read cur
     tim <- do
-      cur' <- read cur
       sig' <- read sig
       pure $ foldr1 minnonan $ tabulate @f (t cur' sig')
-    cur' <- read cur
-    let s = vadd cur' $ direction <&> (* tim)
+    let vadd v w = tabulate @f \i ->
+          -- elementwise error-compensated vector addition
+          add (com' ! i) (v ! i) (w ! i)
+        s = vadd cur' $ direction <&> (* tim)
         newcur = fst <$> s
         newcom = snd <$> s
     write cur newcur

diff --git a/Lib.hs b/Lib.hs
@@ -11,6 +11,8 @@ import Prelude hiding (read)
 -- | march along a line segment, finding all intersections
 -- with grid points
 --
+-- the starting point is NOT included in the output
+--
 -- a compensated sum is used to reduce floating point error
 --
 -- the returned list being infinite, it is recommended to

diff --git a/Lib.hs b/Lib.hs
@@ -0,0 +1,69 @@
+-- | march along a line segment, finding all intersections with grid points
+module Lib (march) where
+
+import Control.Monad.Fix
+import Control.Monad.ST.Lazy
+import Data.Functor
+import Data.Functor.Rep
+import Data.STRef.Lazy
+import Prelude hiding (read)
+
+-- | march along a line segment, finding all intersections
+-- with grid points
+--
+-- a compensated sum is used to reduce floating point error
+--
+-- the returned list being infinite, it is recommended to
+-- use 'take' to limit the number of points to be computed
+march ::
+  forall f a.
+  ( Applicative f,
+    Foldable f,
+    Representable f,
+    RealFloat a
+  ) =>
+  -- | starting point
+  f a ->
+  -- | direction (no need to be normalized)
+  f a ->
+  -- | list of (delta time, point) pairs
+  [(a, f a)]
+march start direction = runST do
+  let fi = fromIntegral :: Int -> a
+      new = newSTRef
+      read = readSTRef
+      write = writeSTRef
+      minnonan a b
+        | isNaN a = b
+        | isNaN b = a
+        | otherwise = min a b -- if both are NaN, then pick either
+  cur <- new start
+  com <- new $ pure 0 -- Kahan sum compensator
+  sig <- new $ signum <$> direction
+  fix \this -> do
+    let (!) = index
+        t cur' sig' i =
+          -- solve for time to next intersection
+          let s = fi (floor $ cur' ! i) + sig' ! i
+           in (s - cur' ! i) / direction ! i
+        add c x y =
+          -- Kahan's compensated sum (x += y)
+          let y' = y - c
+              s = x + y'
+              c' = (s - x) - y'
+           in (s, c')
+    com' <- read com
+    let vadd v w = tabulate @f \i ->
+          -- elementwise compensated vector addition
+          add (com' ! i) (v ! i) (w ! i)
+    tim <- do
+      cur' <- read cur
+      sig' <- read sig
+      pure $ foldr1 minnonan $ tabulate @f (t cur' sig')
+    cur' <- read cur
+    let s = vadd cur' $ direction <&> (* tim)
+        newcur = fst <$> s
+        newcom = snd <$> s
+    write cur newcur
+    write com newcom
+    ((tim, newcur) :) <$> this