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

[Bart]

If you want it for all numbers from 1 to n, you need to use a sieve: https://github.com/benibela/bbutils/blob/master/bbutils.pas#L4725-L4777

Nice one, but: in order to calculate Phi(30030) you have to calculate Pji(1)..Phi(30030).
Mine just factorizes it into 2, 3, 5, 7, 11, 13 and then does (13+1) fractional multiplications.

But very wel suited to create a lookup table of phi(1 to N)!
Great job.

Bart

[rabbit_dance]

i found my program could not be executed properly on windows platform.

[MathMan]

You can take a look here - http://stackoverflow.com/questions/1024640/calculating-phik-for-1kn

But you need to translate the source given (Visual Basic) to Pascal yourself (shouldn't be too hard though). In general the parallel computation of Eulers totient function for 1-n can be made much faster than n individual computations.

Regards,
MathMan

[rabbit_dance]

sorry, i cant read it.

[Eugene Loza]

why so complex? because i want to design a most efficiency program.

Complex is not always efficient.
1. Most efficient program is not the one that calculates something in a most efficient way. I.e. you have some task set, and this one is a sub-task. If you have to calculate just a few Phis then no need to worry of efficiency (you've spent much more time on struggling with bugs already). If your program calculates (or will calculate) something big based on Phis then maybe a pre-calculated data file would do a much more efficient job (and again, you won't really care about efficiency of its generation).
2. You can simply 'compare' your results to correct ones you've got by an inefficient (and robust) method to test for its accuracy/correctness.
3. Moreover, you can 'compare' efficiency of "efficient" and "inefficient" methods ... to find out that "inefficient" method works more efficiently than an "efficient" one, e.g. by simple threading. I've happened a few such problems.
4. Efficient method might be unstable (i.e. provide wrong results under some conditions). E.g. if you solve an equation by direct solution and forget to check whether it has three distinctive real roots, you may get weird results (a real-life example wasted half a year of researchers time).

[MathMan]

There's another way to calculate Phi: https://en.wikipedia.org/wiki/Euler's_totient_function#Computing_Euler.27s_totient_function
But I haven't tested its efficiency, it might be even less efficient (for N<1000000) due to multiple calculation of prime numbers. And definitely the complex code would need a lot of debugging.

It would also give rounding errors:
93296 consists of the following primes: 1, 2, 7, 17
The Result would be 93296 * (1-1/2) * (1-1/7) * (1-1/17) = 37632, but calculating this with floating point numbers will give a wrong result.
Calculating it using Fractions (using integer multipliction, greates common divisor etc) may overflow very easily (e.g. when a number has some primes > 1000).

The multiple computing of prime numbers can be overcome by having a static table of primes (up til Sqrt(High(UInt64))?).

It would be nice to build such a thing though.
The fractions part is easy, I already have a fractions unit.

Bart

Which one wouldn't do of course, as 93296*( 1-1/2 )*( 1-1/7 )*( 1-1/17 ) = 93296*( 2-1 )*( 7-1 )*( 17-1 )/( 2*7*17 ) = 93296/( 2*7*17 )*( 1*6*16 ). The first Division is an exact one - meaning the remainder=0. So this can be done fast, without fractions, in the number magnitude of the requested number (a 32 bit unsigned int in case of 93296 ).

[Bart]

Which one wouldn't do of course, as 93296*( 1-1/2 )*( 1-1/7 )*( 1-1/17 ) = 93296*( 2-1 )*( 7-1 )*( 17-1 )/( 2*7*17 ) = 93296/( 2*7*17 )*( 1*6*16 ). The first Division is an exact one - meaning the remainder=0. So this can be done fast, without fractions, in the number magnitude of the requested number (a 32 bit unsigned int in case of 93296 ).

Yes, you can factor out all this and Yes in the end it will always be an exact divison (by definition, if not, something went wrong in the frist place), but give the prime factors can you do that in Pascal in such a way that no floating point arithmetic is involved?
For prime factors p1,p2,..pn is it alway the case that the denominator is (p1*p2*..*pn)((p1-1)*(p2-1)*...(pn-1)) ?

(I'm just too lazy to work that out by mself, plus I'm at work and not supposed to write here at this point in time.)

Bart

[MathMan]

Yes, you can factor out all this and Yes in the end it will always be an exact divison (by definition, if not, something went wrong in the frist place), but give the prime factors can you do that in Pascal in such a way that no floating point arithmetic is involved?
For prime factors p1,p2,..pn is it alway the case that the denominator is (p1*p2*..*pn)((p1-1)*(p2-1)*...(pn-1)) ?

(I'm just too lazy to work that out by mself, plus I'm at work and not supposed to write here at this point in time.)

Bart

In fact it can - look at the link I provided to rabbit_dance (though as said this is Visual Basic). The calculation shown there is good for computing Phi(1..k) to k=2^63-1 using QWORD.

Cheers,
MathMan

[Bart]

Thought about it.

(1-(1/p)) = (p-1)/p, so the nominator will have (p1-1)*(p2-1)*..*(pn-1) and the denominator the products of the primes p1..pn.

The product of p1*p2*..*pn will always be less than or equal to the original number (N) we factorized.
And this product must be a factor of the original number (N mod (p1*p2*..*pn) = 0).

So first divide N by the product of the primes.
Then mutilply by (p1-1)*(p2-1)*..*(pn-1), and this is guaranteed to be less than N, so the math will be safe from overflowing.

Hope I didn't mess up my maths here.

Bart
(Not a real MathMan, but fond of numbers.)

[MathMan]

Thought about it.

(1-(1/p)) = (p-1)/p, so the nominator will have (p1-1)*(p2-1)*..*(pn-1) and the denominator the products of the primes p1..pn.

The product of p1*p2*..*pn will always be less than or equal to the original number (N) we factorized.
And this product must be a factor of the original number (N mod (p1*p2*..*pn) = 0).

So first divide N by the product of the primes.
Then mutilply by (p1-1)*(p2-1)*..*(pn-1), and this is guaranteed to be less than N, so the math will be safe from overflowing.

Hope I didn't mess up my maths here.

Bart
(Not a real MathMan, but fond of numbers.)

You thought correct 

[rabbit_dance]

Crash!

Code: [Select]

C:\Users\Bart\LazarusProjecten\ConsoleProjecten\bugs\euler>euler
10
  2 0  3 1  4 2  5 2  6 4  7 2  8 6  9 4  10 6  11 4Marked memory at $001A5C20 i
nvalid
Wrong size : 8 allocated 4 freed
  $00409EB9
  $00404C86
  $004014EC  PRINTING,  line 16 of euler.lpr
  $004014EC  PRINTING,  line 16 of euler.lpr
  $004014EC  PRINTING,  line 16 of euler.lpr
  $004014EC  PRINTING,  line 16 of euler.lpr
  $004014EC  PRINTING,  line 16 of euler.lpr
  $004014EC  PRINTING,  line 16 of euler.lpr
  $004014EC  PRINTING,  line 16 of euler.lpr
Call trace for block $001A5C20 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Heap dump by heaptrc unit
13 memory blocks allocated : 307/320
2 memory blocks freed     : 223/240
11 unfreed memory blocks : 84
True heap size : 131072 (128 used in System startup)
True free heap : 130064
Should be : 130160
Call trace for block $001A5C70 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5C20 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5BD0 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5B80 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5B30 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5AE0 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5A90 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5A40 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A59F0 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A59A0 size 8
  $00401570  ASSIGNMENT,  line 33 of euler.lpr
  $00401621  main,  line 50 of euler.lpr
Call trace for block $001A5950 size 8

Bart

could you tell what causes the above error appeared?

[garlar27]

could you tell what causes the above error appeared?

The errors marked are unfreed memory block detected by HeapTrace.

One on line 16: Wrong size : 8 allocated 4 freed

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

1. IF P^.NEXT=NIL THEN EXIT();
2. PRINTING(P); //<===== LINE 16
3. FREEMEM(P,SIZEOF(LINKED_PTR));
4.  

though it states "Line 16" maybe the problem is the following line. I never used Freemem but it seems you are not freeing the right size. My guess is that you are creating a memory space of 8 bytes for a "LINK_NODE" and you are freeing a memory space of 4 Bytes  the size of "LINKED_PTR".

The other on line 33:

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

1. PROCEDURE ASSIGNMENT(PTR:LINKED_PTR;N:INTEGER);
2. VAR K:INTEGER;
3. BEGIN
4.         FOR K:=1 TO N DO
5.         BEGIN
6.                 SUCCESSOR(PTR); //<====== Line 33
7.         END;
8. END;
9.  

Take a look to what happens to "PTR". Without checking the logic I can tell that there might be a problem with how parameters are passed: "procedure ASSIGNMENT" receives the parameter "PTR" by value then you call "SUCCESOR" which receives the parameter "PTR" by VARIABLE and set a new value for ASSIGNMENT's PTR parameter whic value is lost after leaving the procedure.
Check what you wanted to do and what you are actually doing.

Note that if you change ASSIGNMENT's function declaration to:

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

1. PROCEDURE ASSIGNMENT(const PTR:LINKED_PTR;N:INTEGER);

it will no longer compile.

EDIT: if you are using typed pointers then it is better to use "New(PTR)" to create the pointer and space in memory and use "Dispose(PTR)" to release that memory without errors in the memory size.

[Bart]

You thought correct 

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

Code: [Select]

function EulerPhi(N: UInt64): UInt64;
var
  PrimeFactors: IntArray;
  Len, i: Integer;
  Numerator, Denominator: UInt64;
begin
  if (N = 0) then
    raise ERangeError.Create('Phi(0) is undefined.');
  if (N = 1) then
    Exit(1);
  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]);
  if (Len = 1) and (PrimeFactors[0] = N) then //N is a prime, all numbers below it are relative prime to it
    Exit(N-1);
  Result := 0;
  Numerator := 1;
  Denominator := 1;
  // 1 - (1/p) = (p-1) / p
  for i := 0 to Len - 1 do
  begin
    Denominator := Denominator * PrimeFactors[i];
    Numerator := Numerator * (PrimeFactors[i] - 1);
    //writeln('Numerator = ',Numerator, ', Denominator = ',Denominator);
  end;
  Result := N div Denominator;
  //Denominator is the product of all prime PrimeFactors of N, so N must be divisible by it
  if ((N mod Denominator) <> 0) then
    raise Exception.CreateFmt('EulerProduct: Fatal error: cannot resolve %d/%d to an integer',[N, Denominator]);
  Result := Result * Numerator;
end;

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

Bart

[rabbit_dance]

i found the results of execution of the program compiled of the revised source code on linux and windows are different.

[rabbit_dance]

i found the results of execution of the program compiled of the revised source code on linux and windows are different.

FREEMEM(P,SIZEOF(LINKED_PTR));==>FREEMEM(P,SIZEOF(LINK_NODE));
THIS IS THE CHANGE.