# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/compute_modp,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ # Procedures for users # none # Procedures for use only in the rest of the diffop package: # compute_modm, reconstruct, modulus2 (only used in try_factorization in DFactor) # Procedures for use only in this file: # compute_modp, modp_table # Purpose of this file: switch to modular arithmetic, and when done # then repair the damage that was done to the tables that contain the # values of the laurent series. macro( modulus=`DEtools/diffop/modulus`, set_laurents=`DEtools/diffop/set_laurents`, description_laurent=`DEtools/diffop/description_laurent`, accuracy_laurent=`DEtools/diffop/accuracy_laurent`, value_laurent=`DEtools/diffop/value_laurent`, rem_lift=`DEtools/diffop/rem_lift` ): ####################### # mod p computation #----------------------------------------------------- ####################### diffop_modp macro( compute_modm=`DEtools/diffop/compute_modm`, reconstruct=`DEtools/diffop/reconstruct`, modulus2=`DEtools/diffop/modulus2` ): macro( compute_modp=`DEtools/diffop/compute_modp`, modp_table=`DEtools/diffop/modp_table`, reconstruct=`DEtools/diffop/reconstruct` ): compute_modp:=proc() global accuracy_laurent, description_laurent, set_laurents, value_laurent, rem_lift, modp_table; local i,procedure,ac_old,de_old,va_old; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if nargs=0 then # no modular computation accuracy_laurent:=table(modp_table[0][1]); description_laurent:=table(modp_table[0][2]); set_laurents:=modp_table[0][3]; value_laurent:=table(modp_table[0][4]); rem_lift:=table(modp_table[0][5]); return fi; userinfo(6,'diffop',`computing modulo`,modulus); modp_table[0]:=[op(op(accuracy_laurent)),op(op( description_laurent)),set_laurents,op(op(value_laurent)) ,op(op(rem_lift))]; if modp_table[modulus]<>evaln(modp_table[modulus]) then accuracy_laurent:=table(modp_table[modulus][1]); description_laurent:=table(modp_table[modulus][2]); value_laurent:=table(modp_table[modulus][4]); rem_lift:=table(modp_table[modulus][5]); ac_old:=table(modp_table[0][1]); de_old:=table(modp_table[0][2]); va_old:=table(modp_table[0][4]); for i in set_laurents minus modp_table[modulus][3] do accuracy_laurent[i]:=ac_old[i]; description_laurent[i]:=de_old[i]; value_laurent[i]:=va_old[i] od; for i in modp_table[0][5] do if rem_lift[op(i)[1]] =evaln(rem_lift[op(i)[1]]) then rem_lift[op(i)[1]]:=op(i)[2] fi od else compute_modp() fi; procedure:=args[1]; i:=traperror(procedure(args[2..nargs])); if i=lasterror then userinfo(2,'diffop',modulus,`gives error`,i); compute_modp(); return `new prime needed` else modp_table[modulus]:=[op(op(accuracy_laurent)),op(op( description_laurent)),set_laurents,op(op(value_laurent)) ,op(op(rem_lift))] fi; compute_modp(); i end: # Compute modulo prime power. compute_modm:=proc(exponent) global modulus, modulus2; local v; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; modulus:=modulus2^exponent; v:=compute_modp(args[2..nargs]); while v=`new prime needed` do modulus2:=nextprime(modulus2*5); modulus:=modulus2^exponent; v:=compute_modp(args[2..nargs]) od; modulus:=0; v end: # Input: v is a list of polynomials, with coefficients in Z/(modl Z). # Output: a corresponding list with rational number coefficients. reconstruct:=proc(v,modl,N) local w,i,t,r; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if nargs=2 then reconstruct(args,floor(evalf(sqrt(modl/5)))) # 5 instead of 2, because that way we can be almost sure that if # an answer is returned, it will be the right answer. That way we # don't compute a rightdivision in vain. elif type(v,list) then w:=NULL; for i in v do w:=w,reconstruct(i,modl,N); if w[nops([w])]='failed' then return 'failed' fi od; [w] elif indets(v,'name')<>{} then t:=[op(indets(v,'name'))]; t:=t[1]; w:=0; for i from ldegree(v,t) to degree(v,t) do r:=reconstruct(coeff(v,t,i),modl,N); if r='failed' then return r fi; w:=w+r*t^i od; w elif iratrecon(v,modl,N,N,'r','t') then eval(r/t) else 'failed' fi end: #savelib('compute_modp','compute_modm','reconstruct'):