# [seqfan] Re: How to generate A000047?

Harry J. Smith hjsmithh at sbcglobal.net
Sat Jan 24 00:06:26 CET 2009

```This program had a problem for n > 30. I think I have fixed it:

/*==================================================*/
/* Dimensioned to compute values up to 2^50 */
/*==================================================*/

#include<stdio.h>

int D[296000]; /* Odd primes up to 2^22 + 53, was 23001*/

/* Factorize and test */

int Test(long long n)
{
int k,c;
for (k=0; (long long)D[k]*D[k]<=n; k++) // added (long long)
{
c=0; while ((n%D[k])==0){c++; n/=D[k];}
if ((c&1)==1 && ((D[k]&7)==3||(D[k]&7)==5))
return 0;
}
if ((n&7)==3 ||(n&7)==5) return 0;
return 1;
}

int main()
{
int n,k,nD=1,p;
long long x, cntD=1,cntP=0;

/* Init prime numbers */

D[0]=3;
for (n=5; n<(1<<22); n+=2) // 22 WAS 18
{
D[nD++]=n;
for (k=0; D[k]*D[k]<=n; k++)
if ((n%D[k])==0) {nD--; break;}
}

/* compute A47 */

printf("0 1\n");
for (p=1; p<=50; p++) // 50 WAS 36
{
cntP+=cntD;
for (x=(1LL<<(p-1))+1; x<(1LL<<p); x+=2)
cntD+=Test(x);
printf("%d %lld\n",p,cntP+cntD);
}
return 0;
}
/*==================================================*/

-Harry

> -----Original Message-----
> From: seqfan-bounces at list.seqfan.eu [mailto:seqfan-bounces at list.seqfan.eu]
On
> Behalf Of Giovanni Resta
> Sent: Tuesday, January 20, 2009 7:01 AM
> To: Sequence Fanatics Discussion list
> Subject: [seqfan] Re: How to generate A000047?
>
> I used the trivial C program below to compute the values up to 2^30 in
> about 72 minutes. The only "trick" I uses was to test only odd values.
>
> Giovanni
>
> 0 1
> 1 2
> 2 3
> 3 5
> 4 8
> 5 15
> 6 26
> 7 48
> 8 87
> 9 161
> 10 299
> 11 563
> 12 1066
> 13 2030
> 14 3885
> 15 7464
> 16 14384
> 17 27779
> 18 53782
> 19 104359
> 20 202838
> 21 394860
> 22 769777
> 23 1502603
> 24 2936519
> 25 5744932
> 26 11249805
> 27 22048769
> 28 43248623
> 29 84894767
> 30 166758141
>
> /*==================================================*/
> /* Dimensioned to compute values up to 2^36 */
> /*==================================================*/
>
> #include<stdio.h>
>
> int D[23001]; /* Odd primes up to 2^18 */
>
> /* Factorize and test */
>
> int Test(long long n)
> {
>    int k,c;
>    for (k=0; D[k]*D[k]<=n; k++)
>      {
>        c=0; while ((n%D[k])==0){c++; n/=D[k];}
>        if ((c&1)==1 && ((D[k]&7)==3||(D[k]&7)==5))
>           return 0;
>      }
>    if ((n&7)==3 ||(n&7)==5) return 0;
>    return 1;
> }
>
> int main()
> {
>    int n,k,nD=1,p;
>    long long x, cntD=1,cntP=0;
>
> /* Init prime numbers */
>
>    D[0]=3;
>    for (n=5; n<(1<<18); n+=2)
>      {
>        D[nD++]=n;
>        for (k=0; D[k]*D[k]<=n; k++)
> 	if ((n%D[k])==0) {nD--; break;}
>      }
>
> /* compute A47 */
>
>    printf("0 1\n");
>    for (p=1; p<=36; p++)
>      {
>        cntP+=cntD;
>        for (x=(1LL<<(p-1))+1; x<(1LL<<p); x+=2)
> 	cntD+=Test(x);
>        printf("%d %lld\n",p,cntP+cntD);
>      }
>    return 0;
> }
> /*==================================================*/
>
>
>
>
> _______________________________________________
>
> Seqfan Mailing list - http://list.seqfan.eu/

```