/* http://www.research.att.com/~njas/sequences/a036284.c.txt Compute the sequences A036284 - A036286, A037096 - A037097 and A136378 - A136387 in Neil J.A. Sloane's Online Encyclopedia of Integer Sequences. See http://www.research.att.com/~njas/sequences/A036284 This program last edited by Antti Karttunen, 29. December 2007. (mail: his-firstname.his-surname@gmail.com) The original version of this program, written Sunday 29. November 1998, can be found at the address: http://www.iki.fi/karttu/matikka/vertfibo.htm Uses the following libraries: GNU MP library: http://gmplib.org/ (or: http://www.matematik.su.se/~tege/gmp/ ) (version 4.2.2, released 2007-09-11.) Victor Shoup's NTL-library: http://shoup.net/ntl/ (Works at least with: version 5.4.1 -- Release date: 2007.05.09) Compile e.g. as: g++ -I/usr/local/include -L/usr/local/lib A036284.c -o A036284 -lntl -lgmp -lm Identity A: if n > 1 then floor(V(n)/(2^(3*(2^(n-1))))) XOR (V(n)%(2^(3*(2^(n-1))))) = ((6/7)*((2^(3*(2^(n-1))))-1)) (where V(n) is A036284(n), and % means modulo and floor integer truncation towards zero.) or in other words: Take V(n) in octal or binary, halve it from the middle and xor the halves, and the result is a string of 6's (or 110's in binary). Note that this works also other way: V(n) = V(n)%(2^(3*(2^(n-1)))) + 2^(3*(2^(n-1))) * ((V(n)%(2^(3*(2^(n-1))))) XOR ((6/7)*((2^(3*(2^(n-1))))-1))) that is, the whole V(n) can be computed just from the first half, because the second half is the first half XORed with the string of octal 6's. Conjectures: In sequences A036284 and A036286 after n = 3 all or some terms are divisible by 10. Check for the divisibility by 5, by using the "rule of eleven" (5 is "eleven" of the base 4). */ #include #include #include /* Use these only if you use CGAL library's syntactic sugaring wrapper class: #include #include typedef CGAL_Gmpz Z; */ #define twos_power(n) (((unsigned long int) 1) << (n)) #define MPZ_getbit(n,i) (!mpz_sgn((n)) ? 0 : (i == (mpz_scan1 ((n), i)))) /* * void vertbinfibo(,,3) reversed should be in binary: * * 012345678901234567890123 * * 321098765432109876543210 * 000000110001011011101010 * * 131072 + 65536 + 4096 + 1024 + 512 + 128 + 64 + 32 + 8 + 2 = 202474 * = 404948/2 * */ #include NTL_CLIENT typedef unsigned long int ULLI; void initFromBinexp(GF2X *p_z2p,ULLI n) { int i = 0; unsigned char bytes[65]; while(0 != n) { bytes[i] = (n&255); i++; n >>= 8; } *p_z2p = GF2XFromBytes(bytes,(long)i); /* p_z2p->normalize(); Not necessary here. */ } void MPZ_to_GF2X(GF2X *p_z2p,mpz_t src) { int i = 0; unsigned char *bytes; size_t countp = 0; void (*freefunc) (void *, size_t); mp_get_memory_functions (NULL, NULL, &freefunc); bytes = (unsigned char *) mpz_export (NULL, &countp, -1, 1, 0, 0, src); *p_z2p = GF2XFromBytes(bytes,(long)countp); (freefunc)(bytes, 1); /* p_z2p->normalize(); Not necessary here. */ } void GF2X_to_MPZ(mpz_ptr p_dst,GF2X src) { long n_bytes = NumBytes(src); unsigned char *bytes = (unsigned char *) malloc(n_bytes+1); if(NULL == bytes) { fputs("\nAn internal error occurred! Cannot allocate more bytes!\n", stderr); exit(1); } BytesFromGF2X(bytes, src, n_bytes); mpz_import (p_dst, (size_t) n_bytes, (int) -1, (size_t) 1, 0, 0, bytes); free(bytes); } ZZ GF2XtoZZ(GF2X p_z2p) { ZZ z; long int n_bytes = NumBytes(p_z2p); unsigned char *bytes = (unsigned char *) malloc(n_bytes+1); if(NULL == bytes) { fputs("\nAn internal error occurred! Cannot allocate more bytes!\n", stderr); exit(1); } BytesFromGF2X(bytes,p_z2p,n_bytes); z = ZZFromBytes(bytes,n_bytes); free(bytes); return(z); } int number_of_factors_with_multiplicities(mpz_t src) { int totfact,j,u; GF2X gf2x_pol; vec_pair_GF2X_long factors; MPZ_to_GF2X(&gf2x_pol, src); CanZass(factors, gf2x_pol); u=factors.length(); for(totfact=0,j=0; j 0)) { /* printf("i=%lu.\n", i); */ i++; mpz_set (result, item); } while((!max_j || (j < max_j)) && (rule150reverse(item, result) > 0)) { /* printf("j=%lu.\n", j); */ j++; mpz_set (result, item); } *d90 = i; *d150 = j; return(i+j); } /* Maple procedure: A036284:=proc(n) option remember; local a, b, c, i, j, k, l, s, x, y, z; if (0 = n) then (6) else a := 0; b := 0; s := 0; x := 0; y := 0; k := 3*(2^(n-1)); l := 3*(2^n); j := 0; for i from 0 to l do z := bit_i(A036284(n-1),(j)); # j := (i mod k); c := (a + b + (`if`((x = y),x,(z+1))) mod 2); if(c <> 0) then s := s + (2^i); fi; a := b; b := c; x := y; y := z; j := j + 1; if(j = k) then j := 0; fi; od; RETURN(s); fi: end: */ /* * * Here we do binary addition in disguise, so as to keep the * sequence term a(n+1) strictly dependent only of the * previous term a(n). That is, we deduce "what the carry should be" * by examining what the sum bit (m0) of bits m2 and m1 is. * * +-----+ * |Carry| * +-----+-----+ * | n2 | m2 | * +-----+-----+ * | n1 | m1 | * +-----+-----+ * | n0 | m0 | * +-----+-----+ * * Carry = (((m2 XOR m1) AND (NOT m0)) OR (m2 AND m1)) * * m2, m1, m0 Carry * * 0 0 0 0 * * 0 1 1 0 * * 1 0 1 0 * * 1 1 0 1 * * 0 0 1 0 * * 0 1 0 1 * * 1 0 0 1 * * 1 1 1 1 * * * length = the length of cycle for bit-n of Fibonacci numbers. * prevlen = the length of cycle one to the right (it's length/2). * i = index to the bits of the integer to be formed (0 - (length-1)) * j = index to the bits of the prev_vbf, like i, but wraps back to 0 * at prevlen. * */ void vertbinfibo(mpz_t result, mpz_t prev_vbf, unsigned long int n) { if(0 == n) { mpz_set_ui (result, 6); return; } else { int i,j,k; int bit_n_2, bit_n_1, bit_n_0, bit_m_2, bit_m_1, bit_m_0; unsigned long int length = 3*(twos_power(n)); unsigned long int prevlen = length >> 1; /* 3*(twos_power(n-1)) */ for(i=j= bit_n_2 = bit_n_1 = bit_n_0 = bit_m_2 = bit_m_1 = bit_m_0 = 0, mpz_set_ui (result, 0); (i < length); /* The termination test. */ i++, j = ((j < (prevlen-1)) ? (j+1) : 0), bit_n_2 = bit_n_1, bit_n_1 = bit_n_0, bit_m_2 = bit_m_1, bit_m_1 = bit_m_0) { bit_m_0 = (j == (k = mpz_scan1 (prev_vbf, j))); bit_n_0 = bit_n_2 ^ bit_n_1 ^ (((bit_m_2 ^ bit_m_1) & !bit_m_0) | (bit_m_2 & bit_m_1)); if(bit_n_0) { mpz_setbit (result, i); } /* For debugging: printf("i=%lu j=%lu k=%lu bits_m2-m0=%d,%d,%d bits_n2-n0=%d,%d,%d\n", i,j,k, bit_m_2, bit_m_1, bit_m_0, bit_n_2, bit_n_1, bit_n_0); */ } return; } } /* * * Here we do binary addition in disguise, so as to keep the * sequence term a(n+1) strictly dependent only of the * previous terms a(n) and a(n-1). That is, we deduce "what the carry was" * by examining what the sum bit (m0) of bits l1 and m1 is. * * Later we should check what possibilities there are for the simplification of * the boolean expression in this particular case. * * +-----+ * |Carry| * +-----+-----+-----+ * | m1 | l1 | k1 | Previous row shifted once left + * +-----+-----+-----+ * | n1 | m1 | l1 | Previous row * +-----+-----+-----+ * | n0 | m0 | = 3 * previous row. * +-----+-----+ * * Carry = (((l1 XOR m1) AND (NOT m0)) OR (l1 AND m1)) * * l1, m1, m0 Carry * * 0 0 0 0 * * 0 1 1 0 * * 1 0 1 0 * * 1 1 0 1 * * 0 0 1 0 * * 0 1 0 1 * * 1 0 0 1 * * 1 1 1 1 * * * length = the length of the period for bit-n of powers of 3. * prevlen = the length of the period one to the right (length/2). * prev2len = the length of the period two columns to the right (length/4) * i = index to the bits of the integer to be formed (0 - (length-1)) * j = index to the bits of the prev_vbf, like i, but wraps back to 0 * at prevlen. * k = index to the bits of prev2vbf, like j, but wraps back to 0 at * prev2len. * */ void vertbinpow3(mpz_t result, mpz_t prev_vbf, mpz_t prev2vbf, unsigned long int n) { if(0 == n) { mpz_set_ui (result, 1); return; } else if(1 == n) { mpz_set_ui (result, 2); return; } else { int i,j,k; int bit_n_1, bit_n_0, bit_m_1, bit_m_0, bit_l_1, bit_l_0, bit_ind = 0; unsigned long int length = (twos_power(n)); unsigned long int prevlen = length >> 1; /* (twos_power(n-1)) */ unsigned long int prev2len = prevlen >> 1; /* (twos_power(n-2)) */ int prev_is_zero = !mpz_sgn (prev_vbf); int prev2_is_zero = !mpz_sgn (prev2vbf); for(i=j=k=0, bit_n_1 = bit_n_0 = bit_m_1 = bit_m_0 = bit_l_1 = bit_l_0 = 0, mpz_set_ui (result, 0); (i < length); /* The termination test. */ i++, j = ((j+1) & ~prevlen), /* j = ((j+1) mod prevlen) */ k = ((k+1) & ~prev2len), /* k = ((k+1) mod prev2len) */ bit_n_1 = bit_n_0, bit_m_1 = bit_m_0, bit_l_1 = bit_l_0) { bit_l_0 = (prev2_is_zero ? 0 : (k == (mpz_scan1 (prev2vbf, k)))); bit_m_0 = (prev_is_zero ? 0 : (j == (bit_ind = mpz_scan1 (prev_vbf, j)))); /* Carry = (((l1 XOR m1) AND (NOT m0)) OR (l1 AND m1)) */ bit_n_0 = bit_m_1 ^ bit_n_1 ^ (((bit_l_1 ^ bit_m_1) & !bit_m_0) | (bit_l_1 & bit_m_1)); if(bit_n_0) { mpz_setbit (result, i); } /* printf( "i=%lu j=%lu bit_ind=%lu k=%lu bits_m1-m0=%d,%d bits_n1-n0=%d,%d bits_l1-l0=%d,%d\n", i,j,bit_ind, k, bit_m_1, bit_m_0, bit_n_1, bit_n_0, bit_l_1, bit_l_0); */ } return; } } int main(int argc, char **argv) { mpz_t bitseq, prev_bitseq, prev2bitseq, outseq; unsigned long int count = 0, i; unsigned long int A_number = 0; int opt_S_line = 0; int first_printed = 0; int base = 10; int opt_reverse = 0, opt_pow3 = 0; int opt_popcount = 0, opt_binwidth = 0; int opt_start_offset = 0; int opt_distinct_factors = 0, opt_factors_with_multiplicities = 0; int opt_divide_3s_out = 0, opt_divide_9s_out = 0; int opt_divide_extra_2s = 0, opt_divide_extra_3s = 0; int opt_GF2X_factorize = 0; int opt_halfterms = 0, opt_hamdist = 0; char *progname, *s, c; progname = *argv; if(argc < 2) { print_usage: fprintf(stderr, "Usage: %s [-P|-S|-B|-X] [-b base] -A 36284-36286|37096-37097|136378-136387 n_terms\n", progname); exit(1); } while(NULL != (s = *(++argv))) { if(*s && isdigit(*s)) { count = atol(s); } else if('-' == *s) { while(c = *++s) { switch(c) { case 'A': { char *numarg = ++s; if(!*numarg) { numarg = *++argv; } if(!numarg || !*numarg || !isdigit(*numarg)) { fprintf(stderr, "%s: Option '-%c' expects a numeric argument!\n", progname,c); goto print_usage; } A_number = atol(numarg); goto suuret_rattaat; } case 'b': { char *numarg = ++s; if(!*numarg) { numarg = *++argv; } if(!numarg || !*numarg || !isdigit(*numarg)) { fprintf(stderr, "%s: Option '-%c' expects a numeric argument!\n", progname,c); goto print_usage; } base = atol(numarg); goto suuret_rattaat; } case 'h': { opt_halfterms++; break; } case 'H': { opt_hamdist++; break; } case 'r': { opt_reverse++; break; } case '3': { opt_pow3++; break; } case 'P': { opt_popcount++; break; } case 'S': { opt_S_line++; break; } case 'B': { opt_binwidth++; break; } case 'X': { opt_GF2X_factorize++; break; } default: { fprintf(stderr,"%s: An unknown option '%c' !\n", progname,*(s+1)); goto print_usage; } } } } else { goto print_usage; } suuret_rattaat: ; } if(!count) { goto print_usage; } /* Here are the 10 A-numbers you requested: A136378 --- A136387. */ switch(A_number) { case 36284: { break; } case 36285: { base = 8; break; } case 36286: { opt_reverse = 1; break; } case 37096: { opt_pow3 = 1; break; } case 37097: { opt_pow3 = 1; opt_halfterms = 1; opt_start_offset = 2; break; } case 136378: // A091221(A036284(n)) { opt_distinct_factors = 1; break; } case 136379: // A091222(A036284(n)) { opt_factors_with_multiplicities = 1; break; } case 136380: { opt_divide_9s_out = 1; opt_start_offset = 1; break; } case 136381: { opt_divide_9s_out = 1; opt_start_offset = 1; base = 8; break; } case 136382: { opt_divide_9s_out = 1; opt_divide_extra_2s = 2; opt_start_offset = 1; break; } case 136383: { opt_divide_9s_out = 1; opt_divide_extra_2s = 2; opt_start_offset = 1; base = 8; break; } case 136384: { opt_divide_9s_out = 1; opt_divide_extra_2s = 1; opt_divide_extra_3s = 1; opt_start_offset = 1; break; } case 136385: { opt_divide_9s_out = 1; opt_divide_extra_2s = 1; opt_divide_extra_3s = 1; opt_start_offset = 1; base = 8; break; } case 136386: { opt_pow3 = 1; opt_divide_3s_out = 1; opt_start_offset = 3; break; } case 136387: { opt_pow3 = 1; opt_divide_3s_out = 1; opt_start_offset = 3; base = 2; break; } default: { fprintf(stderr,"%s: Sequence with A-number %lu not implemented!\n", progname,A_number); goto print_usage; } } mpz_init (prev2bitseq); mpz_init (prev_bitseq); mpz_init (bitseq); mpz_init (outseq); /* THE MAIN LOOP: */ for(i=0; i <= count; i++, mpz_set (prev2bitseq, prev_bitseq), mpz_set (prev_bitseq, bitseq) ) { GF2X z2p; if(opt_pow3) { vertbinpow3(bitseq, prev_bitseq, prev2bitseq, i); } else { vertbinfibo(bitseq, prev_bitseq, i); } if(i < opt_start_offset) { continue; } if(opt_S_line) { if(!first_printed) { printf("%%S A%06lu ",A_number); first_printed = 1; } } else { printf("%lu ", i); } if(opt_reverse) { mpz_t copy; mpz_init (copy); reverse_bits(outseq, bitseq, 3*(twos_power(i))); reverse_bits(copy, outseq, 3*(twos_power(i))); if(mpz_cmp(copy, bitseq)) { fprintf(stdout,"\nAn internal error occurred! Double-reversion doesn't yield the same number: ("); mpz_out_str(stdout, base, copy); fprintf(stdout," != "); mpz_out_str(stdout, base, bitseq); fprintf(stdout,")\n"); } } else { mpz_set(outseq, bitseq); } if(opt_halfterms) { mpz_t tmpseq; mpz_init (tmpseq); /* Compute the lower half of the term: */ unsigned long int halflen = ((opt_pow3 ? 1 : 3) * twos_power(i))/2; mpz_tdiv_r_2exp (tmpseq, outseq, halflen); mpz_set(outseq, tmpseq); } if(opt_divide_9s_out) { int n_times = twos_power(i-1)-1; mpz_t res; mpz_init(res); divide_by_x_n_times(res,outseq,9,n_times); divide_by_x_n_times(res,res,2,opt_divide_extra_2s); divide_by_x_n_times(outseq,res,3,opt_divide_extra_3s); } else if(opt_divide_3s_out) { int n_times; mpz_t res; mpz_init(res); if(opt_halfterms) { n_times = (twos_power(i-2))-1; } else { n_times = (twos_power(i-1)+twos_power(i-2))-1; } divide_by_x_n_times(res,outseq,3,n_times); divide_by_x_n_times(res,res,2,opt_divide_extra_2s); divide_by_x_n_times(outseq,res,3,opt_divide_extra_3s); } /* Now output. Either a number of factors, or some bit-seq is to be output: */ if(opt_distinct_factors) { fprintf(stdout,"%u",number_of_distinct_factors(outseq)); } else if(opt_factors_with_multiplicities) { fprintf(stdout,"%u",number_of_factors_with_multiplicities(outseq)); } else /* The bit-sequence. */ { mpz_out_str(stdout, base, outseq); } /* Some extra-information about the sequence: */ if(opt_binwidth) { printf(" BW=%lu.", mpz_sizeinbase (outseq, 2)); } if(opt_popcount) { printf(" POP= %lu.", mpz_popcount(outseq)); } if(opt_hamdist) { printf(" HD= %lu.", mpz_hamdist(bitseq, prev_bitseq)); } if(opt_GF2X_factorize && (0 != mpz_sgn(outseq))) { vec_pair_GF2X_long factors; int j,u,totfact,tot3fact,tot7fact; mpz_t mpz_3, mpz_7; mpz_init_set_ui (mpz_3, 3); mpz_init_set_ui (mpz_7, 7); MPZ_to_GF2X(&z2p,outseq); ZZ copy_of = GF2XtoZZ(z2p); cout << ":" << copy_of; // For checking! cout << ":"; CanZass(factors, z2p); u=factors.length(); cout << "," << u << ":("; for(totfact=0,tot3fact=0,tot7fact=0,j=0; j