[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/





More information about the SeqFan mailing list