1 /*-------------- Telecommunications & Signal Processing Lab ---------------
2 McGill University
3
4 Routine:
5 void MSratio (double Val, long int *N, long int *D, double tol,
6 long int MaxN, long int MaxD)
7
8 Purpose:
9 Find a ratio of integers approximation to a floating point value
10
11 Description:
12 A continued fraction expansion approximation to a given floating point
13 value is calculated iteratively. At successive stages in the expansion,
14 the approximation error for the convergents alternates in sign but with
15 diminishing magnitude. The expansion is terminated when either two
16 conditions is satisfied.
17 1: The magnitude of either the numerator or the denominator of the
18 approximation is larger than a given limit. If so, the secondary
19 convergents are tested to see if one provides a smaller error than
20 backing off to the convergent from the previous iteration.
21 2: The absolute value of the difference between the approximation (N/D) and
22 the given number is less than a specified value (tol).
23 If the input value is exactly a ratio of integers, this routine will find
24 those integers if tol is set to zero, and MaxN and MaxD are appropriately
25 large.
26
27 Consider that case that the iteration is stopped by either the numerator
28 exceeding MaxN or the denominator exceeding MaxD. The resulting fraction
29 N/D is the fraction nearest Val amongst those fractions with numerator and
30 denominator not exceeding the given limits.
31
32 References:
33 R. B. J. T. Allenby and E. J. Redfern, "Introduction to Number Theory with
34 Computing", Edward Arnold, 1989.
35
36 I. Niven and H. S. Zuckerman, "An Introduction to the Theory of Numbers",
37 4th ed., John Wiley & Sons, 1980.
38
39 Parameters:
40 -> double Val
41 Value to be approximated. This value should have a magnitude between
42 1/MaxD and MaxN.
43 <- long int N
44 Numerator in the rational approximation to Val. N is negative if Val is
45 negative.
46 <- long int D
47 Denominator in the rational approximation to Val.
48 -> double tol
49 Error tolerance used to terminate the approximation
50 -> long int MaxN
51 Maximum value for N. This value should be at least floor(|Val|).
52 -> long int MaxD
53 Maximum value for D. This value should be at least floor(1/|Val|).
54
55 Author / revision:
56 P. Kabal Copyright (C) 2003
57 $Revision: 1.7 $ $Date: 2003/05/09 02:29:38 $
58
59 -------------------------------------------------------------------------*/
60
61 #include <math.h>
62
63 #include <libtsp.h>
64
65 #define ABSV(x) (((x) < 0) ? -(x) : (x))
66
67
68 void
MSratio(double Val,long int * N,long int * D,double tol,long int MaxN,long int MaxD)69 MSratio (double Val, long int *N, long int *D, double tol, long int MaxN,
70 long int MaxD)
71
72 {
73 double pk, pkm1, pkm2, qk, qkm1, qkm2;
74 double AV, AVS, ak;
75 double rk, sk;
76 double m, mp, mq;
77
78 /*
79 The continued fraction expansion of a value x is of the form
80 x = a[0] + 1
81 -----------------
82 a[1] + 1
83 ----------
84 a[2] + ...
85
86 = {a[0]; a[1], a[2], ... }
87 The a[k] are integers (all positive, except possibly for a[0]). The value
88 a[0] can be obtained from x,
89 a[0] = floor(x).
90 The continued fraction term represents x1 = 1/(x - floor(x)),
91 a[1] = floor(x1).
92 This process can be continued until the remainder is zero (x is exactly a
93 a ratio of integers) or indefinitely (x is irrational).
94
95 The truncated continued fraction expansion approximates x, with the error
96 decreasing as more terms are added. The value of the continued fraction with
97 terms up to a[k] is the convergent C[k],
98 p[k]
99 C[k] = {a[0]; a[1], ... , a[k]} = ---- .
100 q[k]
101 The numerator and denominator satisfy the relationships (k >= 0)
102 p[k] = a[k] p[k-1] + p[k-2], with initial values p[-2] = 0, p[-1] = 1
103 q[k] = a[k] q[k-1] + q[k-2], with initial values q[-2] = 1, q[-1] = 0
104 The even convergents C[0], C[2], C[4], ... form an increasing sequence that
105 approaches the value x from below, while the odd convergents C[1], C[3}, ...
106 form a decreasing sequence which approaches x from above.
107
108 Properties of the convergent C[k] = p[k]/q[k]:
109 1: The convergent C[k] has the property that among fractions with
110 denominator less than or equal to q[k], the C[k] is the one nearest to x.
111 2: The fraction C[k] is in reduced form, i.e. p[k] is relatively prime to
112 q[k].
113
114 The secondary convergents lie between two even or two odd convergents,
115 p[n-1] p[n-1]+p[n] p[n-1]+2p[n]
116 C[n-1] = ------, -----------, ------------, ... ,
117 q[n-1] q[n-1]+q[n] q[n-1]+2q[n]
118
119 p[n-1]+a[n+1]p[n] p[n+1]
120 ----------------- = ------ = C[n+1]
121 q[n-1]+a[n+1]q[n] q[n+1]
122 The secondary convergents are either increasing (n odd) or decreasing
123 (n even). The denominators are increasing in this sequence. No other
124 fraction can be placed into this sequence and maintain these properties.
125
126 The closest fraction to x which has a denominator less than a given value M
127 is either a convergent or a secondary convergent.
128 */
129
130 /* Initialization */
131 pkm1 = 0.0;
132 pk = 1.0;
133 qkm1 = 1.0;
134 qk = 0.0; /* C[-1] = 1/0 */
135
136 AVS = ABSV (Val);
137 AV = AVS;
138
139 /* Main loop for the continued fraction expansion */
140 while (1) {
141 ak = floor (AV);
142
143 /* Recursive calculation of the numerator and denominator */
144 /* (all terms are integer-valued and positive) */
145 pkm2 = pkm1;
146 pkm1 = pk;
147 pk = ak * pkm1 + pkm2;
148 qkm2 = qkm1;
149 qkm1 = qk;
150 qk = ak * qkm1 + qkm2;
151
152 if (pk > MaxN || qk > MaxD) {
153
154 pk = pkm1; /* Back off to the previous convergent */
155 qk = qkm1;
156
157 /* Try a secondary convergent */
158 mp = floor ((MaxN - pkm2) / pkm1);
159 mq = floor ((MaxD - qkm2) / qkm1);
160 m = (mp <= mq) ? mp : mq;
161 if (m > 0) {
162 rk = m * pkm1 + pkm2;
163 sk = m * qkm1 + qkm2;
164 if (ABSV (AVS - rk/sk) < ABSV (AVS - pk/qk)) {
165 pk = rk;
166 qk = sk;
167 }
168 }
169 break;
170 }
171
172 /* Check the approximation error */
173 if (ABSV (AVS - pk / qk) <= tol)
174 break;
175
176 /* Next term
177 (AV should be positive, since the denominator is (AV - floor(AV)).
178 If it is negative, we are being limited by machine precision.)
179 */
180 AV = 1.0 / (AV - ak);
181 if (AV < 0.0)
182 break;
183 }
184
185 if (Val < 0.0)
186 pk = -pk;
187 *N = (long int) pk;
188 *D = (long int) qk;
189
190 return;
191 }
192