andydude Wrote:Gottfried Wrote:If you want, I post the pari/gp routine here for re-implementation/improving. It's quick&dirty, fast only on some range of parameters (seems to cover the slow range of the wikipedia-lambert, for instance) and I didn't investigate the math below the surface.
It seems like you have all the best code! Could you post your pari/gp code for the matrix-exp and LambertW? I also remember that you posted some other pari/gp code somewhere, but I can't find the post...
Andrew Robbins
Here the complete findfixpt-code, where the lambert-w-wikipedia method is selected for some range of b
[update] I also should mention, that I work with default precision of 200 digits in pari/gp. So if you implement, for instance, the matrix-exponential using a lower precision the results may not be as expected...
To check the appropriateness of the parameter-range-separation you should examine the three alternative fixpoint-searcher again. I don't have my earlier ad-hoc-checks at hand and just left this for future investigation, while I'm using mostly "simple" bases b.
[/update]
Code:
\\ finds a fixpoint for (real) base-parameter b.
\\ Depending on its value the best of three methods are chosen
{ findfixpt(b, prec=1e-40, maxiter=1000) = local(r);
if(b>=exp(-exp(1))&b<=exp(exp(-1)),return(LH(b))); \\ 1) use "LH" for some b
if(b<4,return(findfp(b,maxiter,prec))); \\ 2) use "findfp" for some other b
return(findfpLog(b,prec,maxiter)); \\ 3) use "findfplog" for all other b
}
\\ 3) --------------------------------------------------------------
{ findfpLog(b, prec=1e-40, maxiter=1000) = local(lb=1/log(b), x2, x1=b+0.1);
\\ using hint by Henryk
for(k=1,maxiter,
x2=log(x1)*lb;
if(abs(x2-x1)<prec,return(x2));
x1=x2
);
return(0); }
\\ 2) LH(x) - uses wikipedia Lambert-W -----------------------------------
{ LW(x, prec=1E-80, maxiters=200) = local(w, we, w1e);
\\ implementation by translation of wikipedia-python source
w=0;
for(i=1,maxiters
,we=w*exp(w);
w1e=(w+1)*exp(w);
if(prec>abs((x-we)/w1e),return(w),w=w-(we-x)/(w1e-(w+2)*(we-x)/(2*w+2)))
);
print("W doesn't converge fast enough for ",x,"( we=",we);
return(0); }
LH(t11, prec=1E-80, maxiters=200) = exp(-LW(-log(t11),prec,maxiters))
\\ 1) findfp - "Lasso-method" uses fdiff1 ------------------------------------------------
\\ fdiff1 - uses vmean
vmean(v, n) = return(sum(k=1,n,v[k])/n);
\\ iterates a function tmp=(base^tmp + a*tmp)/(a+1) which circles around the fixpoint
\\ "a" is an empirical parameter to enhance convergence, here set to optimal(?) value 2
\\ performs one "round" of the spiral, returns its center
\\ workarrays vectors "list" and "diff" are provided by calling function, so they
\\ need not be created at each function-call
{ fdiff1(base=2.0, start=1+I, a=2.0, list, diff) = local(tmp1, tmp, findex);
tmp1=start;findex=1;
for(k=1,length(list), \\ compute iterations of function(see below) for some length
list[k]=tmp1;
\\ is the current iteration near the first value at k=1 (one round closed)?
diff[k]=abs(tmp1-list[1]);
if(k>1
,if(findex==1
,if(diff[k]<diff[k-1],findex=-1)
,if(diff[k]>diff[k-1],return([k,tmp1,vmean(list,k)])))); \\ return center of round
tmp=tmp1;
tmp1=(base^tmp+a*tmp)/(a+1); \\ this seems to spiral around the fixpoint
);
return([0,0,0]) ; }
{ findfp(s1, maxiter=60, prec=1e-44) = local(diff, list, res, val, oldval);
diff=vectorv(100);list=vectorv(100);
val=0.824678546142+1.56743212385*I; \\ an empirically useful init-point. May be enhanced in another version
oldval=val;
for(k=1,maxiter,
res=fdiff1(s1,val,,list,diff);
val=res[3]; \\ use the center-point as new guess
if(abs(s1^val-val)<prec,return(val));
oldval=val;
);
return(val); }
\\ ----------------------------------------------------------The optimized matrixexponential simply uses a cache, but no improvement over the "naive" exponential-series is made. (to learn about problems see Moler/Van Loan: "19 dubious ways to compute the matrix-exponential")
Code:
\\ performs matrix-exponential using matrix-cache
\\ optimized "naive" implementation of exp-series for matrices
\\ uses a matrix-cache for keeping an array of "hmatrices"-count powers of the argument.
\\ then
\\ exp(X) = I + X/1! + X^2/2! + X^3/3! + ...
\\ is replaced by, for instance cache-size "hmatrices= 3"
\\ cache = [I,X,X^2], X3 = X^3
\\ exp(X) = I *(cache[1]/0! + cache[2]/1! + cache[3]/2!)
\\ + X3 *(cache[1]/3! + cache[2]/4! + cache[3]/5!)
\\ + X3^2*(cache[1]/6! + cache[2]/7! + cache[3]/8!)
\\ + ...
\\ where also the powers of X3 are sequentially computed "on the fly"
\\ Cache should have size at least >8 to benefit from optimization, while many terms are requested
\\ uses MAsum(X) = sum of absolute values of entries in X as terminating-criterion
\\ for nilpotent matrices provide the parameters terms=0, prec=0
{ MExpOpt(m, terms=0, prec=1e-32, hmatrices=8) = local(ff1, ff2, ff3, hm, h, s, k, rs, cs);
rs=rows(m);cs=cols(m);if(rs<>cs,print("MExpOpt: no square matrix");return(m));
hm=vector(hmatrices);
ff1=matid(rs);
for(h=0,hmatrices-1,
hm[1+h]=ff1;ff1=ff1*m);
if(terms==0,terms=rs);
ff2=matid(rs);
s=sum(h=0,hmatrices-1,hm[1+h]/h!);
forstep(k=hmatrices,terms,hmatrices,
ff2=ff2*ff1;
ff3=ff2*sum(h=0,hmatrices-1,hm[1+h]/(k+h)!);
s=s+ff3;
if(prec<>0,if(MASum(ff3)<prec,print(k," terms used");break()));
);
return(s); }
Gottfried Helms, Kassel