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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
-- |
-- Module: Math.NumberTheory.MoebiusInversion
-- Copyright: (c) 2012 Daniel Fischer
-- Licence: MIT
-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>
--
-- Generalised Möbius inversion
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Math.NumberTheory.MoebiusInversion
( generalInversion
, totientSum
) where
import Control.Monad
import Control.Monad.ST
import Data.Proxy
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as MG
import Math.NumberTheory.Roots
import Math.NumberTheory.Utils.FromIntegral
-- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@,
-- computed via generalised Möbius inversion.
-- See <http://mathworld.wolfram.com/TotientSummatoryFunction.html> for the
-- formula used for @totientSum@.
--
-- >>> import Data.Proxy
-- >>> totientSum (Proxy :: Proxy Data.Vector.Unboxed.Vector) 100 :: Int
-- 3044
-- >>> totientSum (Proxy :: Proxy Data.Vector.Vector) 100 :: Integer
-- 3044
totientSum
:: (Integral t, G.Vector v t)
=> Proxy v
-> Word
-> t
totientSum _ 0 = 0
totientSum proxy n = generalInversion proxy (triangle . fromIntegral) n
where
triangle k = (k * (k + 1)) `quot` 2
-- | @generalInversion g n@ evaluates the generalised Möbius inversion of @g@
-- at the argument @n@.
--
-- The generalised Möbius 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 Möbius 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
-- Möbius function are known. A slightly different formula, used here,
-- does not need the values of the Möbius 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@.
--
-- Since the function arguments are used as array indices, the domain of
-- @f@ is restricted to 'Int'.
--
-- 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
:: (Num t, G.Vector v t)
=> Proxy v
-> (Word -> t)
-> Word
-> t
generalInversion proxy fun n = case n of
0 ->error "Möbius inversion only defined on positive domain"
1 -> fun 1
2 -> fun 2 - fun 1
3 -> fun 3 - 2*fun 1
_ -> runST (fastInvertST proxy (fun . intToWord) (wordToInt n))
fastInvertST
:: forall s t v.
(Num t, G.Vector v t)
=> Proxy v
-> (Int -> t)
-> Int
-> ST s t
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 <- MG.unsafeNew (mk0 + 1) :: ST s (G.Mutable v s t)
MG.unsafeWrite small 0 0
MG.unsafeWrite small 1 $! fun 1
when (mk0 >= 2) $
MG.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 <- MG.unsafeRead small m
let nxtk = kmax i (m+1)
mloop (acc - fromIntegral (k-nxtk)*val) nxtk (m+1)
kloop !acc k
| k == 0 = do
MG.unsafeWrite small i $! acc
calcit switch change (i+1)
| otherwise = do
val <- MG.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 <- MG.unsafeNew k0 :: ST s (G.Mutable v s t)
let calcbig :: Int -> Int -> Int -> ST s (G.Mutable v s t)
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 <- MG.unsafeRead small m
let nxtk = kmax i (m+1)
mloop (acc - fromIntegral (k-nxtk)*val) nxtk (m+1)
kloop !acc k
| k == 0 = do
MG.unsafeWrite large (j-1) $! acc
calcbig switch change (j-1)
| otherwise = do
let m = i `quot` (2*k+1)
val <- if m <= mk0
then MG.unsafeRead small m
else MG.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
MG.unsafeRead mvec 0
|
