Last active
July 30, 2025 03:53
-
-
Save axionbuster/66f4c7e4d0ce1d56042cb3c596385076 to your computer and use it in GitHub Desktop.
Fast modular geometric sum
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| -- https://stackoverflow.com/a/9486650 (no relation with author) | |
| 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 | |
| -- 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)) | |
| (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 | |
| {- | |
| === 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) | |
| -} | |
| -- blazingly fast! |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment