Mertens function: Difference between revisions
Not a robot (talk | contribs) (Add APL) |
Not a robot (talk | contribs) (Add Fortran) |
||
Line 711: | Line 711: | ||
59 zero crossings. |
59 zero crossings. |
||
</pre> |
</pre> |
||
=={{header|Fortran}}== |
|||
<lang Fortran> program Mertens |
|||
implicit none |
|||
integer M(1000), n, k, zero, cross |
|||
C Generate Mertens numbers |
|||
M(1) = 1 |
|||
do 10 n=2, 1000 |
|||
M(n) = 1 |
|||
do 10 k=2, n |
|||
M(n) = M(n) - M(n/k) |
|||
10 continue |
|||
C Print table |
|||
write (*,"('The first 99 Mertens numbers are:')") |
|||
write (*,"(' ')",advance='no') |
|||
k = 9 |
|||
do 20 n=1, 99 |
|||
write (*,'(I3)',advance='no') M(n) |
|||
k = k-1 |
|||
if (k .EQ. 0) then |
|||
k=10 |
|||
write (*,*) |
|||
end if |
|||
20 continue |
|||
C Calculate zeroes and crossings |
|||
zero = 0 |
|||
cross = 0 |
|||
do 30 n=2, 1000 |
|||
if (M(n) .EQ. 0) then |
|||
zero = zero + 1 |
|||
if (M(n-1) .NE. 0) cross = cross+1 |
|||
end if |
|||
30 continue |
|||
40 format("M(N) is zero ",I2," times.") |
|||
write (*,40) zero |
|||
50 format("M(N) crosses zero ",I2," times.") |
|||
write (*,50) cross |
|||
end program</lang> |
|||
{{out}} |
|||
<pre>The first 99 Mertens numbers are: |
|||
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) is zero 92 times. |
|||
M(N) crosses zero 59 times.</pre> |
|||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |
Revision as of 16:56, 28 January 2021
You are encouraged to solve this task according to the task description, using any language you may know.
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
8080 Assembly
<lang 8080asm>MAX: equ 1000 ; Amount of numbers to generate org 100h ;;; Generate Mertens numbers lxi b,1 ; Start at place 1; BC = current Mertens number lxi h,MM ; First one is 1 dad b mvi m,1 outer: inx b ; Next Mertens number lxi h,MM dad b mvi m,1 ; Initialize at 1 lxi d,2 ; DE = inner loop counter ('k'), starts at 2 ;;; Now we need to find BC/DE, but there is no hardware divide ;;; We also need to be somewhat clever so it doesn't take forever inner: push d ; Keep both loop counters safe on the stack push b xchg ; Divisor in HL mov d,b ; Dividend in DE mov e,c lxi b,100h ; B = counter, C = zero double: dad h ; Double divisor inr b ; Increment counter call cdehl ; Dividend <= divisor? jnc double ; If so, keep doubling mov a,b ; Keep counter mov b,c ; BC = 0 push b ; Push result variable on stakc (initial 0) mov b,a ; Restore counter xchg ; HL = dividend, DE = doubled divisor subtr: mov a,l ; Try HL -= DE sub e mov l,a mov a,h sbb d mov h,a xthl ; Get result accumulator from stack cmc ; Flip borrow mov a,l ; Rotate into result ral mov l,a mov a,h ral mov h,a mov a,l ; Retrieve flag rar xthl ; Retrieve rest of divisor jc $+4 ; If borrow, dad d ; Add dividend back into divisor xra a ; DE >> 1 ora d rar mov d,a mov a,e rar mov e,a dcr b ; Are we there yet? jnz subtr ; If not, try another subtraction pop h ; HL = quotient ;;; Division is done, do lookup and subraction lxi d,MM ; Look up M[outer/inner] dad d mov e,m ; E = M[BC/DE] pop b ; Restore BC (n) lxi h,MM dad b mov a,m ; A = M[BC] sub e ; A = M[BC] - M[BC/DE] mov m,a ; M[BC] = A pop d ; Restore DE (k) ;;; Update loops inx d ; k++ call cbcde ; DE <= BC? jnc inner lxi h,MAX call chlbc ; BC <= MAX? jnc outer ;;; Print table lxi d,frst99 call puts lxi h,MM+1 ; Start of Merten numbers mvi c,9 ; Column counter table: mov a,m ; Get Merten number ana a ; Set flags mvi b,' ' ; Space jp prtab ; If positive, print space-number-space mvi b,'-' ; Otherwise, print minus sign cma ; And negate the number (make positive) inr a prtab: adi '0' ; Make ASCII digit mov d,a ; Keep number mov a,b ; Print space or minus sign call putc mov a,d ; Restore number call putc ; Print number mvi a,' ' ; Print space call putc dcr c ; Decrement column counter jnz tnext lxi d,nl ; End of columns - print newline call puts mvi c,10 ; Column counter tnext: inx h ; Table done? mov a,l cpi 100 jnz table ; If not, keep going ;;; Find zeroes and crossings lxi b,0 ; B=zeroes, C=crossings lxi d,MAX ; Counter lxi h,MM+1 count: mov a,m ; Get number ana a ; Zero? jnz cnext inr b ; If so, add zero dcx h ; Previous number also zero? mov a,m inx h ana a jz cnext inr c ; If not, add crossiong cnext: inx h dcx d mov a,d ora e jnz count lxi d,zero ; Print zeroes call puts mov a,b call puta lxi d,cross ; Print crossings call puts mov a,c call puta lxi d,tms jmp puts ;;; Print character in A using CP/M, keeping registers putc: push b push d push h mov e,a mvi c,2 call 5 jmp resrgs ;;; Print number in A, keeping registers puta: push b push d push h lxi h,num putad: mvi c,-1 putal: inr c sui 10 jnc putal adi 10+'0' dcx h mov m,a mov a,c ana a jnz putad xchg mvi c,9 call 5 jmp resrgs ;;; Print string in DE using CP/M, keeping registers puts: push b push d push h mvi c,9 call 5 resrgs: pop h pop d pop b ret cdehl: mov a,d cmp h rnz mov a,e cmp l ret cbcde: mov a,b cmp d rnz mov a,c cmp e ret chlbc: mov a,h cmp b rnz mov a,l cmp c ret ;;; Strings db '***' num: db '$' frst99: db 'First 99 Mertens numbers:',13,10,' $' nl: db 13,10,'$' zero: db 'M(N) is zero $' cross: db ' times.',13,10,'M(N) crosses zero $' tms: db ' times.$' ;;; Numbers are stored page-aligned after program MM: equ ($/256)*256+256 </lang>
- Output:
First 99 Mertens numbers: 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) is zero 92 times. M(N) crosses zero 59 times.
8086 Assembly
<lang asm>MAX: equ 1000 ; Amount of Mertens numbers to generate puts: equ 9 ; MS-DOS syscall to print a string putch: equ 2 ; MS-DOS syscall to print a character cpu 8086 org 100h section .text ;;; Generate Mertens numbers mov bx,M ; BX = pointer to start of Mertens numbers mov si,1 ; Current Mertens number mov [si+bx],byte 1 ; First Mertens number is 1 outer: inc si ; Next Mertens number mov [si+bx],byte 1 ; Starts out at 1... mov cx,2 ; CX = from 2 to current number, inner: mov ax,si ; Divide current number, xor dx,dx div cx ; By CX mov di,ax mov al,[di+bx] ; Get value at that location sub [si+bx],al ; Subtract from current number inc cx cmp cx,si jbe inner cmp si,MAX jbe outer ;;; Print the table mov dx,frst99 ; First string call outstr mov si,1 ; Start at index 1 mov dh,9 ; Column count table: mov cl,[si+bx] ; Get item test cl,cl mov dl,' ' jns .print ; Positive? mov dl,'-' ; Otherwise, it is negative, neg cl ; print ' ' and negate .print: call putc ; Print space or minus add cl,'0' ; Add ASCII 0 mov dl,cl call putc ; Print number mov dl,' ' call putc ; Print space dec dh ; One less column left jnz .next mov dx,nl ; Print newline call outstr mov dh,10 .next: inc si ; Done yet? cmp si,100 jb table ; If not, print next item from table ;;; Calculate zeroes and crossings xor cx,cx ; CL = zeroes, CH = crossings mov si,1 mov al,[si+bx] ; AL = current item zc: inc si mov ah,al ; AH = previous item mov al,[si+bx] test al,al ; Zero? jnz .next inc cx ; Then increment zero counter test ah,ah ; Previous one also zero? jz .next inc ch ; Then increment crossing counter .next: cmp si,MAX ; Done yet? jbe zc ;;; Print zeroes and crossings mov dx,zero call outstr mov al,cl call putal mov dx,cross call outstr mov al,ch call putal mov dx,tms jmp outstr putc: mov ah,putch ; Print character int 21h ret ;;; Print AL in decimal format putal: mov di,num .loop: aam ; Extract digit add al,'0' ; Store digit dec di mov [di],al mov al,ah ; Rest of number test al,al ; Done? jnz .loop ; If not, get more digits mov dx,di ; Otherwise, print string outstr: mov ah,puts int 21h ret section .data db '***' ; Number output placeholder num: db '$' frst99: db 'First 99 Mertens numbers:',13,10,' $' nl: db 13,10,'$' zero: db 'M(N) is zero $' cross: db ' times.',13,10,'M(N) crosses zero $' tms: db ' times.$' section .bss mm: resb MAX ; Mertens numbers M: equ mm-1 ; 1-based indexing</lang>
- Output:
First 99 Mertens numbers: 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) is zero 92 times. M(N) crosses zero 59 times.
APL
<lang APL>mertens←{
step ← {⍵,-⍨/⌽1,⍵[⌊n÷1↓⍳n←1+≢⍵]} m1000 ← step⍣999⊢,1 zero ← m1000+.=0 cross ← +/(~∧1⌽⊢)m1000≠0 ⎕←'First 99 Mertens numbers:' ⎕←10 10⍴'∘',m1000 ⎕←'M(N) is zero ',(⍕zero),' times.' ⎕←'M(N) crosses zero ',(⍕cross),' times.'
}</lang>
- Output:
First 99 Mertens numbers: ∘ 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) is zero 92 times. M(N) crosses zero 59 times.
BASIC
Note that if you actually try this on a real 8-bit micro, this will take literal hours to run. Interpreted BASIC is not fast.
For comparison, the 8080 and 8086 assembly versions above run in around 4 minutes and 15 seconds respectively, on real hardware. This took nearly 4 hours on the 8080 (using MBASIC) and a little over 1 hour on the 8086 (using GWBASIC).
<lang BASIC>10 DEFINT C,Z,N,K,M: DIM M(1000) 20 M(1)=1 30 FOR N=2 TO 1000 40 M(N)=1 50 FOR K=2 TO N: M(N) = M(N)-M(INT(N/K)): NEXT 60 NEXT 70 PRINT "First 99 Mertens numbers:" 80 PRINT " "; 90 FOR N=1 TO 99 100 PRINT USING "###";M(N); 110 IF N MOD 10 = 9 THEN PRINT 120 NEXT 130 C=0: Z=0 140 FOR N=1 TO 1000 150 IF M(N)=0 THEN Z=Z+1: IF M(N-1)<>0 THEN C=C+1 160 NEXT 170 PRINT "M(N) is zero";Z;"times." 180 PRINT "M(N) crosses zero";C;"times." 190 END</lang>
- Output:
First 99 Mertens numbers: 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) is zero 92 times. M(N) crosses zero 59 times.
C
<lang c>#include <stdio.h>
- include <stdlib.h>
int* mertens_numbers(int max) {
int* m = malloc((max + 1) * sizeof(int)); if (m == NULL) return m; m[1] = 1; for (int n = 2; n <= max; ++n) { m[n] = 1; for (int k = 2; k <= n; ++k) m[n] -= m[n/k]; } return m;
}
int main() {
const int max = 1000; int* mertens = mertens_numbers(max); if (mertens == NULL) { fprintf(stderr, "Out of memory\n"); return 1; } printf("First 199 Mertens numbers:\n"); const int count = 200; for (int i = 0, column = 0; i < count; ++i) { if (column > 0) printf(" "); if (i == 0) printf(" "); else printf("%2d", mertens[i]); ++column; if (column == 20) { printf("\n"); column = 0; } } int zero = 0, cross = 0, previous = 0; for (int i = 1; i <= max; ++i) { int m = mertens[i]; if (m == 0) { ++zero; if (previous != 0) ++cross; } previous = m; } free(mertens); printf("M(n) is zero %d times for 1 <= n <= %d.\n", zero, max); printf("M(n) crosses zero %d times for 1 <= n <= %d.\n", cross, max); return 0;
}</lang>
- Output:
First 199 Mertens numbers: 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 M(n) is zero 92 times for 1 <= n <= 1000. M(n) crosses zero 59 times for 1 <= n <= 1000.
C++
<lang cpp>#include <iomanip>
- include <iostream>
- include <map>
class mertens_calculator { public:
int mertens_number(int);
private:
std::map<int, int> cache_;
};
int mertens_calculator::mertens_number(int n) {
auto i = cache_.find(n); if (i != cache_.end()) return i->second; int m = 1; for (int k = 2; k <= n; ++k) m -= mertens_number(n/k); cache_.emplace(n, m); return m;
}
void print_mertens_numbers(mertens_calculator& mc, int count) {
int column = 0; for (int i = 0; i < count; ++i) { if (column > 0) std::cout << ' '; if (i == 0) std::cout << " "; else std::cout << std::setw(2) << mc.mertens_number(i); ++column; if (column == 20) { std::cout << '\n'; column = 0; } }
}
int main() {
mertens_calculator mc; std::cout << "First 199 Mertens numbers:\n"; print_mertens_numbers(mc, 200); int zero = 0, cross = 0, previous = 0; for (int i = 1; i <= 1000; ++i) { int m = mc.mertens_number(i); if (m == 0) { ++zero; if (previous != 0) ++cross; } previous = m; } std::cout << "M(n) is zero " << zero << " times for 1 <= n <= 1000.\n"; std::cout << "M(n) crosses zero " << cross << " times for 1 <= n <= 1000.\n"; return 0;
}</lang>
- Output:
First 199 Mertens numbers: 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 M(n) is zero 92 times for 1 <= n <= 1000. M(n) crosses zero 59 times for 1 <= n <= 1000.
Cowgol
<lang cowgol>include "cowgol.coh";
const MAX := 1000;
- Table output
sub printtab(n: int8) is
if n<0 then print_char('-'); n := -n; else print_char(' '); end if; print_char(n as uint8 + '0'); print_char(' ');
end sub;
- Generate Merten numbers
var M: int8[MAX+1]; M[0] := 0; M[1] := 1;
var n: @indexof M := 2; while n < @sizeof M loop
M[n] := 1; var k: @indexof M := 2; while k <= n loop M[n] := M[n] - M[n/k]; k := k + 1; end loop; n := n + 1;
end loop;
- Find zeroes and crossings
var zero: uint8 := 0; var cross: uint8 := 0; n := 1; while n < @sizeof M loop
if M[n] == 0 then zero := zero + 1; if M[n-1] != 0 then cross := cross + 1; end if; end if; n := n + 1;
end loop;
- Print table
print("The first 99 Mertens numbers are:\n"); print(" "); n := 1; var col: uint8 := 9; while n < 100 loop
printtab(M[n]); col := col - 1; if col == 0 then print_nl(); col := 10; end if; n := n + 1;
end loop;
print("M(n) is zero "); print_i8(zero); print(" times\n"); print("M(n) crosses zero "); print_i8(cross); print(" times\n");</lang>
- Output:
The first 99 Mertens numbers are: 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) is zero 92 times M(n) crosses zero 59 times
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.
Fortran
<lang Fortran> program Mertens
implicit none integer M(1000), n, k, zero, cross
C Generate Mertens numbers
M(1) = 1 do 10 n=2, 1000 M(n) = 1 do 10 k=2, n M(n) = M(n) - M(n/k) 10 continue
C Print table
write (*,"('The first 99 Mertens numbers are:')") write (*,"(' ')",advance='no') k = 9 do 20 n=1, 99 write (*,'(I3)',advance='no') M(n) k = k-1 if (k .EQ. 0) then k=10 write (*,*) end if 20 continue
C Calculate zeroes and crossings
zero = 0 cross = 0 do 30 n=2, 1000 if (M(n) .EQ. 0) then zero = zero + 1 if (M(n-1) .NE. 0) cross = cross+1 end if 30 continue 40 format("M(N) is zero ",I2," times.") write (*,40) zero 50 format("M(N) crosses zero ",I2," times.") write (*,50) cross end program</lang>
- Output:
The first 99 Mertens numbers are: 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) is zero 92 times. M(N) crosses zero 59 times.
FreeBASIC
<lang freebasic>function padto( i as ubyte, j as integer ) as string
return wspace(i-len(str(j)))+str(j)
end function
dim as integer M( 1 to 1000 ), n, col, k, psum dim as integer num_zeroes = 0, num_cross = 0 dim as string outstr
M(1) = 1 for n = 2 to 1000
psum = 0 for k = 2 to n psum += M(int(n/k)) next k M(n) = 1 - psum if M(n) = 0 then num_zeroes += 1 if M(n-1)<>0 then num_cross += 1 end if end if
next n
print using "There are ### zeroes in the range 1 to 1000."; num_zeroes print using "There are ### crossings in the range 1 to 1000."; num_cross print "The first 100 Mertens numbers are: "
for n=1 to 100
outstr += padto(3, M(n))+" " if n mod 10 = 0 then print outstr outstr = "" end if
next n</lang>
- Output:
There are 92 zeroes in the range 1 to 1000. There are 59 crossings in the range 1 to 1000. The first 100 Mertens numbers are: 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
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
Haskell
<lang haskell>import Data.List.Split (chunksOf) import qualified Data.MemoCombinators as Memo import Math.NumberTheory.Primes (unPrime, factorise) import Text.Printf (printf)
moebius :: Integer -> Int moebius = product . fmap m . factorise
where m (p, e) | unPrime p == 0 = 0 | e == 1 = -1 | otherwise = 0
mertens :: Integer -> Int mertens = Memo.integral (\n -> sum $ fmap moebius [1..n])
countZeros :: [Integer] -> Int countZeros = length . filter ((==0) . mertens)
crossesZero :: [Integer] -> Int crossesZero = length . go . fmap mertens
where go (x:y:xs) | y == 0 && x /= 0 = y : go (y:xs) | otherwise = go (y:xs) go _ = []
main :: IO () main = do
printf "The first 99 terms for M(1..99):\n\n " mapM_ (printf "%3d" . mertens) [1..9] >> printf "\n" mapM_ (\row -> mapM_ (printf "%3d" . mertens) row >> printf "\n") $ chunksOf 10 [10..99] printf "\nM(n) is zero %d times for 1 <= n <= 1000.\n" $ countZeros [1..1000] printf "M(n) crosses zero %d times for 1 <= n <= 1000.\n" $ crossesZero [1..1000]</lang>
- Output:
The first 99 terms for M(1..99): 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) is zero 92 times for 1 <= n <= 1000. M(n) crosses zero 59 times for 1 <= n <= 1000.
J
<lang J>mu =: 0:`(1 - 2 * 2|#@{.)@.(1: = */@{:)@(2&p:)"0 M =: +/@([: mu 1:+i.)
m1000 =: (M"0) 1+i.1000 zero =: +/ m1000 = 0 cross =: +/ (-.*.1:|]) m1000 ~: 0
echo 'The first 99 Merten numbers are' echo 10 10$ __, 99{.m1000 echo 'M(N) is zero ',(":zero),' times.' echo 'M(N) crosses zero ',(":cross),' times.' exit</lang>
- Output:
The first 99 Merten numbers are __ 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) is zero 92 times. M(N) crosses zero 0 times.
Java
<lang java> public class MertensFunction {
public static void main(String[] args) { System.out.printf("First 199 terms of the merten function are as follows:%n "); for ( int n = 1 ; n < 200 ; n++ ) { System.out.printf("%2d ", mertenFunction(n)); if ( (n+1) % 20 == 0 ) { System.out.printf("%n"); } } for ( int exponent = 3 ; exponent<= 8 ; exponent++ ) { int zeroCount = 0; int zeroCrossingCount = 0; int positiveCount = 0; int negativeCount = 0; int mSum = 0; int mMin = Integer.MAX_VALUE; int mMinIndex = 0; int mMax = Integer.MIN_VALUE; int mMaxIndex = 0; int nMax = (int) Math.pow(10, exponent); for ( int n = 1 ; n <= nMax ; n++ ) { int m = mertenFunction(n); mSum += m; if ( m < mMin ) { mMin = m; mMinIndex = n; } if ( m > mMax ) { mMax = m; mMaxIndex = n; } if ( m > 0 ) { positiveCount++; } if ( m < 0 ) { negativeCount++; } if ( m == 0 ) { zeroCount++; } if ( m == 0 && mertenFunction(n - 1) != 0 ) { zeroCrossingCount++; } } System.out.printf("%nFor M(x) with x from 1 to %,d%n", nMax); System.out.printf("The maximum of M(x) is M(%,d) = %,d.%n", mMaxIndex, mMax); System.out.printf("The minimum of M(x) is M(%,d) = %,d.%n", mMinIndex, mMin); System.out.printf("The sum of M(x) is %,d.%n", mSum); System.out.printf("The count of positive M(x) is %,d, count of negative M(x) is %,d.%n", positiveCount, negativeCount); System.out.printf("M(x) has %,d zeroes in the interval.%n", zeroCount); System.out.printf("M(x) has %,d crossings in the interval.%n", zeroCrossingCount); } } private static int MU_MAX = 100_000_000; private static int[] MU = null; private static int[] MERTEN = null; // Compute mobius and merten function via sieve private static int mertenFunction(int n) { if ( MERTEN != null ) { return MERTEN[n]; } // Populate array MU = new int[MU_MAX+1]; MERTEN = new int[MU_MAX+1]; MERTEN[1] = 1; int sqrt = (int) Math.sqrt(MU_MAX); for ( int i = 0 ; i < MU_MAX ; i++ ) { MU[i] = 1; } for ( int i = 2 ; i <= sqrt ; i++ ) { if ( MU[i] == 1 ) { // for each factor found, swap + and - for ( int j = i ; j <= MU_MAX ; j += i ) { MU[j] *= -i; } // square factor = 0 for ( int j = i*i ; j <= MU_MAX ; j += i*i ) { MU[j] = 0; } } } int sum = 1; for ( int i = 2 ; i <= MU_MAX ; i++ ) { if ( MU[i] == i ) { MU[i] = 1; } else if ( MU[i] == -i ) { MU[i] = -1; } else if ( MU[i] < 0 ) { MU[i] = 1; } else if ( MU[i] > 0 ) { MU[i] = -1; } sum += MU[i]; MERTEN[i] = sum; } return MERTEN[n]; }
} </lang>
- Output:
First 199 terms of the merten function are as follows: 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 For M(x) with x from 1 to 1,000 The maximum of M(x) is M(586) = 7. The minimum of M(x) is M(665) = -12. The sum of M(x) is -1,572. The count of positive M(x) is 254, count of negative M(x) is 654. M(x) has 92 zeroes in the interval. M(x) has 59 crossings in the interval. For M(x) with x from 1 to 10,000 The maximum of M(x) is M(8,511) = 35. The minimum of M(x) is M(9,861) = -43. The sum of M(x) is -20,409. The count of positive M(x) is 3,965, count of negative M(x) is 5,629. M(x) has 406 zeroes in the interval. M(x) has 256 crossings in the interval. For M(x) with x from 1 to 100,000 The maximum of M(x) is M(48,433) = 96. The minimum of M(x) is M(96,014) = -132. The sum of M(x) is -516,879. The count of positive M(x) is 47,830, count of negative M(x) is 50,621. M(x) has 1,549 zeroes in the interval. M(x) has 949 crossings in the interval. For M(x) with x from 1 to 1,000,000 The maximum of M(x) is M(992,998) = 311. The minimum of M(x) is M(926,265) = -368. The sum of M(x) is -14,244,200. The count of positive M(x) is 472,963, count of negative M(x) is 521,676. M(x) has 5,361 zeroes in the interval. M(x) has 3,269 crossings in the interval. For M(x) with x from 1 to 10,000,000 The maximum of M(x) is M(9,993,034) = 1,143. The minimum of M(x) is M(7,109,110) = -1,078. The sum of M(x) is -194,680,528. The count of positive M(x) is 4,938,188, count of negative M(x) is 5,049,266. M(x) has 12,546 zeroes in the interval. M(x) has 7,646 crossings in the interval. For M(x) with x from 1 to 100,000,000 The maximum of M(x) is M(92,418,127) = 3,290. The minimum of M(x) is M(76,015,339) = -3,448. The sum of M(x) is -608,757,258. The count of positive M(x) is 54,659,906, count of negative M(x) is 45,298,186. M(x) has 41,908 zeroes in the interval. M(x) has 25,525 crossings in the interval.
Julia
The OEIS A002321 reference suggests the Mertens function has a negative bias, which it does below 1 million, but this bias seems to switch to a positive bias by 1 billion. There may simply be large swings in the bias overall, which get larger and longer as the sequence continues. <lang julia>using Primes, Formatting
function moebius(n::Integer)
@assert n > 0 m(p, e) = p == 0 ? 0 : e == 1 ? -1 : 0 return 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
function maximinM(N)
z, cros, lastM, maxi, maxM, mini, minM, sumM, pos, neg = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 for i in 1:N m = μ(i) + lastM if m == 0 && lastM != 0 cros += 1 end sumM += m lastM = m if m > maxM maxi = i maxM = m elseif m < minM mini = i minM = m end if m > 0 pos += 1 elseif m < 0 neg += 1 else z += 1 end end println("\nFor M(x) with x from 1 to $(format(N, commas=true)):") println("The maximum of M(x) is M($(format(maxi, commas=true)) = $maxM.") println("The minimum of M(x) is M($(format(mini, commas=true))) = $minM.") println("The sum of M(x) is $(format(sumM, commas=true)).") println("The count of positive M(x) is $(format(pos, commas=true)), count of negative M(x) is $(format(neg, commas=true)).") println("M(x) has $(format(z, commas=true)) zeroes in the interval.") println("M(x) has $(format(cros, commas=true)) crossings in the interval.") diff = pos - neg if diff > 0 println("Positive M(x) exceed negative ones by $(format(diff, commas=true)).") else println("Negative M(x) exceed positive ones by $(format(-diff, commas=true)).") end
end
foreach(maximinM, (1000, 1_000_000, 1_000_000_000))
</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 For M(x) with x from 1 to 1,000: The maximum of M(x) is M(586 = 7. The minimum of M(x) is M(665) = -12. The sum of M(x) is -1,572. The count of positive M(x) is 254, count of negative M(x) is 654. M(x) has 92 zeroes in the interval. M(x) has 59 crossings in the interval. Negative M(x) exceed positive ones by 400. For M(x) with x from 1 to 1,000,000: The maximum of M(x) is M(992,998 = 311. The minimum of M(x) is M(926,265) = -368. The sum of M(x) is -14,244,200. The count of positive M(x) is 472,963, count of negative M(x) is 521,676. M(x) has 5,361 zeroes in the interval. M(x) has 3,269 crossings in the interval. Negative M(x) exceed positive ones by 48,713. For M(x) with x from 1 to 1,000,000,000: The maximum of M(x) is M(903,087,703 = 10246. The minimum of M(x) is M(456,877,618) = -8565. The sum of M(x) is 510,495,361,261. The count of positive M(x) is 510,200,302, count of negative M(x) is 489,658,577. M(x) has 141,121 zeroes in the interval. M(x) has 85,652 crossings in the interval. Positive M(x) exceed negative ones by 20,541,725.
MAD
<lang MAD> NORMAL MODE IS INTEGER
DIMENSION M(1000) M(1) = 1 THROUGH GENMRT, FOR N=2, 1, N.G.1000 M(N) = 1 THROUGH GENMRT, FOR K=2, 1, K.G.N
GENMRT M(N) = M(N) - M(N/K)
PRINT COMMENT $ FIRST 99 MERTEN NUMBERS ARE$ VECTOR VALUES F9 = $S3,9(I2,S1)*$ VECTOR VALUES F10 = $10(I2,S1)*$ PRINT FORMAT F9, M(1), M(2), M(3), M(4), M(5), M(6), 0 M(7), M(8), M(9) THROUGH SHOW, FOR N=10, 10, N.GE.100
SHOW PRINT FORMAT F10, M(N), M(N+1), M(N+2), M(N+3), M(N+4),
0 M(N+5), M(N+6), M(N+7), M(N+8), M(N+9), M(N+10) ZERO = 0 CROSS = 0 THROUGH ZC, FOR N=1, 1, N.G.1000 WHENEVER M(N).E.0, ZERO = ZERO + 1
ZC WHENEVER M(N).E.0 .AND. M(N-1).NE.0, CROSS = CROSS + 1
VECTOR VALUES FZ = $13HM(N) IS ZERO ,I2,S1,5HTIMES*$ PRINT FORMAT FZ, ZERO VECTOR VALUES FC = $18HM(N) CROSSES ZERO ,I2,S1,5HTIMES*$ PRINT FORMAT FC, CROSS END OF PROGRAM </lang>
- Output:
FIRST 99 MERTEN NUMBERS ARE 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) IS ZERO 92 TIMES M(N) CROSSES ZERO 59 TIMES
Pascal
Nearly the same as Square-free_integers#Pascal
Instead here marking all multiples, starting at factor 2, of a prime by incrementing the factor count.
runtime ~log(n)*n
<lang pascal>program Merten;
{$IFDEF FPC}
{$MODE DELPHI} {$Optimization ON,ALL}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF} uses
sysutils;
const
BigLimit = 10*1000*1000*1000;//1e10
type
tSieveElement = Int8; tpSieve = pInt8; tMoebVal = array[-1..1] of Int64;
var
MertensValues : array[-40000..50500] of NativeInt; primes : array of byte; sieve : array of tSieveElement;
procedure CompactPrimes; //searching for needed primes //last primes are marked with -1 var
pSieve : tpSieve; i,lmt,dp:NativeInt;
Begin
setlength(Primes,74500);//suffices for primes to calc square upto 1e12 //extract difference of primes i := 2; lmt := 0; dp := 2; pSieve :=@sieve[0]; repeat IF pSieve[i]= 0 then Begin //mark for Moebius pSieve[i]:= -1; primes[lmt] := dp; dp := 0; inc(lmt); end; inc(dp); inc(i); until i*i >BigLimit; setlength(Primes,lmt+1);
repeat IF pSieve[i]= 0 then //mark for Moebius pSieve[i]:= -1; inc(i); until i >BigLimit;
end;
procedure SieveSquares; //mark all powers >=2 of prime => all powers = 2 is sufficient var
pSieve : tpSieve; i,sq,k,prime : NativeInt;
Begin
pSieve := @sieve[0]; prime := 0; For i := 0 to High(primes) do Begin prime := prime+primes[i]; sq := prime*prime; k := sq; if sq > BigLimit then break; repeat pSieve[k] := 0; inc(k,sq); until k> BigLimit; end;
end;
procedure initPrimes; var
pSieve : tpSieve; fakt, sieveprime : NativeUint;
begin
pSieve := @sieve[0]; sieveprime := 2; repeat if pSieve[sieveprime]=0 then begin fakt := sieveprime+sieveprime; while fakt <=BigLimit do Begin //count divisors inc(pSieve[fakt]); inc(fakt,sieveprime); end; end; inc(sieveprime); until sieveprime>BigLimit DIV 2; //Möbius of 1 pSieve[1] := 1;
//convert to Moebius For fakt := 2 to BigLimit do Begin sieveprime := pSieve[fakt]; IF sieveprime<>0 then pSieve[fakt] := 1-(2*(sieveprime AND 1)) ; end; CompactPrimes; SieveSquares;
end;
procedure OutMerten10(Lmt,ZeroCross:NativeInt;Const MoebVal:tMoebVal); var
i,j: NativeInt;
Begin
Writeln(lmt:11,MoebVal[-1]:11,MoebVal[1]:11,MoebVal[-1]+MoebVal[1]:11, MoebVal[-1]-MoebVal[1]:7,MoebVal[0]:11); i:= low(MertensValues); while MertensValues[i] = 0 do inc(i); j:= High(MertensValues); while MertensValues[j] = 0 do dec(j); write('Merten min ',i:6,' max ',j:6,' zeros ',MertensValues[0]:8); writeln(' zeroCross ',ZeroCross); writeln;
end;
procedure Count_x10; var
MoebCount: tMoebVal; pSieve : tpSieve; i,lmt,Merten,Moebius,LastMert,ZeroCross: NativeInt;
begin
writeln('[1 to limit]'); Writeln('Limit Moeb. odd Moeb.even sqr-free Merten Zeros');
pSieve := @sieve[0]; For i := -1 to 1 do MoebCount[i]:=0; ZeroCross := 0; LastMert :=1; Merten :=0; lmt := 10; i := 1; repeat while i <= lmt do Begin Moebius := pSieve[i]; inc(MoebCount[Moebius]); inc(Merten,Moebius); inc(MertensValues[Merten]);//MoebCount[1]-MoebCount[-1]]); inc(ZeroCross,ORD( (Merten = 0) AND (LastMert <> 0))); LastMert := Merten; inc(i); end; OutMerten10(Lmt,ZeroCross,MoebCount);
IF lmt >= BigLimit then BREAK; lmt := lmt*10; IF lmt >BigLimit then lmt := BigLimit; until false; writeln;
end;
procedure OutMerten(lmt:NativeInt); var
i,k,m : NativeInt;
Begin
iF lmt> BigLimit then lmt := BigLimit; writeln('Mertens numbers from 1 to ',lmt); k := 9; write(:3); m := 0; For i := 1 to lmt do Begin inc(m,sieve[i]); write(m:3); dec(k); IF k = 0 then Begin writeln; k := 10; end; end; writeln;
end;
procedure OutMoebius(lmt:NativeInt); var
i,k : NativeInt;
Begin
iF lmt> BigLimit then lmt := BigLimit; writeln('Möbius numbers from 1 to ',lmt); k := 19; write(:3); For i := 1 to lmt do Begin write(sieve[i]:3); dec(k); IF k = 0 then Begin writeln; k := 20; end; end; writeln;
end;
Begin
setlength(sieve,BigLimit+1); InitPrimes; SieveSquares; Count_x10; OutMoebius(199); OutMerten(99); setlength(primes,0); setlength(sieve,0);
end.</lang>
- Output:
[1 to limit] Limit Moeb. odd Moeb.even sqr-free Merten Zero's 10 4 3 7 1 3 Merten min -2 max 1 zero's 1 zeroCross 1 100 30 31 61 -1 39 Merten min -4 max 2 zero's 6 zeroCross 5 1000 303 305 608 -2 392 Merten min -12 max 7 zero's 92 zeroCross 59 10000 3053 3030 6083 23 3917 Merten min -43 max 35 zero's 406 zeroCross 256 100000 30421 30373 60794 48 39206 Merten min -132 max 96 zero's 1549 zeroCross 949 1000000 303857 304069 607926 -212 392074 Merten min -368 max 311 zero's 5361 zeroCross 3269 10000000 3039127 3040164 6079291 -1037 3920709 Merten min -1078 max 1143 zero's 12546 zeroCross 7646 100000000 30395383 30397311 60792694 -1928 39207306 Merten min -3448 max 3290 zero's 41908 zeroCross 25525 1000000000 303963673 303963451 607927124 222 392072876 Merten min -8565 max 10246 zero's 141121 zeroCross 85652 10000000000 3039652332 3039618610 6079270942 33722 3920729058 Merten min -35517 max 50286 zero's 431822 zeroCross 262605 Möbius numbers from 1 to 199 1 -1 -1 0 -1 1 -1 0 0 1 -1 0 -1 1 1 0 -1 0 -1 0 1 1 -1 0 0 1 0 0 -1 -1 -1 0 1 1 1 0 -1 1 1 0 -1 -1 -1 0 0 1 -1 0 0 0 1 0 -1 0 1 0 1 1 -1 0 -1 1 0 0 1 -1 -1 0 1 -1 -1 0 -1 1 0 0 1 -1 -1 0 0 1 -1 0 1 1 1 0 -1 0 1 0 1 1 1 0 -1 0 0 0 -1 -1 -1 0 -1 1 -1 0 -1 -1 1 0 -1 -1 1 0 0 1 1 0 0 1 1 0 0 0 -1 0 1 -1 -1 0 1 1 0 0 -1 -1 -1 0 1 1 1 0 1 1 0 0 -1 0 -1 0 0 -1 1 0 -1 1 1 0 1 0 -1 0 -1 1 -1 0 0 -1 0 0 -1 -1 0 0 1 1 -1 0 -1 -1 1 0 1 -1 1 0 0 -1 -1 0 -1 1 -1 0 -1 0 -1 Mertens numbers from 1 to 99 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 real 3m54,249s = 234s //BigLimit = 100*1000*1000;takes 2.017s
Perl
<lang perl>use utf8; use strict; use warnings; use feature 'say'; use List::Util 'uniq';
sub prime_factors {
my ($n, $d, @factors) = (shift, 1); while ($n > 1 and $d++) { $n /= $d, push @factors, $d until $n % $d; } @factors
}
sub μ {
my @p = prime_factors(shift); @p == uniq(@p) ? 0 == @p%2 ? 1 : -1 : 0
}
sub progressive_sum {
my @sum = shift @_; push @sum, $sum[-1] + $_ for @_; @sum
}
my($upto, $show, @möebius) = (1000, 199, ()); push @möebius, μ($_) for 1..$upto; my @mertens = progressive_sum @möebius;
say "Mertens sequence - First $show terms:\n" .
(' 'x4 . sprintf "@{['%4d' x $show]}", @mertens[0..$show-1]) =~ s/((.){80})/$1\n/gr . sprintf("\nEquals zero %3d times between 1 and $upto", scalar grep { ! $_ } @mertens) . sprintf "\nCrosses zero%3d times between 1 and $upto", scalar grep { ! $mertens[$_-1] and $mertens[$_] } 1 .. @mertens;</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
Phix
Based on the stackexchange link, short and sweet but not very fast: 1.4s just for the first 1000... <lang Phix>function Mertens(integer n)
integer res = 1 for k=2 to n do res -= Mertens(floor(n/k)) end for return res
end function sequence s = {" ."} for i=1 to 143 do s = append(s,sprintf("%3d",Mertens(i))) end for puts(1,join_by(s,1,12," "))
integer prev = 1, zeroes = 0, crosses = 0 for n=2 to 1000 do
integer m = Mertens(n) if m=0 then zeroes += 1 crosses += prev!=0 end if prev = m
end for printf(1,"\nMertens[1..1000] equals zero %d times and crosses zero %d times\n",{zeroes,crosses})</lang>
- Output:
Matches the wp table:
. 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 Mertens[1..1000] equals zero 92 times and crosses zero 59 times
Prolog
<lang prolog>:- dynamic mertens_number_cache/2.
mertens_number(1, 1):- !. mertens_number(N, M):-
mertens_number_cache(N, M), !.
mertens_number(N, M):-
N >= 2, mertens_number(N, 2, M, 0), assertz(mertens_number_cache(N, M)).
mertens_number(N, N, M, M):- !. mertens_number(N, K, M, S):-
N1 is N // K, mertens_number(N1, M1), K1 is K + 1, S1 is S - M1, mertens_number(N, K1, M, S1).
print_mertens_numbers(Count):-
print_mertens_numbers(Count, 0).
print_mertens_numbers(Count, Count):-!. print_mertens_numbers(Count, N):-
(N == 0 -> write(' ') ; mertens_number(N, M), writef('%3r', [M]) ), N1 is N + 1, Column is N1 mod 20, (N > 0, Column == 0 -> nl ; true ), print_mertens_numbers(Count, N1).
count_zeros(From, To, Z, C):-
count_zeros(From, To, Z, C, 0, 0, 0).
count_zeros(From, To, Z, C, Z, C, _):-
From > To, !.
count_zeros(From, To, Z, C, Z1, C1, P):-
mertens_number(From, M), (M == 0 -> Z2 is Z1 + 1 ; Z2 = Z1), (M == 0, P \= 0 -> C2 is C1 + 1 ; C2 = C1), Next is From + 1, count_zeros(Next, To, Z, C, Z2, C2, M).
main:-
writeln('First 199 Mertens numbers:'), print_mertens_numbers(200), count_zeros(1, 1000, Z, C), writef('M(n) is zero %t times for 1 <= n <= 1000.\n', [Z]), writef('M(n) crosses zero %t times for 1 <= n <= 1000.\n', [C]).</lang>
- Output:
First 199 Mertens numbers: 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 M(n) is zero 92 times for 1 <= n <= 1000. M(n) crosses zero 59 times for 1 <= n <= 1000.
Raku
(formerly Perl 6)
Mertens number is not defined for n == 0. Raku 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) {
return 0 if n %% 4 or n %% 9 or n %% 25 or n %% 49 or n %% 121; my @p = prime-factors(n); +@p == +@p.unique ?? +@p %% 2 ?? 1 !! -1 !! 0
}
my @mertens = lazy [\+] flat , 1, (2..*).hyper.map: -> \n { μ(n) };
put "Mertens sequence - First 199 terms:\n",
@mertens[^200]».fmt('%3s').batch(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\n\nFirst Mertens equal to:";
for 10, 20, 30 … 100 -> $threshold {
printf "%4d: M(%d)\n", -$threshold, @mertens.first: * == -$threshold, :k; printf "%4d: M(%d)\n", $threshold, @mertens.first: * == $threshold, :k;
}</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 First Mertens equal to: -10: M(659) 10: M(1393) -20: M(2791) 20: M(3277) -30: M(9717) 30: M(8503) -40: M(9831) 40: M(11770) -50: M(24018) 50: M(19119) -60: M(24105) 60: M(31841) -70: M(24170) 70: M(31962) -80: M(42789) 80: M(48202) -90: M(59026) 90: M(48405) -100: M(59426) 100: M(114717)
REXX
The Mertens function is named after Franz Mertens.
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).
The above "feature" was added to make the grid to be aligned with other solutions. <lang rexx>/*REXX pgm computes & shows a value grid of the Mertens function for a range of integers*/ parse arg LO HI grp eqZ xZ . /*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 /* " " " " " " */ if eqZ== | eqZ=="," then eqZ= 1000 /* " " " " " " */ if xZ== | xZ=="," then xZ= 1000 /* " " " " " " */ call genP /*generate primes up to max √ HIHI */
call Franz LO, HI
if eqZ>0 then call Franz 1,-eqZ if xZ>0 then call Franz -1, xZ exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ Franz: parse arg a 1 oa,b 1 ob; @Mertens=' The Mertens sequence from ' a= abs(a); b= abs(b); grid= oa>=0 & ob>=0 /*semaphore used to show a grid title. */ if grid then say center(@Mertens LO " ──► " HI" ", max(50, grp*3), '═') /*show title*/
else say
zeros= 0 /*# of 0's found for Mertens function.*/ Xzero= 0 /*number of times that zero was crossed*/ prev= $= /*$ holds output grid of GRP numbers. */
do j=a to b; _= Mertens(j) /*process some numbers from LO ──► HI.*/ if _==0 then zeros= zeros + 1 /*Is Zero? Then bump the zeros counter*/ if _==0 then if prev\==0 then Xzero= Xzero+1 /*prev ¬=0? " " " Xzero " */ prev= _ if grid then $= $ right(_, 2) /*build grid if A & B are non─negative.*/ 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*/ if oa<0 then say @Mertens a " to " b ' has crossed zero ' Xzero " times." if ob<0 then say @Mertens a " to " b ' has ' zeros " zeros." return /*──────────────────────────────────────────────────────────────────────────────────────*/ Mertens: procedure expose @. !. M.; parse arg n /*obtain a integer to be tested for mu.*/
if M.n\==. then return M.n /*was computed before? Then return it.*/ 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*/ M.n= m; 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: !.=.; M.=!.; hihi= max(HI, eqZ, xZ) /*initialize 2 arrays for memoization. */
@.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6= 13; nP=6 /*assign low primes; # primes. */ do lim=nP until lim*lim>=hihi; end /*only keep primes up to the sqrt(HI). */ do j=@.nP+4 by 2 to hihi /*only find odd primes from here on. */ do k=3 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 The Mertens sequence from 1 to 1000 has 92 zeros. The Mertens sequence from 1 to 1000 has crossed zero 59 times.
Sidef
Built-in: <lang ruby>say mertens(123456789) #=> 1170 say mertens(1234567890) #=> 9163</lang>
Algorithm for computing M(n) in sublinear time:
<lang ruby>func mertens(n) is cached {
var lookup_size = (2 * n.iroot(3)**2) var mertens_lookup = [0]
for k in (1..lookup_size) { mertens_lookup[k] = (mertens_lookup[k-1] + k.moebius) }
static cache = Hash()
func (n) {
if (n <= lookup_size) { return mertens_lookup[n] }
if (cache.has(n)) { return cache{n} }
var M = 1 var s = n.isqrt
for k in (2 .. floor(n/(s+1))) { M -= __FUNC__(floor(n/k)) }
for k in (1..s) { M -= (mertens_lookup[k] * (floor(n/k) - floor(n/(k+1)))) }
cache{n} = M }(n)
}</lang>
Task: <lang ruby>with (200) {|n|
say "Mertens function in the range 1..#{n}:" (1..n).map { mertens(_) }.slices(20).each {|line| say line.map{ "%2s" % _ }.join(' ') }
}
with (1000) {|n|
say "\nIn the range 1..#{n}, there are:" say (1..n->count_by { mertens(_)==0 }, " zeros") say (1..n->count_by { mertens(_)==0 && mertens(_-1)!=0 }, " zero crossings")
}</lang>
- Output:
Mertens function in the range 1..200: 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 -8 In the range 1..1000, there are: 92 zeros 59 zero crossings
Wren
<lang ecmascript>import "/fmt" for Fmt import "/math" for Int
var isSquareFree = Fn.new { |n|
var i = 2 while (i * i <= n) { if (n%(i*i) == 0) return false i = (i > 2) ? i + 2 : i + 1 } return true
}
var mu = Fn.new { |n|
if (n < 1) Fiber.abort("Argument must be a positive integer") if (n == 1) return 1 var sqFree = isSquareFree.call(n) var factors = Int.primeFactors(n) if (sqFree && factors.count % 2 == 0) return 1 if (sqFree) return -1 return 0
}
var M = Fn.new { |x| (1..x).reduce { |sum, n| sum + mu.call(n) } }
System.print("The first 199 Mertens numbers are:") for (i in 0..9) {
for (j in 0..19) { if (i == 0 && j == 0) { System.write(" ") } else { System.write("%(Fmt.dm(3, M.call(i*20 + j))) ") } } System.print()
}
// use the recurrence relationship for the last 2 parts rather than calling M directly var count = 0 var mertens = M.call(1) for (i in 2..1000) {
mertens = mertens + mu.call(i) if (mertens == 0) count = count + 1
} System.print("\nThe Mertens function is zero %(count) times in the range [1, 1000].")
count = 0 var prev = M.call(1) for (i in 2..1000) {
var next = prev + mu.call(i) if (next == 0 && prev != 0) count = count + 1 prev = next
} System.print("\nThe Mertens function crosses zero %(count) times in the range [1, 1000].")</lang>
- Output:
The first 199 Mertens numbers are: 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 The Mertens function is zero 92 times in the range [1, 1000]. The Mertens function crosses zero 59 times in the range [1, 1000].
zkl
<lang zkl>fcn mertensW(n){
[1..].tweak(fcn(n,pm){ pm.incN(mobius(n)); pm.value }.fp1(Ref(0)))
} fcn mobius(n){
pf:=primeFactors(n); sq:=pf.filter1('wrap(f){ (n % (f*f))==0 }); // False if square free if(sq==False){ if(pf.len().isEven) 1 else -1 } else 0
} fcn primeFactors(n){ // Return a list of prime factors of n
acc:=fcn(n,k,acc,maxD){ // k is 2,3,5,7,9,... not optimum if(n==1 or k>maxD) acc.close(); else{
q,r:=n.divr(k); // divr-->(quotient,remainder) if(r==0) return(self.fcn(q,k,acc.write(k),q.toFloat().sqrt())); return(self.fcn(n,k+1+k.isOdd,acc,maxD)) # both are tail recursion
} }(n,2,Sink(List),n.toFloat().sqrt()); m:=acc.reduce('*,1); // mulitply factors if(n!=m) acc.append(n/m); // opps, missed last factor else acc;
}</lang> <lang zkl>mertensW().walk(199) .pump(Console.println, T(Void.Read,19,False), fcn{ vm.arglist.pump(String,"%3d".fmt) });
println("\nIn the first 1,000 terms of the Mertens sequence there are:"); otm:=mertensW().pump(1_000,List); otm.reduce(fcn(s,m){ s + (m==0) },0) : println(_," zeros"); otm.reduce(fcn(p,m,rs){ rs.incN(m==0 and p!=0); m }.fp2( s:=Ref(0) )); println(s.value," zero crossings");</lang>
- Output:
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