how to calculate euler phi function of integer from 1 to n effectively?

[MathMan]

OK, so without fractions (and a second shortcut for prime numbers):

Code: [Select]

function EulerPhi(N: UInt64): UInt64;

  ...

  PrimeFactors := Factorize(N);

  //Calculate N * the product of (1 - (1/p)), where p is a prime number listed in PrimeFactors
  Len := Length(PrimeFactors);
  if (Len = 0) then
    raise ERangeError.CreateFmt('EulerPhi: Fatal error: %d has no prime factors???.',[N]);

  ...


I suspect this code is more efficient than the code from rabbit_dance.

Bart

Regarding the function - I don't know what's behind "Factorize()" but would assume that it will not include the value "1" in the list. If that than the exception you raise for an empty list is probably incorrect, as N is a prime then and Phi( p ) := p-1 iff p prime by definition.

Regarding efficiency - hard to say if the above is more efficient than the function from rabbit_dance as nearly all time will be spend in "Factorize()". Depending on the size of N different factorization algorithms are the most efficient. Assuming you'd always select the best fitting one your code is optimal with respect to calculation of a single number. However the question was "calculate Phi for 1..N" and this can be speeded up significantly if one computes the factorization for the numbers 1..N in parallel.

Cheers,
MathMan

[Bart]

Regarding the function - I don't know what's behind "Factorize()" but would assume that it will not include the value "1" in the list.

Correct, since 1 is not a prime.

If that than the exception you raise for an empty list is probably incorrect, as N is a prime then and Phi( p ) := p-1 iff p prime by definition.

No, if there are no elements in it (which b.t.w. should never happen,I'm just insanely paranoid about my math skills), N does not have prime factors.
If N is a prime there will be 1 element, namely N itself.

Regarding efficiency - hard to say if the above is more efficient than the function from rabbit_dance as nearly all time will be spend in "Factorize()". Depending on the size of N different factorization algorithms are the most efficient. Assuming you'd always select the best fitting one your code is optimal with respect to calculation of a single number. However the question was "calculate Phi for 1..N" and this can be speeded up significantly if one computes the factorization for the numbers 1..N in parallel.

Without testing it will be impossible to say (and rabbit_dance's code does not run properly).
My Factorize() function uses a lookputable for primes (at least for the first 10000 primes), but the algoritme could probably be more optimized.
E.g. it uses a series of div and mod, where a DivMod() would probably more efficient, but this is not available for UInt64 types.

Anyway, for me the fun is not in who writes the fastest, or most efficient, of better scale-able code, the fun is in understanding and figuring out the algorithm and in the programming itself.

Chances are I will never ever need this in my real live programs.

Thanks for your input/feedback.

Bart

[MathMan]

Anyway, for me the fun is not in who writes the fastest, or most efficient, of better scale-able code, the fun is in understanding and figuring out the algorithm and in the programming itself.

Chances are I will never ever need this in my real live programs.

Thanks for your input/feedback.

Bart

Yepp - that's exactly why I am implementing an arbitrary precision library as a pet project ...

If you should have questions regarding mathematical algorithms (or closely related stuff) feel free to ask via private mail.

Cheers,
MathMan

[rabbit_dance]

what is maximal integer that the program could calculate?
how many factor affect this?
I have read computer theoretical book before, it mentioned the depth of recursion is limited.
how to change the depth?

[Leledumbo]

what is maximal integer that the program could calculate?
how many factor affect this?

Natively: Min(maximum bitness the computer can compute, compiler limits). On x86[_64] FPC, it's 18446744073709551615 (High(QWord)). I don't know whether the compiler has 128-bit integer support or not.
Unnatively: Arbitrary, only computer memory serves as the limit.

I have read computer theoretical book before, it mentioned the depth of recursion is limited.
how to change the depth?

Give bigger stack or reduce parameter size. Recursion limit is the same as stack limit, it's controlled by CPU but can also be overriden by OS. Language implementations (i.e.: compilers and interpreters) usually doesn't limit itself, but at least with FPC you can control the limit with -Cs. Note however, stack overflow is a fatal error, one shouldn't let the program continues after stack overflow occurs.

[rabbit_dance]

what is maximal integer that the program could calculate?
how many factor affect this?

Natively: Min(maximum bitness the computer can compute, compiler limits). On x86[_64] FPC, it's 18446744073709551615 (High(QWord)). I don't know whether the compiler has 128-bit integer support or not.
Unnatively: Arbitrary, only computer memory serves as the limit.

I have read computer theoretical book before, it mentioned the depth of recursion is limited.
how to change the depth?

Give bigger stack or reduce parameter size. Recursion limit is the same as stack limit, it's controlled by CPU but can also be overriden by OS. Language implementations (i.e.: compilers and interpreters) usually doesn't limit itself, but at least with FPC you can control the limit with -Cs. Note however, stack overflow is a fatal error, one shouldn't let the program continues after stack overflow occurs.

i have ever heard high precision arithmetic can be implemented by programming, but i want to know is there any implementation, and how to use it.

[Leledumbo]

i have ever heard high precision arithmetic can be implemented by programming, but i want to know is there any implementation, and how to use it.

A lot:

• http://rosettacode.org/wiki/Arbitrary-precision_integers_(included)#Pascal
• http://www.delphiforfun.org/programs/Library/big_integers.htm
• https://bitbucket.org/sergworks/tforge/downloads (use the numerics package)

[Thaddy]

Don't forget NUMLIB. It comes with fpc as standard and is arbitrary precision as well.
Docs are cumbersome (to me).

It is in packages\numlib

[BeniBela]

Or my bigdecimalmath, although there is no reason to use it, if all you need are integers

[rabbit_dance]

PROGRAM EULER_PHI;
TYPE
   LIST=^PTR;
   PTR=RECORD
   NUMBER:INTEGER;
   NEXT:LIST;
END;
VAR P,P1,P2:LIST; VAR N,I,J:INTEGER;
PROCEDURE OUTPUT(P:LIST);
VAR K:LIST;
BEGIN
   WHILE(NOT (P=NIL)) DO
   BEGIN
      WRITE(P^.NUMBER,' ');
      K:=P^.NEXT;
      FREEMEM(P,SIZEOF(PTR));
      P:=K;
   END;
END;
PROCEDURE EXTRACT_ELEMENT(N:INTEGER;VAR P:LIST);
VAR I:INTEGER;
BEGIN
        I:=1;
        WHILE (I<=N) DO
        BEGIN
      P:=P^.NEXT;
                I:=I+1;
        END;

END;
PROCEDURE INIT_LIST(N:INTEGER;VAR P:LIST);
VAR I:INTEGER;VAR K:LIST;
BEGIN
   P:=ALLOCMEM(SIZEOF(PTR));
   K:=P;I:=1;
   WHILE (I<N) DO
   BEGIN
      K^.NEXT:=ALLOCMEM(SIZEOF(PTR));
      K:=K^.NEXT;
      I:=I+1;
   END;
END;
BEGIN
   READLN(N);
   INIT_LIST(N,P);
   P^.NUMBER:=1;
   I:=2;
        P1:=P;
   WHILE(I<=N) DO
   BEGIN
           EXTRACT_ELEMENT(1,P1);
      IF P1^.NUMBER=0 THEN
      BEGIN
         J:=I;
         P2:=P1;
         IF P2^.NUMBER=0 THEN P2^.NUMBER:=J;
         P2^.NUMBER:=(P2^.NUMBER DIV I)*(I-1);
                        J:=J+I;
         WHILE(J<=N) DO
         BEGIN
            EXTRACT_ELEMENT(I,P2);
            IF P2^.NUMBER=0 THEN P2^.NUMBER:=J;
            P2^.NUMBER:=(P2^.NUMBER DIV I)*(I-1);
                                J:=J+I;
                END;
      END;
      I:=I+1;
   END;
   OUTPUT(P);
END.

[rabbit_dance]

Code: Pascal  [Select][+][-]

1. PROGRAM EULER_PHI;
2. TYPE
3.         LIST=^PTR;
4.         PTR=RECORD
5.         NUMBER:INTEGER;
6.         NEXT:LIST;
7. END;
8. VAR P,P1,P2:LIST; VAR N,I,J:INTEGER;
9. PROCEDURE OUTPUT(P:LIST);
10. VAR K:LIST;
11. BEGIN
12.         WHILE(NOT (P=NIL)) DO
13.         BEGIN
14.                 WRITE(P^.NUMBER,' ');
15.                 K:=P^.NEXT;
16.                 FREEMEM(P,SIZEOF(PTR));
17.                 P:=K;
18.         END;
19. END;
20. PROCEDURE EXTRACT_ELEMENT(N:INTEGER;VAR P:LIST);
21. VAR I:INTEGER;
22. BEGIN
23.         I:=1;
24.         WHILE (I<=N) DO
25.         BEGIN
26.                 P:=P^.NEXT;
27.                 I:=I+1;
28.         END;
29.  
30. END;
31. PROCEDURE INIT_LIST(N:INTEGER;VAR P:LIST);
32. VAR I:INTEGER;VAR K:LIST;
33. BEGIN
34.         P:=ALLOCMEM(SIZEOF(PTR));
35.         K:=P;I:=1;
36.         WHILE (I<N) DO
37.         BEGIN
38.                 K^.NEXT:=ALLOCMEM(SIZEOF(PTR));
39.                 K:=K^.NEXT;
40.                 I:=I+1;
41.         END;
42. END;
43. BEGIN
44.         READLN(N);
45.         INIT_LIST(N,P);
46.         P^.NUMBER:=1;
47.         I:=2;
48.         P1:=P;
49.         WHILE(I<=N) DO
50.         BEGIN
51.                 EXTRACT_ELEMENT(1,P1);
52.                 IF P1^.NUMBER=0 THEN
53.                 BEGIN
54.                         J:=I;
55.                         P2:=P1;
56.                         IF P2^.NUMBER=0 THEN P2^.NUMBER:=J;
57.                         P2^.NUMBER:=(P2^.NUMBER DIV I)*(I-1);
58.                         J:=J+I;
59.                         WHILE(J<=N) DO
60.                         BEGIN
61.                                 EXTRACT_ELEMENT(I,P2);
62.                                 IF P2^.NUMBER=0 THEN P2^.NUMBER:=J;
63.                                 P2^.NUMBER:=(P2^.NUMBER DIV I)*(I-1);
64.                                 J:=J+I;
65.                         END;
66.                 END;
67.                 I:=I+1;
68.         END;
69.         OUTPUT(P);
70. END.
71.  

[rabbit_dance]

Code: Pascal  [Select][+][-]

1. PROGRAM EULER_PHI;
2. TYPE
3.         LIST=^PTR;
4.         PTR=RECORD
5.                 NUMBER:INTEGER;
6.                 NEXT:LIST;
7. END;
8. VAR P,P1,P2:LIST; VAR N,I,J:INTEGER;
9. PROCEDURE OUTPUT(P:LIST);
10. VAR K:LIST;
11. BEGIN
12.         WHILE(NOT (P=NIL)) DO
13.         BEGIN
14.                 WRITE(P^.NUMBER,' ');
15.                 K:=P^.NEXT;
16.                 FREEMEM(P,SIZEOF(PTR));
17.                 P:=K;
18.         END;
19. END;
20. PROCEDURE EXTRACT_ELEMENT(N:INTEGER;VAR P:LIST);
21. VAR I:INTEGER;
22. BEGIN
23.         I:=1;
24.         WHILE (I<=N) DO
25.         BEGIN
26.                 P:=P^.NEXT;
27.                 I:=I+1;
28.         END;
29.  
30. END;
31. PROCEDURE INIT_LIST(N:INTEGER;VAR P:LIST);
32. VAR I:INTEGER;VAR K:LIST;
33. BEGIN
34.         P:=ALLOCMEM(SIZEOF(PTR));
35.         K:=P;I:=1;
36.         WHILE (I<N) DO
37.         BEGIN
38.                 K^.NEXT:=ALLOCMEM(SIZEOF(PTR));
39.                 K:=K^.NEXT;
40.                 I:=I+1;
41.         END;
42. END;
43. PROCEDURE COMPUTE(N:INTEGER;VAR P:LIST);
44. VAR P1,P2:LIST;VAR I,J:INTEGER;
45. BEGIN
46.         INIT_LIST(N,P);
47.         P^.NUMBER:=1;
48.         I:=2;
49.         P1:=P;
50.         WHILE(I<=N) DO
51.         BEGIN
52.                 EXTRACT_ELEMENT(1,P1);
53.                 IF P1^.NUMBER=0 THEN
54.                 BEGIN
55.                         P2:=P1;
56.                         IF P2^.NUMBER=0 THEN P2^.NUMBER:=I-1;
57.                         J:=2*I;
58.                         WHILE(J<=N) DO
59.                         BEGIN
60.                                 EXTRACT_ELEMENT(I,P2);
61.                                 IF P2^.NUMBER=0 THEN P2^.NUMBER:=J;
62.                                 P2^.NUMBER:=(P2^.NUMBER DIV I)*(I-1);
63.                                 J:=J+I;
64.                         END;
65.                 END;
66.                 I:=I+1;
67.         END;
68. END;
69. BEGIN
70.         READLN(N);
71.         COMPUTE(N,P);
72.         OUTPUT(P);
73.         READLN();
74. END.
75.