[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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)