/* ------------------------------------------------------------------------- */ /* */ /* REVPOW.C - Powers Reversed with Schroeppel Method. */ /* An example application for gnu multiple precision library (gnu mp/GMP) */ /* Coded by Antti Karttunen 29.10.1998, and placed in Public Domain */ /* */ /* References: */ /* */ /* The multiplier (/2) and mask sequences, as well as an example sequence */ /* of three's powers reversed can now be found from Neil J.A. Sloane's */ /* Encyclopedia of Integer Sequences as A036213, A036214 and A036215 */ /* and accessible through the following URLs: */ /* */ /* http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=036213 */ /* http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=036214 */ /* http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=036215 */ /* */ /* GNU MP library: */ /* http://www.matematik.su.se/~tege/gmp/ */ /* */ /* */ /* CGAL C++ library for geometric calculations containing CGAL_Gmpz */ /* wrapper class (not needed unless you want some syntactic sugaring): */ /* http://www.risc.uni-linz.ac.at/projects/basic/cgal/doc_html-0.9/reference/ref-manual3/Chapter_Numbertype.html */ /* */ /* In Item 167 (Parity, count, reverse bits) of HAKMEM */ /* at URL: http://www.inwap.com/pdp10/hbaker/hakmem/hacks.html#item167 */ /* R. Schroeppel [rcs@cs.arizona.edu] gives the values for n=3 and 4 */ /* of multiplier, mask and divisor for reversing 6 and 8 bits respectively. */ /* This is also mentioned in DECsystem-10/20 Processor Reference Manual */ /* AA-H391A-TK, Chapter 2, User Operations, section 2.15: */ /* Programming Examples: Reversing Order of Digits */ /* */ /* ------------------------------------------------------------------------- */ #include #include #include /* Use these only if you use CGAL library's syntactic sugaring wrapper class: #include #include typedef CGAL_Gmpz Z; */ /* bit reversing multiplier for 2*n bits: 2 2n +2n 2 - 1 2 * -------------- 2n 2 - 1 bit reversing mask for 2*n bits: 2 2 2n +3n+1 2n -n-1 n 2 - 1 (4n+1) 2 - 1 2 * -------------- + 2 * -------------- 2n+1 2n+1 2 - 1 2 - 1 The first term is for the leftmost (n+1) bits (in the original) and the second term is for the rest (n-1) bits that are "wrapped over". n Taking the 2 as a common factor we get: 2 2 2n +3n+1 2n +2n 3n+1 n 2 + 2 - 2 - 1 2 * -------------------------------- 2n+1 2 - 1 bit reversing divisor for 2*n bits: 2n +2 2 - 1 */ /* bit reversing multiplier for 2*n bits: 2 2n +2n 2 - 1 2 * -------------- 2n 2 - 1 */ #define twos_power(n) (((unsigned long int) 1) << (n)) void bit_reversing_multiplier(mpz_t result, unsigned long int n) /* For 2*n bits */ { mpz_t numerator, denominator; mpz_init_set_ui (numerator,0); mpz_init_set_ui (denominator,0); mpz_setbit (numerator, 2*((n+1)*(n))); mpz_sub_ui (numerator, numerator, 1); mpz_setbit (denominator, 2*(n)); mpz_sub_ui (denominator, denominator, 1); mpz_tdiv_q (result, numerator, denominator); mpz_mul_2exp (result, result, 1); /* Multiply by 2**1, i.e. 2 */ mpz_clear (numerator); mpz_clear (denominator); } /* 2 2 2n +3n+1 2n +2n 3n+1 n 2 + 2 - 2 - 1 2 * -------------------------------- 2n+1 2 - 1 */ void bit_reversing_mask(mpz_t result, unsigned long int n) /* For 2*n bits */ { mpz_t numerator, denominator, auxiliary; mpz_init_set_ui (numerator,0); mpz_init_set_ui (denominator,0); mpz_setbit (numerator, ( (2*(n*n)) + (3*n) + 1 ) ); mpz_setbit (numerator, ( (2*(n*n)) + (2*n) ) ); /* This shortcut works only upto n=10 (as 3*10+1 = 31) in 32-bit machines: mpz_sub_ui (numerator, numerator, (twos_power((3*n)+1)+1)); So we use these three calls instead: */ mpz_init_set_ui (auxiliary, 0); mpz_setbit (auxiliary, ((3*n)+1) ); mpz_sub (numerator, numerator, auxiliary); mpz_sub_ui (numerator, numerator, 1); mpz_setbit (denominator, ( (2*(n)) + 1 ) ); mpz_sub_ui (denominator, denominator, 1); mpz_tdiv_q (result, numerator, denominator); mpz_mul_2exp (result, result, n); /* Multiply by 2**n */ mpz_clear (numerator); mpz_clear (denominator); mpz_clear (auxiliary); } /* bit reversing divisor for 2*n bits: 2n +2 2 - 1 */ void bit_reversing_divisor(mpz_t result, unsigned long int n) /* For 2*n bits */ { mpz_set_ui (result, 0); mpz_setbit (result, ( (2*n) + 2 ) ); mpz_sub_ui (result, result, 1); } int main(int argc, char **argv) { mpz_t nth_power, multiplier, mask, divisor, work, work2; unsigned long int n, count, i, j, binlength; int base = 10; if(argc < 3) { fprintf(stderr, "Usage: %s n count [base]\n", *argv); exit(1); } n = atol(*(argv+1)); if(mpz_init_set_str (nth_power, *(argv+1), 10) < 0) { fprintf(stderr, "%s: Please give a valid integer (in decimal system), not: \"%s\"\n", *argv,*(argv+1)); exit(1); } count = atol(*(argv+2)); if(argc > 3) { base = atoi(*(argv+3)); } mpz_init (multiplier); /* .ptr()->mpZ */ mpz_init (mask); mpz_init (divisor); mpz_init (work); mpz_init (work2); /* Output as: i nth_power (length) *multiplier &mask %divisor /[1|2] =nth_power_reversed */ for(i=1; i <= count; i++, mpz_mul_ui (nth_power, nth_power, n)) /* nth_power *= n */ { printf("%3lu. ", i); mpz_out_str(stdout, base, nth_power); binlength = mpz_sizeinbase (nth_power, 2); j = (binlength+1)/2; printf(" (%lu. %lu.) *", binlength, j); bit_reversing_multiplier(multiplier, j); mpz_out_str(stdout, base, multiplier); bit_reversing_mask(mask, j); printf(" &"); mpz_out_str(stdout, base, mask); bit_reversing_divisor(divisor, j); printf(" %%"); mpz_out_str(stdout, base, divisor); /* work = (((nth_power * multiplier) & mask) % divisor); */ mpz_mul (work, nth_power, multiplier); /* work = nth_power*multiplier; */ mpz_and (work, work, mask); /* work &= mask; */ mpz_mod (work, work, divisor); /* work %= divisor; */ if(binlength & 1) /* An odd length, so we have to shift right once */ { printf(" /2 "); mpz_tdiv_q_2exp (work, work, 1); /* work >>= 1; */ } else { printf(" /1 "); } printf(" ="); mpz_out_str(stdout, base, work); if(n) { mpz_tdiv_qr_ui (work, work2, work, n); printf(" /%lu =", n); mpz_out_str(stdout, base, work); printf(" %%%lu =", n); mpz_out_str(stdout, base, work2); } printf("\n"); } }