module facmisc; % Miscellaneous routines used from several sections.
% Authors: A. C. Norman and P. M. A. Moore, 1979.
fluid '(current!-modulus
image!-set!-modulus
modulus!/2
othervars
polyzero
% pt
save!-zset
zerovarset);
global '(largest!-small!-modulus pseudo!-primes teeny!-primes);
% (1) Investigate variables in polynomial.
symbolic procedure multivariatep(a,v);
if domainp a then nil
else if not(mvar a eq v) then t
else if multivariatep(lc a,v) then t
else multivariatep(red a,v);
symbolic procedure variables!-in!-form a;
% collect variables that occur in the form a;
variables!.in!.form(a,nil);
symbolic procedure get!.coefft!.bound(poly,degbd);
% Calculates a coefft bound for the factors of poly. This simple
% bound is that suggested by Paul Wang and Linda p. Rothschild in
% Math. Comp. Vol 29 July 75, p.940, due to Gel'fond.
% Note that for tiny polynomials the bound is forced up to be
% larger than any prime that will get used in the mod-p splitting;
max(get!-height poly * fixexpfloat sumof degbd,110);
symbolic procedure sumof degbd;
if null degbd then 0 else cdar degbd + sumof cdr degbd;
symbolic procedure fixexpfloat n;
% Compute exponential function e**n for potentially large N,
% rounding result up somewhat. Note that exp(10)=22027 or so,
% so if the basic floating point exponential function is accurate
% to 6 or so digits we are protected here against roundoff.
% This could be replaced by ceiling exp n, but is written this
% way to avoid floating point overflow.
% if n>10 then (fixexpfloat(n2)*fixexpfloat(n-n2) where n2 = n/2)
% else 1+fix exp n;
if n>10 then 22027*fixexpfloat(n-10) else ceiling exp float n;
% (2) Minor variations on ordinary algebraic operations.
symbolic procedure quotfail(a,b);
% version of quotf that fails if the division does;
if polyzerop a then polyzero
else begin scalar w;
w:=quotf(a,b);
if didntgo w then errorf list("UNEXPECTED DIVISION FAILURE",a,b)
else return w
end;
symbolic procedure quotfail1(a,b,msg);
% version of quotf that fails if the division does, and gives
% custom message;
if polyzerop a then polyzero
else begin scalar w;
w:=quotf(a,b);
if didntgo w then errorf msg
else return w
end;
% (3) pseudo-random prime numbers - small and large.
symbolic procedure set!-teeny!-primes();
begin scalar i;
i:=-1;
teeny!-primes:=mkvect 9;
putv(teeny!-primes,i:=iadd1 i,3);
putv(teeny!-primes,i:=iadd1 i,5);
putv(teeny!-primes,i:=iadd1 i,7);
putv(teeny!-primes,i:=iadd1 i,11);
putv(teeny!-primes,i:=iadd1 i,13);
putv(teeny!-primes,i:=iadd1 i,17);
putv(teeny!-primes,i:=iadd1 i,19);
putv(teeny!-primes,i:=iadd1 i,23);
putv(teeny!-primes,i:=iadd1 i,29);
putv(teeny!-primes,i:=iadd1 i,31)
end;
set!-teeny!-primes();
symbolic procedure random!-small!-prime();
begin
scalar p;
repeat <<p:=small!-random!-number(); if evenp p then p := iadd1 p>>
until primep p;
return p
end;
symbolic procedure small!-random!-number();
% Returns a smallish number from a distribution strongly favouring
% smaller numbers;
begin scalar w;
% The next lines generate a random value in the range 0 to 1000000.
w := remainder(next!-random!-number(), 1000);
w := remainder(next!-random!-number(),1000) + 1000*w;
if w < 0 then w := w + 1000000;
w:=1.0+1.5*float w/1000000.0; % 1.0 to 2.5
w:=times(w,w); % In range 1.0 to 6.25
return fix exp w; % Should be in range 3 to 518,
% < 21 about half the time;
end;
% symbolic procedure fac!-exp u;
% % Simple exp routine. Assumes that Lisp has a routine for
% % exponentiation of floats by integers. Relative accuracy 4.e-5.
% begin scalar x; integer n;
% n := fix u;
% if (x := (u - float n)) > 0.5 then <<x := x - 1.0; n := n + 1>>;
% u := ee***n;
% return u*((x+6.0)*x+12.0)/((x-6.0)*x+12.0)
% end;
symbolic procedure random!-teeny!-prime l;
% get one of the first 10 primes at random providing it is
% not in the list L or that L says we have tried them all;
if l='all or (length l = 10) then nil
else begin scalar p;
repeat
p:=getv(teeny!-primes,remainder(next!-random!-number(),10))
until not member(p,l);
return p
end;
% symbolic procedure primep n;
% Test if prime. Only for use on small integers.
% n=2 or
% (n>2 and not evenp n and primetest(n,3));
% symbolic procedure primetest(n,trial);
% if igreaterp(itimes(trial,trial),n) then t
% else if iremainder(n,trial)=0 then nil
% else primetest(n,iplus2(trial,2));
% PSEUDO-PRIMES will be a list of all composite numbers which are
% less than 2^24 and where 2926^(n-1) = 3315^(n-1) = 1 mod n.
pseudo!-primes:=mkvect 87;
begin
scalar i,l;
i:=0;
l:= '(2047 4033 33227 38503 56033
137149 145351 146611 188191 226801
252601 294409 328021 399001 410041
488881 512461 556421 597871 636641
665281 722261 742813 873181 950797
1047619 1084201 1141141 1152271 1193221
1373653 1398101 1461241 1584133 1615681
1627921 1755001 1857241 1909001 2327041
2508013 3057601 3363121 3542533 3581761
3828001 4069297 4209661 4335241 4510507
4588033 4650049 4877641 5049001 5148001
5176153 5444489 5481451 5892511 5968873
6186403 6189121 6733693 6868261 6955541
7398151 7519441 8086231 8134561 8140513
8333333 8725753 8927101 9439201 9494101
10024561 10185841 10267951 10606681 11972017
13390081 14063281 14469841 14676481 14913991
15247621 15829633 16253551);
while l do <<
putv(pseudo!-primes,i,car l);
i:=i+1;
l:=cdr l >>
end;
symbolic procedure random!-prime();
begin
% I want a random prime that is smaller than largest-small-modulus.
% I do this by generating random odd integers in the range lsm/2 to
% lsm and filtering them for primality. Prime testing is done using
% a Fermat test followed by lookup in an exception table that was
% laboriously precomputed. This process should be distinctly faster
% than trial-division testing of candidate primes, but the exception
% table is tedious to compute, so I limit lsm to 2**24 here. This is
% both the value that Cambridge Lisp can support directly, an indication
% of how large an exception table I computed using 48 hours of CPU time
% and large enough that primes selected this way will hardly ever
% be unlucky just through being too small.
scalar p,w,oldmod,lsm, lsm2;
lsm := largest!-small!-modulus;
if lsm > 2**24 then lsm := 2**24;
lsm2 := lsm/2;
% W will become 1 when P is prime;
oldmod := current!-modulus;
while not (w=1) do <<
p := remainder(next!-random!-number(), lsm);
if p < lsm2 then p := p + lsm2;
if evenp p then p := p + 1;
set!-modulus p;
w:=modular!-expt(modular!-number 2926,isub1 p);
if w=1
and (modular!-expt(modular!-number 3315,isub1 p) neq 1
or pseudo!-prime!-p p)
then w:=0>>;
set!-modulus oldmod;
return p
end;
symbolic procedure pseudo!-prime!-p n;
begin
scalar low,mid,high,v;
low:=0;
high:=87; % Size of vector of pseudo-primes;
while not (high=low) do << % Binary search in table;
mid:=iquotient(iplus2(iadd1 high,low),2);
% Mid point of (low,high);
v:=getv(pseudo!-primes,mid);
if igreaterp(v,n) then high:=isub1 mid else low:=mid >>;
return (getv(pseudo!-primes,low)=n)
end;
% (4) useful routines for vectors.
symbolic procedure form!-sum!-and!-product!-mod!-p(avec,fvec,r);
% sum over i (avec(i) * fvec(i));
begin scalar s;
s:=polyzero;
for i:=1:r do
s:=plus!-mod!-p(times!-mod!-p(getv(avec,i),getv(fvec,i)),
s);
return s
end;
symbolic procedure form!-sum!-and!-product!-mod!-m(avec,fvec,r);
% Same as above but AVEC holds alphas mod p and want to work
% mod m (m > p) so minor difference to change AVEC to AVEC mod m;
begin scalar s;
s:=polyzero;
for i:=1:r do
s:=plus!-mod!-p(times!-mod!-p(
!*f2mod !*mod2f getv(avec,i),getv(fvec,i)),s);
return s
end;
symbolic procedure reduce!-vec!-by!-one!-var!-mod!-p(v,pt,n);
% Substitute for the given variable in all elements creating a
% new vector for the result. (All arithmetic is mod p).
begin scalar newv;
newv:=mkvect n;
for i:=1:n do
putv(newv,i,evaluate!-mod!-p(getv(v,i),car pt,cdr pt));
return newv
end;
symbolic procedure make!-bivariate!-vec!-mod!-p(v,imset,var,n);
begin scalar newv;
newv:=mkvect n;
for i:=1:n do
putv(newv,i,make!-bivariate!-mod!-p(getv(v,i),imset,var));
return newv
end;
symbolic procedure times!-vector!-mod!-p(v,n);
% product of all the elements in the vector mod p;
begin scalar w;
w:=1;
for i:=1:n do w:=times!-mod!-p(getv(v,i),w);
return w
end;
symbolic procedure make!-vec!-modular!-symmetric(v,n);
% fold each elt of V which is current a modular poly in the
% range 0->(p-1) onto the symmetric range (-p/2)->(p/2);
for i:=1:n do putv(v,i,make!-modular!-symmetric getv(v,i));
% (5) Combinatorial fns used in finding values for the variables.
symbolic procedure make!-zerovarset vlist;
% vlist is a list of pairs (v . tag) where v is a variable name and
% tag is a boolean tag. The procedure splits the list into two
% according to the tags: Zerovarset is set to a list of variables
% whose tag is false and othervars contains the rest;
for each w in vlist do
if cdr w then othervars:= car w . othervars
else zerovarset:= car w . zerovarset;
symbolic procedure make!-zeroset!-list n;
% Produces a list of lists each of length n with all combinations of
% ones and zeroes;
begin scalar w;
for k:=0:n do w:=append(w,kcombns(k,n));
return w
end;
symbolic procedure kcombns(k,m);
% produces a list of all combinations of ones and zeroes with k ones
% in each;
if k=0 or k=m then begin scalar w;
if k=m then k:=1;
for i:=1:m do w:=k.w;
return list w
end
else if k=1 or k=isub1 m then <<
if k=isub1 m then k:=0;
list!-with!-one!-a(k,1 #- k,m) >>
else append(
for each x in kcombns(isub1 k,isub1 m) collect (1 . x),
for each x in kcombns(k,isub1 m) collect (0 . x) );
symbolic procedure list!-with!-one!-a(a,b,m);
% Creates list of all lists with one a and m-1 b's in;
begin scalar w,x,r;
for i:=1:isub1 m do w:=b . w;
r:=list(a . w);
for i:=1:isub1 m do <<
x:=(car w) . x; w:=cdr w;
r:=append(x,(a . w)) . r >>;
return r
end;
symbolic procedure make!-next!-zset l;
begin scalar k,w;
image!-set!-modulus:=iadd1 image!-set!-modulus;
set!-modulus image!-set!-modulus;
w:=for each ll in cdr l collect
for each n in ll collect
if n=0 then n
else <<
k:=modular!-number next!-random!-number();
while (zerop k) or (onep k) do
k:=modular!-number next!-random!-number();
if k>modulus!/2 then k:=k-current!-modulus;
k >>;
save!-zset:=nil;
return w
end;
endmodule;
end;