Artifact 850d677e1d5dc2a1f5c769debe59c5a684cbb5e78366db64451de52f7b094bff:
- Executable file
r37/packages/alg/genmod.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 6132) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/alg/genmod.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 6132) [annotate] [blame] [check-ins using]
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;