module genmod; % Modular arithmetic where the modulus may be any size.
% Authors: A. C. Norman and P. M. A. Moore, 1981.
% Modifications by: John Abbott.
% Note: when balanced_mod is on, the results here are not always in
% range. *modular2f is used to correct this. However, these routines
% should be updated to give balanced results.
fluid '(!*balanced_mod current!-modulus modulus!/2);
global '(largest!-small!-modulus);
symbolic procedure set!-general!-modulus p;
if not numberp p or p=0 then current!-modulus
else begin
scalar previous!-modulus;
previous!-modulus:=current!-modulus;
current!-modulus:=p;
modulus!/2 := p/2;
% Allow for use of small moduli where appropriate.
if p <= largest!-small!-modulus then set!-small!-modulus p;
return previous!-modulus
end;
symbolic procedure general!-modular!-plus(a,b);
begin scalar result;
result:=a+b;
if result >= current!-modulus then result:=result-current!-modulus;
return result
end;
symbolic procedure general!-modular!-difference(a,b);
begin scalar result;
result := a - b;
if result < 0 then result:=result+current!-modulus;
return result
end;
symbolic procedure general!-modular!-number a;
begin
a:=remainder(a,current!-modulus);
if a < 0 then a:=a+current!-modulus;
return a
end;
symbolic procedure general!-modular!-times(a,b);
begin scalar result;
result:=remainder(a*b,current!-modulus);
if result<0
then result := result+current!-modulus; %can this happen?
return result
end;
symbolic procedure general!-modular!-reciprocal a;
% Note this returns a positive result.
if !*balanced_mod and a<0
then general!-reciprocal!-by!-gcd(current!-modulus,
a + current!-modulus,0,1)
else general!-reciprocal!-by!-gcd(current!-modulus,a,0,1);
symbolic procedure general!-modular!-quotient(a,b);
general!-modular!-times(a,general!-modular!-reciprocal b);
symbolic procedure general!-modular!-minus a;
if a=0 then a else current!-modulus - a;
symbolic procedure general!-reciprocal!-by!-gcd(a,b,x,y);
%On input A and B should be coprime. This routine then
%finds X and Y such that A*X+B*Y=1, and returns the value Y
%on input A > B;
if b=0 then rerror(alg,8,"Invalid modular division")
else if b=1 then if y < 0 then y+current!-modulus else y
else begin scalar w;
%N.B. Invalid modular division is either:
% a) attempt to divide by zero directly
% b) modulus is not prime, and input is not
% coprime with it;
w:=quotient(a,b); %Truncated integer division;
return general!-reciprocal!-by!-gcd(b,a-b*w,y,x-y*w)
end;
% The next two functions compute the "reverse" of a binary number.
% This is the number obtained when writing down the binary expansion
% in reverse order. If 2^r divides n (but 2^(r+1) does not) then
% reverse-num(reverse-num(n)) = abs(n)/2^r. r can be computed using
% height2.
symbolic procedure reverse!-num(n);
if n = 0 then n
else if n<0 then -reverse!-num1(-n,ilog2(-n)+1)
else reverse!-num1(n,ilog2(n)+1);
global '(reverse!-num!-table!*);
reverse!-num!-table!* := mkvect 16;
putv(reverse!-num!-table!*,1,8);
putv(reverse!-num!-table!*,2,4);
putv(reverse!-num!-table!*,3,12);
putv(reverse!-num!-table!*,4,2);
putv(reverse!-num!-table!*,5,10);
putv(reverse!-num!-table!*,6,6);
putv(reverse!-num!-table!*,7,14);
putv(reverse!-num!-table!*,8,1);
putv(reverse!-num!-table!*,9,9);
putv(reverse!-num!-table!*,10,5);
putv(reverse!-num!-table!*,11,13);
putv(reverse!-num!-table!*,12,3);
putv(reverse!-num!-table!*,13,11);
putv(reverse!-num!-table!*,14,7);
putv(reverse!-num!-table!*,15,15);
symbolic procedure reverse!-num1(n,bits);
if n = 0 then 0
else if bits = 1 then n
else if bits = 2 then getv(reverse!-num!-table!*,4*n)
else if bits = 3 then getv(reverse!-num!-table!*,2*n)
else if bits = 4 then getv(reverse!-num!-table!*,n)
else begin scalar shift,qr;
shift := 2**(bits/2);
qr := divide(n,shift);
if not evenp bits then shift := shift*2;
return reverse!-num1(cdr qr,bits/2)*shift +
reverse!-num1(car qr,(bits+1)/2)
end;
% Interface to algebraic mode.
flag('(reverse!-num),'integer);
deflist('((reverse!-num rnreverse!-num!*)),'!:rn!:);
%put('fibonacci,'!:rn!:,'rnfibonacci!*);
put('reverse!-num,'number!-of!-args,1);
put('reverse!-num,'simpfn,'simpiden);
symbolic procedure rnreverse!-num!*(x);
(if fixp y then reverse!-num(y)
else !*p2f mksp(list('reverse!-num,y),1))
where y=rnfixchk x;
% Interface to algebraic mode.
put('reverse!-num, 'simpfn, 'simpreverse!-num);
symbolic procedure simpreverse!-num(u);
begin scalar arg;
if length(u) neq 1 then typerr(u,"integer");
arg := simpcar u;
if denr(arg) neq 1 or not fixp(numr(arg))
then rederr("reverse!-num: argument should be an integer");
return reverse!-num(numr(arg)) ./ 1
end;
% This is an iterative version of general!-modular!-expt.
% Its principal advantage over the (simpler) recursive implementation
% is that it avoids excessive memory consumption when both n and the
% modulus are quite large -- try primep(2^10007-1) if you don't believe
% it!
symbolic procedure general!-modular!-expt(a,n);
% Computes a**n modulo current-modulus. Uses Fermat's Little
% Theorem where appropriate for a prime modulus.
if a=0 then if n=0 then rerror(alg,101,"0^0 formed") else 0
else if n=0 then 1
else if n=1 then a
else if n>=current!-modulus-1 and primep current!-modulus
then general!-modular!-expt(a,remainder(n,current!-modulus-1))
else begin scalar x, revn;
while evenp n do <<n := n/2; a := general!-modular!-times(a,a)>>;
revn := reverse!-num n;
x := 1;
while revn>0 do
<<x := general!-modular!-times(x,x);
if not evenp(revn) then x := general!-modular!-times(x,a);
revn := revn/2>>;
return x
end;
endmodule;
end;