# arithmoi / Math / NumberTheory / MoebiusInversion / Int.hs

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137``` ```-- | -- Module: Math.NumberTheory.MoebiusInversion -- Copyright: (c) 2012 Daniel Fischer -- Licence: MIT -- Maintainer: Daniel Fischer -- Stability: Provisional -- Portability: Non-portable (GHC extensions) -- -- Generalised Moebius inversion for 'Int' valued functions. -- {-# LANGUAGE BangPatterns #-} {-# OPTIONS_GHC -fspec-constr-count=8 #-} module Math.NumberTheory.MoebiusInversion.Int ( generalInversion , totientSum ) where import Data.Array.ST import Data.Array.Base import Control.Monad.ST import Math.NumberTheory.Powers.Squares -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@, -- computed via generalised Moebius inversion. -- Arguments less than 1 cause an error to be raised. totientSum :: Int -> Int totientSum = (+1) . generalInversion triangle where triangle n = (n*(n-1)) `quot` 2 -- | @generalInversion g n@ evaluates the generalised Moebius inversion of @g@ -- at the argument @n@. -- -- The generalised Moebius inversion implemented here allows an efficient -- calculation of isolated values of the function @f : N -> Z@ if the function -- @g@ defined by -- -- > -- > g n = sum [f (n `quot` k) | k <- [1 .. n]] -- > -- -- can be cheaply computed. By the generalised Moebius inversion formula, then -- -- > -- > f n = sum [moebius k * g (n `quot` k) | k <- [1 .. n]] -- > -- -- which allows the computation in /O/(n) steps, if the values of the -- Moebius function are known. A slightly different formula, used here, -- does not need the values of the Moebius function and allows the -- computation in /O/(n^0.75) steps, using /O/(n^0.5) memory. -- -- An example of a pair of such functions where the inversion allows a -- more efficient computation than the direct approach is -- -- > -- > f n = number of reduced proper fractions with denominator <= n -- > g n = number of proper fractions with denominator <= n -- > -- -- (a /proper fraction/ is a fraction @0 < n/d < 1@). Then @f n@ is the -- cardinality of the Farey sequence of order @n@ (minus 1 or 2 if 0 and/or -- 1 are included in the Farey sequence), or the sum of the totients of -- the numbers @2 <= k <= n@. @f n@ is not easily directly computable, -- but then @g n = n*(n-1)/2@ is very easy to compute, and hence the inversion -- gives an efficient method of computing @f n@. -- -- For 'Int' valued functions, unboxed arrays can be used for greater efficiency. -- That bears the risk of overflow, however, so be sure to use it only when it's -- safe. -- -- The value @f n@ is then computed by @generalInversion g n@). Note that when -- many values of @f@ are needed, there are far more efficient methods, this -- method is only appropriate to compute isolated values of @f@. generalInversion :: (Int -> Int) -> Int -> Int generalInversion fun n | n == 1 = fun 1 | n == 2 = fun 2 - fun 1 | n == 3 = fun 3 - 2*fun 1 | otherwise = fastInvert fun n fastInvert :: (Int -> Int) -> Int -> Int fastInvert fun n = big `unsafeAt` 0 where !k0 = integerSquareRoot (n `quot` 2) !mk0 = n `quot` (2*k0+1) kmax a m = (a `quot` m - 1) `quot` 2 big = runSTUArray \$ do small <- newArray_ (0,mk0) :: ST s (STUArray s Int Int) unsafeWrite small 0 0 unsafeWrite small 1 (fun 1) unsafeWrite small 2 (fun 2 - fun 1) let 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 <- unsafeRead small m let nxtk = kmax i (m+1) mloop (acc - (k-nxtk)*val) nxtk (m+1) kloop !acc k | k == 0 = do unsafeWrite small i acc calcit switch change (i+1) | otherwise = do val <- 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 <- newArray_ (0,k0-1) let 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 <- unsafeRead small m let nxtk = kmax i (m+1) mloop (acc - (k-nxtk)*val) nxtk (m+1) kloop !acc k | k == 0 = do unsafeWrite large (j-1) acc calcbig switch change (j-1) | otherwise = do let m = i `quot` (2*k+1) val <- if m <= mk0 then unsafeRead small m else 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 calcbig sw ch k0 ```