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