reduce-historical/r38/packages/alg/gcdchk.red

MODULE GCDCHK;   % Check for a unit gcd using modular arithmetic.

% Author: Arthur C. Norman and Mary Ann Moore.

% Modifications by: Anthony C. Hearn.

FLUID '(!*BACKTRACE LIST!-OF!-LARGE!-PRIMES MODULAR!-VALUES);

% LIST!-OF!-LARGE!-PRIMES is a list of the largest pair of adjacent
% primes that can fit in the inum range of the implementation.
% This should be set here in an implementation dependent manner.
% For the time begin, a maximum inum value of 2^23 is assumed.

LIST!-OF!-LARGE!-PRIMES := '(8388449 8388451);

SYMBOLIC PROCEDURE MONIC!-MOD!-P A; 
   IF NULL A THEN NIL 
    ELSE IF DOMAINP A THEN 1
    ELSE IF LC A = 1 THEN A 
    ELSE MULTIPLY!-BY!-CONSTANT!-MOD!-P(A,MODULAR!-RECIPROCAL LC A);
 
SMACRO PROCEDURE BEFORE!-IN!-ORDER(A,B);
%Predicate taking the value true if the polynomial
%a has a leading term which comes strictly before that
%of b in canonical order;
   NULL DOMAINP A AND (DOMAINP B OR LDEG A>LDEG B);
 
SYMBOLIC PROCEDURE UNI!-PLUS!-MOD!-P(A,B);
% Form the sum of the two univariate polynomials a and b
% working over the ground domain defined by the routines
% modular!-plus, modular!-times etc. The inputs to this
% routine are assumed to have coefficients already
% in the required domain;
  IF NULL A THEN B 
   ELSE IF NULL B THEN A 
   ELSE IF BEFORE!-IN!-ORDER(A,B)
    THEN (LT A) .+ UNI!-PLUS!-MOD!-P(RED A,B)
   ELSE IF BEFORE!-IN!-ORDER(B,A)
    THEN (LT B) .+ UNI!-PLUS!-MOD!-P(A,RED B)
   ELSE IF DOMAINP A
    THEN <<A:=MODULAR!-PLUS(A,B); IF A=0 THEN NIL ELSE A>>
   ELSE BEGIN SCALAR W; 
      W:=UNI!-PLUS!-MOD!-P(RED A,RED B);
      B:=MODULAR!-PLUS(LC A,LC B);
      IF B=0 THEN RETURN W; 
      RETURN (LPOW A .* B) .+ W 
   END; 
 
%symbolic procedure uni!-times!-mod!-p(a,b);
%   if (null a) or (null b) then nil
%   else if domainp a then multiply!-by!-constant!-mod!-p(b,a)
%   else if domainp b then multiply!-by!-constant!-mod!-p(a,b)
%   else uni!-plus!-mod!-p(
%         uni!-plus!-mod!-p(uni!-times!-mod!-p(red a,red b),
%                                 uni!-times!-term!-mod!-p(lt a,b)),
%                    uni!-times!-term!-mod!-p(lt b,red a));
 
SYMBOLIC PROCEDURE UNI!-TIMES!-TERM!-MOD!-P(TERM,B);
%Multiply the given polynomial by the given term;
    IF NULL B THEN NIL 
    ELSE IF DOMAINP B THEN <<
       B:=MODULAR!-TIMES(TC TERM,B);
       IF B=0 THEN NIL 
       ELSE (TPOW TERM .* B) .+ NIL >> 
    ELSE BEGIN SCALAR W; 
	W:=MODULAR!-TIMES(TC TERM,LC B);
	IF W=0 THEN RETURN UNI!-TIMES!-TERM!-MOD!-P(TERM,RED B);
        W:= (TVAR TERM TO (TDEG TERM+LDEG B)) .* W; 
	RETURN W .+ UNI!-TIMES!-TERM!-MOD!-P(TERM,RED B)
    END; 
 
SYMBOLIC PROCEDURE UNI!-REMAINDER!-MOD!-P(A,B);
% Remainder when a is divided by b;
    IF NULL B THEN REDERR "B=0 IN REMAINDER-MOD-P" 
    ELSE IF DOMAINP B THEN NIL
    ELSE XUNI!-REMAINDER!-MOD!-P(A,B);
 
SYMBOLIC PROCEDURE XUNI!-REMAINDER!-MOD!-P(A,B);
% Remainder when the univariate modular polynomial a is
% divided by b, given that b is non degenerate;
   IF DOMAINP A OR LDEG A < LDEG B THEN A
   ELSE BEGIN 
    SCALAR Q,W; 
    Q:=MODULAR!-QUOTIENT(MODULAR!-MINUS LC A,LC B); 
% compute -lc of quotient;
    W:= LDEG A - LDEG B; %ldeg of quotient;
    IF W=0 THEN A:=UNI!-PLUS!-MOD!-P(RED A,
      MULTIPLY!-BY!-CONSTANT!-MOD!-P(RED B,Q)) 
    ELSE 
      A:=UNI!-PLUS!-MOD!-P(RED A,UNI!-TIMES!-TERM!-MOD!-P(
            (MVAR B TO W) .* Q,RED B)); 
% the above lines of code use red a and red b because
% by construction the leading terms of the required
% answers will cancel out;
     RETURN XUNI!-REMAINDER!-MOD!-P(A,B)
   END; 
 
SYMBOLIC PROCEDURE MULTIPLY!-BY!-CONSTANT!-MOD!-P(A,N); 
% Multiply the polynomial a by the constant n
% assumes that a is univariate, and that n is coprime with
% the current modulus so that modular!-times(xxx,n) neq 0
% for all xxx;
   IF NULL A THEN NIL 
   ELSE IF N=1 THEN A 
   ELSE IF DOMAINP A THEN MODULAR!-TIMES(A,N)
   ELSE (LPOW A .* MODULAR!-TIMES(LC A,N)) .+
     MULTIPLY!-BY!-CONSTANT!-MOD!-P(RED A,N); 
 
SYMBOLIC PROCEDURE UNI!-GCD!-MOD!-P(A,B);
%Return the monic gcd of the two modular univariate
%polynomials a and b;
    IF NULL A THEN MONIC!-MOD!-P B 
    ELSE IF NULL B THEN MONIC!-MOD!-P A 
    ELSE IF DOMAINP A THEN 1
    ELSE IF DOMAINP B THEN 1
    ELSE IF LDEG A > LDEG B THEN 
      ORDERED!-UNI!-GCD!-MOD!-P(A,B)
    ELSE ORDERED!-UNI!-GCD!-MOD!-P(B,A);
 
SYMBOLIC PROCEDURE ORDERED!-UNI!-GCD!-MOD!-P(A,B);
% As above, but degr a > degr b;
   IF NULL B THEN MONIC!-MOD!-P A 
   ELSE ORDERED!-UNI!-GCD!-MOD!-P(B,UNI!-REMAINDER!-MOD!-P(A,B));
 
SYMBOLIC MACRO PROCEDURE MYERR U;
   LIST('ERRORSET,
	'LIST . 
	   MKQUOTE CAADR U .
	      FOR EACH J IN CDADR U COLLECT LIST('MKQUOTE,J),
	T,'!*BACKTRACE);

SYMBOLIC PROCEDURE MODULAR!-MULTICHECK(U,V,VAR);
   IF ERRORP (U := MYERR MODULAR!-MULTICHECK1(U,V,VAR)) THEN NIL
    ELSE CAR U;

SYMBOLIC PROCEDURE MODULAR!-MULTICHECK1(U,V,VAR);
% TRUE if a modular check tells me that U and V are coprime;
  BEGIN
    SCALAR OLDP,P,MODULAR!-VALUES,UMODP,VMODP;
    P:=LIST!-OF!-LARGE!-PRIMES;
    OLDP:=SETMOD NIL;
TRY!-NEXT!-PRIME:
    MODULAR!-VALUES:=NIL;
    IF NULL P THEN GOTO UNCERTAIN;
    SETMOD CAR P;
    P:=CDR P;
    IF NULL MODULAR!-IMAGE(LC U,VAR) OR NULL MODULAR!-IMAGE(LC V,VAR)
      THEN GO TO TRY!-NEXT!-PRIME;
    UMODP:=MODULAR!-IMAGE(U,VAR);
    VMODP:=MODULAR!-IMAGE(V,VAR);
    P := DOMAINP UNI!-GCD!-MOD!-P(UMODP,VMODP);
UNCERTAIN:
    SETMOD OLDP;
    RETURN P
  END;

SYMBOLIC PROCEDURE MODULAR!-IMAGE(P,VAR);
    IF DOMAINP P
      THEN IF NULL P THEN NIL
	    ELSE IF NOT ATOM P THEN ERROR1()
	    ELSE <<P := MODULAR!-NUMBER P; IF P=0 THEN NIL ELSE P>>
    ELSE BEGIN
      SCALAR V,X,W;
      V:=MVAR P;
      IF V=VAR THEN <<
          X:=MODULAR!-IMAGE(LC P,VAR);
          IF NULL X THEN RETURN MODULAR!-IMAGE(RED P,VAR)
          ELSE RETURN (LPOW P .* X) .+ MODULAR!-IMAGE(RED P,VAR) >>;
      X:=ATSOC(V,MODULAR!-VALUES);
      IF NULL X THEN <<
	  X:=MODULAR!-NUMBER RANDOM CAR LIST!-OF!-LARGE!-PRIMES;
          MODULAR!-VALUES:=(V . X) . MODULAR!-VALUES >>
      ELSE X:=CDR X;
      X:=MODULAR!-EXPT(X,LDEG P);
      W:=MODULAR!-IMAGE(RED P,VAR);
      V:=MODULAR!-IMAGE(LC P,VAR);
      IF NULL V THEN X:=NIL
      ELSE X:=MODULAR!-TIMES(V,X);
      IF W THEN X:=MODULAR!-PLUS(X,W);
      RETURN IF X=0 THEN NIL ELSE X
    END;

ENDMODULE;

END;