定义

以周为单位测量的角度
倍增映射

角度倍增映射^[1]^[2]

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

输入：真分数 $r$ 的形式

十进制比率（如果 r 不是无理数）
十进制数字序列和基数
二进制数字序列和基数

输出：真分数

二进制分数上的倍增映射

倍增映射 d 对分数 x 的二进制表示的影响是简单地移位 x 的位到左侧，丢弃移位到一位的位 = 左移。^[3]

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

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

示例

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

轨道

映射的迭代给出轨道（角度序列）

周期性（输入是有奇数分母的有理数 $d=2^{p}-1$ ）
前周期性（输入是有偶数分母的有理数）
非周期性（当输入是无理数时）

以上所有类型都是无限的。有限二进制展开具有相等的无限前周期表示。

周期性

角度在倍增映射下的周期轨道

注意，这里连接单位圆上 2 个点 z1 和 z2 的弦意味着 $z_{2}=z_{1}^{2}$ 。它并不意味着这些点是同一条射线的着陆点。

有些轨道不会交叉

周期为 2 的轨道（角度在倍增映射下）
周期为 3 的轨道（角度在倍增映射下）
周期为 4 的轨道（角度在倍增映射下）
周期为 5 的轨道（角度在倍增映射下）
周期为 6 的轨道（角度在倍增映射下）
周期 9

但有些轨道会交叉

11/63 在倍增映射下的周期为 6 的轨道
74/511 在倍增映射下的周期为 9 的轨道

角度 $r={\frac {n}{d}}$ 的周期轨道，其中分母 d 为奇数

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

其中

n 是分数的分子 = 从 0 到 (d-1) 的整数
d 是分数的分母 = 正整数
p 是周期 = 正整数

周期为 1（1 个轨道）的角度 = ${\frac {n}{2^{1}-1}}={\frac {n}{1}}$

{0/1 = 1/1 }

周期 2 (1 轨道) 对应角度 = ${\frac {n}{2^{2}-1}}={\frac {n}{3}}$

{1/3, 2/3}

周期 3 (2 轨道) 对应角度 = ${\frac {n}{2^{3}-1}}={\frac {n}{7}}$

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

周期 4 (3 轨道) 对应角度 = ${\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 }

周期 5 (6 轨道) 对应角度 = ${\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}

周期 6 (9 轨道) 对应角度 = ${\frac {n}{2^{6}-1}}={\frac {n}{63}}$ , 仅分子

{ 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]

周期 = 7 对应角度 = ${\frac {n}{2^{7}-1}}={\frac {n}{127}}$

周期 = 8 对应角度 = ${\frac {n}{2^{8}-1}}={\frac {n}{255}}$

周期 = 9 对应角度 = ${\frac {n}{2^{9}-1}}={\frac {n}{511}}$

参见

参数平面上的根点（键点）的外部射线
动态平面
- 落在抛物固定点上的外部射线。（第一个角度是尾流的角度）

前周期

这里分数 r 的分母 d

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

是偶数

d=2^{t}q

并且 q 是奇数。

这里

t 是前周期

示例轨道

前周期 1 和周期 1 对应角度 ${\frac {n}{2}}$

{1/2 , (0/2)}

参见

Misiurewicz 点在参数平面上
动态平面上的树枝状朱利亚集

非周期

比较

实现

整数
- 2 个整数（分子和分母）
- 有理数
  - GMP 库中的 mpq 类型
浮点数
- double

向前

使用 32 位有符号整数将最大前周期限制在约 30
使用 double 表示分数将最大精确位数限制在约 53 位

C 使用整数

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

使用双精度数可以得到精度为 53 位二进制数扩展（双精度浮点数的精度位数）。

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

使用 GMP 库中的有理数类型 mpq 的 C 语言（任意精度）

C 函数（使用 GMP 库）^[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;
}

编译

gcc r.c -lgmp

运行

./a.out

结果

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 函数

doubling_map(n,d):=mod(2*n,d)/d $

或使用有理数作为输入

DoublingMap(r):=
  block([d,n],
        n:ratnumer(r),
        d:ratdenom(r),
        mod(2*n,d)/d)$

Common Lisp 函数

(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 函数^[5]

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

倍增映射的反函数

每个角度 α ∈ R/Z（以转数表示）都具有

一个像 = 2α mod 1 在倍增映射下
“在倍增映射下有两个原像：α/2 和 (α + 1)/2。”^[6] 倍增映射的反函数是一个多值函数。

在 Maxima CAS 中

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

请注意这两个原像之间的差

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

是半圈 = 180 度 = Pi 弧度。

倍增映射 d 下的像和原像
$\alpha$	$d^{1}(\alpha )$	$d^{-1}(\alpha )$
${\frac {1}{2}}$	${\frac {1}{1}}$	$\left\{{\frac {1}{4}},{\frac {3}{4}}\right\}$
${\frac {1}{3}}$	${\frac {2}{3}}$	$\left\{{\frac {1}{6}},{\frac {4}{6}}\right\}$
${\frac {1}{4}}$	${\frac {1}{2}}$	$\left\{{\frac {1}{8}},{\frac {5}{8}}\right\}$
${\frac {1}{5}}$	${\frac {2}{5}}$	$\left\{{\frac {1}{10}},{\frac {6}{10}}\right\}$
${\frac {1}{6}}$	${\frac {1}{3}}$	$\left\{{\frac {1}{12}},{\frac {7}{12}}\right\}$
${\frac {1}{7}}$	${\frac {2}{7}}$	$\left\{{\frac {1}{14}},{\frac {4}{7}}\right\}$

问答

分析

查找扩展类型
查找前周期和周期
转换十进制小数为二进制展开式
格式化展开式

如何查找扩展类型

十进制小数的二进制展开式

/*
 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://wikibooks.cn/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;
}

如何查找前周期和周期？

周期和倍增映射

如何查找周期角度下的倍增映射

视觉
数值
- 从十进制小数的分母（约化有理分数 m/n）中读取周期
- 在二进制展开式（二进制序列）中查找周期/前周期
  - 字符串匹配算法，例如，朴素或暴力搜索，Knuth-Morris-Pratt 算法^[7]
  - 离散傅里叶变换
- 从角度在倍增映射下的行程中读取它

约化有理分数 m/n 的二进制展开式的周期等于2 模n 的乘法阶

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

角度 5/31 在倍增映射下的轨道

C 版本

双精度：正向和逆向倍增映射

/*


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 [email protected]: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;}

任意精度

// 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++ 版本

/*
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 版本

(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 版本

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 版本

Haskell 版本^[8]

从整数类型（Int 或 Integer）到任何其他类型的转换通过“fromIntegral”完成。目标类型会自动推断

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

转换

所有十进制小数都能精确地转换为二进制吗？

并非所有小数都能精确地转换为二进制。只有分母为 2 的幂 (有限) 的那些小数才具有精确的十进制表示。“在其他所有情况下，表示中都会存在误差。误差的大小取决于用于表示它的位数。” ^[9]

十进制转换为二进制

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

精度

我需要什么类型精度的数字来计算倍增映射？

“倍增映射每次迭代都会损失一位精度。”^[10]
将角度存储为精确的有理数 (例如，使用 libgmp 的 mpq_t) 使其能够精确到超过 53 次迭代 (双精度浮点数的精度位数)。仅在需要使用 sin 和 cos 来寻找目标点时才转换为 double。double 对目标点来说已经足够精确。Claude Heiland-Allen^[11]

参见

二进制小数

参考文献

[1] 关于映射 x↦kx (mod Z) 周期轨道的组合类型。第一部分：实现作者：Carsten L. Petersen，Saeed Zakeri

[2] Mark McClure 的 chaos_and_doubling 映射

[3] Mark McClure 的倍增映射

[4] 波兰维基教科书中的 UNIX 系统编程：GMP

[5] Claude Heiland-Allen 使用 Haskell 和 SVG 输出的 lavaurs 算法

[6] 符号动力学和自相似群，作者：VOLODYMYR NEKRASHEVYCH

[7] th stackexchange 问题：有限二进制序列的周期

[8] Claude Heiland-Allen 使用 Haskell 和 SVG 输出的 lavaurs 算法

[9] Davide Borchia 的二进制小数计算器

[10] Mark McClure 的混沌与分形。2.11 关于计算的一些说明

[11] ractalforums.org : 使用牛顿法求解参数外部射线

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]