axionbuster · July 30, 2025 03:53 · Jul 30, 2025 · Jul 30, 2025 · Jul 30, 2025
diff --git a/geometric.hs b/geometric.hs
@@ -101,4 +101,4 @@ main = do
 another test. 1 + 285 + ... + 285 ^ 4294967296 = 3364 (mod 7731)
 -}
 
--- very fast!
+-- blazingly fast!
diff --git a/geometric.hs b/geometric.hs
@@ -1,3 +1,5 @@
+-- https://stackoverflow.com/a/9486650 (no relation with author)
+
 import Control.Monad
 import Data.Foldable
 import GHC.Natural
@@ -15,6 +17,7 @@ geo modulus = go
   squared = (`power` 2)
   plus x y = (x + y) `mod` modulus
   times x y = (x * y) `mod` modulus
+  -- you can trust me, it's OK to right-fold.
   go r 0 = 1
   go r n = case n `divMod` 2 of
    (_, 0) -> 1 `plus` (r `times` go r (n - 1))
@@ -38,3 +41,64 @@ main = do
      modulo2 = 7731
  printf "another test. 1 + %v + ... + %v ^ %v = %v (mod %v)\n"
         base2 base2 count2 (geo modulo2 base2 count2) modulo2
+
+{-
+=== Modulo 2 ===
+1 + 2 + ... + 2 ^ 2 = 1 (mod 2)
+1 + 2 + ... + 2 ^ 3 = 1 (mod 2)
+1 + 2 + ... + 2 ^ 4 = 1 (mod 2)
+
+=== Modulo 3 ===
+1 + 2 + ... + 2 ^ 2 = 1 (mod 3)
+1 + 2 + ... + 2 ^ 3 = 0 (mod 3)
+1 + 2 + ... + 2 ^ 4 = 1 (mod 3)
+
+=== Modulo 4 ===
+1 + 2 + ... + 2 ^ 2 = 3 (mod 4)
+1 + 2 + ... + 2 ^ 3 = 3 (mod 4)
+1 + 2 + ... + 2 ^ 4 = 3 (mod 4)
+
+=== Modulo 5 ===
+1 + 2 + ... + 2 ^ 2 = 2 (mod 5)
+1 + 2 + ... + 2 ^ 3 = 0 (mod 5)
+1 + 2 + ... + 2 ^ 4 = 1 (mod 5)
+
+=== Modulo 6 ===
+1 + 2 + ... + 2 ^ 2 = 1 (mod 6)
+1 + 2 + ... + 2 ^ 3 = 3 (mod 6)
+1 + 2 + ... + 2 ^ 4 = 1 (mod 6)
+
+=== Modulo 7 ===
+1 + 2 + ... + 2 ^ 2 = 0 (mod 7)
+1 + 2 + ... + 2 ^ 3 = 1 (mod 7)
+1 + 2 + ... + 2 ^ 4 = 3 (mod 7)
+
+=== Modulo 8 ===
+1 + 2 + ... + 2 ^ 2 = 7 (mod 8)
+1 + 2 + ... + 2 ^ 3 = 7 (mod 8)
+1 + 2 + ... + 2 ^ 4 = 7 (mod 8)
+
+=== Modulo 9 ===
+1 + 2 + ... + 2 ^ 2 = 7 (mod 9)
+1 + 2 + ... + 2 ^ 3 = 6 (mod 9)
+1 + 2 + ... + 2 ^ 4 = 4 (mod 9)
+
+=== Modulo 10 ===
+1 + 2 + ... + 2 ^ 2 = 7 (mod 10)
+1 + 2 + ... + 2 ^ 3 = 5 (mod 10)
+1 + 2 + ... + 2 ^ 4 = 1 (mod 10)
+
+=== Modulo 11 ===
+1 + 2 + ... + 2 ^ 2 = 7 (mod 11)
+1 + 2 + ... + 2 ^ 3 = 4 (mod 11)
+1 + 2 + ... + 2 ^ 4 = 9 (mod 11)
+
+=== Modulo 12 ===
+1 + 2 + ... + 2 ^ 2 = 7 (mod 12)
+1 + 2 + ... + 2 ^ 3 = 3 (mod 12)
+1 + 2 + ... + 2 ^ 4 = 7 (mod 12)
+
+another test. 1 + 285 + ... + 285 ^ 4294967296 = 3364 (mod 7731)
+-}
+
+-- very fast!
diff --git a/geometric.hs b/geometric.hs
@@ -0,0 +1,40 @@
+import Control.Monad
+import Data.Foldable
+import GHC.Natural
+import Prelude hiding (exp)
+import Text.Printf
+
+type N = Natural
+
+-- geometric sum (1 + r + r^2 + ... + r ^ n).
+-- order: geo [modulus] [base] [last exponent, inclusive].
+geo :: N -> N -> N -> N
+geo modulus = go
+ where
+  power x y = powModNatural x y modulus
+  squared = (`power` 2)
+  plus x y = (x + y) `mod` modulus
+  times x y = (x * y) `mod` modulus
+  go r 0 = 1
+  go r n = case n `divMod` 2 of
+   (_, 0) -> 1 `plus` (r `times` go r (n - 1))
+   (q, _) -> (1 `plus` r) `times` go (squared r) q
+
+main :: IO ()
+main = do
+ let base = 2
+ forM_ [2..12] $ \mo -> do
+  printf "=== Modulo %v ===\n" mo
+  forM_ [2..4] $ \po -> do
+   let exp = sum (map (\i -> powModNatural base i mo) [0..po]) `mod` mo
+       get = geo mo base po
+   printf "1 + %v + ... + %v ^ %v = %v (mod %v)" base base po (geo mo 2 po) mo
+   unless (exp == get) $ do
+    printf " (BAD! exp = %v)" exp
+   putStrLn ""
+  putStrLn ""
+ let base2 = 285
+     count2 = 2 ^ (32 :: Int)
+     modulo2 = 7731
+ printf "another test. 1 + %v + ... + %v ^ %v = %v (mod %v)\n"
+        base2 base2 count2 (geo modulo2 base2 count2) modulo2