Post Reply 
[VA] SRC #013 - Pi Day 2023 Special
03-18-2023, 10:21 AM
Post: #6
RE: [VA] SRC #013 - Pi Day 2023 Special
My initial attempt was, instead of counting the square-free numbers, to count the numbers that have a square prime factor, and to subtract this amount from N.
The count of numbers <=N having the 2^2 factor is IP(N/2^2), those having the factor 3^2 is IP(N/3^2), and so on,
so the count of square-free numbers should be S = N - IP(N/2^2) - IP(N/3^2) - IP(N/5^2) ...
In that way, I only had to explore primes up to SQR(N), so it was pretty fast.
Unfortunately, numbers that are multiple of several square primes such as (2*3)^2 are counted several times, and the result is underestimated:

20 N=12345 ! test case
30 S=N @ P=1
40 P=FPRIM(P+1) @ S=S-IP(N/P^2) @ IF P<SQR(N) THEN 40
50 DISP N;S

>RUN
12345 6793 (true result=7503)

I didn't find a way to easily manage the numbers with multiple square prime factors.
An improvement was to use the current sum S instead of N: S = N ; S=S-IP(S/2^2); S=S-IP(S/3^2) ... but it's still an approximation at most:

40 P=FPRIM(P+1) @ S=S-IP(S/P^2) @ IF P<SQR(N) THEN 40

>RUN
12345 7529 (true result=7503)


So, I resorted to the formula using the Möbius function as disclosed in the post #2 above.
Matter of fact, this formula starts with the same terms as my first attempt:
S = N - IP(N/2^2) - IP(N/3^2) - IP(N/5^2) ...
but "magically" manages the numbers with multiple square prime factors with terms such as +IP(N/(2*3)^2) !

Here is my implementation, (correct) results and timings on my Emu71/DOSBox, at about 150x speed:

10 ! SRC13
20 INPUT N
25 T=TIME
30 S=N @ FOR I=2 TO SQR(N) @ S=S+FNM(I)*IP(N/I^2) @ NEXT I
40 T=TIME-T
50 DISP N;S;T
80 !
90 DEF FNM(N) ! Moebius function
110 C=-1 @ Q=1
120 P=PRIM(N) @ IF P=0 THEN P=N
130 IF P=Q THEN FNM=0 @ END
140 IF P<N THEN C=-C @ N=N/P @ Q=P @ GOTO 120
150 FNM=C @ END DEF

>RUN
12345 7503 .23
10000 60794 .65
567890 345237 1.37
1000000 607926 1.8

However, I'm not fully satisfied by applying a formula without understanding how and why it works.

Now, I'm curious to read Valentin's solution and explanations !

J-F
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: [VA] SRC #013 - Pi Day 2023 Special - J-F Garnier - 03-18-2023 10:21 AM



User(s) browsing this thread: 1 Guest(s)