A078843: A better (best) PARI code

[Max Alekseyev]

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
>