/* Prime-sieve part based on John Moyer's 64bit sieve program (11.December 2002 version), which is available at ftp://ftp.rsok.com/pub/source/sieve2310_64bit.c and is: */ /* Copyright 1999-2000 (C) John Moyer */ /* assume 32 bit ints */ /* mailto:jrm@rsok.com */ /* http://www.rsok.com/~jrm/ */ #include #include #include #include #include /**********************************************************************/ /* */ /* */ /* This part by Antti Karttunen. Collect statistics for certain */ /* OEIS-sequences. */ /* */ /* */ /**********************************************************************/ typedef unsigned long long int ULLI; /* 64 bits at our disposal. */ typedef unsigned long int ULI; /* Plain 32 with no frills. */ int glob_dyckness_checked_only_up_to_n = 16; #define power_of_2(i) (((ULLI)1) << (i)) /* See the section "Number Conversion" at the end of the excerpt: http://www.iki.fi/kartturi/matikka/kl10exmp.txt */ int fprint_ulli(FILE *fp,ULLI x) { int s = 0; if(x >= 10) { s = fprint_ulli(fp,(x/((ULLI)10))); } fputc(('0' + (x%((ULLI)10))),fp); return(s+1); } /* Max exp-value 63 is surely enough. Nobody will compute up to that many terms in the foreseeable future. */ #define vec65zeros { 0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0 }; #define vec128zeros { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }; /* Like above, but we have prefilled the position 1 of the vector with the first prime 2, so that the prime_found continues filling it from the position 2. (Because out starting point is 3, as to avoid handling of the unique even prime 2. */ #define vec128_with_initial_2 \ { 1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }; /* (fibo 93) = 12200160415121876738 is the last Fibonacci less than 18446744073709551616 as (fibo 94) = 19740274219868223167. */ #define COMPUTE_FIBOS_UP_TO 93 ULLI vecA000045[128] = vec128zeros; /* Fibonacci numbers. Precompute them! */ #define A000045(n) (vecA000045[n]) ULLI vecA003714[128] = vec128zeros; /* Fibbinary numbers. Compute a sample. */ void precompute_fibos_to_vector(ULLI *vec,int upto_n) { int i; vec[0] = 0; vec[1] = 1; for(i=2; i <= upto_n; i++) { vec[i] = vec[i-1]+vec[i-2]; } } int index_of_largest_fibo_present(ULLI n) { int ind = COMPUTE_FIBOS_UP_TO; /* Slightly dumb, yes. */ /* Note that when n=0, the loop ends when ind=0, as A000045(0)=0 */ while(n < A000045(ind)) { ind--; } return(ind); } ULLI A003714(ULLI orig_n) { ULLI n = orig_n; ULLI z = 0; while(n != 0) { int i = index_of_largest_fibo_present(n); n -= A000045(i); if(i > (63+2)) { /* We are interested only about the three least significant fibits, thus we can safely ignore the higher fibits, especially if they would corrupt the result because of index wrap-over: */ fprintf(stderr,"A003714("); fprint_ulli(stderr,n); fprintf(stderr,") resulted a fib-index %u > 65, ignored.\n",i); } else { z |= power_of_2(i-2); } } return(z); } /* precompute_fibos_to_vector should have been called first! */ void precompute_fibbinaries_to_vector(ULLI *vec,int upto_n) { int i; for(i=0; i <= upto_n; i++) { vec[i] = A003714((ULLI)i); } } ULLI vecA036378[65] = vec65zeros; /* Here are the 20 A-numbers you requested: 95005 --- 95024. */ /* Here are the 60 A-numbers you requested: 95052 --- 95111. */ ULLI vecA095005[65] = vec65zeros; /* # of Odious primes (A027697). */ ULLI vecA095006[65] = vec65zeros; /* # of Evil primes (A027699). */ ULLI vecA027697[128] = vec128_with_initial_2; /* Odious primes. */ ULLI vecA027699[128] = vec128zeros; /* Evil primes. */ ULLI vecA095007[65] = vec65zeros; /* Primes of the form 4k+1. */ ULLI vecA095008[65] = vec65zeros; /* Primes of the form 4k+3. */ ULLI vecA095009[65] = vec65zeros; /* Primes of the form 8k+1. */ ULLI vecA095010[65] = vec65zeros; /* Primes of the form 8k+3. */ ULLI vecA095011[65] = vec65zeros; /* Primes of the form 8k+5. */ ULLI vecA095012[65] = vec65zeros; /* Primes of the form 8k+7. */ ULLI vecA095013[65] = vec65zeros; /* Primes of the form 8k+-1. */ ULLI vecA095014[65] = vec65zeros; /* Primes of the form 8k+-3. */ ULLI vecA095015[65] = vec65zeros; /* Primes of the form 6k+1. */ ULLI vecA095016[65] = vec65zeros; /* Primes of the form 6k+5. */ ULLI vecA095017[65] = vec65zeros; /* Lesser twin primes (less than A095016). */ ULLI vecA095018[65] = vec65zeros; /* Binarily balanced primes (A066196). */ ULLI vecA095019[65] = vec65zeros; /* # Zero-bit dominant primes (A095071). */ ULLI vecA095020[65] = vec65zeros; /* # One-bit dominant primes (A095070). */ ULLI vecA095052[65] = vec65zeros; /* # Primes with one more 0- than 1-bits. */ ULLI vecA095053[65] = vec65zeros; /* # Primes with one more 1- than 0-bits. */ ULLI vecA095054[65] = vec65zeros; /* # Primes that are not zero-dominant. */ ULLI vecA095055[65] = vec65zeros; /* # Primes that are not one-bit-dominant. */ ULLI vecA095056[65] = vec65zeros; /* Primes with three 1-bits (A081091). */ ULLI vecA095057[65] = vec65zeros; /* Primes with four 1-bits (A095077). */ ULLI vecA095058[65] = vec65zeros; /* # Primes with just one 0-bit. (A095078) */ ULLI vecA095059[65] = vec65zeros; /* # Primes with just 2 0-bits. (A095079) */ ULLI vecA095060[65] = vec65zeros; /* # fibeven primes (A095080). */ ULLI vecA095080[128] = vec128_with_initial_2; /* fibeven primes. */ ULLI vecA095061[65] = vec65zeros; /* # fibodd primes (A095081). */ ULLI vecA095081[128] = vec128zeros; /* fibodd primes. */ ULLI vecA095062[65] = vec65zeros; /* # fib00 primes (A095082). */ ULLI vecA095082[128] = vec128zeros; /* fib00 primes. */ ULLI vecA095063[65] = vec65zeros; /* # fibodious primes (A095083). */ ULLI vecA095083[128] = vec128_with_initial_2; /* fibodious primes. */ ULLI vecA095064[65] = vec65zeros; /* # fibevil primes (A095084). */ ULLI vecA095084[128] = vec128zeros; /* fibevil primes. */ ULLI vecA095065[65] = vec65zeros; /* # fib000 primes (A095085). */ ULLI vecA095085[128] = vec128zeros; /* fib000 primes. */ ULLI vecA095066[65] = vec65zeros; /* # fib001 primes (A095086). */ ULLI vecA095086[128] = vec128zeros; /* fib001 primes. */ ULLI vecA095067[65] = vec65zeros; /* # fib010 primes (A095087). */ ULLI vecA095087[128] = vec128zeros; /* fib010 primes. */ ULLI vecA095068[65] = vec65zeros; /* # fib100 primes (A095088). */ ULLI vecA095088[128] = vec128zeros; /* fib100 primes. */ ULLI vecA095069[65] = vec65zeros; /* # fib101 primes (A095089). */ ULLI vecA095089[128] = vec128zeros; /* fib101 primes. */ ULLI vecA095021[65] = vec65zeros; /* Primes of the form 5k+1. */ ULLI vecA095022[65] = vec65zeros; /* Primes of the form 5k+2. */ ULLI vecA095023[65] = vec65zeros; /* Primes of the form 5k+3. */ ULLI vecA095024[65] = vec65zeros; /* Primes of the form 5k+4. */ ULLI vecA095071[128] = vec128zeros; /* 0-bit dominant primes. */ ULLI vecA095070[128] = vec128zeros; /* 1-bit dominant primes. */ ULLI vecA095072[128] = vec128zeros; /* Primes with one more 0- than 1-bits. */ ULLI vecA095073[128] = vec128zeros; /* Primes with one more 1- than 0-bits. */ ULLI vecA095074[128] = vec128_with_initial_2; /* Not 0-bit-dominant primes. */ ULLI vecA095075[128] = vec128_with_initial_2; /* Not 1-bit-dominant primes. */ ULLI vecA095077[128] = vec128zeros; /* Primes with four 1-bits. */ ULLI vecA095078[128] = vec128_with_initial_2; /* Primes with a single 0-bit. */ ULLI vecA095079[128] = vec128zeros; /* Primes with two 0-bits. */ /* A095100 is for all odd numbers of form 4k+3 whose jacobi-vector = valid Motzkin path. and A095101 is for its complement (among 4k+3 integers). (compute also sum of their diving indices, and its average for each ]2^n,2^(n+1)] range. Moments, etc.) A095090 and A095091 are for the respective counts. */ ULLI vecA095092[65] = vec65zeros; /* Number of A095102 primes. */ ULLI vecA095102[128] = vec128zeros; /* 4k+3 pr. w/ Legendre-vec = Dyck path */ ULLI vecA095093[65] = vec65zeros; /* Number of A095103 primes. */ ULLI vecA095103[128] = vec128zeros; /* 4k+3 pr. w/ Legendre-vec != Dyck path */ ULLI vecA095094[65] = vec65zeros; /* Number of A080114 primes. */ ULLI vecA080114[128] = vec128zeros; /* Primes with valid Dyck path. */ ULLI vecA095095[65] = vec65zeros; /* Number of A080115 primes. */ ULLI vecA080115[128] = vec128_with_initial_2; /* Primes not in A080114. */ void fprint_ULLI_vector(FILE *fp,ULLI *vec,int start,int end) { int i; for(i=start; i <= end; i++) { if(i>start) { fprintf(fp,","); } fprint_ulli(fp,*(vec+i)); } } /* Returns the number of terms printed if finished because no more fits, zero otherwise, when everything has been printed. */ int fprint_ULLI_vector_in_pieces(FILE *fp,ULLI *vec, int start,int end,int max_linelen,int islast) { int i=start; /* Number of terms printed this time. */ int pl=0; /* Print length. */ for(;;) { pl += fprint_ulli(fp,vec[i++]); if(i > end) { return(0); } /* Finished the vector */ if(islast && ((pl+1) >= max_linelen)) { return(i-start); } /* No trailing commas on %U-line */ fprintf(fp,","); pl += 1; if(pl >= max_linelen) { return(i-start); } /* Return non-zero to indicate that more terms should be printed. */ } } int ULLIvec_find_pos_of_1st_larger_than_one(ULLI *vec,int veclen) { int pos_of_1st_term_gte_2 = 1; /* For one-based seqs. only. */ while((pos_of_1st_term_gte_2 <= veclen) && (vec[pos_of_1st_term_gte_2] < 2)) { pos_of_1st_term_gte_2++; } if(pos_of_1st_term_gte_2 > veclen) { pos_of_1st_term_gte_2 = 1; } /* Not found, use 1. */ return(pos_of_1st_term_gte_2); } /* ;; Modified from prime:jacobi-symbol of module factor.scm in A. Jaffer's SLIB ;; (version slib3a1) and available at: ;; http://swiss.csail.mit.edu/~jaffer/slib_toc ;;; Solovay-Strassen Prime Test ;;; if n is prime, then J(a,n) is congruent mod n to a**((n-1)/2) ;;; (modulo p 16) is because we care only about the low order bits. ;;; The odd? tests are inline of (expt -1 ...) ;;; Here's the original: (define (prime:jacobi-symbol p q) (cond ((zero? p) 0) ((= 1 p) 1) ((odd? p) (if (odd? (quotient (* (- (modulo p 16) 1) (- q 1)) 4)) (- (prime:jacobi-symbol (modulo q p) p)) (prime:jacobi-symbol (modulo q p) p))) (else (let ((qq (modulo q 16))) (if (odd? (quotient (- (* qq qq) 1) 8)) (- (prime:jacobi-symbol (quotient p 2) q)) (prime:jacobi-symbol (quotient p 2) q)))))) ;;; Here's my slightly faster fix-version, but still quite stupid: (define (ofix:jacobi-symbol p q) (cond ((zero? p) 0) ((= 1 p) 1) ((odd? p) (if (odd? (fix:lsh (* (- (fix:and p 15) 1) (- (fix:and q 15) 1)) -2)) (- (ofix:jacobi-symbol (modulo q p) p)) (ofix:jacobi-symbol (modulo q p) p))) (else (let ((qq (fix:and q 15))) (if (odd? (fix:lsh (- (* qq qq) 1) -3)) (- (ofix:jacobi-symbol (fix:lsh p -1) q)) (ofix:jacobi-symbol (fix:lsh p -1) q)))))) */ /* Here's our C-version. The teaching: Never implement mathematician's or logician's elegant formulae directly, without thinking! Working in GF(2) with XOR & AND is much more natural for computers than working in isomorphic multiplicative group of {+1, -1} with MUL. We also convert the recursion in our formula to a simple loop. q should be always odd! */ int js_ULI(ULI p,ULI q) { int s = 0; /* 0 in bit-2 stands for +1, 1 in bit-2 for -1. */ ULI new_p; loop: if(0 == p) { return(p); } if(1 == p) { return(1-(s&2)); } /* Convert 1 in bit-1 to -1, 0 to +1. */ if(p&1) /* We have an odd p. */ { /* If both p & q are 3 mod 4, then the sign changes, otherwise stays same: */ s ^= (p&q); /* Only the bit-1 is significant, others are ignored. */ new_p = q % p; /* Could we have a simple minus here as with Euclid? */ q = p; p = new_p; goto loop; } else /* We have an even p. So (2k/q) = (2/q)*(k/q) */ { /* where (2/q) = 1 if q is +-1 mod 8 and -1 if q is +-3 mod 8. */ /* I.e. sign changes only if q's lower bits are (011) or (101), i.e. if the bit-1 and bit-2 xored yield 1. */ s ^= (q^(q>>1)); /* Thus, this does it. */ p >>= 1; goto loop; } } /* Here it is in MIT/GNU Scheme: ;; I hope the compiler is clever enough to see that it is ;; required that p and q are fixnums, and compiles nothing ;; unnecessary. (define (fix:jacobi-symbol p q) (if (not (and (fix:fixnum? p) (fix:fixnum? q) (fix:= 1 (fix:and q 1)))) (error "fix:jacobi-symbol: args must be fixnums, and 2. arg should be odd: " p q ) (let loop ((p p) (q q) (s 0)) ;; 0 in bit-2 stands for +1, 1 in bit-2 for -1. (cond ((fix:zero? p) 0) ((fix:= 1 p) (fix:- 1 (fix:and s 2))) ((fix:= 1 (fix:and p 1)) ;; Odd p ? (loop (fix:remainder q p) p (fix:xor s (fix:and p q))) ) (else ;; It's even. (loop (fix:lsh p -1) q (fix:xor s (fix:xor q (fix:lsh q -1)))) ) ) ) ) ) */ /* Returns the index i of the first point where Sum_{j=1..i} JS(i,n) goes negative, and zero if it doesn't go when checked up to i=k. */ int js_diving_index(ULI n,ULI k) { ULI i; long int s; for(s=0,i=1; i <= k; i++) { s += js_ULI(i,n); /* printf("js_diving_index(%lu,%lu): i=%lu,s=%ld\n",n,k,i,s); */ if(s < 0) { return(i); } } if((0 != s) && (i == n)) { fprintf(stderr, "js_diving_index: Sum of J(1,%lu)..J(n-1,%lu) is not zero: %ld", n,n,s); exit(1); } return(0); } int binwidth(ULLI n) { int i=0; while(0 != n) { i++; n >>= 1; } return(i); } int A000120(ULLI n) { int i=0; while(0 != n) { i += (n&1); n >>= 1; } return(i); } /* This is the characteristic function of A000069, i.e. A000120(n) mod 2 */ int A010060(ULLI n) { int i=0; while(0 != n) { i ^= (n&1); n >>= 1; } return(i); } char *w6d(char *tb,int n) { sprintf(tb,"%06u",n); return(tb); } void vec_sum_checker(FILE *fp,int veclen,int offset, int Anum_sum,ULLI *sumvec, int Anum1summand,ULLI *summand1vec, int Anum2summand,ULLI *summand2vec) { int i,matched; char tb1[81] = { 'A' }, tb2[81] = { 'A' }, tb3[81] = { 'A' }; char *Astr1summand = (w6d(tb1+1,Anum1summand)-1); char *Astr2summand = (w6d(tb2+1,Anum2summand)-1); char *Astr_sum = (w6d(tb3+1,Anum_sum)-1); for(i=offset, matched=0; i <= veclen; i++) { if((summand1vec[i] + summand2vec[i]) == sumvec[i]) { matched++; } else { fprintf(fp,"%s[%u] != (%s[%u]+%s[%u]), i.e. ", Astr_sum,i, Astr1summand,i, Astr2summand,i); fprint_ulli(fp,sumvec[i]); fprintf(fp," != ("); fprint_ulli(fp,summand1vec[i]); fprintf(fp,"+"); fprint_ulli(fp,summand2vec[i]); fprintf(fp,")\n"); } } fprintf(fp,"%s matched %s+%s in %u positions.\n", Astr_sum,Astr1summand,Astr2summand,matched); } #define OEIS_S_T_U_LINE_MAXLEN 64 void output_OEIS_sequence(FILE *fp,int Anum, ULLI *vec,int veclen, int offset, char *name, char *author_info, char *datestr, char *Y_line) { char tb[81] = { 'A' }; char *Astr = (w6d(tb+1,Anum)-1); int pos_of_1st_term_gte_2 = ULLIvec_find_pos_of_1st_larger_than_one(vec,veclen); int printed_only_up_to_nth_term = 0; int sec_printed_only_up_to_nth_term = 0; int upto_n = veclen; fprintf(fp,"%%I %s\n",Astr); /* Old way: all the terms to %S-line: { fprintf(fp,"%%S %s ",Astr); fprint_ULLI_vector(fp,vec,offset,veclen); } */ { fprintf(fp,"%%S %s ",Astr); printed_only_up_to_nth_term = fprint_ULLI_vector_in_pieces(fp,vec,offset,veclen, OEIS_S_T_U_LINE_MAXLEN,0); fprintf(fp,"\n"); } if(printed_only_up_to_nth_term > 0) /* Continue onto the %T -line ? */ { fprintf(fp,"%%T %s ",Astr); sec_printed_only_up_to_nth_term = fprint_ULLI_vector_in_pieces(fp,vec+printed_only_up_to_nth_term, offset, upto_n-printed_only_up_to_nth_term, OEIS_S_T_U_LINE_MAXLEN,0); fprintf(fp,"\n"); } if(sec_printed_only_up_to_nth_term > 0) /* Continue onto the %U -line ? */ { printed_only_up_to_nth_term += sec_printed_only_up_to_nth_term ; fprintf(fp,"%%U %s ",Astr); fprint_ULLI_vector_in_pieces(fp,vec+printed_only_up_to_nth_term,offset, upto_n-printed_only_up_to_nth_term, OEIS_S_T_U_LINE_MAXLEN,1); fprintf(fp,"\n"); } fprintf(fp,"%%N %s %s\n",Astr,name); /* Name should end with period. */ fprintf(fp,"%%Y %s %s\n",Astr,Y_line); /* fprintf(fp,"%%H %s A. Karttunen, C-program for computing the initial terms of this sequence\n", Astr); */ fprintf(fp,"%%K %s nonn\n",Astr); fprintf(fp,"%%O %s %u,%u\n",Astr,offset,pos_of_1st_term_gte_2); fprintf(fp,"%%A %s %s, %s\n", Astr, author_info, datestr); fprintf(fp,"\n"); fflush(fp); } /**********************************************************************/ /* */ /**********************************************************************/ #define STEP_CONST 128 int small_primes[343] = { 2,3,5,7,11,13,17,19,23,29, 31,37,41,43,47,53,59,61,67,71, 73,79,83,89,97,101,103,107,109,113, 127,131,137,139,149,151,157,163,167,173, 179,181,191,193,197,199,211,223,227,229, 233,239,241,251,257,263,269,271,277,281, 283,293,307,311,313,317,331,337,347,349, 353,359,367,373,379,383,389,397,401,409, 419,421,431,433,439,443,449,457,461,463, 467,479,487,491,499,503,509,521,523,541, 547,557,563,569,571,577,587,593,599,601, 607,613,617,619,631,641,643,647,653,659, 661,673,677,683,691,701,709,719,727,733, 739,743,751,757,761,769,773,787,797,809, 811,821,823,827,829,839,853,857,859,863, 877,881,883,887,907,911,919,929,937,941, 947,953,967,971,977,983,991,997,1009,1013, 1019,1021,1031,1033,1039,1049,1051,1061,1063,1069, 1087,1091,1093,1097,1103,1109,1117,1123,1129,1151, 1153,1163,1171,1181,1187,1193,1201,1213,1217,1223, 1229,1231,1237,1249,1259,1277,1279,1283,1289,1291, 1297,1301,1303,1307,1319,1321,1327,1361,1367,1373, 1381,1399,1409,1423,1427,1429,1433,1439,1447,1451, 1453,1459,1471,1481,1483,1487,1489,1493,1499,1511, 1523,1531,1543,1549,1553,1559,1567,1571,1579,1583, 1597,1601,1607,1609,1613,1619,1621,1627,1637,1657, 1663,1667,1669,1693,1697,1699,1709,1721,1723,1733, 1741,1747,1753,1759,1777,1783,1787,1789,1801,1811, 1823,1831,1847,1861,1867,1871,1873,1877,1879,1889, 1901,1907,1913,1931,1933,1949,1951,1973,1979,1987, 1993,1997,1999,2003,2011,2017,2027,2029,2039,2053, 2063,2069,2081,2083,2087,2089,2099,2111,2113,2129, 2131,2137,2141,2143,2153,2161,2179,2203,2207,2213, 2221,2237,2239,2243,2251,2267,2269,2273,2281,2287, 2293,2297,2309 }; /* expand scheme from 30 to 2310 30 == 2*3*5 and 8 == (2-1)*(3-1)*(5-1) 2*3*5*7*11 == 2310 1*2*4*6*10 == 480 bits == 60 bytes ==> prime or not prime for 38.5 integers per byte In each 2310 integers, for N >= 1, the numbers that might be prime are: N*2310+1, N*2310+13, N*2310+17, . . . N*2310+13*13, N*2310+13*17, . . . N*2310+2309 */ unsigned long maskbits[32] = { 0x00000001, 0x00000002, 0x00000004, 0x00000008, 0x00000010, 0x00000020, 0x00000040, 0x00000080, 0x00000100, 0x00000200, 0x00000400, 0x00000800, 0x00001000, 0x00002000, 0x00004000, 0x00008000, 0x00010000, 0x00020000, 0x00040000, 0x00080000, 0x00100000, 0x00200000, 0x00400000, 0x00800000, 0x01000000, 0x02000000, 0x04000000, 0x08000000, 0x10000000, 0x20000000, 0x40000000, 0x80000000 }; /* value of mask_bit_index is bit corresponding to integer modulo 2310 mask_bit_index[i] is zero if this is integer that could not possibly be prime or bit number from 1 to 480 if it could be prime. */ unsigned int mask_bit_index[2310]; /* bit is one to 480 */ void clear_one_mask_bit(unsigned long *sieve_array, int bit) { int k; k = (bit - 1) >> 5 ; /* divide by 32 */ sieve_array[k] &= ~(maskbits[(bit - 1) & 31]); /* modulo 32 */ #ifdef DEBUG printf("sieve_array[%d]=0x%08x, ~(maskbits[(%d - 1) & 31]=0x%08x\n", k, sieve_array[k], bit, ~(maskbits[(bit - 1) & 31])); #endif } /* bit is one to 480 */ int test_one_mask_bit(unsigned long *sieve_array, int bit) { int k; k = (bit - 1) >> 5 ; /* divide by 32 */ return (sieve_array[k] & (maskbits[(bit - 1) & 31])); /* modulo 32 */ } extern char *optarg; extern int optind, opterr, optopt; int getopt (int argc, char *const argv[], const char *opts); /* Max exp-value 63 is surely enough. Nobody will compute up to that many terms in the foreseeable future. */ int main(int argc, char *argv[]) { int expi; int max_primes_collected; ULLI start, stop; char *progname = argv[0]; char *s; /* One might wonder why I (AK) don't use getopt as jrm did, but the reason is that I never have... */ arg_loop: if(NULL != (s = *++argv)) { if('-' == *s) { switch(*(s+1)) { case 'd': { if(!*(s+2)) { ++argv; if(!(*argv) || !isdigit(**(argv))) { numeric_required: fprintf(stderr, "%s: Option -%c requires a numeric parameter!\n", progname,*(s+1)); goto usage; } else { glob_dyckness_checked_only_up_to_n = atoi(*argv); } } else /* Something following directly after -d */ { if(!isdigit(*(s+2))) { goto numeric_required; } glob_dyckness_checked_only_up_to_n = atoi(s+2); } goto arg_loop; } default: { fprintf(stderr,"%s: Unknown option %s !\n",progname,s); goto usage; break; } } } expi = atoi(*argv); if((expi < 1) || (expi > 63)) { usage: fprintf(stderr, "Usage: %s [-d num] exponent, where the exponent is in range [1,63]\n", progname); exit(1); } start = 3; /* stop = power_of_2(expi+1)-1; */ stop = power_of_2(expi+1)+1; /* Because we search twin primes also. */ } else { goto usage; } precompute_fibos_to_vector(vecA000045,COMPUTE_FIBOS_UP_TO); precompute_fibbinaries_to_vector(vecA003714,127); max_primes_collected = 101; compute_primes(start,stop,max_primes_collected); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95005,vecA095005,95006,vecA095006); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95007,vecA095007,95008,vecA095008); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95013,vecA095013,95014,vecA095014); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95015,vecA095015,95016,vecA095016); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95019,vecA095019,95054,vecA095054); vec_sum_checker(stdout,expi,1, 95013,vecA095013,95009,vecA095009,95012,vecA095012); vec_sum_checker(stdout,expi,1, 95014,vecA095014,95010,vecA095010,95011,vecA095011); vec_sum_checker(stdout,expi,1, 95054,vecA095054,95018,vecA095018,95020,vecA095020); vec_sum_checker(stdout,expi,1, 95055,vecA095055,95018,vecA095018,95019,vecA095019); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95063,vecA095063,95064,vecA095064); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95060,vecA095060,95061,vecA095061); vec_sum_checker(stdout,expi,1, 95060,vecA095060,95062,vecA095062,95067,vecA095067); vec_sum_checker(stdout,expi,1, 95061,vecA095061,95066,vecA095066,95069,vecA095069); vec_sum_checker(stdout,expi,1, 95062,vecA095062,95065,vecA095065,95068,vecA095068); vec_sum_checker(stdout,expi,1, 95008,vecA095008,95092,vecA095092,95093,vecA095093); vec_sum_checker(stdout,expi,1, 36378,vecA036378,95094,vecA095094,95095,vecA095095); output_OEIS_sequence(stdout, 45,vecA000045, max_primes_collected, 1, "Fibonacci numbers: F(n) = F(n-1) + F(n-2), F(0) = 0, F(1) = 1, F(2) = 1, ...", "njas", "", "HERE JUST FOR CHECKING!" ); output_OEIS_sequence(stdout, 3714,vecA003714, max_primes_collected, 1, "Fibbinary numbers", "njas", "", "HERE JUST FOR CHECKING!" ); output_OEIS_sequence(stdout, 36378,vecA036378, expi, 1, "Number of primes p such that 2^n < p < 2^(n+1).", "Labos E. (labos(AT)ana1.sote.hu)", "May 13 2004", "a(n) = A095005(n)+A095006(n) = A095007(n) + A095008(n)" " = A095013(n) + A095014(n)" " = A095015(n) + A095016(n) (for n > 1)" " = A095021(n)+A095022(n)+A095023(n)+A095024(n)" /* " = A095018(n)+A095019(n)+A095020(n)" */ " = A095019(n)+A095054(n)" " = A095020(n)+A095055(n)" " = A095060(n)+A095061(n)" " = A095063(n)+A095064(n)." ); output_OEIS_sequence(stdout, 95005,vecA095005, expi, 1, "Number of odious primes (A027697) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095006(n)." ); output_OEIS_sequence(stdout, 95006,vecA095006, expi, 1, "Number of evil primes (A027699) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095005(n)." ); output_OEIS_sequence(stdout, 27697,vecA027697,vecA027697[0], 1, "Odious primes: primes with odd number of 1's in binary expansion.", "njas", "", "Complement of A027699 in A000040." " Union of A091206\{3} and odious members of A091209." " Cf. A095005." ); output_OEIS_sequence(stdout, 27699,vecA027699,vecA027699[0], 1, "Evil primes: primes with even number of 1's in binary expansion.", "njas", "", "Complement of A027697 in A000040." " Cf. A095006." ); output_OEIS_sequence(stdout, 95007,vecA095007, expi, 1, "Number of 4k+1 primes (A002144) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095008(n) = A095009(n)+A095011(n)." ); output_OEIS_sequence(stdout, 95008,vecA095008, expi, 1, "Number of 4k+3 primes (A002145) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095007(n) = A095010(n)+A095012(n)." ); /* Cf. A091126, A091127, A091128, A091129. */ output_OEIS_sequence(stdout, 95009,vecA095009, expi, 1, "Number of 8k+1 primes (A007519) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095013(n)-A095012(n) = A095007(n)-A095011(n). " "Cf. A091126." ); output_OEIS_sequence(stdout, 95010,vecA095010, expi, 1, "Number of 8k+3 primes (A007520) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095014(n)-A095011(n) = A095008(n)-A095012(n). " "Cf. A091127." ); output_OEIS_sequence(stdout, 95011,vecA095011, expi, 1, "Number of 8k+5 primes (A007521) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095014(n)-A095010(n). Cf. A091128." ); output_OEIS_sequence(stdout, 95012,vecA095012, expi, 1, "Number of 8k+7 primes (A007522) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095013(n)-A095009(n). Cf. A091129." ); output_OEIS_sequence(stdout, 95013,vecA095013, expi, 1, "Number of 8k+-1 primes (A001132) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095014(n) = A095009(n)+A095012(n)." ); output_OEIS_sequence(stdout, 95014,vecA095014, expi, 1, "Number of 8k+-3 primes (A003629) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095013(n) = A095010(n)+A095011(n)." ); output_OEIS_sequence(stdout, 95015,vecA095015, expi, 1, "Number of 6k+1 primes (A002476) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095016(n) (apart the initial term)." ); output_OEIS_sequence(stdout, 95016,vecA095016, expi, 1, "Number of 6k+5 primes (A007528) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095015(n) (apart the initial term)." ); output_OEIS_sequence(stdout, 95017,vecA095017, expi, 1, "Number of lesser twin primes (A001359) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095016, A036378." ); output_OEIS_sequence(stdout, 95021,vecA095021, expi, 1, "Number of 5k+1 primes (A030430) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A036378." ); output_OEIS_sequence(stdout, 95022,vecA095022, expi, 1, "Number of 5k+2 primes (A030432) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A036378." ); output_OEIS_sequence(stdout, 95023,vecA095023, expi, 1, "Number of 5k+3 primes (A030431) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A036378." ); output_OEIS_sequence(stdout, 95024,vecA095024, expi, 1, "Number of 5k+4 primes (A030433) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A036378." ); output_OEIS_sequence(stdout, 95018,vecA095018, expi, 1, "Number of binarily balanced primes (A066196) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095005-A095006." ); output_OEIS_sequence(stdout, 95019,vecA095019, expi, 1, "Number of zero-bit dominant primes (A095071) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); output_OEIS_sequence(stdout, 95020,vecA095020, expi, 1, "Number of one-bit dominant primes (A095070) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); output_OEIS_sequence(stdout, 95052,vecA095052, expi, 1, "Number of primes with #0-bits = 1+#1-bits (A095072) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); output_OEIS_sequence(stdout, 95053,vecA095053, expi, 1, "Number of primes with #1-bits = 1+#0-bits (A095073) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); output_OEIS_sequence(stdout, 95054,vecA095054, expi, 1, "Number of primes with #0-bits <= #1-bits (A095074) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095018(n)+A095020(n)." ); output_OEIS_sequence(stdout, 95055,vecA095055, expi, 1, "Number of primes with #1-bits <= #0-bits (A095075) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095018(n)+A095019(n)." ); output_OEIS_sequence(stdout, 95056,vecA095056, expi, 1, "Number of primes with three 1-bits (A081091) in range ]2^n,2^(n+1)].", "Labos E. (labos(AT)ana1.sote.hu) & " "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); output_OEIS_sequence(stdout, 95057,vecA095057, expi, 1, "Number of primes with four 1-bits (A095077) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); output_OEIS_sequence(stdout, 95058,vecA095058, expi, 1, "Number of primes with a single 0-bit (A095078) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); output_OEIS_sequence(stdout, 95059,vecA095059, expi, 1, "Number of primes with two 0-bits (A095079) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095018." ); /* Then a few prime-sequences themselves: */ output_OEIS_sequence(stdout, 95070,vecA095070, vecA095070[0], /* Contains the count of collected primes. */ 1, "One-bit dominant primes, i.e. primes whose binary expansion contains" " more 1's than 0's.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Intersect of A000040 & A072600. Complement of A095075 in A000040." " Cf. A095020." ); output_OEIS_sequence(stdout, 95071,vecA095071, vecA095071[0], /* Contains the count of collected primes. */ 1, "Zero-bit dominant primes, i.e. primes whose binary expansion contains" " more 0's than 1's.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Complement of A095074 in A000040." " Subsets: A095072, ... Cf. A095019." ); output_OEIS_sequence(stdout, 95072,vecA095072, vecA095072[0], /* Contains the count of collected primes. */ 1, "Primes in whose binary expansion the number of 0-bits is one more" " than the number of 1-bits.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Intersect of A000040 & A031444. Subset of A095071." " Complement of A0950xx in A000040." "Cf. A095052." ); output_OEIS_sequence(stdout, 95073,vecA095073, vecA095073[0], /* Contains the count of collected primes. */ 1, "Primes in whose binary expansion the number of 1-bits is one more" " than the number of 0-bits.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Intersect of A000040 & A031448. Subset of A095070." "Cf. A095053." ); output_OEIS_sequence(stdout, 95074,vecA095074, vecA095074[0], /* Contains the count of collected primes. */ 1, "Primes in whose binary expansion the number of 0-bits is less than" " or equal to number of 1-bits.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Complement of A095071 in A000040." " Differs from A057447 first time at n=18, where a(n)=71," " while A057447(18)=67." " Cf. A095054." ); output_OEIS_sequence(stdout, 95075,vecA095075, vecA095075[0], /* Contains the count of collected primes. */ 1, "Primes in whose binary expansion the number of 1-bits is less than" " or equal to number of 0-bits.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", " Complement of A095070 in A000040." "Cf. A095055." ); output_OEIS_sequence(stdout, 95077,vecA095077, vecA095077[0], /* Contains the count of collected primes. */ 1, "Primes with four 1-bits in their binary expansion.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Subset of A027699. " "Differs from A085448 first time at n=19, where a(n)=337, " "while A085448 continues from there with 311, whose binary expansion " "has six 1-bits, not four. Cf. A095057." ); output_OEIS_sequence(stdout, 95078,vecA095078, vecA095078[0], /* Contains the count of collected primes. */ 1, "Primes with a single 0-bit in their binary expansion.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Intersect of A000040 & A030130. Cf. A095058." ); output_OEIS_sequence(stdout, 95079,vecA095079, vecA095079[0], /* Contains the count of collected primes. */ 1, "Primes with two 0-bits in their binary expansion.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Intersect of A000040 & A0xxxxx. Cf. A095059." ); output_OEIS_sequence(stdout, 95060,vecA095060, expi, 1, "Number of fibeven primes (A095080) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095061(n) = A095062(n)+A095067(n)." ); output_OEIS_sequence(stdout, 95080,vecA095080,vecA095080[0], 1, "Fibeven primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends with zero.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Intersect of A000040 & A022342. Union of A095082 & A095087." " Cf. A095060, A095081." ); output_OEIS_sequence(stdout, 95061,vecA095061, expi, 1, "Number of fibodd primes (A095081) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095060(n) = A095066(n)+A095069(n)." ); output_OEIS_sequence(stdout, 95081,vecA095081,vecA095081[0], 1, "Fibodd primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends with one.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Intersect of A000040 & A003622. Union of A095086 & A095089." " Cf. A095061, A095080." ); output_OEIS_sequence(stdout, 95062,vecA095062, expi, 1, "Number of fib00 primes (A095082) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095060(n)-A095067(n) = A095065(n)+A095068(n)." ); output_OEIS_sequence(stdout, 95082,vecA095082,vecA095082[0], 1, "Fib00 primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends with two zeros.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095062. Union of A095085 & A095088." ); output_OEIS_sequence(stdout, 95065,vecA095065, expi, 1, "Number of fib000 primes (A095085) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095062(n)-A095068(n). Cf. A095066-A095067." ); output_OEIS_sequence(stdout, 95085,vecA095085,vecA095085[0], 1, "Fib000 primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends with three zeros.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095065." ); output_OEIS_sequence(stdout, 95066,vecA095066, expi, 1, "Number of fib001 primes (A095086) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095061(n)-A095069(n). Cf. A095065 & A095067." ); output_OEIS_sequence(stdout, 95086,vecA095086,vecA095086[0], 1, "Fib001 primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends with two zeros and final 1.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095066." ); output_OEIS_sequence(stdout, 95067,vecA095067, expi, 1, "Number of fib010 primes (A095087) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095060(n)-A095062(n)." ); output_OEIS_sequence(stdout, 95087,vecA095087,vecA095087[0], 1, "Fib010 primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends with zero, one and zero.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095067." ); output_OEIS_sequence(stdout, 95068,vecA095068, expi, 1, "Number of fib100 primes (A095088) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095062(n)-A095065(n)." ); output_OEIS_sequence(stdout, 95088,vecA095088,vecA095088[0], 1, "Fib100 primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends with one and two final zeros.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095068." ); output_OEIS_sequence(stdout, 95069,vecA095069, expi, 1, "Number of fib101 primes (A095089) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A095061(n)-A095066(n)." ); output_OEIS_sequence(stdout, 95089,vecA095089,vecA095089[0], 1, "Fib101 primes, i.e. primes whose Zeckendorf-expansion A014417(p)" " ends as one, zero, one.", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "Cf. A095069." ); output_OEIS_sequence(stdout, 95063,vecA095063, expi, 1, "Number of fibodious primes (A095083) in range ]2^n,2^(n+1)].", "Antti Karttunen (his-firstname.his-surname(AT)iki.fi)", "May 28 2004", "a(n) = A036378(n)-A095064(n)." ); form 4k+1. */ else { vecA095008[w]++; } /* Primes of the form 4k+3. */ switch(p & 7) /* I.e. modulo 8. */ { case 1: { vecA095009[w]++; /* Primes of the form 8k+1. */ vecA095013[w]++; /* Primes of the form 8k+-1. */ break; } case 3: { vecA095010[w]++; /* Primes of the form 8k+3. */ vecA095014[w]++; /* Primes of the form 8k+-3. */ break; } case 5: { vecA095011[w]++; /* Primes of the form 8k+5. */ vecA095014[w]++; /* Primes of the form 8k+-3. */ break; } case 7: { vecA095012[w]++; /* Primes of the form 8k+5. */ vecA095013[w]++; /* Primes of the form 8k+-1. */ break; } default: { fprintf(stderr,"Prime "); fprint_ulli(stderr,p); fprintf(stderr, " with weight %d and odiousness %d,", w,is_odious); fprintf(stderr," has impossible congruence modulo 8: %d\n", ((int)(p & 7))); exit(1); } } mod6 = (p % 6); if(1 == mod6) { vecA095015[w]++; } /* Primes of the form 6k+1. */ else if(5 == mod6) { vecA095016[w]++; } /* Primes of the form 6k+5. */ switch(p % 5) { case 0: { break; } /* Five itself. Ignore. */ case 1: { vecA095021[w]++; break; } case 2: { vecA095022[w]++; break; } case 3: { vecA095023[w]++; break; } case 4: { vecA095024[w]++; break; } } if(1) /* Check the Zeckendorf-expansion as well. */ { ULLI ze = A003714(p); int is_fibodious = A010060(ze); if(is_fibodious) { vecA095063[w]++; add_p_to_vector(p,vecA095083); } else /* No, it's fibevil. */ { vecA095064[w]++; add_p_to_vector(p,vecA095084); } switch(ze & 7) /* We are interested about the three rightmost fibits. */ { case 0: /* 000 */ { vecA095065[w]++; add_p_to_vector(p,vecA095085); vecA095062[w]++; add_p_to_vector(p,vecA095082); /* 00 */ vecA095060[w]++; add_p_to_vector(p,vecA095080); /* 0 */ break; } case 1: /* 001 */ { vecA095066[w]++; add_p_to_vector(p,vecA095086); vecA095061[w]++; add_p_to_vector(p,vecA095081); /* 1 */ break; } case 2: /* 010 */ { vecA095067[w]++; add_p_to_vector(p,vecA095087); vecA095060[w]++; add_p_to_vector(p,vecA095080); /* 0 */ break; } case 4: /* 100 */ { vecA095068[w]++; add_p_to_vector(p,vecA095088); vecA095062[w]++; add_p_to_vector(p,vecA095082); /* 00 */ vecA095060[w]++; add_p_to_vector(p,vecA095080); /* 0 */ break; } case 5: /* 101 */ { vecA095069[w]++; add_p_to_vector(p,vecA095089); vecA095061[w]++; add_p_to_vector(p,vecA095081); /* 1 */ break; } default: { fprintf(stderr,"Prime "); fprint_ulli(stderr,p); fprintf(stderr, " with weight %d and odiousness %d,", w,is_odious); fprintf(stderr," has impossible Zeckendorf-expansion: "); fprint_ulli(stderr,ze); fprintf(stderr," as modulo 8 it results %u.\n", ((int)(ze&7))); exit(1); } } } if((prev_prime + 2) == p) /* Found a pair of primes. (twin primes). */ { int w2 = binwidth(prev_prime)-1; vecA095017[w2]++; /* Lesser twin primes (subset of 6k+5 primes). */ } prev_prime = p; } int compute_primes(ULLI start, ULLI stop,int max_primes_collected) { unsigned int b; unsigned long s; ULLI k, j, ii; ULLI i; unsigned int bit_index = 0; ULLI indx; int c; unsigned long *list; FILE *fp; ULLI stop_here, t_stop_here; char ifnam[2048]; int errflag = 0; int write_flag = 0; s = (stop+2309)/2310 * 60 +60; /* bytes required */ fprintf(stderr,"Computing from %Lu to %Lu, %lu bytes required.\n\n\n", start,stop,s); /* print the first few small primes if they were requested */ for ( i = 0 ; i < 344 ; i ++) { if ( small_primes[i] > stop ) return 0; /* return to OS if nothing more to do */ if ( small_primes[i] >= start ) { prime_found(small_primes[i],max_primes_collected); } } fprintf(stderr,"attempting to malloc %lu bytes\n", s); list = malloc(s); if ( list == NULL ) { fprintf(stderr, "Could not allocate %lu bytes of memory\n", s); exit(1); } memset(list,0xff,s); j = 5; /* create array mapping 2310 integers to 480 bits */ mask_bit_index[0] = bit_index++; mask_bit_index[1] = bit_index++; for ( ii = 2; ii < 2310 ; ii++ ) { while (small_primes[j] < ii ) j++; /* if modulo any of the prime factors of 2310, then could not be prime */ if ( ii == small_primes[j] || ((ii%2)!=0 && (ii%3)!=0 && (ii%5)!=0 && (ii%7)!=0 && (ii%11)!=0)) { mask_bit_index[ii] = bit_index++; } else mask_bit_index[ii] = 0; #ifdef DEBUG printf("mask_bit_index[%d]=%u, j=%d\n", ii, mask_bit_index[ii], j); #endif } indx = 0L; stop_here = (unsigned long) (sqrt((stop+2309.0)/2310.0 *2310.0) + 0.5); for ( k = 0 ; k*2310 < stop ; k+=STEP_CONST ) { /* no need to sieve numbers larger than this range */ t_stop_here = (k+STEP_CONST)*2310UL; if ( stop_here < t_stop_here ) t_stop_here = stop_here; for( i = 13 ; i <= t_stop_here ; i+=2 ) { /* start with 13 since multiples of 2,3,5,7,11 are handled by the storage method */ /* increment by 2 for odd numbers only */ /* if ( (i >= ((k+STEP_CONST)* 2310))) break; /* no need to sieve numbers larger than this range, do next k */ b = i % 2310; /* i could not possibly be prime if remainder is 2,3,4,7,11 */ if ( mask_bit_index[b] == 0 || (i < k*2310 && (test_one_mask_bit(&list[i/2310 *15], mask_bit_index[b]))==0 ) ) continue; /* or this one already marked so it is not a prime */ /* */ if ( k == 0 ) indx = i*i; else indx = (k*2310) - (k*2310)%i +i; if ( (indx & 1) == 0 ) indx += i; /* start with i*i since any integer < i has already been sieved */ /* add 2 * i to avoid even numbers and mark all multiples of this prime */ for ( ; indx < (k+STEP_CONST)*2310 && indx <= (stop+2309)/2310*2310 ; indx +=(i+i)) { b = indx % 2310; /* modulo 2310 */ if ( mask_bit_index[b] != 0 ) { clear_one_mask_bit(&list[indx/2310 *15],mask_bit_index[b]); #ifdef DEBUG printf("indx = %lu, mask_bit_index[%d]=%d\n", indx, b, mask_bit_index[b]); #endif } } /* for indx */ } /* for i */ if ( start < (k+STEP_CONST)*2310 ) /* are there some to print now? */ { if ( k*2310 < start ) i = start; else i = k*2310; if ( (i & 1) == 0 ) i++; /* force it to be odd */ if ( 2311 >= i ) i = 2311; for ( ; i <= stop && i < (k+STEP_CONST)*2310; i+=2 ) { b = i % 2310; if ( mask_bit_index[b] != 0 && i >= start && (test_one_mask_bit(&list[i/2310 *15], mask_bit_index[b]))!=0 ) { prime_found(i,max_primes_collected); } } /* for i */ } } /* for k */ if(write_flag) { if(NULL == (fp = fopen(ifnam,"wb"))) { perror(ifnam); return 1; } i = fwrite( list, 1L, s, fp ); if(i != s) { fprintf(stderr,"write error: i=%Lu, s=%lu\n",i,s); } } return 0; }