Phil Carmody on Tue, 13 Nov 2007 01:23:21 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: wishlist: evaluating polcyclo() and others |
--- Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr> wrote: > [ Posted on behalf of Joerg, who is being persecuted by qsecretary -- K.B.] I find that qsecretary is more often persecuted by yahoo's spam filters! > ----- Forwarded message from Joerg Arndt <arndt@jjj.de> ----- > For the Chebyshev polynomials (and all types > of linear recurrences) use the function frec() from: > http://www.jjj.de/pari/fastrec.inc.gp > > Both Chebyshev T and U can also be computed via > http://www.jjj.de/pari/cheby.inc.gp > Slightly faster methods are described on > pp.648-649 of the fxtbook > ( http://www.jjj.de/fxt/#fxtbook ) If we're contributing simple polynomial algorithms, then here's what I use for generalised Lucas sequences (and primitive parts thereof). It's heavily based on a seed fed to me by Dr. David Broadhurst, and has evolved to a state of 'local optimality'. --- 8< ------ uv(p,q,n)=2*lift(Mod((p+x)/2,x^2+4*q-p^2)^n) u(p,q,n)=2*polcoeff(lift(Mod((p+x)/2,x^2+4*q-p^2)^n),1) v(p,q,n)=2*polcoeff(lift(Mod((p+x)/2,x^2+4*q-p^2)^n),0) puv(p,q,n)=local(dn,nl,uvt,put=1,pvt=1,mt);dn=divisors(n);nl=length(dn);for(k=1, nl,uvt=uv(p,q,n/dn[k]);mt=moebius(dn[k]);put*=polcoeff(uvt,1)^mt;if(dn[k]%2==1,p vt*=polcoeff(uvt,0)^mt));[put,pvt] pu(p,q,n)=local(dn,nl);dn=divisors(n);nl=length(dn);prod(k=1,nl,u(p,q,n/dn[k])^m oebius(dn[k])) pv(p,q,n)=local(dn,nl);dn=divisors(n);nl=length(dn);prod(k=1,nl,if(dn[k]%2==1,v( p,q,n/dn[k])^moebius(dn[k]),1)) ------------- The '/2' perturbs me a little, but back when it was evolving, I couldn't get anything without it as quick, as the numbers grew too large (my n's were big). The primitive parts might easily be highly sub-optimal, there's almost certainly some Nash-Mihailescu tree evaluation method that's faster. Hmmm, changing subject entirely - chinese remainders! Here's something that I threw together a long time ago: --- 8< -------- naivemultichin(ls) = local(ans);ans=Mod(0,1);for(k=1,length(ls),ans=chinese(ans,ls[k]));ans multichin2(ls, s, l) = local(ans); ans=Mod(0,1);for(k=s,s+l-1,ans=chinese(ans,ls [k]));ans recchin(ls, s, l) = if(l>16,chinese(recchin(ls,s,l\2),recchin(ls,s+l\2,(l+1)\2)) ,multichin2(ls,s,l)) --------------- recchin is a simple divide and conquer which takes a vector of Mods, a start position in that vector, and a length. It relies on lazy (fast) vector copies. If vector copies aren't fast, then it needs a re-write. However, it's way faster than a loop of chineses (naivemultichin in the above). I have a feeling that Bill's latest lambda additions might mean that there are simpler ways to perform the above, for my shame I've not been keeping up with the bleeding edge. Everything in this mail is hereby placed in the public domain. Name changes of all functions heartily supported! :-) Phil () ASCII ribbon campaign () Hopeless ribbon campaign /\ against HTML mail /\ against gratuitous bloodshed [stolen with permission from Daniel B. Cristofani] ____________________________________________________________________________________ Never miss a thing. Make Yahoo your home page. http://www.yahoo.com/r/hs