{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Math.NumberTheory.MoebiusInversion
( generalInversion
, totientSum
) where
import Control.Monad
import Control.Monad.ST
import qualified Data.Vector.Mutable as MV
import Math.NumberTheory.Powers.Squares
totientSum :: Int -> Integer
totientSum n
| n < 1 = 0
| otherwise = generalInversion (triangle . fromIntegral) n
where
triangle k = (k*(k+1)) `quot` 2
generalInversion :: (Int -> Integer) -> Int -> Integer
generalInversion fun n
| n < 1 = error "Möbius inversion only defined on positive domain"
| n == 1 = fun 1
| n == 2 = fun 2 - fun 1
| n == 3 = fun 3 - 2*fun 1
| otherwise = fastInvert fun n
fastInvert :: (Int -> Integer) -> Int -> Integer
fastInvert fun n = runST (fastInvertST fun n)
fastInvertST :: forall s. (Int -> Integer) -> Int -> ST s Integer
fastInvertST fun n = do
let !k0 = integerSquareRoot (n `quot` 2)
!mk0 = n `quot` (2*k0+1)
kmax a m = (a `quot` m - 1) `quot` 2
small <- MV.unsafeNew (mk0 + 1) :: ST s (MV.MVector s Integer)
MV.unsafeWrite small 0 0
MV.unsafeWrite small 1 $! (fun 1)
when (mk0 >= 2) $
MV.unsafeWrite small 2 $! (fun 2 - fun 1)
let calcit :: Int -> Int -> Int -> ST s (Int, Int)
calcit switch change i
| mk0 < i = return (switch,change)
| i == change = calcit (switch+1) (change + 4*switch+6) i
| otherwise = do
let mloop !acc k !m
| k < switch = kloop acc k
| otherwise = do
val <- MV.unsafeRead small m
let nxtk = kmax i (m+1)
mloop (acc - fromIntegral (k-nxtk)*val) nxtk (m+1)
kloop !acc k
| k == 0 = do
MV.unsafeWrite small i $! acc
calcit switch change (i+1)
| otherwise = do
val <- MV.unsafeRead small (i `quot` (2*k+1))
kloop (acc-val) (k-1)
mloop (fun i - fun (i `quot` 2)) ((i-1) `quot` 2) 1
(sw, ch) <- calcit 1 8 3
large <- MV.unsafeNew k0 :: ST s (MV.MVector s Integer)
let calcbig :: Int -> Int -> Int -> ST s (MV.MVector s Integer)
calcbig switch change j
| j == 0 = return large
| (2*j-1)*change <= n = calcbig (switch+1) (change + 4*switch+6) j
| otherwise = do
let i = n `quot` (2*j-1)
mloop !acc k m
| k < switch = kloop acc k
| otherwise = do
val <- MV.unsafeRead small m
let nxtk = kmax i (m+1)
mloop (acc - fromIntegral (k-nxtk)*val) nxtk (m+1)
kloop !acc k
| k == 0 = do
MV.unsafeWrite large (j-1) $! acc
calcbig switch change (j-1)
| otherwise = do
let m = i `quot` (2*k+1)
val <- if m <= mk0
then MV.unsafeRead small m
else MV.unsafeRead large (k*(2*j-1)+j-1)
kloop (acc-val) (k-1)
mloop (fun i - fun (i `quot` 2)) ((i-1) `quot` 2) 1
mvec <- calcbig sw ch k0
MV.unsafeRead mvec 0