A078843: A better (best) PARI code

[Max Alekseyev]

Peter Pein has pointed out a bug in my appi() code that showed up in some cases.
This is corrected code:

{ appi(k,n,m=2) = local(r=0);
 if(k==0,return(1));
 if(k==1,return(primepi(n)));
 forprime(p=m, floor(sqrtn(n,k)+1e-20),
r+=appi(k-1,n\p,p)-(k==2)*(primepi(p)-1));
 r }

Regards,
Max

On Jan 6, 2008 6:28 AM, 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
> >
>