A078843: A better (best) PARI code

[zak seidov]

Max,

My lastthree  values agree with yours!
Only times are different:
Mmca code spent about 6 h (?!) -
much more than your 5 min(?!)

{48,2650785552}
{49,4045250962}
{50,6171386419}

I'll try to use your PARI code,
hopefully it'l work.

Thanks a lot,
zak

--- Max Alekseyev <maxale at gmail.com> wrote:

> Zak,
> 
> This is a straightforward PARI/GP implementation of
> AlmostPrimePi function:
> 
> { appi(k,n,m=2) =
> if(k==0,return(1));
> if(k==1,return(primepi(n)));
> local(r=0);
> forprime(p=m, ceil(sqrtn(n,k)),
> r+=appi(k-1,n\p,p)-(k==2)*(primepi(p)-1));
> r }
> 
> appi(k,n) gives the number of k-almostprimes not
> exceeding n.
> 
> I also suggest the following incremental formula for
> A078843:
> 
> A078843(n) = A078843(n-1) + appi3(n-k,[3^n/2^k])
> 
> where k = ceil(n*c) with c = log(5/3)/log(5/2) =
> 0.55749295065024006729857073190835923443
> and appi3(k,n) is the number of k-almostprimes not
> divisible by 3 and
> not exceeding n.
> This is its PARI/GP implementation:
> 
> { appi3(k,n) = appi(k,n) - if(k>=1,appi(k-1,n\3)) }
> 
> The following code computes and prints out first 50
> elements of
> A078843, using the above incremental formula:
> 
> ? a=1; for(n=1,50, k=ceil(n*c);
> a+=appi3(n-k,3^n\2^k); print1(a,", "); )
> 2, 3, 5, 8, 14, 23, 39, 64, 103, 169, 269, 427, 676,
> 1065, 1669, 2628,
> 4104, 6414, 10023, 15608, 24281, 37733, 58503,
> 90616, 140187, 216625,
> 334527, 516126, 795632, 1225641, 1886570, 2901796,
> 4460359, 6851532,
> 10518476, 16138642, 24748319, 37932129, 58110457,
> 88981343, 136192537,
> 208364721, 318653143, 487128905, 744398307,
> 1137129971, 1736461477,
> 2650785552, 4045250962, 6171386419,
> 
> In particular, it confirms the values computed by
> Zak.
> 
> Computing 50 first terms took about 5 minutes here.
> As the PARI/GP
> built-in function primepi() is not that fast,
> replacing it with an
> efficient implementation (e.g., array of precomputed
> values) may speed
> up computation further and help to obtain many more
> new terms. I'd be
> interested to see how fast this approach is in
> Mathematica which has
> an efficient PrimePi[] implementation.
> 
> Regards,
> Max
> 
> On Jan 6, 2008 4:44 AM, zak seidov
> <zakseidov at yahoo.com> wrote:
> > Dear PARI gurus,
> >
> > can anyone please
> >
> > write PARI code for EWW/RGWV formulas and
> MATHEMATICA
> > code in  A078843?
> >
> > My newbie PARI level allows me only to suggest
> i.o.f.
> >
> > %o A078843 (PARI)
> > a(n)=sum(i=1,3^n,if(bigomega(i)-n,0,1))
> >
> > a shorter (faster?) one:
> >
> > a(n) = sum(i=2^n,3^n,bigomega(i))
> >
> > but for larger n's it is yet not fast.
> >
> > Thanks a lot,
> > zak
> >
> > PS
> > My current additional terms in A078843:
> > {n,A078843(n)}
> > {38,37932129}
> > {39,58110457}
> > {40,88981343}
> > {41,136192537}
> > {42,208364721}
> > {43,318653143}
> > {44,487128905}
> > {45,744398307}
> > {46,1137129971}
> >
> >
> >      
>
____________________________________________________________________________________
> > Never miss a thing.  Make Yahoo your home page.
> > http://www.yahoo.com/r/hs
> >
> 



      ____________________________________________________________________________________
Never miss a thing.  Make Yahoo your home page. 
http://www.yahoo.com/r/hs