######################################################################### # decimal.g # defines evalf(), Exp() and decimal numbers. (C) Jean MICHEL, may 1997 # # The main function defined in this file is a function evalf similar to # the Maple function with the same name. # # gap> evalf(1/3); # 0.3333333333 # # The precision (number of digits after the decimal point) should be set # once and for all by calling the function SetDecimalPrecision (previously # computed numbers become invalid --- and print wrong --- after a call to # SetDecimalPrecision so it should not be called in the middle of a computation) # # gap> SetDecimalPrecision(20); # gap> evalf(1/3); # 0.33333333333333333333 # # The result is a ``decimal number'', which is a fixed-precision real # number for which the operations .+,-,*,/,^,< are defined. # # gap> evalf(1/3)+1; # 1.3333333333 # gap> last^3; # 2.3703703704 # # The functions Exp and Pi() are also defined by this file. Their # code provides an example of how one can use decimal numbers. # # gap> Exp(1); # 2.7182818285 # # It is also possible to evalf() cyclotomic numbers. For this reason # this file needs also the file complex.g. # # gap> evalf(ER(2)); # 1.4142135623 # gap> evalf(E(3)); # -0.5+0.8660254038I # gap> last^3; # 1 # gap> Exp(Pi()*evalf(E(4))); # 1 # # To make evaluation of cyclotomics relatively fast, a cache of previously # computed primitive roots of unity is maintained (this cache is of course # reinitialised when SetDecimalPrecision() is called). # # There is a last function defined in this file we should mention. # IsDecimal(x) returns true iff x is a decimal number. ######################################################################### DecimalOps:=OperationsRecord("DecimalOps"); SetDecimalPrecision:=function(n) DecimalOps.cachedRootsOfOne:=[]; DecimalOps.precision:=n+1; Unbind(DecimalOps.pi); end; IsDecimal:=n->IsRec(n) and IsBound(n.operations) and n.operations=DecimalOps; evalf:=x->x; #to stop complaints until defined DecimalOps.\<:=function(a,b)return evalf(a).mantissax*b);fi; a:=evalf(a);b:=evalf(b);res:=ShallowCopy(a); res.mantissa:=DecimalOps.roundquotient(a.mantissa*b.mantissa, 10^DecimalOps.precision); return res; end; DecimalOps.\/:=function(a,b)local res; a:=evalf(a);b:=evalf(b);res:=ShallowCopy(a); res.mantissa:=DecimalOps.roundquotient(a.mantissa*10^DecimalOps.precision, b.mantissa); return res; end; DecimalOps.\^:=function(h,n)local p; p:=1; if n<0 then h:=1/h;n:=-n;fi; while n>0 do if n mod 2 <> 0 then p:=p*h;fi; h:=h*h; n:=QuoInt(n,2); od; return p; end; DecimalOps.String:=function(a)local res,as,i; res:=""; if a.mantissa<0 then Add(res,'-');fi; as:=DecimalOps.roundquotient(AbsInt(a.mantissa),10); Append(res,String(QuoInt(as,10^(DecimalOps.precision-1)))); as:=String(as mod(10^(DecimalOps.precision-1))); as:=Concatenation(".",List([1..DecimalOps.precision-1-Length(as)],x->'0'),as); i:=First([Length(as),Length(as)-1..1],x->as[x]<>'0'); as:=as{[1..i]}; if as="." then as:="";fi; Append(res,as); return String(res); end; DecimalOps.Print:=function(a)Print(String(a));end; Exp:=function(x)local res,i,p,nres; res:=-1;nres:=0;p:=evalf(1);i:=1; while res<>nres do res:=nres;nres:=nres+p;p:=p*evalf(x/i);i:=i+1; od; return res; end; Pi:=function()local n,a,i,r; if not IsBound(DecimalOps.pi) then n:=DecimalOps.precision; DecimalOps.precision:=n+1; DecimalOps.pi:=0; a:=[0,2,2,1];i:=1; repeat r:=evalf(a[1+i mod 4]/((-4)^QuoInt(i,4)*i)); DecimalOps.pi:=DecimalOps.pi+r; i:=i+1; until i mod 4<>1 and r=evalf(0); DecimalOps.precision:=n; DecimalOps.pi.mantissa:=DecimalOps.roundquotient(DecimalOps.pi.mantissa,10); fi; return DecimalOps.pi; end; # if you don't want to use complex numbers comment out lines # with a # at the end #Read("complex.g"); # evalf:=function(x)local res,N,v,zeta,zetai,i; if IsDecimal(x) then return x; elif IsInt(x) then return rec(mantissa:=x*10^DecimalOps.precision,operations:=DecimalOps); elif IsRat(x) then return evalf(Numerator(x))/evalf(Denominator(x)); elif IsCyc(x) then # N:=NofCyc(x); # v:=List(CoeffsCyc(x,N),evalf); # res:=v[1]; # for i in [2..N] do # if i=2 then # if IsBound(DecimalOps.cachedRootsOfOne[N]) then zeta:=DecimalOps.cachedRootsOfOne[N]; else zeta:=Exp(2*Pi()*Complex(0,1)/N); # DecimalOps.cachedRootsOfOne[N]:=zeta; fi; zetai:=zeta; # fi; # res:=res+v[i]*zetai; # zetai:=zeta*zetai; # od; # if GaloisCyc(x,-1)=x then res:=res.r;fi; # elif IsComplex(x) then # return Complex(evalf(x.r),evalf(x.i)); # elif IsList(x) then return List(x,evalf); else Error("does not know how to evalf(",x,")"); fi; return res; end; SetDecimalPrecision(10);