Essays/Mertens Function
Basics
The Mertens function on positive integer n is the sum of the Möbius function on the integers from 1 to n .
mobius =: */ @: - @: ~: @: q: mertens0=: +/ @: (mobius"0) @: >: @: i. mertens0"0 >:i.5 10 1 0 _1 _1 _2 _1 _2 _2 _2 _1 _2 _2 _3 _2 _1 _1 _2 _2 _3 _3 _2 _1 _2 _2 _2 _1 _1 _1 _2 _3 _4 _4 _3 _2 _1 _1 _2 _1 0 0 _1 _2 _3 _3 _3 _2 _3 _3 _3 _3
Faster and leaner computations of the Mertens function are possible.
Avoiding 0s
mertens1=: 3 : 0
p=. i.@(1&>.)&.(p:^:_1) 1+<.%:y
b=. -. > ([ +. #@[ $ {.&1@])~&.>/(<"0 *:p),<(1+y)$0
+/ _1 ^ #@q: I. b
)
(mertens1"0 -: mertens1"0) >: i.5 10
1
mertens1"0 ]10^i.8
1 _1 1 2 _23 _48 212 1037
mertens1 avoids applying the Möbius function to exactly those numbers having a 0 value; that is, all numbers containing a square; that is, all multiples of the squares of 2 3 5 7 11 13 ... up to the largest prime less than or equal to the square root of n . The Möbius function on the remaining numbers (I. b in mertens1) is just the parity of the length of the prime factorization, 1 if even and _1 if odd.
The second line in mertens1 is equivalent to b=. -. +./ (1+y) ($ {.&1)"0 *:p but is much leaner by not creating a large intermediate matrix.
Avoiding Factoring
mq=: 0&$: : (4 : 0)
p=. i.&.(p:^:_1) 10>.1+<.y%2
q=. ,. p {.~ p (<:i.0:) <.%:y
z=. (2<.y) $ (1,p:^:_1)`((i.1 0) ; [: ,. i.&.(p:^:_1))@.x >:y
for. i. 0 >. <: (*/\p:i.20) I. >:y do.
p=. p {.~ p (<:i.0:) <.y%*/{.q
q=. q #~ ({:"1 q)<y%*/"1 q
j=. 1 + p i. {:"1 q
k=. j -~ p I. <.1+y%*/"1 q
q=. (k#q) ,. (j,:"0 k) ;@:(<;.0) p
z=. z, #`<@.x q
end.
)
mertens2=: -/ @ mq
mertens2 avoids factoring altogether, but instead enumerates lists of distinct primes of each length. The sub-function mq provides counts if the left argument is 0 and the actual tables of primes if it is 1 ; the monad mq is equivalent to 0&mq .
For example:
mq 1e6 1 78498 209867 206964 92966 18387 1235 8 -/ mq 1e6 212 mertens2 1e6 212 1 mq 45 ┌┬──┬────┬─────┐ ││ 2│2 3│2 3 5│ ││ 3│2 5│2 3 7│ ││ 5│2 7│ │ ││ 7│2 11│ │ ││11│2 13│ │ ││13│2 17│ │ ││17│2 19│ │ ││19│3 5│ │ ││23│3 7│ │ ││29│3 11│ │ ││31│3 13│ │ ││37│5 7│ │ ││41│ │ │ ││43│ │ │ └┴──┴────┴─────┘ #&> 1 mq 45 1 14 12 2 mq 45 1 14 12 2
Given integer n and q , the c-column table of prime factorizations of square-free numbers less than or equal to n , the dyad n nextq q generates the (1+c)-column table of such prime factorizations.
nextq=: 4 : 0
p=. i.&.(p:^:_1) 1+<.x%*/{.y
y=. y #~ ({:"1 y)<x%*/"1 y
j=. (**/$y) * 1 + p i. {:"1 y
k=. j -~ p I. <.1+x%*/"1 y
(k#y) ,. (j,:"0 k) ;@:(<;.0) p
)
45 nextq i.1 0
2
3
5
7
11
13
17
19
23
29
31
37
41
43
45 nextq 45 nextq i.1 0
2 3
2 5
2 7
2 11
2 13
2 17
2 19
3 5
3 7
3 11
3 13
5 7
45 nextq 45 nextq 45 nextq i.1 0
2 3 5
2 3 7
45 nextq&.>^:(<4) <i.1 0
┌┬──┬────┬─────┐
││ 2│2 3│2 3 5│
││ 3│2 5│2 3 7│
││ 5│2 7│ │
││ 7│2 11│ │
││11│2 13│ │
││13│2 17│ │
││17│2 19│ │
││19│3 5│ │
││23│3 7│ │
││29│3 11│ │
││31│3 13│ │
││37│5 7│ │
││41│ │ │
││43│ │ │
└┴──┴────┴─────┘
(1 mq 45) -: 45 nextq&.>^:(<4) <i.1 0
1
y=: ?1e6
p: i.20
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
*/\ p: i.20
2 6 30 210 2310 30030 510510 9.69969e6 2.23093e8 6.46969e9 2.0056e11 7.42074e12 3.0425e14 ...
1 + (*/\p:i.20) I. >:y
8
(1 mq y) -: y nextq&.>^:(<8) <i.1 0
1
The last equivalence suggests an alternative computation of 1&mq :
qn =: 1 + (*/\p:i.20) I. >: mq1=: 3 : 'y nextq&.>^:(<qn y) <i.1 0' (1&mq -: mq1) ?1e6 1
mqcheck is like mq but incorporates checks.
mqcheck=: 0&$: : (4 : 0)
z=. x mq y
assert. ($z) = ,1 + (*/\p:i.20) I. >:y
if. 0=x do.
assert. (0 < z) *. z = <.z
assert. (mertens0 y) = -/ z
assert. 1={.z
assert. (((1<y)$1){z) = p:^:_1 >:y
else.
assert. (32=3!:0 z) *. 2=#@$&>z
assert. (i.#z) -: {:@$&> z
assert. * #&> z
assert. 1 p: ; ,&.> z
assert. z -: ~. &.> z
assert. z -: ~."1&.> z
assert. z -: /:~ &.> z
assert. z -: /:~"1&.> z
assert. y >: ; */"1&.> z
assert. (0 mq y) -: #&> z
assert. (mertens0 y) = -/ #&> z
assert. ({.z) = <i.1 0
assert. (((1<y)$1){z) = <,.i.&.(p:^:_1) >:y
assert. ({.&.>}.z) -: <\p:i.<:#z
end.
)
Avoiding Factoring (2)
A leaner variation of mertens2 derives by observing that for present purposes, the prime factorization can be replaced by the product and the largest prime in a factorization. Thus:
mr=: 0&$: : (4 : 0)
p=. i.&.(p:^:_1) 1+<.y%2
r=. ,. i.&.(p:^:_1) 1+<.%:y
z=. (2<.y) $ (1,p:^:_1)`((,1) ; i.&.(p:^:_1))@.x >:y
for. i. 0 >. <: (*/\p:i.20) I. >:y do.
p=. p {.~ p (<:i.0:) <.y%{.{.r
r=. r #~ ({:"1 r)<y%{."1 r
j=. 1 + p i. {:"1 r
k=. j -~ p I. <.1+y%{."1 r
r=. (k#{."1 r) (* ,. ]) (j,:"0 k) ;@:(<;.0) p
z=. z, #`(<@:({."1))@.x r
end.
)
mertens3=: -/ @ mr
1 mr y are the square-free integers less than or equal to y grouped by the lengths of the prime factorizations. 0 mr y (also mr y) are the number of integers in each group (#&> 1 mr y). For example:
mr 45 1 14 12 2 -/ mr 45 _3 mertens3 45 _3 1 mr 45 ┌─┬─────────────────────────────────────┬──────────────────────────────────┬─────┐ │1│2 3 5 7 11 13 17 19 23 29 31 37 41 43│6 10 14 22 26 34 38 15 21 33 39 35│30 42│ └─┴─────────────────────────────────────┴──────────────────────────────────┴─────┘ (1 mr 45) -: */"1&.> 1 mq 45 1 (1&mr -: */"1&.>@(1&mq))"0 >: 2 10 ?@$ 1e6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
mv y exploits 1 mr y to calculate the Mertens function from 1 to y . It is much more efficient than the equivalent mertens3"0 >:i.y . Thus:
mv=: 3 : '+/\ ((#&>t)#($t)$1 _1) (<:;t=. 1 mr y)} y$0' mv 30 1 0 _1 _1 _2 _1 _2 _2 _2 _1 _2 _2 _3 _2 _1 _1 _2 _2 _3 _3 _2 _1 _2 _2 _2 _1 _1 _1 _2 _3 mertens3"0 >: i.30 1 0 _1 _1 _2 _1 _2 _2 _2 _1 _2 _2 _3 _2 _1 _1 _2 _2 _3 _3 _2 _1 _2 _2 _2 _1 _1 _1 _2 _3 (mv -: mertens3"0@:>:@:i.) >:?1e4 1 load 'plot' plot (>:i.1e6);mv 1e6
Benchmark
expression time space mertens0 1e6 5.60969 8.59231e7 mertens1 1e6 2.89146 1.78297e7 mertens2 1e6 0.14057 1.09080e7 mertens3 1e6 0.13749 9.02624e6 mertens0 1e7 mertens1 1e7 60.78100 1.51000e8 mertens2 1e7 1.64406 9.17554e7 mertens3 1e7 1.59134 7.49949e7
Collected Definitions
mobius =: */ @: - @: ~: @: q:
mertens0=: +/ @: (mobius"0) @: >: @: i.
mertens1=: 3 : 0
p=. i.@(1&>.)&.(p:^:_1) 1+<.%:y
b=. -. > ([ +. #@[ $ {.&1@])~&.>/(<"0 *:p),<(1+y)$0
+/ _1 ^ #@q: I. b
)
mq=: 0&$: : (4 : 0)
p=. i.&.(p:^:_1) 10>.1+<.y%2
q=. ,. p {.~ p (<:i.0:) <.%:y
z=. (2<.y) $ (1,p:^:_1)`((i.1 0) ; [: ,. i.&.(p:^:_1))@.x >:y
for. i. 0 >. <: (*/\p:i.20) I. >:y do.
p=. p {.~ p (<:i.0:) <.y%*/{.q
q=. q #~ ({:"1 q)<y%*/"1 q
j=. 1 + p i. {:"1 q
k=. j -~ p I. <.1+y%*/"1 q
q=. (k#q) ,. (j,:"0 k) ;@:(<;.0) p
z=. z, #`<@.x q
end.
)
mertens2=: -/ @ mq
nextq=: 4 : 0
p=. i.&.(p:^:_1) 1+<.x%*/{.y
y=. y #~ ({:"1 y)<x%*/"1 y
j=. (**/$y) * 1 + p i. {:"1 y
k=. j -~ p I. <.1+x%*/"1 y
(k#y) ,. (j,:"0 k) ;@:(<;.0) p
)
qn =: 1 + (*/\p:i.20) I. >:
mq1=: 3 : 'y nextq&.>^:(<qn y) <i.1 0'
mqcheck=: 0&$: : (4 : 0)
z=. x mq y
assert. ($z) = ,1 + (*/\p:i.20) I. >:y
if. 0=x do.
assert. (0 < z) *. z = <.z
assert. (mertens0 y) = -/ z
assert. 1={.z
assert. (((1<y)$1){z) = p:^:_1 >:y
else.
assert. (32=3!:0 z) *. 2=#@$&>z
assert. (i.#z) -: {:@$&> z
assert. * #&> z
assert. 1 p: ; ,&.> z
assert. z -: ~. &.> z
assert. z -: ~."1&.> z
assert. z -: /:~ &.> z
assert. z -: /:~"1&.> z
assert. y >: ; */"1&.> z
assert. (0 mq y) -: #&> z
assert. (mertens0 y) = -/ #&> z
assert. ({.z) = <i.1 0
assert. (((1<y)$1){z) = <,.i.&.(p:^:_1) >:y
assert. ({.&.>}.z) -: <\p:i.<:#z
end.
)
mr=: 0&$: : (4 : 0)
p=. i.&.(p:^:_1) 1+<.y%2
r=. ,. i.&.(p:^:_1) 1+<.%:y
z=. (2<.y) $ (1,p:^:_1)`((,1) ; i.&.(p:^:_1))@.x >:y
for. i. 0 >. <: (*/\p:i.20) I. >:y do.
p=. p {.~ p (<:i.0:) <.y%{.{.r
r=. r #~ ({:"1 r)<y%{."1 r
j=. 1 + p i. {:"1 r
k=. j -~ p I. <.1+y%{."1 r
r=. (k#{."1 r) (* ,. ]) (j,:"0 k) ;@:(<;.0) p
z=. z, #`(<@:({."1))@.x r
end.
)
mertens3=: -/ @ mr
mv=: 3 : '+/\ ((#&>t)#($t)$1 _1) (<:;t=. 1 mr y)} y$0'
See also
Contributed by Roger Hui.