Mertens function: Difference between revisions
m (→{{header|Perl 6}}: added zkl header) |
|||
| Line 240: | Line 240: | ||
Crosses zero 59 times between 1 and 1000</pre> |
Crosses zero 59 times between 1 and 1000</pre> |
||
=={{header|zkl}}== |
|||
<lang zkl></lang> |
|||
<lang zkl></lang> |
|||
{{out}} |
|||
<pre> |
|||
</pre> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
Revision as of 01:11, 26 January 2020
The Mertens function M(x) is the count of square-free integers up to x that have an even number of prime factors, minus the count of those that have an odd number.
It is an extension of the Möbius function. Given the Möbius function μ(n), the Mertens function M(x) is the sum of the Möbius numbers from n == 1 through n == x.
- Task
- Write a routine (function, procedure, whatever) to find the Mertens number for any positive integer x.
- Use that routine to find and display here, on this page, at least the first 99 terms in a grid layout. (Not just one long line or column of numbers.)
- Use that routine to find and display here, on this page, the number of times the Mertens function sequence is equal to zero in the range M(1) through M(1000).
- Use that routine to find and display here, on this page, the number of times the Mertens function sequence crosses zero in the range M(1) through M(1000). (Crossing defined as this term equal to zero but preceding term not.)
- See also
This is not code golf. The stackexchange link is provided as an algorithm reference, not as a guide.
- Related tasks
Factor
<lang factor>USING: formatting grouping io kernel math math.extras math.ranges math.statistics prettyprint sequences ;
! Take the cumulative sum of the mobius sequence to avoid ! summing lower terms over and over.
- mertens-upto ( n -- seq ) [1,b] [ mobius ] map cum-sum ;
"The first 199 terms of the Mertens sequence:" print 199 mertens-upto " " prefix 20 group [ [ "%3s" printf ] each nl ] each nl
"In the first 1,000 terms of the Mertens sequence there are:" print 1000 mertens-upto [ [ zero? ] count bl pprint bl "zeros." print ] [
2 <clumps> [ first2 [ 0 = not ] [ zero? ] bi* and ] count bl pprint bl "zero crossings." print
] bi</lang>
- Output:
The first 199 terms of the Mertens sequence:
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 -2 -2 -3 -3 -2 -2 -1 0 -1
-1 -2 -1 -1 -1 0 -1 -2 -2 -1 -2 -3 -3 -4 -3 -3 -3 -2 -3 -4
-4 -4 -3 -4 -4 -3 -2 -1 -1 -2 -2 -1 -1 0 1 2 2 1 1 1
1 0 -1 -2 -2 -3 -2 -3 -3 -4 -5 -4 -4 -5 -6 -5 -5 -5 -4 -3
-3 -3 -2 -1 -1 -1 -1 -2 -2 -1 -2 -3 -3 -2 -1 -1 -1 -2 -3 -4
-4 -3 -2 -1 -1 0 1 1 1 0 0 -1 -1 -1 -2 -1 -1 -2 -1 0
0 1 1 0 0 -1 0 -1 -1 -1 -2 -2 -2 -3 -4 -4 -4 -3 -2 -3
-3 -4 -5 -4 -4 -3 -4 -3 -3 -3 -4 -5 -5 -6 -5 -6 -6 -7 -7 -8
In the first 1,000 terms of the Mertens sequence there are:
92 zeros.
59 zero crossings.
Go
<lang go>package main
import "fmt"
func mertens(to int) ([]int, int, int) {
if to < 1 {
to = 1
}
merts := make([]int, to+1)
primes := []int{2}
var sum, zeros, crosses int
for i := 1; i <= to; i++ {
j := i
cp := 0 // counts prime factors
spf := false // true if there is a square prime factor
for _, p := range primes {
if p > j {
break
}
if j%p == 0 {
j /= p
cp++
}
if j%p == 0 {
spf = true
break
}
}
if cp == 0 && i > 2 {
cp = 1
primes = append(primes, i)
}
if !spf {
if cp%2 == 0 {
sum++
} else {
sum--
}
}
merts[i] = sum
if sum == 0 {
zeros++
if i > 1 && merts[i-1] != 0 {
crosses++
}
}
}
return merts, zeros, crosses
}
func main() {
merts, zeros, crosses := mertens(1000)
fmt.Println("Mertens sequence - First 199 terms:")
for i := 0; i < 200; i++ {
if i == 0 {
fmt.Print(" ")
continue
}
if i%20 == 0 {
fmt.Println()
}
fmt.Printf(" % d", merts[i])
}
fmt.Println("\n\nEquals zero", zeros, "times between 1 and 1000")
fmt.Println("\nCrosses zero", crosses, "times between 1 and 1000")
}</lang>
- Output:
Mertens sequence - First 199 terms:
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 -2 -2 -3 -3 -2 -2 -1 0 -1
-1 -2 -1 -1 -1 0 -1 -2 -2 -1 -2 -3 -3 -4 -3 -3 -3 -2 -3 -4
-4 -4 -3 -4 -4 -3 -2 -1 -1 -2 -2 -1 -1 0 1 2 2 1 1 1
1 0 -1 -2 -2 -3 -2 -3 -3 -4 -5 -4 -4 -5 -6 -5 -5 -5 -4 -3
-3 -3 -2 -1 -1 -1 -1 -2 -2 -1 -2 -3 -3 -2 -1 -1 -1 -2 -3 -4
-4 -3 -2 -1 -1 0 1 1 1 0 0 -1 -1 -1 -2 -1 -1 -2 -1 0
0 1 1 0 0 -1 0 -1 -1 -1 -2 -2 -2 -3 -4 -4 -4 -3 -2 -3
-3 -4 -5 -4 -4 -3 -4 -3 -3 -3 -4 -5 -5 -6 -5 -6 -6 -7 -7 -8
Equals zero 92 times between 1 and 1000
Crosses zero 59 times between 1 and 1000
Julia
<lang julia>using Primes
function moebius(n::Integer)
@assert n > 0 m(p, e) = p == 0 ? 0 : e == 1 ? -1 : 0 reduce(*, m(p, e) for (p, e) in factor(n) if p ≥ 0; init=1)
end μ(n) = moebius(n)
mertens(x) = sum(n -> μ(n), 1:x) M(x) = mertens(x)
print("First 99 terms of the Mertens function for positive integers:\n ") for n in 1:99
print(lpad(M(n), 3), n % 10 == 9 ? "\n" : "")
end println("\nM(n) equals 0: ", sum(n -> M(n) == 0, 1:1000), " times from 1 to 1000.") println("\nM(n) \"crosses\" 0: ",
sum(n -> M(n) == 0 && M(n-1) != 0, 1:1000), " times from 1 to 1000.")
</lang>
- Output:
First 99 terms of the Mertens function for positive integers:
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 -2 -2 -3 -3 -2 -2 -1 0 -1
-1 -2 -1 -1 -1 0 -1 -2 -2 -1
-2 -3 -3 -4 -3 -3 -3 -2 -3 -4
-4 -4 -3 -4 -4 -3 -2 -1 -1 -2
-2 -1 -1 0 1 2 2 1 1 1
M(n) equals 0: 92 times from 1 to 1000.
M(n) "crosses" 0: 59 times from 1 to 1000.
Perl 6
Mertens number is not defined for n == 0. Perl 6 arrays are indexed from 0 so store a blank value at position zero to keep x and M(x) aligned.
<lang perl6>use Prime::Factor;
sub μ (Int \n) {
my @p = prime-factors(n); +@p == +Bag(@p).keys ?? +@p %% 2 ?? 1 !! -1 !! 0
}
my @mertens = lazy [\+] flat ' ', 1, (2..*).map: -> \n { μ(n) };
put "Mertens sequence - First 199 terms:\n",
@mertens[^200]».fmt('%3s').rotor(20).join("\n"),
"\n\nEquals zero ", +@mertens[1..1000].grep( !* ),
' times between 1 and 1000', "\n\nCrosses zero ",
+@mertens[1..1000].kv.grep( {!$^v and @mertens[$^k]} ),
' times between 1 and 1000';</lang>
- Output:
Mertens sequence - First 199 terms:
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 -2 -2 -3 -3 -2 -2 -1 0 -1
-1 -2 -1 -1 -1 0 -1 -2 -2 -1 -2 -3 -3 -4 -3 -3 -3 -2 -3 -4
-4 -4 -3 -4 -4 -3 -2 -1 -1 -2 -2 -1 -1 0 1 2 2 1 1 1
1 0 -1 -2 -2 -3 -2 -3 -3 -4 -5 -4 -4 -5 -6 -5 -5 -5 -4 -3
-3 -3 -2 -1 -1 -1 -1 -2 -2 -1 -2 -3 -3 -2 -1 -1 -1 -2 -3 -4
-4 -3 -2 -1 -1 0 1 1 1 0 0 -1 -1 -1 -2 -1 -1 -2 -1 0
0 1 1 0 0 -1 0 -1 -1 -1 -2 -2 -2 -3 -4 -4 -4 -3 -2 -3
-3 -4 -5 -4 -4 -3 -4 -3 -3 -3 -4 -5 -5 -6 -5 -6 -6 -7 -7 -8
Equals zero 92 times between 1 and 1000
Crosses zero 59 times between 1 and 1000
REXX
(Currently working on the 2nd and 3rd task's requirements.)
Programming note: This REXX version supports the specifying of the low and high values to be generated,
as well as the group size for the grid (it can be specified as 1 which will show a vertical list).
A null value will be shown as a bullet (•) when showing the Möbius and/or Mertens value of for zero (which can be changed easily). <lang rexx>/*REXX pgm computes & shows a value grid of the Mertens function for a range of integers*/ /*───────────────────────────────────────────────{function is named after Franz Mertens}*/ parse arg LO HI grp . /*obtain optional arguments from the CL*/ if LO== | LO=="," then LO= 0 /*Not specified? Then use the default.*/ if HI== | HI=="," then HI= 199 /* " " " " " " */ if grp== | grp=="," then grp= 20 /* " " " " " " */ !.=. /* ______ */ call genP HI /*generate primes up to the √ HI */ say center(' The Mertens sequence from ' LO " ──► " HI" ", max(50, grp*3), '═') /*title*/ $= /*variable holds output grid of GRP #s.*/
do j=LO to HI; $= $ right( Mertens(j), 2) /*process some numbers from LO ──► HI.*/
if words($)==grp then do; say substr($, 2); $= /*show grid if fully populated,*/
end /* and nullify it for more #s.*/
end /*j*/ /*for small grids, using wordCnt is OK.*/
if $\== then say substr($, 2) /*handle any residual numbers not shown*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ Mertens: procedure expose @. !.; parse arg n /*obtain a integer to be tested for mu.*/
if n<1 then return '∙' /*handle special cases of non─positive#*/
m= 0 /*the sum of all the MU's (so far). */
do k=1 for n; m= m + mobius(k) /*sum the MU's up to N. */
end /*k*/ /* [+] mobius function uses memoization*/
return m /*return the sum of all the MU's. */
/*──────────────────────────────────────────────────────────────────────────────────────*/ mobius: procedure expose @. !.; parse arg x 1 ox /*obtain a integer to be tested for mu.*/
if !.x\==. then return !.x /*X computed before? Return that value*/
if x<1 then return '∙' /*handle special case of non-positive #*/
#= 0 /*start with a mu value of zero. */
do k=1; p= @.k /*get the Kth (pre─generated) prime.*/
if p>x then leave /*prime (P) > X? Then we're done. */
if p*p>x then do; #= #+1; leave /*prime (P**2 > X? Bump # and leave.*/
end
if x//p==0 then do; #= #+1 /*X divisible by P? Bump mu number. */
x= x % p /* Divide by prime. */
if x//p==0 then return 0 /*X÷by P? Then return zero*/
end
end /*k*/ /*# (below) is almost always small, <9*/
!.ox= -1 ** #; return !.ox /*raise -1 to the mu power, memoize it.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/ genP: @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6= 13; nP=6 /*assign low primes; # primes. */
do lim=nP until lim*lim>=HI; end /*only keep primes up to the sqrt(HI). */
do j=@.nP+4 by 2 to HI /*only find odd primes from here on. */
if j// 3==0 then iterate /*is J divisible by #3 Then not prime*/
parse var j -1 _;if _==5 then iterate /*Is last digit a "5"? " " " */
if j// 7==0 then iterate /*is J divisible by 7? " " " */
if j//11==0 then iterate /* " " " " 11? " " " */
if j//13==0 then iterate /*is " " " 13? " " " */
do k=7 while k*k<=j /*divide by some generated odd primes. */
if j // @.k==0 then iterate j /*Is J divisible by P? Then not prime*/
end /*k*/ /* [↓] a prime (J) has been found. */
nP= nP+1; if nP<=HI then @.nP=j /*bump prime count; assign prime to @.*/
end /*j*/; return</lang>
- output when using the default inputs:
Output note: note the use of a bullet (•) to signify that a "null" is being shown (for the 0th entry).
══════════ The Mertens sequence from 0 ──► 199 ══════════ ∙ 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 -2 -2 -3 -3 -2 -2 -1 0 -1 -1 -2 -1 -1 -1 0 -1 -2 -2 -1 -2 -3 -3 -4 -3 -3 -3 -2 -3 -4 -4 -4 -3 -4 -4 -3 -2 -1 -1 -2 -2 -1 -1 0 1 2 2 1 1 1 1 0 -1 -2 -2 -3 -2 -3 -3 -4 -5 -4 -4 -5 -6 -5 -5 -5 -4 -3 -3 -3 -2 -1 -1 -1 -1 -2 -2 -1 -2 -3 -3 -2 -1 -1 -1 -2 -3 -4 -4 -3 -2 -1 -1 0 1 1 1 0 0 -1 -1 -1 -2 -1 -1 -2 -1 0 0 1 1 0 0 -1 0 -1 -1 -1 -2 -2 -2 -3 -4 -4 -4 -3 -2 -3 -3 -4 -5 -4 -4 -3 -4 -3 -3 -3 -4 -5 -5 -6 -5 -6 -6 -7 -7 -8