Fractals/Mathematics/doubling

Definition Angle doubling mapOn Combinatorial Types of Periodic Orbits of the Map x↦kx (mod Z). Part I: Realization Authors: Carsten L. Petersen, Saeed Zakerichaos_and_doubling map by Mark McClure Input: proper fraction in forms decimal ratio ( if r is not irrational ) decimal sequence of digits and radix binary sequence of digits and radix Output: proper fraction Doubling map on binary fraction Effect of doubling map d on binary representation of fraction x is to simply shift the bits of x to the left, discarding the bit that shifted into the ones place = left shift.doubling map by Mark McClure Example Orbits Iteration of the map gives orbit ( sequence of angles): periodic ( input is a rational number with odd denominator ) preperiodic ( input is rational number with even denominator ) aperiodic ( when input is irrational number ) All above types are infinite. Finite binary expansion has equal infinite preperiodic representation. periodic Periodic orbits of angles under doubling map Note that here chord joining 2 points z1 and z2 on unit circle means that . It does not mean that these points are landing points of the same ray. Some orbits do not cross : but some do cross: Periodic orbits for angles with odd denominator d Where: n is numerator of the fraction = integer from 0 to (d-1) d is denominator of the fraction = positive integer p is a period = positive integer Period 1 ( 1 orbit) for angles = {0/1 = 1/1 } Period 2 ( 1 orbit) for angles = {1/3, 2/3} Period 3 ( 2 orbits) for angles = {1/7, 2/7, 4/7 } {3/7, 6/7, 5/7 } Period 4 ( 3 orbits) for angles = {1/15, 2/15, 4/15, 8/15 } {3/15, 6/15, 12/15, 9/15 } {7/15, 14/15, 13/15, 11/15 } Period 5 ( 6 orbits) for angles = { 1/31, 2/31, 4/31, 8/31, 16/31} { 3/31, 6/31, 12/31, 24/31, 17/31} { 5/31, 10/31, 20/31, 9/31, 18/31} { 7/31, 14/31, 28/31, 25/31, 19/31} {11/31, 22/31, 13/31, 26/31, 21/31} {15/31, 30/31, 29/31, 27/31, 23/31} Period 6 ( 9 orbits ) for angles = , only numerators { 1, 2, 4, 8,16,32} { 3, 6,12,24,48,33} { 5,10,20,40,17,34} { 7,14,28,56,49,35} [11,22,44,25,50,37} [13,26,52,41,19,38} [15,30,60,57,51,39} [23,46,29,58,53,43] [31,62,61,59,55,47] Period = 7 for angles = Period = 8 for angles = Period = 9 for angles = See also external ray for root points ( bond point) on the parameter plane: wake subwake island = midget = mini Mandelbrot set dynamic plane external rays that land on the parabolic fixed point . ( first angle is angle of the wake) preperiodic thumb|right|Preperiodic orbit Here denominator d of fraction r is even: and q is odd. Here t is preperiod Example orbits: Preperiod 1 and period 1 for angles {1/2 , (0/2)} See Misiurewicz points on the parameter plane Dendrite Julia sets on the dynamic plane aperiodic thumb|right|Irrational rotation Compare Siegel disc Herman rings ( rational maps) Implementations integer: 2 integers ( numerator and denominator) rational number mpq type from GMP library floating point double Forward Using 32bit signed int limits maximum preperiod to about 30 Using double for fractions limits maximum number of accurate bits to about 53 C using integers #include <stdio.h> #include <stdlib.h> // malloc #include <limits.h> // INT_MAX #include <assert.h> #include <math.h> /* gcc d.c -Wall -Wextra -lm example input fractions wikibooks Fractals/Mathematics/binary */ int nMax; /* https://stackoverflow.com/questions/19738919/gcd-function-for-c The GCD function uses Euclid's Algorithm. It computes A mod B, then swaps A and B with an XOR swap. */ int gcd(int a, int b) { int temp; while (b != 0) { temp = a % b; a = b; b = temp; } return a; } /* n/d -> n_new/d_new */ int give_reduced_fraction(const int n, const int d, int * n_new, int *d_new){ int divisor = gcd(n,d); if (divisor != 1) { *n_new = n / divisor; *d_new = d / divisor;} else { *n_new = n; *d_new = d; } return 0; } /* Binary expansion can be : * finite - decimal ratio number with even denominator which is a power of 2. Note that it has also equal infinite preperiodic representation * infinite ** periodic : preperiod = 0, period > 0, rational number with odd denominator which is not a power of 2 ** preperiodic ( = eventually periodic) : preperiod > 0, period > 0, rational number with even denominator ** aperiodic : preperiod = 0, period = 0 ( or infinity), irrational number input is a rational number t = n/d so here are only 2 types of result: * 0 = preperiodic ( = eventually periodic) : preperiod > 0, period > 0, rational number with even denominator * 1 = periodic : preperiod = 0, period > 0, rational number with odd denominator which is not a power of 2 */ int give_dynamic_type(const int n, const int d){ // boolean flag int is_odd = 0; int is_POT = 0; if (d % 2) { // d % 2 > 0 fprintf(stdout, "denominator %d is odd\n", d); // parzysty is_odd = 1; } else { // d % 2 = 0 fprintf(stdout, "denominator %d is even\n", d); // nieparzysty is_odd = 0; } if (d>0 && ((d & (d-1)) == 0)) { fprintf(stdout, "denominator %d is power of 2 ( POT) = 2^%d\n", d, (int) log2(d)); is_POT = 1;} else { fprintf(stdout, "denominator %d is not the power of 2\n", d); is_POT = 0;} // wikibooks Fractals/Mathematics/binary#preperiod_and_period if ( is_odd == 0 && is_POT == 1) { fprintf(stdout, "Binary expansion of %d/%d is finite and has equal infinite preperiodic representation\n", n,d); //fprintf(stdout, "denominator is even and POT\n"); return 0; } // preperiodic ( = eventually periodic) : preperiod > 0, period > 0, rational number with even denominator if (is_odd == 0 && is_POT == 0) { fprintf(stdout, "Binary expansion of %d/%d is preperiodic\n", n,d); //fprintf(stdout, "denominator is even and not POT\n"); return 0; } // periodic : preperiod = 0, period > 0, rational number with odd denominator which is not a power of 2 if (is_odd == 1 && is_POT == 0) { fprintf(stdout, "Binary expansion of %d/%d is periodic\n", n,d); //fprintf(stdout, "denominator is not even (odd) and not POT\n"); return 1; } return 0; } int give_period(const int n0, const int d0){ int n = n0; int d = d0; int i; int iMax = 20; // period_max int period = 0; // printf(" i\t numerator/denominator \n"); // header for ( i=0 ; i< iMax; i++){ // printf("%3d:\t %3d / %3d \n", i, n, d); // check signed integer overflow before it will happen if ( n > nMax ) {printf(" signed integer overflow will happen = wrap; use extended precision \n"); return 1;} // angle doubling map modulo 1 n = 2*n; n = n % d; // if (n == n0) { period = i+1; // printf(" preperiod = 0 and period of fraction under doubling map is %d\n", period); return period;} } return 0; // period not found, maybe period > iMax } int give_periodic_orbit(const int period, int n, int d, int orbit[]){ int i; int iMax = period; // for ( i=0 ; i< iMax; i++){ orbit[i] = n; // check signed integer overflow before it will happen if ( n > nMax ) {printf(" signed integer overflow will happen = wrap; use extended precision \n"); return 1;} // angle doubling map modulo 1 n = 2*n; n = n % d; } return 0; } int give_preperiod(const int period, const int n0, const int d0, int orbit[]){ int n = n0; int d = d0; int i; int iMax = period; // for ( i=0 ; i< iMax; i++){ if (orbit[i] == n) return i; // angle doubling map modulo 1 n = 2*n; n = n % d; } return 0; } /* input : reduced proper rational fraction in t = n0/d0 output * perperiod as a result * period as a pointer */ int give_preperiod_and_period(const int n0, const int d0, int *period){ int n = n0; int d = d0; *period = 0; int preperiod = 0; int preperiod_max = 1000; if ( preperiod_max > INT_MAX - *period ){printf("give_preperiod_and_period error: preperiod_max > INT_MAX - period, signed integer overflow will happen = ; use extended precision \n"); return -1;} int i; int iMax = preperiod_max; // preperiod_max // go thru preperiodic part to periodic part for ( i=0 ; i< iMax; i++){ //printf("%3d:\t %3d / %3d \n", i, n, d); // check signed integer overflow before it will happen if ( n > nMax ) {printf("give_preperiod_and_period error: signed integer overflow will happen = wrap; use extended precision \n"); return 1;} // angle doubling map modulo 1 n = 2*n; n = n % d; } // now it should be periodic part *period = give_period (n,d); // periodic orbit int *orbit; // only numerators orbit = (int *) malloc( *period * sizeof(* orbit)); give_periodic_orbit(*period, n, d, orbit); preperiod = give_preperiod( *period, n0, d0, orbit); free(orbit); return preperiod; } void give_orbit(const int n0, const int d0, const int preperiod, const int period){ int n = n0; int d = d0; int i; int iMax = preperiod+period; // preperiod_max for ( i=0 ; i< iMax; i++){ if ( i < preperiod) { printf("%3d:\t %3d / %3d \t preperiodic part \n", i, n, d);} else {printf("%3d:\t %3d / %3d \t periodic part \n", i, n, d);} // angle doubling map modulo 1 n = 2*n; n = n % d; } } /* Algorithm:[36] Multiply the input decimal fraction by two from above result take integer part as the binary digit take the fractional part as the starting point for the next step repeat until you either get to 0 or a periodic number read the number starting from the top - the first binary digit is the first digit after the comma */ void print_binary_expansion(const int n0, const int d0, const int preperiod, const int period){ int n = n0; int d = d0; int i_max = preperiod+period; int i; double t = (double) n/d; double t_fractional; double t_integer; int binary_digit; printf("formated infinite binary expansion of %d/%d is 0.", n0,d0); for (i = 0; i < i_max; ++i){ // doubling map t *= 2.0; // split t_fractional = modf(t, &t_integer); // binary_digit = (int) t_integer; if (i== preperiod) {printf("(");} printf("%d", binary_digit); // take the fractional part as the starting point for the next step t = t_fractional; } printf(")\n"); } int describe_fraction(const int n0, const int d0){ // proper decimal fraction // t = n/d //int n0 = 1; // //int d0 = 66; // // tr = n0r/d0r = t after reduction int n0r ; // int d0r ; // int n; int d; int period = 0; int preperiod = 0; assert(n0 > 0); assert(d0 > 1); assert(n0 < d0); // proper fraction // ------------ Reducing Fractions to Lowest Terms ------------------------------- // ----------- wikipedia Irreducible_fraction ---------------------------- give_reduced_fraction(n0, d0, &n0r, &d0r); if (n0 == n0r && d0 ==d0r ) {printf(" fraction = %d/%d\t is irreducible = in lowest terms \n", n0, d0 ); } else {printf(" fraction = %d/%d\t after reduction is %d/%d \n", n0, d0, n0r,d0r); } n = n0r; d = d0r; assert(n > 0); assert(d > 1); assert(n < d); int type = give_dynamic_type(n,d); // ------------------compute preperiod and period under doubling map ------------------------- printf("fraction %d/%d under doubling map has: \n", n0r, d0r); if (type ==0 ){ preperiod = give_preperiod_and_period( n0r, d0r, &period); if (preperiod > 0) { printf("\t preperiod = %d and period = %d\n", preperiod, period); if (period == 0 ) {printf("period = 0 means period NOT FOUND !!!\n\t maybe period > iMax in give_period \n\tor maybe preperiod_max in give_preperiod_and_period is to big \n");}}} // -------------- else { period = give_period(n,d); if (period > 0) {printf(" preperiod = 0 and period = %d\n", period); } else {printf(" preperiod = 0 and period of fraction under doubling map > 0 but period NOT FOUND !!! maybe period > iMax in give_period \n");}} // ------------------------------------------------- give_orbit( n0, d0, preperiod, period); // ----------formated infinite binary expansion --------------------- print_binary_expansion( n0r, d0r, preperiod, period); return 0; } int main(void) { nMax = INT_MAX / 2; // signed integer describe_fraction(1,22); return 0; } C using double numbers using double gives precision of 53 bits of binary number expansion (the number of bits of accuracy of double precision floating point) t *= 2.0; // t = 2*t if (t > 1.0) t--; // t = t modulo 1 /* */ #include <stdio.h> // printf #include <math.h> // fabs #define iMax 8 // int main(){ double t0 ; double t ; double ti; // final t after iMax iterations double tr; // double dt; int itinerary[iMax]= {0}; int i; t0 = (double) 1/7; t = t0; // check the input : it should be 0.0 <= t < 1.0 if (t>1.0) {printf("t is > 1.0\n"); return 1;} if (t<0.0) {printf("t is < 0.0\n"); return 1;} printf("forward iteration of doubling map\n"); for(i=0; i<iMax; i++){ printf("t%d = %f", i, t); t = t*2.0; // doubling if (t>1.0) { itinerary[i]= 1; t = t - 1.0; printf(" wrap\n");} // modulo 1 else printf("\n"); } printf("t%d = %f\n", i, t); // ti = t; printf("\nbackward iteration of doubling map = halving map \n"); // for(i=iMax; i>0; i--){ // reverse counting printf("t%d = %f", i, t); if (itinerary[i-1]==1) { // i-1 !!! t = t + 1.0; printf(" unwrap\n");} // modulo 1 else printf("\n"); t = t/2.0; // halving } printf("t%d = %f\n", i, t); tr = t; // printf("\n\nresults \n"); printf("t0 = %f\n", t0); printf("t%d = %f\n", iMax, ti); dt = fabs(t0- tr); printf("tr = %f\n", tr); printf("dt = fabs(t0- tr) = %f\n", dt ); printf("\nitinerary:\n"); for(i=0; i<iMax; i++) printf("itinerary[%d] = %d \n", i, itinerary[i]); printf("\ndecimal %f has binary expansion = 0.", t0); for(i=0; i<iMax; i++) printf("%d", itinerary[i]); printf("\n"); if (dt < 0.0000000001) printf("program works good !\n"); else printf("program fails !\n"); return 0;} C using rational type mpq from GMP library ( arbitrary precision) C function (using GMP library)Programowanie w systemie UNIX: GMP in polish wikibooks: // rop = (2*op ) mod 1 void mpq_doubling(mpq_t rop, const mpq_t op) { mpz_t n; // numerator mpz_t d; // denominator mpz_inits(n, d, NULL); // mpq_get_num (n, op); // mpq_get_den (d, op); // n = (n * 2 ) % d mpz_mul_ui(n, n, 2); mpz_mod( n, n, d); // output mpq_set_num(rop, n); mpq_set_den(rop, d); mpz_clears(n, d, NULL); } /* C programme using gmp gcc r.c -lgmp http://gmplib.org/manual/Rational-Number-Functions.html#Rational-Number-Functions */ #include <stdio.h> #include <gmp.h> // input = num/den unsigned int num = 1; unsigned int den = 6; int main () { mpq_t q; // rational number; int b =2 ; // base of numeral system mpz_t n ; mpz_t d ; mpf_t f ; char *sr; char *sf; mp_exp_t exponent ; // holds the exponent for the result string // init and set variables mpq_init (q); // Initialize r and set it to 0/1. mpf_init (f); mpz_init_set_ui(n,num); mpz_init_set_ui(d,den); mpq_set_num(q,n); mpq_set_den(q,d); mpq_canonicalize (q); // It is the responsibility of the user to canonicalize the assigned variable before any arithmetic operations are performed on that variable. mpf_set_q(f,q); // convertions sr = mpq_get_str (NULL, b, q); // rational number = fraction : from decimal to binary // If str is NULL, the result string is allocated using the current allocation function sf = mpf_get_str (NULL, &exponent, b, 0, f); // floating point number : from decimal to binary // If n_digits is 0 then that accurate maximum number of digits are generated. // print gmp_printf ("a rational number %Zd / %Zd is in a canonical form ( decimal fraction) : %Qd\n",n,d, q); // gmp_printf(" numerator of rational number = %Zd \n", mpq_numref(q)); // gmp_printf(" denominator of rational number = %Zd \n", mpq_denref(q)); // gmp_printf (" binary rational ( string ) : %s \n", sr); // gmp_printf (" decimal floating point number : %.Ff \n", f); // // The output of the GMP programme is in the form mantissa,'e',exponent where the value is // 0.mantissa * 2^exponent // GMP represents the floating point numbers (for base 2) as a pair of exponent and mantissa (and sign). //The generated string is a fraction, with an implicit radix point immediately to the left of the first digit. // Trailing zeros are not returned. // For example, the number 3.1416 would be returned as string "31416" and exponent 1 so : // 3.1416 = 0.31416*10^1 // another example : 1/6 as a exponential form of binary number will be : // mantissa = 101010101010101010101010101010101010101010101010101010101010101011 // exponet = -2 // so 1/6 = 0.101010101010101010101010101010101010101010101010101010101010101011 * 2 ^(-2) = 0.0(01) gmp_printf (" mantissa of binary floating ( string ) : %s \n", sf); // gmp_printf (" exponent : %ld \n", exponent); // printf ("binary floating number in exponential form =0.%s*%d^%ld\n", sf,b, exponent); // clear memory mpq_clear (q); mpz_clear (n); mpz_clear (d); mpf_clear (f); return 0; } Compile: gcc r.c -lgmp run ./a.out result: a rational number 1 / 6 is in a canonical form ( decimal fraction) : 1/6 numerator of rational number = 1 denominator of rational number = 6 binary rational ( string ) : 1/110 decimal floating point number : 0.166666666666666666667 mantissa of binary floating ( string ) : 101010101010101010101010101010101010101010101010101010101010101011 exponent : -2 binary floatin in exponential form =0.101010101010101010101010101010101010101010101010101010101010101011*2^-2 Maxima CAS Maxima CAS function using numerator and denominator as an input doubling_map(n,d):=mod(2*n,d)/d $ or using rational number as an input DoublingMap(r):= block([d,n], n:ratnumer(r), d:ratdenom(r), mod(2*n,d)/d)$ Common Lisp function (defun doubling-map (ratio-angle) " period doubling map = The dyadic transformation (also known as the dyadic map, bit shift map, 2x mod 1 map, Bernoulli map, doubling map or sawtooth map " (let* ((n (numerator ratio-angle)) (d (denominator ratio-angle))) (setq n (mod (* n 2) d)) ; (2 * n) modulo d (/ n d))) ; result = n/d Haskell Haskell functionlavaurs' algorithm in Haskell with SVG output by Claude Heiland-Allen -- by Claude Heiland-Allen -- type Q = Rational double :: Q -> Q double p | q >= 1 = q - 1 | otherwise = q where q = 2 * p C++ // mndcombi.cpp by Wolf Jung (C) 2010. // http://mndynamics.com/indexp.html // n is a numerator // d is a denominator // f = n/d is a rational fraction (angle in turns) // twice is doubling map = (2*f) mod 1 // n and d are changed (arguments passed to function by reference) void twice(unsigned long long int &n, unsigned long long int &d) { if (n >= d) return; if (!(d & 1)) { d >>= 1; if (n >= d) n -= d; return; } unsigned long long int large = 1LL; large <<= 63; //avoid overflow: if (n < large) { n <<= 1; if (n >= d) n -= d; return; } n -= large; n <<= 1; large -= (d - large); n += large; } Inverse function of doubling map Every angle α ∈ R/Z measured in turns has: one image = 2α mod 1 under doubling map "two preimages under the doubling map: α/2 and (α + 1)/2.".SYMBOLIC DYNAMICS AND SELF-SIMILAR GROUPS by VOLODYMYR NEKRASHEVYCH Inverse of doubling map is multivalued function. In Maxima CAS: InvDoublingMap(r):= [r/2, (r+1)/2]; Note that difference between these 2 preimages is half a turn = 180 degrees = Pi radians. + Images and preimages under doubling map d Q&A analysis find expansion type find preperiod and period convert decimal fraction to binary expansion format the expansion How to find expansion type Binary expansion of decimal fraction /* gcc i.c -Wall -Wextra -lm https://stackoverflow.com/questions/108318/how-can-i-test-whether-a-number-is-a-power-of-2 (n>0 && ((n & (n-1)) == 0)) https://stackoverflow.com/questions/160930/how-do-i-check-if-an-integer-is-even-or-odd if (x % 2) { } // x is odd An integer is even if it is a multiple of two, power of 2, true if num is perfectly divisible by 2 : mod(24,2) = 0 and odd if it is not. */ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <assert.h> int main(void) { int n = 1; int d = 4; // // boolean flag int is_odd = 0; int is_POT = 0; //int type; // assert( n > 0); assert( d > 0); assert( n < d); // proper fraction if (d % 2) { // d % 2 > 0 fprintf(stdout, "denominator %d is odd\n", d); // parzysty is_odd = 1; } else { // d % 2 = 0 fprintf(stdout, "denominator %d is even\n", d); // nieparzysty is_odd = 0; } if (d>0 && ((d & (d-1)) == 0)) { fprintf(stdout, "denominator %d is power of 2 ( POT) = 2^%d\n", d, (int) log2(d)); is_POT = 1;} else { fprintf(stdout, "denominator %d is not the power of 2\n", d); is_POT = 0;} // https://en.wikibooks.org/wiki/Fractals/Mathematics/binary#preperiod_and_period if ( is_odd == 0 && is_POT == 1) { fprintf(stdout, "Binary expansion of %d/%d is finite and has equal infinite representation\n", n,d); fprintf(stdout, "denominator is even and POT\n"); } // preperiodic ( = eventually periodic) : preperiod > 0, period > 0, rational number with even denominator if (is_odd == 0 && is_POT == 0) { fprintf(stdout, "Binary expansion of %d/%d is preperiodic\n", n,d); fprintf(stdout, "denominator is even and not POT\n"); } // periodic : preperiod = 0, period > 0, rational number with odd denominator which is not a power of 2 if (is_odd == 1 && is_POT == 0) { fprintf(stdout, "Binary expansion of %d/%d is periodic\n", n,d); fprintf(stdout, "denominator is not even (odd) and not POT\n"); } return 0; } How to find preperiod and period ? Period and the doubling map How to find the period of angle under doubling map visual numerical read period from denominator of decimal fraction ( reduced rational fraction m/n ) find period/preperiod in the binary expansion ( binary sequence) string matching algorithms, for instance, naive or brute-force search, Knuth-Morris-Pratt algorithmmath stackexchange question: period-of-a-finite-binary-sequence discrete Fourier transform read it from the itinerary of angle under doubling map Period of binary expansion of reduced rational fraction m/n is equal to the multiplicative order of 2 modulo n: C version double precision: forward and inverse doubling map /* doubling map 2*t mod 1 how to invert doubling map Inverse of doubling map is multivalued function: 2 preimages t/2 and (t+1)/2 to choose proper preimage one needs extra information from forward iteration = itinerary itinerary : list of symbols for coding the orbits of a given dynamical system by partitioning the space X and forming an itinerary http://www.maths.qmul.ac.uk/~sb/cf_chapter4.pdf see also how to convert proper decimal fraction to binary fraction commons :Binary_decomposition_of_dynamic_plane_for_f0(z)_%3D_z%5E2.png ---------- git -------------------- cd existing_folder git init git remote add origin git@gitlab.com:adammajewski/doubling_map.git git add . git commit -m "Initial commit" git push -u origin master */ #include <stdio.h> // printf #include <math.h> // fabs #define iMax 8 // int main(){ double t0 ; double t ; double ti; // final t after iMax iterations double tr; // double dt; int itinerary[iMax]= {0}; int i; t0 = (double) 1/7; t = t0; // check the input : it should be 0.0 <= t < 1.0 if (t>1.0) {printf("t is > 1.0\n"); return 1;} if (t<0.0) {printf("t is < 0.0\n"); return 1;} printf("forward iteration of doubling map\n"); for(i=0; i<iMax; i++){ printf("t%d = %f", i, t); // wikipedia Dyadic_transformation t = t*2.0; // doubling if (t>1.0) { itinerary[i]= 1; t = t - 1.0; printf(" wrap\n");} // modulo 1 else printf("\n"); } printf("t%d = %f\n", i, t); // ti = t; printf("\nbackward iteration of doubling map = halving map \n"); // for(i=iMax; i>0; i--){ // reverse counting printf("t%d = %f", i, t); if (itinerary[i-1]==1) { // i-1 !!! t = t + 1.0; printf(" unwrap\n");} // modulo 1 else printf("\n"); t = t/2.0; // halving } printf("t%d = %f\n", i, t); tr = t; // printf("\n\nresults \n"); printf("t0 = %f\n", t0); printf("t%d = %f\n", iMax, ti); dt = fabs(t0- tr); printf("tr = %f\n", tr); printf("dt = fabs(t0- tr) = %f\n", dt ); printf("\nitinerary:\n"); for(i=0; i<iMax; i++) printf("itinerary[%d] = %d \n", i, itinerary[i]); printf("\ndecimal %f has binary expansion = 0.", t0); for(i=0; i<iMax; i++) printf("%d", itinerary[i]); printf("\n"); if (dt < 0.0000000001) printf("program works good !\n"); else printf("program fails !\n"); return 0;} arbitrary precision // gcc d.c -lgmp -Wall #include <stdio.h> #include <gmp.h> // a multiple precision integer, as defined by the GMP library. The C data type for such integers is mpz_t int print_z(mpz_t z, int base, char *s){ printf("%s= ", s); mpz_out_str (stdout, 10, z); printf (" for base = %d\n", base); return 0; } // rop = (2*op) mod 1 // wikipedia : dyadic_transformation or doubling map void mpq_doubling(mpq_t rop, const mpq_t op) { mpz_t n; // numerator mpz_t d; // denominator mpz_inits(n, d, NULL); // mpq_get_num (n, op); // mpq_get_den (d, op); // n = (n * 2 ) % d mpz_mul_ui(n, n, 2); mpz_mod( n, n, d); // output mpq_set_num(rop, n); mpq_set_den(rop, d); mpz_clears(n, d, NULL); } int main () { int i; // unsigned long int e = 89; // exponent is also a period of doubling map unsigned long int b = 2; // arbitrary precision variables from GMP library mpz_t n ; // numerator of q mpz_t d ; // denominator of q mpq_t q; // rational number q = n/d // init and set variables mpz_init_set_ui(n, 1); // d = (2^e) -1 // http://fraktal.republika.pl/mset_external_ray.html mpz_init(d); mpz_ui_pow_ui(d, b, e) ; // d = b^e mpz_sub_ui(d, d, 1); // d = d-1 // q = n/d mpq_init (q); // mpq_set_num(q,n); mpq_set_den(q,d); mpq_canonicalize (q); // It is the responsibility of the user to canonicalize the assigned variable before any arithmetic operations are performed on that variable. // print //print_z(d, 10, "d "); //print_z(n, 10, "n "); gmp_printf ("q = %Qd\n",q); // // for (i=0; i<(1+2*e) ; i++){ mpq_doubling(q, q); gmp_printf ("q = %Qd\n",q); // } // clear memory mpq_clear (q); mpz_clears(n, d, NULL); return 0; } C++ version /* based on : mndcombi.cpp by Wolf Jung (C) 2010. http://mndynamics.com/indexp.html which is the part of Mandel 5.5 multiplatform C++ GUI program using QT on the same licence as above "The function is computing the preperiod and period (of n/d under doubling map) and setting the denominator to 2^preperiod*(2^period - 1) if possible. So 1/5 becomes 3/15 and 2/10 becomes 3/15 as well. The period is returned as the value of the function, n and d are changed ( Arguments passed to function by reference) and the preperiod is returned in k." (Wolf Jung) Question : if result is >=0 why do not use unsigneg char or unsigned int for type of result ??? */ int normalize(unsigned long long int &n, unsigned long long int &d, int &k) { if (!d) return 0; // d==0 error n %= d; while (!(n & 1) && !(d & 1)) { n >>= 1; d >>= 1; } int p; unsigned long long int n0, n1 = n, d1 = d, np; k = 0; while (!(d1 & 1)) { k++; d1 >>= 1; if (n1 >= d1) n1 -= d1; } n0 = n1; for (p = 1; p <= 65 - k; p++) { twice(n1, d1); if (n1 == n0) break; } if (k + p > 64) return 0; // more then max unsigned long long int np = 1LL; np <<= (p - 1); np--; np <<= 1; np++; //2^p - 1 for p <= 64 n0 = np; d >>= k; n1 = d; if (n1 > n0) { n1 = n0; n0 = d; } while (1) { d1 = n0 % n1; if (!d1) break; n0 = n1; n1 = d1; } //gcd n1 n /= d/n1; n *= np/n1; d = np << k; return p; } Lisp version (defun give-period (ratio-angle) "gives period of angle in turns (ratio) under doubling map" (let* ((n (numerator ratio-angle)) (d (denominator ratio-angle)) (temp n)) ; temporary numerator (loop for p from 1 to 100 do (setq temp (mod (* temp 2) d)) ; (2 x n) modulo d = doubling) when ( or (= temp n) (= temp 0)) return p ))) Maxima CAS version DoublingMap(r):= block([d,n], n:ratnumer(r), d:ratdenom(r), mod(2*n,d)/d)$ /* Tests : GivePeriod (1/7) 3 GivePeriod (1/14) 0 GivePeriod (1/32767) 15 GivePeriod (65533/65535) 16 Gives 0 if : * not periodic ( preperiodic ) * period >pMax */ GivePeriod (r):= block([rNew, rOld, period, pMax, p], pMax:100, period:0, p:1, rNew:DoublingMap(r), while ((p<pMax) and notequal(rNew,r)) do (rOld:rNew, rNew:DoublingMap(rOld), p:p+1 ), if equal(rNew,r) then period:p, period ); Haskell version Haskell versionlavaurs' algorithm in Haskell with SVG output by Claude Heiland-Allen Conversion from an integer type (Int or Integer) to anything else is done by "fromIntegral". The target type is inferred automatically -- by Claude Heiland-Allen -- import Data.List (findIndex, groupBy) -- type N = Integer -- type Q = Rational period :: Q -> N period p = let Just i = (p ==) `findIndex` drop 1 (iterate double p) in fromIntegral i + 1 conversion Convert binary fraction to decimal ratio convert decimal fraction to binary knowledgedoor calculators Can all decimal fractions be converted exactly to binary? Not all. Only those for which denominator is a power of 2 ( finite ) have exact decimal representation. "In every other case, there will be an error in the representation. The error's magnitude depends on the number of digits used to represent it." binary-fraction calculator by Davide Borchia decimal to binary #include <stdio.h> #include <stdlib.h> // malloc #include <limits.h> // INT_MAX #include <assert.h> #include <math.h> /* gcc d.c -Wall -Wextra -lm example input fractions in wikibooks /Fractals/Mathematics/binary */ int nMax; /* https://stackoverflow.com/questions/19738919/gcd-function-for-c The GCD function uses Euclid's Algorithm. It computes A mod B, then swaps A and B with an XOR swap. */ int gcd(int a, int b) { int temp; while (b != 0) { temp = a % b; a = b; b = temp; } return a; } /* n/d -> n_new/d_new */ int give_reduced_fraction(const int n, const int d, int * n_new, int *d_new){ int divisor = gcd(n,d); if (divisor != 1) { *n_new = n / divisor; *d_new = d / divisor;} else { *n_new = n; *d_new = d; } return 0; } /* Binary expansion can be : * finite - decimal ratio number with even denominator which is a power of 2. Note that it has also equal infinite preperiodic representation * infinite ** periodic : preperiod = 0, period > 0, rational number with odd denominator which is not a power of 2 ** preperiodic ( = eventually periodic) : preperiod > 0, period > 0, rational number with even denominator ** aperiodic : preperiod = 0, period = 0 ( or infinity), irrational number input is a rational number t = n/d so here are only 2 types of result: * 0 = preperiodic ( = eventually periodic) : preperiod > 0, period > 0, rational number with even denominator * 1 = periodic : preperiod = 0, period > 0, rational number with odd denominator which is not a power of 2 */ int give_dynamic_type(const int n, const int d){ // boolean flag int is_odd = 0; int is_POT = 0; if (d % 2) { // d % 2 > 0 fprintf(stdout, "denominator %d is odd\n", d); // parzysty is_odd = 1; } else { // d % 2 = 0 fprintf(stdout, "denominator %d is even\n", d); // nieparzysty is_odd = 0; } if (d>0 && ((d & (d-1)) == 0)) { fprintf(stdout, "denominator %d is power of 2 ( POT) = 2^%d\n", d, (int) log2(d)); is_POT = 1;} else { fprintf(stdout, "denominator %d is not the power of 2\n", d); is_POT = 0;} // wikibooks: Fractals/Mathematics/binary#preperiod_and_period if ( is_odd == 0 && is_POT == 1) { fprintf(stdout, "Binary expansion of %d/%d is finite and has equal infinite preperiodic representation\n", n,d); //fprintf(stdout, "denominator is even and POT\n"); return 0; } // preperiodic ( = eventually periodic) : preperiod > 0, period > 0, rational number with even denominator if (is_odd == 0 && is_POT == 0) { fprintf(stdout, "Binary expansion of %d/%d is preperiodic\n", n,d); //fprintf(stdout, "denominator is even and not POT\n"); return 0; } // periodic : preperiod = 0, period > 0, rational number with odd denominator which is not a power of 2 if (is_odd == 1 && is_POT == 0) { fprintf(stdout, "Binary expansion of %d/%d is periodic\n", n,d); //fprintf(stdout, "denominator is not even (odd) and not POT\n"); return 1; } return 0; } int give_period(const int n0, const int d0){ int n = n0; int d = d0; int i; int iMax = 20; // period_max int period = 0; // printf(" i\t numerator/denominator \n"); // header for ( i=0 ; i< iMax; i++){ // printf("%3d:\t %3d / %3d \n", i, n, d); // check signed integer overflow before it will happen if ( n > nMax ) {printf(" signed integer overflow will happen = wrap; use extended precision \n"); return 1;} // angle doubling map modulo 1 n = 2*n; n = n % d; // if (n == n0) { period = i+1; // printf(" preperiod = 0 and period of fraction under doubling map is %d\n", period); return period;} } return 0; // period not found, maybe period > iMax } int give_periodic_orbit(const int period, int n, int d, int orbit[]){ int i; int iMax = period; // for ( i=0 ; i< iMax; i++){ orbit[i] = n; // check signed integer overflow before it will happen if ( n > nMax ) {printf(" signed integer overflow will happen = wrap; use extended precision \n"); return 1;} // angle doubling map modulo 1 n = 2*n; n = n % d; } return 0; } int give_preperiod(const int period, const int n0, const int d0, int orbit[]){ int n = n0; int d = d0; int i; int iMax = period; // for ( i=0 ; i< iMax; i++){ if (orbit[i] == n) return i; // angle doubling map modulo 1 n = 2*n; n = n % d; } return 0; } /* input : reduced proper rational fraction in t = n0/d0 output * perperiod as a result * period as a pointer */ int give_preperiod_and_period(const int n0, const int d0, int *period){ int n = n0; int d = d0; *period = 0; int preperiod = 0; int preperiod_max = 1000; if ( preperiod_max > INT_MAX - *period ){printf("give_preperiod_and_period error: preperiod_max > INT_MAX - period, signed integer overflow will happen = ; use extended precision \n"); return -1;} int i; int iMax = preperiod_max; // preperiod_max // go thru preperiodic part to periodic part for ( i=0 ; i< iMax; i++){ //printf("%3d:\t %3d / %3d \n", i, n, d); // check signed integer overflow before it will happen if ( n > nMax ) {printf("give_preperiod_and_period error: signed integer overflow will happen = wrap; use extended precision \n"); return 1;} // angle doubling map modulo 1 n = 2*n; n = n % d; } // now it should be periodic part *period = give_period (n,d); // periodic orbit int *orbit; // only numerators orbit = (int *) malloc( *period * sizeof(* orbit)); give_periodic_orbit(*period, n, d, orbit); preperiod = give_preperiod( *period, n0, d0, orbit); free(orbit); return preperiod; } void give_orbit(const int n0, const int d0, const int preperiod, const int period){ int n = n0; int d = d0; int i; int iMax = preperiod+period; // preperiod_max for ( i=0 ; i< iMax; i++){ if ( i < preperiod) { printf("%3d:\t %3d / %3d \t preperiodic part \n", i, n, d);} else {printf("%3d:\t %3d / %3d \t periodic part \n", i, n, d);} // angle doubling map modulo 1 n = 2*n; n = n % d; } } /* Algorithm:[36] Multiply the input decimal fraction by two from above result take integer part as the binary digit take the fractional part as the starting point for the next step repeat until you either get to 0 or a periodic number read the number starting from the top - the first binary digit is the first digit after the comma */ void print_binary_expansion(const int n0, const int d0, const int preperiod, const int period){ int n = n0; int d = d0; int i_max = preperiod+period; int i; double t = (double) n/d; double t_fractional; double t_integer; int binary_digit; printf("formated infinite binary expansion of %d/%d is 0.", n0,d0); for (i = 0; i < i_max; ++i){ // doubling map t *= 2.0; // split t_fractional = modf(t, &t_integer); // binary_digit = (int) t_integer; if (i== preperiod) {printf("(");} printf("%d", binary_digit); // take the fractional part as the starting point for the next step t = t_fractional; } printf(")\n"); } int describe_fraction(const int n0, const int d0){ // proper decimal fraction // t = n/d //int n0 = 1; // //int d0 = 66; // // tr = n0r/d0r = t after reduction int n0r ; // int d0r ; // int n; int d; int period = 0; int preperiod = 0; assert(n0 > 0); assert(d0 > 1); assert(n0 < d0); // proper fraction // ------------ Reducing Fractions to Lowest Terms ------------------------------- // ----------- wikipedia Irreducible_fraction ---------------------------- give_reduced_fraction(n0, d0, &n0r, &d0r); if (n0 == n0r && d0 ==d0r ) {printf(" fraction = %d/%d\t is irreducible = in lowest terms \n", n0, d0 ); } else {printf(" fraction = %d/%d\t after reduction is %d/%d \n", n0, d0, n0r,d0r); } n = n0r; d = d0r; assert(n > 0); assert(d > 1); assert(n < d); int type = give_dynamic_type(n,d); // ------------------compute preperiod and period under doubling map ------------------------- printf("fraction %d/%d under doubling map has: \n", n0r, d0r); if (type ==0 ){ preperiod = give_preperiod_and_period( n0r, d0r, &period); if (preperiod > 0) { printf("\t preperiod = %d and period = %d\n", preperiod, period); if (period == 0 ) {printf("period = 0 means period NOT FOUND !!!\n\t maybe period > iMax in give_period \n\tor maybe preperiod_max in give_preperiod_and_period is to big \n");}}} // -------------- else { period = give_period(n,d); if (period > 0) {printf(" preperiod = 0 and period = %d\n", period); } else {printf(" preperiod = 0 and period of fraction under doubling map > 0 but period NOT FOUND !!! maybe period > iMax in give_period \n");}} // ------------------------------------------------- give_orbit( n0, d0, preperiod, period); // ----------formated infinite binary expansion --------------------- print_binary_expansion( n0r, d0r, preperiod, period); return 0; } int main(void) { nMax = INT_MAX / 2; // signed integer describe_fraction(1,22); return 0; } precision What number type precision do I need for computing doubling map? "the doubling map loses one bit of precision with every iterate."Chaos and Fractals by Mark McClure. 2.11 A few notes on computation storing the angle as an exact rational (e.g. mpq_t using libgmp), allows it to be accurate for more than 53 iterations (the number of bits of accuracy of double precision floating point). Convert to double only when needing to do sin and cos for finding the target point. Double is plenty accurate for target point. Claude Heiland-Allenfractalforums.org : parameter-external-ray-using-newton-method See also binary fraction References