(* ----------------------------------------------------------------------------
 * $Id: SACPRIM.mi,v 1.4 1995/11/05 08:45:56 kredel Exp $
 * ----------------------------------------------------------------------------
 * This file is part of MAS.
 * ----------------------------------------------------------------------------
 * Copyright (c) 1989 - 1992 Universitaet Passau
 * ----------------------------------------------------------------------------
 * $Log: SACPRIM.mi,v $
 * Revision 1.4  1995/11/05 08:45:56  kredel
 * Cosmetic Changes.
 *
 * Revision 1.2  1994/10/06  11:34:52  kredel
 * Cosmetic changes
 *
 * Revision 1.1.1.1  1993/06/11  12:54:12  kredel
 * Initial Version 0.7 of MAS from University of Passau
 *
 * Revision 1.3  1992/10/15  16:28:20  kredel
 * Changed rcsid variable
 *
 * Revision 1.2  1992/02/12  13:19:17  pesch
 * Moved CONST Definition to the right place.
 *
 * Revision 1.1  1992/01/22  15:08:44  kredel
 * Initial revision
 *
 * ----------------------------------------------------------------------------
 *)

IMPLEMENTATION MODULE SACPRIM;

(* SAC Factorization and Prime Number Implementation Module. *)



(* Import Lists and Declarations. *)

FROM MASELEM IMPORT MASEVEN, MASREM, MASQREM;

FROM MASBIOS IMPORT SWRITE, BLINES;

FROM MASSTOR IMPORT LIST, SIL, BETA, LISTVAR,
                    COMP, SFIRST, LIST1, 
                    LENGTH, ADV, INV, FIRST, RED;

FROM MASERR IMPORT ERROR, severe;

FROM SACLIST IMPORT LIST2, LIST3, LSRCH;

FROM SACSET IMPORT LBIBMS;             

FROM SACD IMPORT DGCD;

FROM SACI IMPORT IDREM, ISQRT, IDQ, IMAX, ISUM, IDQR,
                 ICOMP, IDIF, IPROD, IQR, IQ;

FROM SACM IMPORT MDLCRA, MDPROD, MIEXP, MDEXP, MDDIF, MDNEG;

CONST rcsidi = "$Id: SACPRIM.mi,v 1.4 1995/11/05 08:45:56 kredel Exp $";
CONST copyrighti = "Copyright (c) 1989 - 1992 Universitaet Passau";



PROCEDURE DPGEN(ML, K: LIST): LIST;
(*Digit prime generator. K and m are positive beta-integers.  
L is the list (p(1),...,p(r)) of all prime numbers p such that 
m le p lt m+2*K, with p(1) lt p(2) lt ... lt p(r).  
A local array is used.*)
VAR  P: ARRAY[1..1000] OF BOOLEAN;
     L, PL, ML1, ML2, HL, DL, SL, QL, RL: LIST;  
     i, k: INTEGER;
BEGIN
(*5*) (*initialize.*) k:=INTEGER(K); 
      IF k > 1000 THEN k:=1000;
         ERROR(severe,"in DPGEN, only 2*1000 numbers checked.");
         END;
      IF MASEVEN(ML) THEN ML1:=ML+1 ELSE ML1:=ML END; 
      HL:=LIST(2*k)-2; ML2:=ML1+HL;
      FOR i:=1 TO k DO P[i]:=TRUE END;
(*2*) (*mark proper odd multiples of dl for dl=3 and dl=6*nl+1 or
      dl=6*nl-1 with dl**2 le ml2.*) 
      DL:=3; SL:=0;
      LOOP   QL:=ML2 DIV DL;
             IF QL < DL THEN (*go to 3;*) EXIT; END;
             RL:=ML1 MOD DL;
             IF (RL+HL >= DL) OR (RL = 0) THEN
                IF RL = 0 THEN i:=1; ELSE
                   IF RL MOD 2 = 0  
                      THEN i:=INTEGER(DL-(RL DIV 2)+1); 
                      ELSE i:=INTEGER(((DL-RL) DIV 2)+1); END;
                   END;
                IF ML1 <= DL THEN i:=i+INTEGER(DL); END;
                WHILE i <= k DO P[i]:=FALSE; i:=i+INTEGER(DL); END;
                END;
             IF SL = 1 THEN DL:=DL+4; SL:=2; ELSE
                IF SL = 2 THEN DL:=DL+2; SL:=1; ELSE DL:=5; SL:=2;
                   END;
                END;
             END; (*loop*)
(*3*) (*construct prime list.*) L:=SIL; PL:=ML2; i:=k;
      REPEAT IF P[i] THEN L:=COMP(PL,L); END;
             PL:=PL-2; i:=i-1;
             UNTIL i = 0;
      IF ML = 1 THEN SFIRST(L,2); END;
      IF ML = 2 THEN L:=COMP(2,L); END;
      RETURN(L);
(*6*) END DPGEN;


PROCEDURE FRESL(NL: LIST; VAR ML,L: LIST);
(*Fermat residue list.  n is a positive integer with no prime divisors
less than 17.  m is a positive beta-integer and L is an ordered list
of the elements of Z(m) such that if x**2-n is a square then x is
congruent to a (modulo m) for some a in L.*)
VAR  AL1, AL2, AL3, AL4, BL1, H, HL, J1Y, J2Y, KL, KL1, L1, M, ML1:
     LIST;
BEGIN
(*1*) (*modulus 2**5.*) AL1:=IDREM(NL,32); AL2:=MASREM(AL1,16);
      AL3:=MASREM(AL2,8); AL4:=MASREM(AL3,4);
      IF AL4 = 3 THEN ML:=4;
         IF AL3 = 3 THEN BL1:=2; ELSE BL1:=0; END;
         ELSE
         IF AL3 = 1 THEN ML:=8;
            IF AL2 = 1 THEN BL1:=1; ELSE BL1:=3; END;
            ELSE ML:=16;  
            CASE INTEGER(AL1 DIV 8) OF
                 0 : BL1:=3; |
                 1 : BL1:=7; |
                 2 : BL1:=5; |
                 3 : BL1:=1; 
                 END;
            END;
         END;
      IF ML = 4 THEN L:=LIST1(BL1); ELSE J1Y:=ML-BL1;
         L:=LIST2(BL1,J1Y); END;
      KL:=LENGTH(L);
(*2*) (*modulus 3**3.*) AL1:=IDREM(NL,27); AL2:=MASREM(AL1,3);
      IF AL2 = 2 THEN ML1:=3; KL1:=1; L1:=LIST1(0); ELSE ML1:=27;
         KL1:=4; L1:=FRLSM(ML1,AL1); END;
(*3*) (*combine.*) L:=MDLCRA(ML,ML1,L,L1); ML:=ML*ML1; KL:=KL*KL1;
(*4*) (*modulus 5**2.*) AL1:=IDREM(NL,25); AL2:=MASREM(AL1,5);
      IF (AL2 = 2) OR (AL2 = 3) THEN ML1:=5; J1Y:=AL2-1; J2Y:=6-AL2;
         L1:=LIST2(J1Y,J2Y); KL1:=2; ELSE ML1:=25; L1:=FRLSM(ML1,AL1);
         KL1:=7; END;
(*5*) (*combine.*)
      IF ML1 >= (BETA DIV ML) THEN RETURN; END;
      KL:=KL*KL1; L:=MDLCRA(ML,ML1,L,L1); ML:=ML*ML1;
(*6*) (*moduli 7,11,13.*) M:=LIST3(7,11,13); H:=LIST3(64,48,0);
      LOOP   ADV(M,ML1,M);
             IF ML1 >= (BETA DIV ML) THEN RETURN; END;
             AL1:=IDREM(NL,ML1); L1:=FRLSM(ML1,AL1); KL1:=LENGTH(L1);
             L:=MDLCRA(ML,ML1,L,L1); ML:=ML*ML1; KL:=KL*KL1;
             ADV(H,HL,H);
             IF KL > HL THEN RETURN; END;
             END;
(*9*) RETURN; END FRESL;


PROCEDURE FRLSM(ML,AL: LIST): LIST;
(*Fermat residue list, single modulus.  m is a positive beta-integer.
a belongs to Z(m).  L is a list of the distinct b in Z(m) such
that b**2-a is a square in Z(m).*)
VAR  IL, ILP, JL, L, MLP, S, SL, SLP, SP: LIST;
BEGIN
(*1*) (*compute list of squares of z(ml).*) MLP:=ML DIV 2; S:=BETA;
      FOR IL:=0 TO MLP DO SL:=MDPROD(ML,IL,IL); S:=COMP(SL,S); END;
(*2*) (*compute l.*) L:=BETA; SP:=S;
      FOR IL:=MLP TO 0 BY -1 DO ADV(SP,SL,SP); SLP:=MDDIF(ML,SL,AL);
          JL:=LSRCH(SLP,S);
          IF JL <> 0 THEN L:=COMP(IL,L); ILP:=MDNEG(ML,IL);
             IF ILP <> IL THEN L:=COMP(ILP,L); END;
             END;
          END;
      RETURN(L);
(*5*) END FRLSM;


PROCEDURE GDPGEN(ML: LIST; KL: INTEGER): LIST;
(*Gaussian digit prime generator. m and k are positive beta-integers. 
L is the list (p(1),...,p(r)) of all prime numbers p such that m is 
less than or equal to p, p is less than m+4*k and p is congruent to 
3 mod 4, with p(1) lt p(2) lt ... lt p(r).  A local array is used.*)
VAR  P: ARRAY[1..1000] OF BOOLEAN;
     RL4, RL6, QLB, MLB, JL, MLS,  
     L, PL, J1Y, ML1, ML2, HL, DL, SL, QL, RL: LIST;  
     i: INTEGER;
BEGIN
(*5*) (*initialize.*) RL:=MASREM(ML,4); J1Y:=ML+3; ML1:=J1Y-RL; 
      J1Y:=LIST(4*KL);
      J1Y:=ML1+J1Y; ML2:=J1Y-4;
      IF KL > 1000 THEN KL:=1000;
         ERROR(severe,"in GDPGEN, only 4*1000 numbers checked.");
         END;
      FOR i:=1 TO KL DO P[i]:=TRUE; END;
(*2*) (*mark proper multiples of dl=3, dl=6*nl+1 and dl=6*nl-1 with
            dl**2 le ml2.*) DL:=3; RL4:=3; RL6:=3;
      LOOP   QLB:=ML2 DIV DL;
             IF DL > QLB THEN (*go to 3;*) EXIT; END;
             MLB:=QLB*DL; RL:=MASREM(MLB,4);
             IF RL = 3 THEN JL:=0; ELSE
                IF RL4 = 1 THEN JL:=RL+1; ELSE JL:=3-RL; END;
                END;
             J1Y:=JL*DL; MLS:=MLB-J1Y; J1Y:=MLS-ML1; J1Y:=J1Y DIV 4;
             i:=INTEGER(J1Y)+1; QL:=QLB-JL;
             WHILE (i > 0) AND (QL > 1) DO P[i]:=FALSE; i:=i-INTEGER(DL);
                   QL:=QL-4; END;
             IF RL6 = 1 THEN DL:=DL+4; RL6:=5; ELSE DL:=DL+2;
                IF RL4 = 3 THEN RL4:=1; ELSE RL4:=3; END;
                IF RL6 = 5 THEN RL6:=1; ELSE RL6:=RL6+2; END;
                END;
             END;
(*3*) (*construct prime list.*) L:=SIL; PL:=ML2;
      FOR i:=KL TO 1 BY -1 DO
          IF P[i] THEN L:=COMP(PL,L); END;
          PL:=PL-4; END;
      RETURN(L);
(*6*) END GDPGEN;


PROCEDURE IFACT(NL: LIST): LIST;
(*Integer factorization.  n is a positive integer. F is a list 
(q(1), q(2),...,q(h)) of the prime factors of n, q(1) le q(2) le ...
le q(h), with n equal to the product of the q(i).*)
VAR  AL, BL, CL, F, FP, J1Y, ML, MLP, PL, RL, SL, TL: LIST;
BEGIN
(*1*) (*find small factors of nl.*) ISPD(NL,F,ML);
      IF ML = 1 THEN RETURN(F); END;
      F:=INV(F); AL:=1000;
REPEAT
(*2*) (*test for primality.*)
      IF ML < BETA THEN MLP:=ML-1; RL:=MDEXP(ML,3,MLP); ELSE
         MLP:=IDIF(ML,1); RL:=MIEXP(ML,3,MLP); END;
      IF RL = 1 THEN (*go to 5;*) 
         (*5*) (*selfridge primality test.*) FP:=IFACT(MLP);
         SL:=ISPT(ML,MLP,FP);
         IF SL = 1 THEN F:=COMP(ML,F); F:=INV(F); RETURN(F); END;
         (*go to 3;*)         
         END;
(*3*) (*search for a medium divisor.*) ISQRT(ML,CL,TL); J1Y:=IDQ(CL,3);
      BL:=IMAX(5000,J1Y);
      IF ICOMP(AL,BL) > 0 THEN (*go to 4;*) PL:=1; 
         ELSE
         IMPDS(ML,AL,BL,PL,ML);
         IF PL <> 1 THEN AL:=PL; F:=COMP(PL,F); (*go to 2;*) END;
         END;
UNTIL PL = 1;
      AL:=BL;
(*4*) (*search for large divisor.*) BL:=CL; ILPDS(ML,AL,BL,PL,ML);
      IF PL <> 1 THEN F:=COMP(PL,F); END;
      F:=COMP(ML,F); F:=INV(F); RETURN(F);
(*8*) END IFACT;


PROCEDURE ILPDS(NL,AL,BL: LIST; VAR PL,NLP: LIST);
(*Integer large prime divisor search.  n is a positive integer with
no prime divisors less than 17.  1 le a le b le n.  A search is made
for a divisor p of the integer n, with a le p le b.  If such a p
is found then np=n/p, otherwise p=1 and np=n.  A modular version
of fermats method is used, and the search goes from a to b.*)
VAR  J1Y, L, LP, ML, QL, RL, RL1, RL2, SL, TL, XL, XL1, XL2, YL,
     YLP: LIST;
BEGIN
(*1*) (*compute boundary values of xl.*) IQR(NL,BL,QL,RL);
      XL1:=ISUM(BL,QL); IDQR(XL1,2,XL1,SL);
      IF (RL <> 0) OR (SL <> 0) THEN XL1:=ISUM(XL1,1); END;
      QL:=IQ(NL,AL); XL2:=ISUM(AL,QL); XL2:=IDQ(XL2,2);
(*2*) (*compute and sort residue list.*) FRESL(NL,ML,L); L:=LBIBMS(L);
      L:=INV(L);
(*3*) (*find starting residue.*) RL:=IDREM(XL2,ML); LP:=L;
      WHILE (LP <> SIL) AND (RL < FIRST(LP)) DO LP:=RED(LP); END;
      IF LP = SIL THEN LP:=L; SL:=ML; ELSE SL:=0; END;
      ADV(LP,RL1,LP); J1Y:=SL+RL; SL:=J1Y-RL1; XL:=IDIF(XL2,SL);
(*4*) (*test successive values of xl.*)
      WHILE ICOMP(XL,XL1) >= 0 DO J1Y:=IPROD(XL,XL);
            YLP:=IDIF(J1Y,NL); ISQRT(YLP,YL,TL);
            IF TL = 0 THEN PL:=IDIF(XL,YL); NLP:=ISUM(XL,YL); RETURN;
               END;
            IF LP <> SIL THEN ADV(LP,RL2,LP); SL:=RL1-RL2; ELSE
               ADV(L,RL2,LP); J1Y:=ML+RL1; SL:=J1Y-RL2; END;
            RL1:=RL2; XL:=IDIF(XL,SL); END;
(*5*) (*no divisor found.*) PL:=1; NLP:=NL; RETURN;
(*8*) END ILPDS;


PROCEDURE IMPDS(NL,AL,BL: LIST; VAR PL,QL: LIST);
(*Integer medium prime divisor search.  n, a and b are positive
integers such that a le b le n and n has no
positive divisors less than a.  If n has a prime
divisor in the closed interval from a to b then p is the least
such prime and q=n/p.  Otherwise p=1 and q=n.*)
VAR  J1Y, LP, PLP, RL, RL1, RL2: LIST;
BEGIN
(*1*) (*determine first divisor.*) RL:=IDREM(AL,210); LP:=UZ210;
      WHILE RL > FIRST(LP) DO LP:=RED(LP); END;
      RL1:=FIRST(LP); J1Y:=RL1-RL; PL:=ISUM(AL,J1Y);
(*2*) (*repeated trial division.*)
      WHILE ICOMP(PL,BL) <= 0 DO
            IF PL < BETA THEN IDQR(NL,PL,QL,RL); ELSE
               IQR(NL,PL,QL,RL); END;
            IF RL = 0 THEN RETURN; END;
            LP:=RED(LP);
            IF LP = SIL THEN LP:=UZ210; RL2:=RL1-210; ELSE
               RL2:=RL1; END;
            RL1:=FIRST(LP);
            IF PL < BETA THEN J1Y:=PL+RL1; PLP:=J1Y-RL2;
               IF PLP >= BETA THEN J1Y:=PLP-BETA; PL:=LIST2(J1Y,1);
                  ELSE PL:=PLP; END;
               ELSE J1Y:=RL1-RL2; PL:=ISUM(PL,J1Y); END;
            END;
(*3*) (*no divisors found.*) PL:=1; QL:=NL; RETURN;
(*6*) END IMPDS;


PROCEDURE ISPD(NL: LIST; VAR F,ML: LIST);
(*Integer small prime divisors.  n is a positive integer.
F is a list of primes (q(1),q(2),...,q(h)), h non-negative,
q(1) le q(2) le ... lt q(h), such that n is equal to m times the
product of the q(i) and m is not divisible by any prime in SMPRM.
Either m=1 or m gt 1,000,000.*)
VAR  LP, PL, QL, RL: LIST;
BEGIN
(*1*) F:=BETA; ML:=NL; LP:=SMPRM;
      REPEAT PL:=FIRST(LP);
             IF ML < BETA THEN MASQREM(ML,PL,QL,RL); ELSE
                IDQR(ML,PL,QL,RL); END;
             IF RL = 0 THEN F:=COMP(PL,F); ML:=QL; ELSE LP:=RED(LP);
                END;
             UNTIL (QL <= PL) OR (LP = SIL);
      IF (QL <= PL) AND (ML <> 1) THEN F:=COMP(ML,F); ML:=1; END;
      F:=INV(F); RETURN;
(*4*) END ISPD;


PROCEDURE ISPT(ML,MLP,F: LIST): LIST;
(*Integer selfridge primality test.  m is an integer greater than or
equal to 3.  mp=m-1.  F is a list (q(1),q(2),...,q(k)),
q(1) le q(2) le ... le q(k), of the prime factors of mp, with
mp equal to the product of the q(i). An attempt is made to find a 
root of unity modulo m of order m-1.  If the existence of such a root 
is discovered then m is prime and s=1.  If it is discovered that no such
root exists then m is not a prime and s=-1.  Otherwise the primality
of m remains uncertain and s=0.*)
VAR  AL, BL, FP, MLPP, PL, PL1, PP, QL, QL1, SL: LIST;
BEGIN
(*1*) (*initialize outer loop.*) FP:=F; QL1:=1; PL1:=1;
      LOOP
(*2*) (*get next divisor of mlp, if any.*)
      REPEAT IF FP = SIL THEN SL:=1; RETURN(SL); END;
             ADV(FP,QL,FP);
             UNTIL ICOMP(QL,QL1) > 0;
      QL1:=QL;
(*3*) (*try successive small primes.*) PP:=SMPRM;
      REPEAT IF PP = SIL THEN SL:=0; RETURN(SL); END;
             ADV(PP,PL,PP);
             IF PL > PL1 THEN PL1:=PL; AL:=MIEXP(ML,PL,MLP);
                IF AL <> 1 THEN SL:=-1; RETURN(SL); END;
                END;
             MLPP:=IQ(MLP,QL); BL:=MIEXP(ML,PL,MLPP);
             UNTIL BL <> 1;
(*4*) (*return for next divisor.*) (*go to 2;*)
      END;
(*7*) RETURN(SL); END ISPT;


(* initialization of SACPRIM *)

VAR IL: LIST;

BEGIN

LISTVAR(SMPRM); LISTVAR(UZ210);

(*compute small primes list.*) SMPRM:=DPGEN(1,500);

(*compute units of z sub 210.*) UZ210:=SIL;
FOR IL:=209 TO 1 BY -2 DO
    IF DGCD(210,IL) = 1 THEN
       UZ210:=COMP(IL,UZ210); END;
       END;

END SACPRIM.


(* -EOF- *)