module MCPowser where

import Data.Bits
import Data.List
import Data.Ratio
import Data.Tuple

{-
This is a draft of an implementation of the algorithm described in Section 5 of
https://maths-people.anu.edu.au/~brent/pd/rpb242.pdf .

The MC is not 100% sure it really works in all cases, but it delivered correct
results for all of his test cases.
-}

ld = go 0
  where go acc n = if n <= 1 then acc else go (acc + 1) (n `div` 2)

negateIfOdd n x = if even n then x else -x

fact n = product [2..n]

{-# INLINE timesTwoPow #-}
timesTwoPow :: Integer -> Integer -> Integer
timesTwoPow x n = x `shiftL` fromIntegral n

{-# INLINE divTwoPow #-}
divTwoPow :: Integer -> Integer -> Integer
divTwoPow x n = x `shiftR` fromIntegral n

{-# INLINE extractBits #-}
extractBits :: Integer -> Integer -> (Integer, Integer)
extractBits n x = (x .&. (shiftL 1 (fromIntegral n) - 1), shiftR x (fromIntegral n))

everyOther :: [a] -> [a]
everyOther (x : _ : xs) = x : everyOther xs
everyOther _            = []

tangents :: Integer -> [Integer]
tangents n = ts
  where -- p is the precision (i.e. number of bits per "block") to use
        p = n * (ld n + 1)
        zs = reverse $ scanl' (*) 1 [2*n,2*n-1..1]
        -- s and c approximate sin(2^(-p)) and cos(2^(-p)), shifted with a power of two to get an integer
        !s = sum $ zipWith (\k x -> negateIfOdd k $ x `timesTwoPow` (2*(n-k-1)*p))   [0..n-1] (everyOther (tail zs))
        !c = sum $ zipWith (\k x -> negateIfOdd k $ x `timesTwoPow` ((2*(n-k)-1)*p)) [0..n-1] (everyOther zs)
        -- v is then a (shifted) approximation for tan(2^(-p))
        !v = ((fact (2*n-1) `timesTwoPow` ((2*n-1)*p)) * s) `div` c
        -- we carve out the modified Tangent numbers from the bit blocks of v
        ts' = genericTake n $ unfoldr (Just . extractBits (2*p)) v
        -- and compute the actual tangent numbers
        ts = zipWith (\t k -> t `div` product [2*k..2*n-1]) (reverse ts') [1..]

tangent :: Integer -> Integer
tangent = last . tangents . succ

bernoulli :: Integer -> Rational
bernoulli 0 = 1
bernoulli 1 = -0.5
bernoulli n | odd n = 0
bernoulli n = negateIfOdd (n `div` 2) $ (n * tangent (n `div` 2 - 1)) % (2 ^ n - 4 ^ n)
  where negateIfOdd n = if even n then id else negate