Essays/Mertens Function

From J Wiki
Jump to: navigation, search

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

Mertens.jpg

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.