[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]

(*
                Daniel Doubrovkine - dblock@infomaniak.ch
   part of the Expression Calculator 2.0 Multithread, destribute freely
              http://www.infomaniak.ch/~dblock/express.htm
    (ref. Haussman, Simple Algebra course at the University of Geneva)

                            "Euler's Indice"
       euler's indice is a insomorphe function phi: Zm->Z1*Z2...Zn
    means the decomposition of a number into a primes product is unique
    with phi(n):=n*(1-1/p1)*...*(1-1/pn) where pi are primes taken once
    from the decomposition.
    ex: (following algorith steps):
       168=2*84
       168=2*2*42
       168=2*2*2*21
       168=2*2*2*3*7 where all are primes!
  now, phi(168)=168*(1-1/2)(1-1/3)(1-1/7)=48 (this is always an integer!)
       thus, 168 has 48 primes to it.
*)

(*defining a dynamic primes array, so this function will never be limited*)
type
   PPrimesArray=^TPrimesArray;
   TPrimesArray=record
      Prime:longint;
      NextPrime:PPrimesArray;
      end;


   (*adding a prime to the array*)
   procedure PrimesAdd(var AllPrimes: PPrimesArray;Prime:longint);
   var
      TempPrimes:PPrimesArray;
   begin
      TempPrimes:=AllPrimes;
      while (TempPrimes<>nil) do begin
         if TempPrimes^.Prime=Prime then exit;
         TempPrimes:=TempPrimes^.NextPrime;
         end;
      new(TempPrimes);
      TempPrimes^.Prime:=Prime;
      TempPrimes^.NextPrime:=AllPrimes;
      AllPrimes:=TempPrimes;
      end;

   (*removing the primes array*)
   procedure PrimesRemove(AllPrimes: PPrimesArray);
   var
      TempPrimes:PPrimesArray;
   begin
      while AllPrimes<>nil do begin
         TempPrimes:=AllPrimes;
         AllPrimes:=AllPrimes^.NextPrime;
         Dispose(TempPrimes);
         end;
      end;

   (*by the way, this routine uses some code from swag*)
   function GeneratePrimes(var AllPrimes: PPrimesArray;n: longint):longint;
    procedure CalculateTPrime(n: longint);
    var
      TempPrime: PPrimesArray;
      CurrentPrime: LongInt;
      prime:boolean;
      number,max_div,divisor,lastprime:longint;
    begin
      CurrentPrime:=AllPrimes^.Prime;
      for number:=CurrentPrime to MaxLongInt do begin
         max_div:=Round(Sqrt(number)+0.5);
         prime:=number mod 2 <> 0;
         divisor:=3;
         while prime and (divisor<max_div) do begin
               prime:=number mod divisor <>0;
               divisor:=divisor+2;
               end;
         if prime then begin
         if AllPrimes^.Prime<CurrentPrime then PrimesAdd(AllPrimes,CurrentPrime);
            CurrentPrime:=number;
         if n<CurrentPrime then exit;
         end;
         end;
         end;

   var
      TempPrime: PPrimesArray;
      CurrentPrime: LongInt;
      prime:boolean;
      number,max_div,divisor,lastprime:longint;
   begin
     PrimesAdd(AllPrimes,1);
     PrimesAdd(AllPrimes,2);
      if (n>-1) and (n<1) then begin
         writeln('Prime requires values |x|>1');
         halt;
         end;
      if n>maxlongint then begin
         writeln('prime limit too large');
         halt;
         end;

      n:=abs(n);

      if (n<=AllPrimes^.Prime) then begin
      TempPrime:=AllPrimes;
         while TempPrime<>nil do begin
               CurrentPrime:=TempPrime^.Prime;
               if n>=TempPrime^.Prime then begin
                  GeneratePrimes:=TempPrime^.Prime;
                  exit;
                  end;
           TempPrime:=TempPrime^.NextPrime;
           end;
           end;
      CalculateTPrime(n);
      GeneratePrimes:=AllPrimes^.Prime;
      end;


function phi(var AllPrimes:PPrimesArray;Value:longint):longint;
var
   UsedPrimes:PPrimesArray;
   (*this is partailly from swag...*)
   procedure Factoring(lin:longint);
   var
      lcnt:longint;
   begin
      lcnt:=2;
      if GeneratePrimes(AllPrimes,lin)=lin then begin
         write(' ',lin);
         PrimesAdd(UsedPrimes,lin);
         end else
         while(lcnt*lcnt<=lin) do begin
            if (lin mod lcnt) = 0 then begin
               if GeneratePrimes(AllPrimes,lcnt)<>lcnt then factoring(lcnt)
                  else begin
                     write(' ',lcnt);
                     primesadd(UsedPrimes,lcnt);
                     end;
               if GeneratePrimes(AllPrimes,lin div lcnt)<>(lin div lcnt) then factoring(lin div lcnt)
                  else begin
                     write(' ',lin div lcnt);
                     primesadd(UsedPrimes,lin div lcnt);
                     end;
               exit;
               end;
         lcnt:=lcnt+1;
         end;
      end;

   var
      FinalResult:real;
      TempPrimes:PPrimesArray;
   begin
     if Value=0 then begin
        Phi:=0;
        exit;
        end;
     Value:=Abs(Value);
     FinalResult:=Value;
     UsedPrimes:=nil;
     write('Decomposition of ',Value:0,':');
     Factoring(Value);
     writeln;
     TempPrimes:=UsedPrimes;
     while TempPrimes<>nil do begin
        FinalResult:=FinalResult*(1-1/TempPrimes^.Prime);
        TempPrimes:=TempPrimes^.NextPrime;
        end;
     phi:=trunc(FinalResult);
     PrimesRemove(UsedPrimes);
     writeln('Euler''s indice for ',Value:0,' is ',trunc(FinalResult));
     writeln('(there are ',trunc(FinalResult),' primes to ',Value,')');
   end;

var
   AllPrimes: PPrimesArray;
begin
     Phi(AllPrimes,168);
     PrimesRemove(AllPrimes);
   end.

[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]