Added faster exponentiation modulo logic

README.md

@@ -20,5 +20,5 @@ Run `stack exec pi-digits-exe`, or execute the binary found in `./.stack-work/in
 - [X] ~~Replace the mess in `prompt` with a monad.~~ Issue corrected by using `optparse-applicative`.
 - [X] Separate into different files. Possibly parsing and logic files.
 - [X] ~~Find a better way to transfer "globals", like `delim` and `printFun`.~~ `parseIndex` handles all parsing and calling logic now.
-- [ ] Implement a faster mod operation, to allow for larger numbers (like 12345678901234567890). It will likely be implemented with the algorithm explained in the paper.
+- [X] Implement a faster mod operation, to allow for larger numbers (like 12345678901234567890). It will likely be implemented with the algorithm explained in the paper.
 - [ ] Customize print behavior and frequency. Flush output every N digits, or something similar.

src/Logic.hs

@@ -26,7 +26,7 @@ hexPi n =
 sumPi :: Integer -> Integer -> Double
 sumPi n x =
   let
-    summation1 = sum [(fromIntegral ((16^(n-k) `mod` ((8*k)+x)))) / (fromIntegral ((8*k)+x)) | k <- [0..n]]
+    summation1 = sum [(fromIntegral (fastMod 16 (n-k) ((8*k)+x))) / (fromIntegral ((8*k)+x)) | k <- [0..n]]
     summation2 = sum [16^^(n-k) / (fromIntegral ((8*k)+x)) | k <- [(n+1)..5000]]
   in
     summation1 + summation2
@@ -35,3 +35,36 @@ sumPi n x =
 -- The list of answers.
 hexDigits :: [Integer]
 hexDigits = [hexPi x | x <- [0..]]
+
+
+fastMod :: Integer -> Integer -> Integer -> Integer
+fastMod b n k =
+  let
+    t = largestT 0 n
+  in
+    a b (fromIntegral n) k 1 t
+
+-- Get the largest t such that t^2 <= n
+largestT :: Integer -> Integer -> Double
+largestT t n
+  | (2 ^ (t + 1)) <= n = largestT (t + 1) n
+  | otherwise = 2^^t
+
+a :: Integer -> Double -> Integer -> Integer -> Double -> Integer
+a b n k r t =
+  if n >= t then
+    let
+      r' = (b * r) `mod` k
+      n' = n - t
+      t' = t / 2
+    in
+      if t' >= 1 then
+          a b n' k ((r' ^ 2) `mod` k) t'
+      else
+        r'
+  else
+    let t' = t / 2 in
+      if t' >= 1 then
+          a b n k ((r ^ 2) `mod` k) t'
+      else
+        r