{PROGRAMS FROM "PRIMES AND PROGRAMMING" BY P.J.GIBLIN,
PUBLISHED BY CAMBRIDGE UNIVERSITY PRESS.
NOTE THAT THE DRIVE AND DIRECTORY MAY NEED TO BE SPECIFIED
WHEN USING THE FILE PRIMES.DAT


PLEASE REPORT ANY ERRORS TO pjgiblin@liv.ac.uk}



PROGRAM P1_3_5;
 VAR
   a,b,r:extended;
 BEGIN
   writeln ('type in 2 positive numbers separated by a space ');
   read (a,b);
   write ('GCD ',a:0:0,' and ',b:0:0, ' is ');
   WHILE b > 0 DO
      BEGIN
        r:=a - b *  INT(a/b);
        a:=b;
        b:=r;
      END;
   writeln (a:0:0);
 END.

PROGRAM P1_3_16;     {Finds linear combination of a and b such that
                 as + bt = gcd(a,b)}
 VAR
   a,b,a0,b0,q,r,s,s1,s2,t,t1,t2:extended;
 BEGIN
   writeln ('type in 2 numbers separated by a space ');
   read (a,b);
   a0:=a;
   b0:=b;
   s2:=1;
   s1:=0;
   t2:=0;
   t1:=1;
   WHILE b > 0 DO
      BEGIN
        q:=INT(a/b);
        IF (q > a/b) THEN q:=q-1;
        r:=a - b * q;
        s:=s2 - q * s1;
        t:=t2 - q * t1;
        a:=b;
        b:=r;
        s2:=s1;  s1:=s;
        t2:=t1;  t1:=t;
              END;
   writeln ('gcd of   ',a0:0:0,' and   ',b0:0:0,' is ',   a:0:0);
   writeln ('linear combination is ',a0:0:0,'*',s2:0:0,'+',b0:0:0,'*',
                             t2:0:0,'=',a:0:0);

 END.


PROGRAM P2_1_3;

  VAR
    n,q,max:extended;
    qdividesn:boolean;

  BEGIN
    writeln;
    write(' Max =  ? ');
    readln(max);
    n:=5;
    writeln('prime numbers are  ');
    write('2,3');
    WHILE n<=max do
        BEGIN
            q:=3;
                 REPEAT
                       BEGIN
                            qdividesn:=(n=q*INT(n/q));
                            q:=q+2;
                       END;
                 UNTIL (q > sqrt(n)) OR (qdividesn);
            IF NOT(qdividesn) THEN write(',',n:0:0);

            {It is here that the new prime is found: this n
            is prime.  So if you want to count the primes,
            or compare one prime with the next, etc., then
            you do it here.  Remember that if you end up writing
            several instructions then you must
            replace THEN by THEN BEGIN...
            END; with the instructions in the ... }

            n:=n+2;
        END;
    END.



PROGRAM P2_1_4;
  {prints out the first no_primes, starting with 2, and stores them
      in an array }
  CONST
        size = 1000;
  VAR
    n:extended;
    i,j,no_primes:integer;
    pass:boolean;
    p:array[1..size] of extended;

  BEGIN
    writeln;
    write(' How many primes ( < 1000) ? ');
    readln(no_primes);
    p[1]:=2;
    p[2]:=3;
    writeln(2,',',3);
    n:=3;
    FOR i:=3 TO no_primes DO
      BEGIN
        pass:=FALSE;
        WHILE NOT(pass) DO
          BEGIN
            n:=n+2;
            j:=2;
            pass:=TRUE;
            WHILE(pass AND (j<i) AND (p[j]*p[j]<=n)) DO
              BEGIN
                pass:=NOT(n=p[j]*INT(n/p[j]));
                j:=j+1;
              END;
            IF pass THEN
            BEGIN
              p[i]:=n;
              write(n:0:0,',');
            END;
          END
      END;
  END.

PROGRAM P2_1_10;
 VAR
      n,q,max: extended;
      qdividesn: boolean;
      f:text;

 BEGIN
      assign(f,'primes.dat');    {If  your  disk  is in
                                    drive    A   you   may   need   to   put
                                    'a:primes.dat' here}
      rewrite(f);            {Here, f is the name of  the  file  within  the
                        program, and primes.dat is its name on the disk}
      write ('max = ?');
      readln(max);
      n:=3; 
      writeln('prime      numbers  are   '); {You might as well write primes 
                          on  the  screen  as  well  so  you can see what is 
                          happening} 
      write ('2,3'); 
      write(f,2:7);                 {These  lines  write  2  and  3  to  the  
                                    file,  leaving  7 spaces for each} 
      write(f,3:7); 
      WHILE n <= max DO
           BEGIN 
                q:=3;
                     REPEAT 
                          BEGIN 
                          qdividesn:=(n = q * INT(n/q));
                          q:= q + 2; 
                          END; 
                     UNTIL (q > SQRT(n)) OR (qdividesn);
                IF NOT(qdividesn) THEN 
                     BEGIN 
                     write(',',n:0:0); 
                     write(f,n:7:0);  {This writes n into the file,  leaving 
                                    7 spaces for it and no decimal places} 
                     END; 
                n:= n+2; 
           END; 

 close(f);   {This MUST be included, to close the file} 
 END.





 PROGRAM P2_1_11;
 VAR
      f : text;
      i,m: integer;
      prime : extended;
 BEGIN
      assign(PRIMES.DAT');    {If this is in disk drive A
                               then  you  may need to put a:PRIMES.DAT}
      reset(f);
      writeln('How many primes?  '); 
      read(m); 
      FOR i:= 1 TO m DO 
           BEGIN 
           read(f,prime); 
           write(prime:0:0,','); 
           END; 
      close(f); 
 END. 



PROGRAM P2_1_13;
VAR f:text;
    n,prime:extended;

BEGIN
    writeln('Type n');
    read(n);
    assign(f,'primes.dat');  {You may need to specify drive here}
    reset(f);
    WHILE n>1 DO
    BEGIN
        read(f,prime);
            WHILE (n/prime = INT(n/prime)) DO
            BEGIN
                write(prime:0:0,',');
                n:=n/prime;
            END;
    END;
    close(f);
    readln;
END.

PROGRAM P2_1_15;
{Trial division with one factor > 10^5 but < 10^10}

VAR f:text;
    n,prime:extended;

BEGIN
    writeln;
    writeln('Type n');
    read(n);
    assign(f,'primes.dat');  {You may need to specify drive here}
    reset(f);
    prime:=1;
    WHILE ((n>1) AND (not (eof(f))) AND (prime*prime <=n)) DO
    BEGIN
        read(f,prime);
            WHILE ((n/prime = INT(n/prime)) AND (prime*prime <=n ) ) DO
            BEGIN
                write(prime:0:0,',');
                n:=n/prime;
            END;
    END;
    IF n>1 THEN
        write(n:0:0);
    close(f);
    readln;
END.



PROGRAM P2_1_17;

 {Special factorization of r=ux, s=vy with (u,v)=(x,y)=1 and
 uv=lcm(r,s), xy=gcd(r,s).  This finds x and y from r and s
 using the file of primes}

 VAR 
      f: text; 
      r,s,x,y,p,xtemp,ytemp,ri,si: extended; 
 BEGIN 
      assign(f, 'primes.dat');  {you may need to specify drive here}
      reset(f);
      writeln('input r and s');
      readln(r,s);
      x:=1; y:=1;
      WHILE ((r>1) OR (s>1)) DO
      BEGIN
           read(f,p);
           xtemp:=1;  ytemp:=1; ri:=0; si:=0;
 {Here, ri and si correspond to the power of the prime p  in the prime-power
 expansions of r and s}
           WHILE ((r/p) = INT(r/p)) DO
           BEGIN
                ri:=ri+1;
                r:=r/p;
                xtemp:=xtemp*p;
           END;
           WHILE ((s/p) = INT(s/p)) DO
           BEGIN
                si:=si+1;
                s:=s/p;
                ytemp:=ytemp*p;
           END;
           IF (ri>=si) THEN
              y:=y*ytemp         {Note no semicolon here!}
           ELSE
                x:=x*xtemp;
           writeln('p=',p:0:0,'  x=',x:0:0,'  y=',y:0:0);
      END;
 END.


PROGRAM P2_1_21;
 VAR
     x,y,n: extended;
     ydividesx: boolean;
     f: file of extended;
 BEGIN
     assign(f, 'smallprimes');  {Just to change the name}
     rewrite(f);
     writeln('n = ?');
     readln(n);
     writeln('prime numbers are ');
     x:= 2;
     write(f,x);  {This is the way to write 2 into the file}
     write(2,',');  {So 2 is written to the screen too}
     x:= 3;
     write(f,x);
     write(3);
     WHILE x <= n DO
         BEGIN 
             y:= 3; 
             REPEAT
                 BEGIN 
                     ydividesx:= (x = y* INT(x/y)); 
                     y:= y+2;
                 END; 
             UNTIL (y > SQRT(x)) OR (ydividesx); 
             IF NOT (ydividesx) THEN 
                 BEGIN 
                 write(f,x); 
                 write(',',x:0:0); 
                 END; 
              x:=x+2; 
           END; 
     close(f);
 END.

PROGRAM P2_1_21A;
 VAR
     k,m:longint;
     n:extended;   {Or integer and real will do}
     f: file of extended;
 BEGIN
     assign(f, 'smallprimes');
     reset(f);
     writeln('type the value of k');  {The k'th element will be pulled out}
     readln(k);
         seek(f,k-1);   {Why k-1? Because files are numbered from 0 !}
         read(f,n);     {So the k'th element is called n}
         write(n:0:0);
     close(f);
 END.





PROGRAM P2_2_1;  {primes between r and r + 2n -1}
VAR
     r,c,m,s,i,j: extended;
     n,k: integer;
     key:array[1..10000] of integer; {So n should not be larger than 10000}
BEGIN
     writeln;
     writeln('What is the value of r (even)?  ');  readln(r);
     writeln('What is the value of n (< 10000)');  readln(n);
     c:= r + 2*n + 1;
     m:= INT(SQRT(c));
     FOR k:= 1 TO n DO
          key[k]:= 0;
     i:= 3;
     WHILE i <= m DO
          BEGIN
               s:= 2*INT((r+i)/(2*i)) + 1;
               IF s = 1 THEN s:= 3;
               j:= s;
               WHILE j <= INT(c/i) DO
                    BEGIN
                         k:=TRUNC((i*j - r  +  1)/2);   {This  makes  k  an 
                                                            integer}
                         key[k]:= 1;
                         j:= j + 2;
                    END;
               i:= i + 2;
          END;
     FOR k:= 1 TO n DO
          BEGIN
               IF key[k] = 0 THEN
               write((r + 2*k - 1):0:0,',');
          END;
END.

PROGRAM P2_4_1;
  {uses the sieve of Eratosthenes to print out the primes up to 2000 }
  CONST
       f:array[0..12] of integer = (3,5,7,11,13,17,19,23,29,31,37,41,43);
  VAR
    i,j,prime:integer;
    key:array[1..2000] of boolean;


  BEGIN
    writeln;
    FOR i:=1 TO 2000 DO
      key[i]:=TRUE;
    j:=0;
    REPEAT
      BEGIN
        prime:=f[j];
        i:=prime*prime;
        WHILE i<= 2000 DO
          BEGIN
            key[i]:=FALSE;
            i:=i+2*prime;
          END;
          j:=j+1;
      END;
    UNTIL (prime=43);
    writeln(' primes less than 2000 are :-');
       i:=3;
    WHILE i<= 2000 DO
      BEGIN
        IF key[i] THEN
        write(',',i);
        i:=i+2
      END;
    END.


 PROGRAM P2_4_3;
 CONST
     n=1000;
 VAR
     p,q,i,j,k: longint;  {Better than integer because q gets big}
     a: array[1..n] of boolean;
 BEGIN 
     FOR j:=1 TO n DO 
         a[j]:= true; 
     p:= 3;   q:= 4;    j:=1; 
     WHILE j<=n DO 
         BEGIN 
             IF a[j]= true THEN 
                 BEGIN 
                     write(p, ','); 
                     k:=q; 
                     WHILE k<=n DO 
                         BEGIN
                             a[k]:=false;
                             k:=k+p;
                         END;
                 END;
             j:= j+1; p:= p+2; q:=q+2*p-2; {Think of 2q+1=p }
         END;
 END.


 PROGRAM P2_4_8;
  {computes the number of primes <= x, where x is at most 9,699,690,
       using Legendre's sieve method}
  CONST
       max=443;
  VAR
    f:text;
    i,p:integer;
    s0,s1,x,rootx,total,ans:extended;
    primarr:array[1..max] of extended;
    finish:boolean;

    FUNCTION sieve(no_terms:integer):extended;
      {computes the terms for the Legendre's sieve method }
    VAR
          cond,cond4,exit:boolean; prod,s:extended; i,j:integer;
          count:array[1..10] of integer;

        PROCEDURE test(VAR prod:extended; VAR cond,exit:boolean);
        VAR
            i:integer;  cond1,cond2,cond3:boolean;
        BEGIN
            prod:=1;
            FOR i:=1 TO no_terms DO
                BEGIN
                   prod:=prod*primarr[count[i]];
                   cond1:=prod > x;
                   cond2:=primarr[count[no_terms]] > rootx;
                   cond:=cond1 OR cond2
                END;
            cond3:=count[no_terms]=count[1]+(no_terms-1);
            IF cond AND cond3 THEN
                exit:=TRUE
        END;         {of test}

    BEGIN  {sieve}
      exit:=FALSE;  prod:=1;  s:=0;
      FOR i:=1 TO no_terms DO
        count[i]:=i;
      test(prod,cond,exit);
      WHILE NOT(exit) DO
        BEGIN
           WHILE cond AND NOT(exit) DO
              BEGIN
              j:=2;
              WHILE (j<= no_terms-1) AND NOT(exit)  DO
                 BEGIN
                    cond4:=count[no_terms]=count[j]+(no_terms-j);
                    IF cond AND cond4 THEN
                       BEGIN
                       count[j-1]:=count[j-1]+1;
                       FOR i:=j TO no_terms DO
                          count[i]:=count[i-1]+1;
                       test(prod,cond,exit);
                       j:=no_terms
                       END
                    ELSE
                    j:=j+1
                 END;
              END;
          WHILE NOT(cond) DO
              BEGIN
                 s:=s+INT(x/prod);
                 count[no_terms]:=count[no_terms]+1;
                 test(prod,cond,exit);
                 IF cond THEN
                     BEGIN
                        count[no_terms-1]:=count[no_terms-1]+1;
                        count[no_terms]:=count[no_terms-1]+1;
                        test(prod,cond,exit);
                     END;
              END;
        END;
      sieve:=s
    END;   {of sieve}

  BEGIN   {main program}
    assign(f,'primes.dat');  {You may need to specify drive here}
    reset(f);
    writeln('Program to compute the number of primes <= x=9,699,690');
    writeln('using the legendre''s sieve technique');
    write('x= ? ');   readln(x);
    rootx:=sqrt(x);
    FOR i:=1 TO max DO
      read(f,primarr[i]);
    total:=0;  s0:=INT(x)-1;  writeln('s0= ',s0:8:0);
    i:=1; s1:=0;
    WHILE primarr[i] <= rootx DO
       BEGIN
          s1:=s1+INT(x/primarr[i]);
          i:=i+1;
       END;
    writeln('s1= ',s1:8:0);
    total:=s0-s1;
    finish:=FALSE;
    WHILE NOT(finish) DO
       BEGIN
          i:=2;
          WHILE (i<= 7) or (NOT finish) DO
             BEGIN
                ans:=sieve(i);
                IF ans=0 THEN
                finish:=TRUE
                ELSE
                   BEGIN
                      writeln('s',i:1,'= ',ans:8:0);
                      IF ODD(i) THEN
                         ans:=-ans;
                      total:=total+ans
                   END;
                i:=i+1;
             END;
       END;
    p:=0; i:=1;
    WHILE primarr[i]<= rootx DO
       BEGIN
          p:=p+1;
          i:=i+1;
       END;
    writeln(' p= ',p:8);
    write('the number of primes up to ',x:8:0);
    writeln(' is ',total+p:8:0);
    close(f);
    writeln
  END.



 PROGRAM P2_4_12;
 {Computes the function phi(x,a)}
 VAR
      b,k   : integer;
      y     : real;
      f     : text;
      p     : array[1..100] of integer;  {This will be an array  of  primes,
 drawn from the file PRIMES.DAT.  Note that this array is available  to  the
 function phi without further declaration}

 FUNCTION phi(x:real; a:integer):integer;
 BEGIN
      IF a>1 THEN
           phi:= phi(x, a-1)-phi(x/p[a], a-1)
      ELSE
           phi:= trunc((x+1)/2);
 END;

 BEGIN {main program}
      writeln;
      assign(f, 'primes.dat');  {You may need to specify drive here}
      reset(f);
      for k:= 1 to 100 do
           read(f, p[k]);
      close(f);

      writeln('Type the values of x and a separated by a space');
      readln(y,b); 
      writeln(phi(y,b)); 
 END. 


PROGRAM P4_2_3;
{PowerRule}

VAR
 a, asaved, n, nsaved, r, m : extended;


PROCEDURE PowerRule (base, power, modulus:extended;
                            VAR        res:extended);
VAR
     d: extended;

BEGIN
res:=1;
    WHILE power>0 DO
      BEGIN
         d:=power-2*INT(power/2);
         IF d=1 THEN
           BEGIN
              res:=base*res;
              res:=res-modulus*INT(res/modulus);
           END;
         base:=base*base;
         base:=base-modulus*INT(base/modulus);
         power:=(power-d)/2;
      END;

END; {Of PowerRule}


BEGIN {of main program}

  writeln;
  writeln ('Type the values of a, n and m, separated by spaces');
  readln (asaved, nsaved, m);

  n:=nsaved;
  a:=asaved;

  writeln;
  PowerRule(a,n,m,r);

  writeln;
  writeln ('The value of (', asaved:0:0, ' to the power ', nsaved:0:0,
           ') mod ', m:0:0, ' is ', r:0:0)

END. { of main program }




PROGRAM P4_3_4;   {PowerRule using Head's algorithm}
VAR
 a, asaved, n, nsaved, r, m, capT, t:extended;


PROCEDURE MultiplyModM (x,y,modulus : extended;
                      VAR prod  : extended);
{Note that in most examples it will be the  same  modulus  m  throughout  a
calculation, so it is more efficient to calculate T (here called capT)  and
t (see 3.1) once during the main program than to calculate  it  within  the
procedure each time it is used.  The result, prod, is xy mod m}

VAR
   a,b,c,d,z,e,f,v,v1,g,h,j,j1,k: extended;   {v1 and j1     are      extra
           variables used in the calculation when allowing for t to be < 0}
BEGIN
   a:= INT(x/capT);  b:= x - a * capT;
   c:= INT(y/capT);  d:= y - c * capT;
   z:= a * d + b * c;   z:= z - modulus * INT(z/modulus);
   e:= INT(a*c/capT);  f:= a * c - e * capT;
   v:= z + e * t;
   v1:= v/modulus;
   IF (v>0) OR (v1=INT(v1)) THEN v:= v - modulus * INT(v1)
      ELSE  v:=v - modulus * (INT(v1) - 1);  {This allows for v < 0}
   g:= INT(v/capT);  h:= v - g * capT;
   j:= (f + g) * t;
   j1:=j/modulus;
   IF (j>0) OR (j1=INT(j1)) THEN j:= j - modulus * INT(j1)
      ELSE j:= j - modulus * (INT(j1) - 1);
   k:= j + b * d;  k:= k - modulus * INT(k/modulus);
   prod:= h * capT + k; prod:= prod - modulus * INT(prod/modulus);
END;  {of the PROCEDURE MultiplyModM}



PROCEDURE PowerRule (base, power, modulus:extended;
                            VAR        res:extended);
VAR
     d: extended;

BEGIN
res:=1;
    WHILE power>0 DO
    BEGIN
         d:=power-2*INT(power/2);
         IF d=1 THEN
              MultiplyModM(base,res,modulus,res);  {Here we use Head's
                                                       algorithm}
         MultiplyModM(base,base,modulus,base);        {also here}
         power:=(power-d)/2;
    END;

END; {Of PowerRule}



BEGIN {of main program}

  writeln;
  writeln ('Type the values of a, n and m separated by spaces');
  readln (asaved, nsaved, m);

  capT:=int(sqrt(m)+0.5);
  t:=capT*capT-m;
  n:=nsaved;
  a:=asaved;

  writeln;
  PowerRule(a,n,m,r);

  writeln;
  writeln ('The value of (', asaved:0:0, ' to the power ', nsaved:0:0,
           ') mod ', m:0:0, ' is ', r:0:0);

END { of main program }.



PROGRAM P5_1_5; {Miller's Test as in the book, with improved
                                power algorithm}

VAR
 n, b, n1, r:extended;


PROCEDURE PowerRule (base, power, modulus:extended;
                            VAR        res:extended);
VAR
     d: extended;

BEGIN
res:=1;
    WHILE power>0 DO
    BEGIN
         d:=power-2*INT(power/2);
         IF d=1 THEN
         BEGIN
              res:=base*res;
              res:=res-modulus*INT(res/modulus);
         END;
         base:=base*base;
         base:=base-modulus*INT(base/modulus);
         power:=(power-d)/2;
    END;

END; {Of Powerrule}



BEGIN {of main program}

  writeln;
  writeln ('Type the values of n and b');
  readln (n,b);

  n1:=n-1;
  r:=1;

  WHILE (n1=INT(n1)) AND (r=1) DO
        BEGIN
             write('Power = ',n1:0:0);
             PowerRule(b,n1,n,r);
             n1:=n1/2;
             writeln('   Residue is  ',r:0:0);
        END;
  IF (r=1) OR (r = n-1) THEN
        writeln(n:0:0, ' has passed Millers test to base ',b:0:0)
  ELSE
        writeln(n:0:0, ' has failed Millers test to base ',b:0:0);
  readln;
  END.  {of main program}




PROGRAM P5_1_6a; {Miller's Test as in the book, with improved
                                power algorithm}


{THIS VERSION USES HEAD'S ALGORITHM}


VAR
 n, b, n1, r, capT, t:extended;

PROCEDURE MultiplyModM ({inputs}  x,y,modulus : extended;
                        {output}  VAR prod  : extended);

  { Multiplies x and y modulo m storing the result in prod. }

  VAR
    a, b, c, d, z, e, f, v, v1, g, h, j, j1, k : extended;

  BEGIN

    a:=int(x/capT); b:=x - a* capT;
    c:=int(y/capT); d:=y - c* capT;
    z:=a*d + b*c;
    z:=z - modulus*int(z/modulus);
    e:=int(a*c/capT); f:=a*c-e*capT;
    v:=z + e*t;
    v1:=v/modulus;
    if (v>0) or (v1=int(v1)) then v:=v-modulus*int(v1)
        else v:=v-modulus * (int(v1)-1);
    g:=int(v/capT);    h:=v - g*capT;
    j:=(f + g)*t;
    j1:=j/modulus;
    if (j>0) or (j1=int(j1)) then j:=j-modulus*int(j1)
        else j:= j-modulus*(int(j1)-1);
    k:=j + b*d;
    k:=k-modulus*int(k/modulus);
    prod:=h*capT+k;
    prod:=prod-modulus*int(prod/modulus);

  END { of MultiplyModM };


PROCEDURE PowerRule (base, power, modulus:extended;
                            VAR        res:extended);
VAR
     d: extended;

BEGIN
res:=1;
    WHILE power>0 DO
    BEGIN
         d:=power-2*INT(power/2);
         IF d=1 THEN
         BEGIN
              MultiplyModM(base,res,modulus,res);
         END;
         MultiplyModM(base,base,modulus,base);
         power:=(power-d)/2;
    END;

END; {Of Powerrule}



BEGIN {of main program}

  writeln;
  writeln ('Type the values of n and b');
  readln (n,b);

  capT:=int(sqrt(n)+0.5);
  t:=capT*capT-n;
  n1:=n-1;
  r:=1;

  WHILE (n1=INT(n1)) AND (r=1) DO
        BEGIN
             write('Power = ',n1:0:0);
             PowerRule(b,n1,n,r);
             n1:=n1/2;
             writeln('   Residue is  ',r:0:0);
        END;
  IF (r=1) OR (r = n-1) THEN
        writeln(n:0:0, ' has passed Millers test to base ',b:0:0)
  ELSE
        writeln(n:0:0, ' has failed Millers test to base ',b:0:0);
  END.  {of main program}



PROGRAM P6_3_5;   {n-1 over q test for primality}
VAR
 a, asaved, n, nsaved, r, m, capT, t:extended;
 newn,q:extended;


PROCEDURE MultiplyModM (x,y,modulus : extended;
                      VAR prod  : extended);
{Note that in most examples it will be the  same  modulus  m  throughout  a
calculation, so it is more efficient to calculate T (here called capT)  and
t (see 3.1) once during the main program than to calculate  it  within  the
procedure each time it is used.  The result, prod, is xy mod m}

VAR
   a,b,c,d,z,e,f,v,v1,g,h,j,j1,k: extended;   {v1 and j1     are      extra
           variables used in the calculation when allowing for t to be < 0}
BEGIN
   a:= INT(x/capT);  b:= x - a * capT;
   c:= INT(y/capT);  d:= y - c * capT;
   z:= a * d + b * c;   z:= z - modulus * INT(z/modulus);
   e:= INT(a*c/capT);  f:= a * c - e * capT;
   v:= z + e * t;
   v1:= v/modulus;
   IF (v>0) OR (v1=INT(v1)) THEN v:= v - modulus * INT(v1)
      ELSE  v:=v - modulus * (INT(v1) - 1);  {This allows for v < 0}
   g:= INT(v/capT);  h:= v - g * capT;
   j:= (f + g) * t;
   j1:=j/modulus;
   IF (j>0) OR (j1=INT(j1)) THEN j:= j - modulus * INT(j1)
      ELSE j:= j - modulus * (INT(j1) - 1);
   k:= j + b * d;  k:= k - modulus * INT(k/modulus);
   prod:= h * capT + k; prod:= prod - modulus * INT(prod/modulus);
END;  {of the PROCEDURE MultiplyModM}



PROCEDURE PowerRule (base, power, modulus:extended;
                            VAR        res:extended);
VAR
     d: extended;

BEGIN
res:=1;
    WHILE power>0 DO
    BEGIN
         d:=power-2*INT(power/2);
         IF d=1 THEN
              MultiplyModM(base,res,modulus,res);  {Here we use Head's
                                                       algorithm}
         MultiplyModM(base,base,modulus,base);        {also here}
         power:=(power-d)/2;
    END;

END; {Of PowerRule}



BEGIN {of main program}


writeln ('Type the value of n ');
readln (nsaved);
m:=nsaved;   {We always work modulo the original n}
capT:=INT(SQRT(m)+0.5);
t:=capT*capT-m;
asaved:=2;  {The first base we try is 2}
writeln('a=2; input q=0 to go to another a and a=0 to finish');
         WHILE (asaved <> 0) DO
            BEGIN
            q:=1;    {We first test n-1 over q for q=1}
            WHILE (q <> 0) DO
                BEGIN
                   n:=(nsaved-1)/q;
                   newn:=n;
                   a:=asaved;
                   PowerRule (a,n,m,r);
                   writeln ('The value of (', asaved:0:0, ' to the power ',
                         newn:0:0, ') mod ', m:0:0, ' is ',r:0:0);
                   writeln('Type q');  {Enter q=0 to end this loop}
                            readln(q);
                   END;
              writeln('Type a');    {Enter a=0 to end this loop}
              readln(asaved);
              END;
      END. { of main program }



PROGRAM P7_2_3;  {RSA public key cipher}

VAR
    message_length,num_blocks,block_length,i,j : integer;
    s,numstring : string;
    x: array[1..32] of char;
    a,b,c: array[1..32] of extended;
    power_of_ten, max_power: extended;
    t, capT: extended;

    n,p,q,phi:extended;   {n is the product of two primes p and q}
    enc,dec:extended;       {enc is the encrypting exponent, and
                               must be coprime to (p-1)(q-1).
                             dec is the decrypting exponent.}
    u:integer;
    success:boolean;      {This is used to signal when an encrypting
                           exponent has been found which is coprime
                           to (p-1)(q-1)}


PROCEDURE MultiplyModN (x,y,modulus : extended; VAR prod  : extended);

  { Multiplies x and y modulo "modulus" (here this is n) using Head's
  algorithm and stores the result in prod. }

  VAR
    a, b, c, d, z, e, f, v, v1, g, h, j, j1, k : extended;

  BEGIN

    a:=int(x/capT); b:=x - a* capT;
    c:=int(y/capT); d:=y - c* capT;
    z:=a*d + b*c;
    z:=z - modulus*int(z/modulus);
    e:=int(a*c/capT); f:=a*c-e*capT;
    v:=z + e*t;
    v1:=v/modulus;
    if (v>0) or (v1=int(v1)) then v:=v-modulus*int(v1)
        else v:=v-modulus * (int(v1)-1);
    g:=int(v/capT);    h:=v - g*capT;
    j:=(f + g)*t;
    j1:=j/modulus;
    if (j>0) or (j1=int(j1)) then j:=j-modulus*int(j1)
        else j:= j-modulus*(int(j1)-1);
    k:=j + b*d;
    k:=k-modulus*int(k/modulus);
    prod:=h*capT+k;
    prod:=prod-modulus*int(prod/modulus);

  END { of MultiplyModN };


PROCEDURE Powerrule (base, power, modulus:extended;
                            VAR        res:extended);
VAR
     e: extended;

BEGIN
res:=1;
    WHILE power>0 DO
    BEGIN
         e:=power-2*INT(power/2);
         IF e=1 THEN
              MultiplyModN(base,res,modulus,res);
         MultiplyModN(base,base,modulus,base);
         power:=(power-e)/2;
    END;

END; {Of Powerrule}



PROCEDURE Inverse(a:extended; VAR inv:extended);
{This finds the inverse of a modulo phi (=(p-1)(q-1), or p-1 when q=1)
when they are coprime, and stores the result in inv}

VAR
   b,a0,b0,r,s,s1,s2,t,t1,t2,q,quotient:extended;

BEGIN
   b:= phi; a0:=a;  b0:=b;
   s2:=1;   s1:=0;
   t2:=0;   t1:=1;
   WHILE b > 0 DO
       BEGIN
           quotient:=a/b;
               IF ((quotient>0) OR (quotient = INT(quotient))) THEN
                   q:=INT(quotient)
               ELSE
                   q:=INT(quotient)-1;
            r:=a - b * q;
            s:=s2 - q * s1;
            t:=t2 - q * t1;
            a:=b;
            b:=r;
            s2:=s1;
            t2:=t1;
            s1:=s;
            t1:=t;
        END;
   IF (a=1) THEN
        BEGIN
            writeln(a0:0:0,' and ',b0:0:0, ' are coprime');
            quotient:=s2/b0;
               IF ((quotient>0) or (quotient = int(quotient))) THEN
                   s2:=s2-int(quotient)*b0
               ELSE
                   s2:=s2-(int(quotient)-1)*b0;
            inv:=s2;
            success:=true;
        END
    ELSE
        writeln('Your choice of e is not coprime to (p-1)(q-1)');
END;  {of Inverse}



BEGIN  {main program}
    writeln('Type primes p and q (or p prime and q = 1) ');
    readln(p,q);
    n:=p*q;

    capT:=INT(SQRT(n)+0.5);
    t:=capT*capT-n;         {These are needed in MultiplyModN}

    IF (q>1) THEN                  {Using just a single prime}
        phi:=(p-1)*(q-1)               {as in exponentiation ciphers}
    ELSE                           {we can put q=1}
        phi:=p-1;

    writeln('n = ',n:0:0);
    success:=false;
    WHILE (success=false) DO
        BEGIN
        writeln('Type e, the encryption exponent');
        readln(enc);
        Inverse(enc,dec);
        END;

        writeln('Inverse is ',dec:0:0);

    writeln('Enter the block length');
        {This refers to the number of CHARACTERS in a block.  Using
        just the letters A..Z and the space, all ACSII codes are <=90.
        Hence 9090..90 with block_length 90's must be < n}
    readln(block_length);

    writeln('Enter the message length');
    readln(message_length);

    writeln('Type the string in CAPITALS');
    readln(s);

    num_blocks:=message_length div block_length;  {Integer division!}
    FOR i:= 1 to num_blocks DO
        BEGIN
        a[i]:=0;
            FOR j:=1 to block_length DO
                BEGIN
                    x[j]:=s[j+block_length*(i-1)];
                    a[i]:=a[i]*100 + ord(x[j]);
                END;
         writeln('This block becomes the number ',a[i]:0:0);
            {This is the (2*block_length)-digit number
            associated to the i-th block of characters}

         Powerrule(a[i],enc,n,b[i]);
         writeln('Encrypted equivalent of this block is ',b[i]:0:0);
             {b[i] is the encrypted equivalent
              of the i-th block of characters}

         Powerrule(b[i],dec,n,c[i]);
         writeln('Applying the inverse gives ',c[i]:0:0);
            {This should give the same number as a[i]!}
         writeln;

    END ;

    {Now we need to turn the numbers c[i] back into characters so we
    need to read the number c[i] in groups of two digits, from
    left to right}

    writeln('The decrypted message is: ');
    writeln;

    power_of_ten:=1;
    FOR j:=2 TO block_length DO
        power_of_ten:=power_of_ten*100;
        max_power:=power_of_ten;

    FOR i:=1 TO num_blocks DO
        BEGIN
             power_of_ten:=max_power;
                 FOR j:=1 TO block_length DO
                     BEGIN
                         u:=TRUNC(c[i]/power_of_ten);
                         c[i]:=c[i] - u*power_of_ten;
                         power_of_ten:=power_of_ten/100;
                         Write(chr(u));
                     END;
        END;
        writeln;
END.



PROGRAM P9_1_3;

{  Calculates the function d(n), where n is input by the user,
   using prime numbers obtained from a file. }


VAR
  prime                   : extended;
  n, nsaved, total, power : extended;
  f                       : text;


BEGIN

  writeln ('Type the value of n');
  readln (n);
  nsaved:=n;

  Assign (f,'primes.dat');  {You may need to specify drive here}
  reset (f);

  total:=1;
  WHILE (n>1) DO
  BEGIN
    read(f,prime);
    power:=1;
    WHILE n/prime = Int(n/prime) DO
    BEGIN
      n:=Int(n/prime);
      power:=power + 1
    END;
    total:=total*power
  END;

    writeln ('d(', nsaved:0:0, ') = ', total:0:0)

END.


PROGRAM P9_2_5;

{  Calculates the function sigma(n), where n is input by the user,
   using prime numbers obtained from a file. }


VAR
  prime                   : extended;
  n, nsaved, sigma, power_of_prime, sum_of_powers: extended;
  f                       : text;


BEGIN

  writeln ('Type the value of n');
  readln(n);
  nsaved:=n;

  Assign (f,'primes.dat');  {You may need to specify drive here}
  reset (f);

  sigma:=1;
  WHILE (n>1) DO
  BEGIN
    read(f,prime);
    power_of_prime:=1;
    sum_of_powers:=1;
    WHILE (n/prime = INT(n/prime)) DO
    BEGIN
      n:=INT(n/prime);
      power_of_prime:=power_of_prime*prime;
      sum_of_powers:=sum_of_powers + power_of_prime;
    END;
    sigma:=sigma * sum_of_powers;
  END;

    writeln ('sigma(', nsaved:0:0, ') = ', sigma:0:0)

END.



PROGRAM P10_2_6;
{Continued fraction of root n}


VAR
    n,capP,capQ,rootn,x,a: extended;
    k,flag:integer;

BEGIN
    writeln('input the value of n');
    readln(n);
    rootn:=INT(SQRT(n));
    capP:=0; capQ:=1;
    k:=0; flag:=0;
    WHILE (flag<2) DO
        BEGIN
            IF (flag=1) THEN flag:=2;
            x:=(capP+rootn)/capQ;
            a:=INT(x);
            {These are xk and ak}
            writeln(k,'     ',capP:0:0,'       ',capQ:0:0,'     ',a:0:0);
            capP:=a*capQ-capP;
            capQ:=(n-capP*capP)/capQ;
            {These are Pk+1 and and Qk+1}
            k:=k+1;
            IF (capQ=1) THEN flag:=1;
        END;
END.



PROGRAM P10_5_5;


CONST
  arraysize = 100;

VAR
  n, rootn,  capP, capQ, Qstar, p, pp, ppp, p1: extended;
  maxprime, a, b, lhs, rhs, diff, sum, gcd, number : extended;
  p_array,   Qstar_array                   : ARRAY[1..arraysize] OF extended;
  maxk, i, k, rownum, terms                : integer;
  yorn                                     : char;   {A yes/no answer}
  f                                        : text;   {The file of primes}


PROCEDURE Factors(number, maxprime: extended);
VAR
    prime:extended;

BEGIN
    IF (number < 0) THEN
    BEGIN
        write('-');
        number:= -number;
    END;

    reset(f);
    WHILE number>1 DO
    BEGIN
        read(f,prime);
        IF(prime <= maxprime) THEN
        BEGIN
            WHILE (number/prime = INT(number/prime)) DO
            BEGIN
                write(prime:0:0,'*');
                number:=number/prime;
            END;
        END
        ELSE
            BEGIN
               write('***'); {This indicates that the prime exceeds
                                     maxprime}
               number:=1;
            END;
    END;
    writeln;
    close(f);
END;   {of Factors}


PROCEDURE Euclid(a,b     : extended;
                  VAR gcd: extended);
VAR
    r : extended;

BEGIN
    WHILE (b>0) DO
        BEGIN
  r:= a - b*INT(a/b);
            a:= b;
            b:= r;
        END;
    gcd:= a;
END; {Of Euclid}



BEGIN {********  main program  ********}

    Assign (f,'primes.dat');  {You may need to specify the drive here}

    writeln ('Type the value of n followed by the maximum k+1');
    readln (n,maxk);
    writeln ('Type the value of a bound for the primes used');
    writeln ('in factorising by trial division');
    readln (maxprime);
    writeln (' n = ', n:0:0);
    writeln (' k + 1   p[k] mod n   Q[k+1]* ');

    rootn:=INT(SQRT(n));
    capP:=0;
    capQ:=1;
    pp:=1;
    ppp:=0;
    
    i:=1;
        WHILE (i<=maxk)  DO
            BEGIN
            a:=INT((capP + rootn)/capQ);
            p:=a*pp + ppp;
            p:=p - INT(p/n)*n;
            capP:=a*capQ - capP;
            capQ:=(n - capP*capP)/capQ;
                IF (i<>Int(i/2)*2) THEN
                   Qstar:=-capQ
                ELSE
                   Qstar:=capQ;
                IF (p>n/2) THEN
                   p1:=p - n
                ELSE
                   p1:=p;
     {So p1 is the smallest residue in absolute value of p}
            ppp:=pp;
            pp:=p;
            write (i:4, p1:11:0, Qstar:11:0, ' = ');
            Factors (Qstar,maxprime);
            p_array[i]:=ABS(p1);   {We don't need to bother with
                                        the sign of p}
            Qstar_array[i]:=Qstar;
            i:= i + 1
        END {of WHILE i<=k};

{*****   The rest of the program tries the factoring method on n *****}

      writeln ('Do you wish to try a factor? (y/n)');
      readln (yorn);

      WHILE (yorn='y') OR (yorn='Y') DO
          BEGIN
              lhs:=1;
              rhs:=1;
              writeln ('How many terms to make a square?');
              readln (terms);
    
              FOR i:=1 TO terms DO
                  BEGIN
                      writeln ('What is the number of the row?');
                      readln (rownum);  
                      lhs:=lhs*p_array[rownum];
                      lhs:=lhs - n*INT(lhs/n);
                      rhs:=rhs*Qstar_array[rownum]
                  END;

              rhs:=sqrt(rhs);  
                
        { At this point, (lhs) contains the product ,modulo n, of the p
          and (rhs) contains the square root of the product of the Q*.
          Next the difference (diff) will be computed. Its g.c.d.
          with n is a factor (hopefully proper) of n. }
                   
        writeln ('lhs = ', lhs:0:0, ' rhs = ', rhs:0:0);   
        diff:=lhs - rhs;
        Euclid(diff,n,gcd);
        writeln ('A factor of ', n:0:0, ' is ', gcd:0:0);
        sum:=lhs + rhs;
        Euclid(sum,n,gcd);
        writeln('A factor of ', n:0:0, ' is ', gcd:0:0);
        writeln ('Do you wish to try another factor? (y/n)');
        readln (yorn);
        END  { of yorn is 'yes' }

END { of main program }.

      
      
PROGRAM P11_4_6;

{ Calculates the Jacobi symbol (n/k) where
  n and k are > 0 and n and k are coprime.
  Successive values of the ratio are displayed, with a minus sign as
  appropriate, as the calculation proceeds.
  The user is first prompted for a value of n and then for as many values
  of k as are required, a value of two or less signifying the end of
  the run. }
                  
                      
VAR
  count                             : integer; 
  x, y, n, nsaved, k, ksaved, tempk : extended;
  sign                              : boolean;
                  

              
PROCEDURE WriteRatio ({input} n, k : extended; sign : boolean);

  { Displays the intermediate values of (n/k) with appropriate sign. }
          
  BEGIN
                   
    writeln;
    IF NOT sign THEN
      write (' -(', n:0:0, '/', k:0:0,')')
    ELSE
      write ('  (', n:0:0, '/', k:0:0,')')
        
  END { of WriteRatio };
        

BEGIN {main program}

  writeln ('Type the value of n');
  readln (n);
  nsaved:=n;
  writeln ('Now type the first value of k');
  readln (k);   
 ksaved:=k;

  WHILE ksaved>2 DO
  BEGIN
    sign:=TRUE;
    WHILE n>1 DO
    BEGIN
      n:=n - k*int(n/k);
      WriteRatio (n,k,sign);
      count:=0;       
   
      WHILE n/2=Int(n/2) DO
      BEGIN
        n:=n/2;
        count:=count + 1
      END;
              
      IF count/2 <> Int(count/2) THEN
      BEGIN
        x:=k - 8*int(k/8);
        IF (x=3) OR (x=5) THEN
        sign:= NOT sign;
      END;         

      IF count>0 THEN
        WriteRatio (n,k,sign);

      IF n>1 THEN
      BEGIN
        tempk:=k;
        k:=n;
        n:=tempk;
        x:=k - 4*int(k/4);
        y:=n - 4*int(n/4);
        IF (x=3) AND (y=3) THEN   
          sign:= NOT sign;
        WriteRatio (n,k,sign);
      END;
    END { of WHILE n>1 };
    writeln;
    IF (sign = TRUE) THEN   
        writeln('Hence (',nsaved:0:0,'/',ksaved:0:0,') = 1')
        ELSE
        writeln('Hence (',nsaved:0:0,'/',ksaved:0:0,') = -1');
    
    writeln ('Type the next value of k (a value less than 3 ends the run)');
    readln (k);
    ksaved:=k;
    n:=nsaved         
  END { of WHILE ksaved>2 loop };
      
END { of main program }.
        
        
PROGRAM P11_5_5;   {Implementing Shanks's method for solving x^2=a (modp)
                     using a quadratic non-residue z (mod p)}
      

CONST
  maxbit = 100;
        
VAR
 x,p,a,p1,odd,z,c,n,b,y,k1,x2 : extended;
 s,i,looplength,index: integer;
        

PROCEDURE PowerRule (base, power, modulus:extended;
                            VAR        res:extended);
VAR
     d: extended;

BEGIN
res:=1;
    WHILE power>0 DO
      BEGIN
         d:=power-2*INT(power/2);
         IF d=1 THEN
           BEGIN
              res:=base*res;
              res:=res-modulus*INT(res/modulus);

      END;
         base:=base*base;
         base:=base-modulus*INT(base/modulus);
         power:=(power-d)/2;
      END;
    
END; {Of PowerRule}

  
BEGIN {of main program}
    writeln('type p,a and z');
    readln(p,a,z);
    p1:=p-1;
    index:=0;
    WHILE (p1/2 = INT(p1/2)) DO
        BEGIN
            p1:=p1/2;
            index:=index+1;
        END;   
        s:=index;
        odd:=p1;
        PowerRule(z,odd,p,c);
        PowerRule(a,odd,p,n);  
        writeln('n:=',n:0:0);
        k1:=(odd+1)/2;
        PowerRule(a,k1,p,x);
        writeln('x:=',x:0:0);
   
        WHILE (n>1) DO
            BEGIN
                looplength:=s;
                y:=n;
                FOR i:=1 TO looplength DO
                    BEGIN
                        IF (y=1) THEN
                            BEGIN
                                y:=c;
                                s:=i-1;
                            END
                            ELSE
                            BEGIN
                                y:=y*y;
                                y:=y-p*INT(y/p);
                            END;
                    END;
                b:=y;
                c:=b*b;  c:=c-p*INT(c/p);
                x:=b*x;  x:=x-p*INT(x/p);
                writeln('x=',x:0:0);
                n:=c*n;  n:=n-p*INT(n/p);
                writeln('n=',n:0:0);
            END;
writeln('x:=',x:0:0);
x2:=x*x;                            {These three lines just for checking}
x2:=x2-p*INT(x2/p);
writeln('x*x=',x2:0:0);
readln;