[seqfan] Re: How to generate A000047?

[Giovanni Resta]

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;
}
/*==================================================*/