# Definition

Angle doubling map[1][2]

${\displaystyle T(r)=(2r){\bmod {1}}}$

Input: proper fraction ${\displaystyle r}$  in forms

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.[3]

${\displaystyle x=0.b_{1}b_{2}b_{3}\dots .}$

${\displaystyle d(x)=0.b_{2}b_{3}\dots .}$

Example

• ${\displaystyle 2*{\frac {11}{31}}=2*0.01011=0.10110={\frac {22}{31}}}$

# Orbits

Iteration of the map gives orbit ( sequence of angles):

• periodic ( input is a rational number with odd denominator ${\displaystyle d=2^{p}-1}$ )
• 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 ${\displaystyle z_{2}=z_{1}^{2}}$ . 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 ${\displaystyle r={\frac {n}{d}}}$  with odd denominator d

${\displaystyle r={\frac {n}{d}}={\frac {n}{2^{p}-1}}}$



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 = ${\displaystyle {\frac {n}{2^{1}-1}}={\frac {n}{1}}}$

{0/1 = 1/1 }



Period 2 ( 1 orbit) for angles = ${\displaystyle {\frac {n}{2^{2}-1}}={\frac {n}{3}}}$

{1/3, 2/3}


Period 3 ( 2 orbits) for angles = ${\displaystyle {\frac {n}{2^{3}-1}}={\frac {n}{7}}}$

{1/7, 2/7, 4/7 }
{3/7, 6/7, 5/7 }


Period 4 ( 3 orbits) for angles = ${\displaystyle {\frac {n}{2^{4}-1}}={\frac {n}{15}}}$

{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 = ${\displaystyle {\frac {n}{2^{5}-1}}={\frac {n}{31}}}$

{ 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 = ${\displaystyle {\frac {n}{2^{6}-1}}={\frac {n}{63}}}$ , 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 = ${\displaystyle {\frac {n}{2^{7}-1}}={\frac {n}{127}}}$

Period = 8 for angles = ${\displaystyle {\frac {n}{2^{8}-1}}={\frac {n}{255}}}$

Period = 9 for angles = ${\displaystyle {\frac {n}{2^{9}-1}}={\frac {n}{511}}}$

## preperiodic

Here denominator d of fraction r

${\displaystyle r={\frac {n}{d}}}$

is even:

${\displaystyle d=2^{t}q}$

and q is odd.

Here

• t is preperiod

Example orbits:

Preperiod 1 and period 1 for angles ${\displaystyle {\frac {n}{2}}}$

{1/2 , (0/2)}


See

Compare

# 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;}

}
}

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)[4]:

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


-- 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:

In Maxima CAS:

InvDoublingMap(r):= [r/2, (r+1)/2];


Note that difference between these 2 preimages

${\displaystyle {\frac {\alpha }{2}}-{\frac {\alpha +1}{2}}={\frac {1}{2}}}$

is half a turn = 180 degrees = Pi radians.

 ${\displaystyle \alpha }$ ${\displaystyle d^{1}(\alpha )}$ ${\displaystyle d^{-1}(\alpha )}$ ${\displaystyle {\frac {1}{2}}}$ ${\displaystyle {\frac {1}{1}}}$ ${\displaystyle \left\{{\frac {1}{4}},{\frac {3}{4}}\right\}}$ ${\displaystyle {\frac {1}{3}}}$ ${\displaystyle {\frac {2}{3}}}$ ${\displaystyle \left\{{\frac {1}{6}},{\frac {4}{6}}\right\}}$ ${\displaystyle {\frac {1}{4}}}$ ${\displaystyle {\frac {1}{2}}}$ ${\displaystyle \left\{{\frac {1}{8}},{\frac {5}{8}}\right\}}$ ${\displaystyle {\frac {1}{5}}}$ ${\displaystyle {\frac {2}{5}}}$ ${\displaystyle \left\{{\frac {1}{10}},{\frac {6}{10}}\right\}}$ ${\displaystyle {\frac {1}{6}}}$ ${\displaystyle {\frac {1}{3}}}$ ${\displaystyle \left\{{\frac {1}{12}},{\frac {7}{12}}\right\}}$ ${\displaystyle {\frac {1}{7}}}$ ${\displaystyle {\frac {2}{7}}}$ ${\displaystyle \left\{{\frac {1}{14}},{\frac {4}{7}}\right\}}$

# 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 ?

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)
• 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:

${\displaystyle Period_{2}(m/n)={ord}_{n}(2)}$

### 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

commons :Binary_decomposition_of_dynamic_plane_for_f0(z)_%3D_z%5E2.png

---------- git --------------------
cd existing_folder
git init
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);
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
);



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

### 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." [9]

### 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;}

}
}

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."[10]
• 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-Allen[11]