From 70b4c2177eb9ffa9dcfff0ebeac47a5ec1622cee Mon Sep 17 00:00:00 2001 From: SheetJS Date: Tue, 24 Dec 2013 23:06:06 -0500 Subject: [PATCH] version bump 0.2.0: continued fraction method --- README.md | 7 +++++ frac.js | 31 +++++++++++++++--- frac.md | 88 ++++++++++++++++++++++++++++++++++++++++++++++++---- package.json | 2 +- 4 files changed, 116 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index 070bbe9..e14baf8 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,10 @@ Rational approximation to a floating point number with bounded denominator. Uses the Mediant Method +This module also provides an implementation of the continued fraction method as +described by Aberth in "A method for exact computation with rational numbers", +which appears to be used by spreadsheet programs for displaying fractions + ## JS Installation and Usage In node: @@ -30,3 +34,6 @@ For example: > frac(Math.PI,100) // [ 0, 22, 7 ] > frac(Math.PI,100,true) // [ 3, 1, 7 ] ``` + +`frac.cont` implements the Aberth algorithm (input and output specifications +match the original `frac` function) diff --git a/frac.js b/frac.js index 0451e66..65c5e3b 100644 --- a/frac.js +++ b/frac.js @@ -4,17 +4,38 @@ var frac = function(x, D, mixed) { if(x !== n1) while(d1 <= D && d2 <= D) { var m = (n1 + n2) / (d1 + d2); if(x === m) { - if(d1 + d2 <= D) d1+=d2, n1+=n2, d2=D+1; + if(d1 + d2 <= D) { d1+=d2; n1+=n2; d2=D+1; } else if(d1 > d2) d2=D+1; else d1=D+1; break; } - else if(x < m) n2 = n1+n2, d2 = d1+d2; - else n1 = n1+n2, d1 = d1+d2; + else if(x < m) { n2 = n1+n2; d2 = d1+d2; } + else { n1 = n1+n2; d1 = d1+d2; } } - if(d1 > D) d1 = d2, n1 = n2; + if(d1 > D) { d1 = d2; n1 = n2; } if(!mixed) return [0, n1, d1]; var q = Math.floor(n1/d1); return [q, n1 - q*d1, d1]; }; -if(typeof module !== undefined) module.exports = frac; +frac.cont = function cont(x, D, mixed) { + var sgn = x < 0 ? -1 : 1; + var B = x * sgn; + var P_2 = 0, P_1 = 1, P = 0; + var Q_2 = 1, Q_1 = 0, Q = 0; + var A = B|0; + while(Q_1 < D) { + A = B|0; + P = A * P_1 + P_2; + Q = A * Q_1 + Q_2; + if((B - A) < 0.00000001) break; + B = 1 / (B - A); + P_2 = P_1; P_1 = P; + Q_2 = Q_1; Q_1 = Q; + } + if(Q > D) { Q = Q_1; P = P_1; } + if(Q > D) { Q = Q_2; P = P_2; } + if(!mixed) return [0, sgn * P, Q]; + var q = Math.floor(sgn * P/Q); + return [q, sgn*P - q*Q, Q]; +}; +if(typeof module !== 'undefined') module.exports = frac; diff --git a/frac.md b/frac.md index 2d8f76f..f6d145f 100644 --- a/frac.md +++ b/frac.md @@ -48,7 +48,7 @@ denominator) ``` if(x === m) { - if(d1 + d2 <= D) d1+=d2, n1+=n2, d2=D+1; + if(d1 + d2 <= D) { d1+=d2; n1+=n2; d2=D+1; } else if(d1 > d2) d2=D+1; else d1=D+1; break; @@ -58,8 +58,8 @@ denominator) Otherwise shrink the range: ``` - else if(x < m) n2 = n1+n2, d2 = d1+d2; - else n1 = n1+n2, d1 = d1+d2; + else if(x < m) { n2 = n1+n2; d2 = d1+d2; } + else { n1 = n1+n2; d1 = d1+d2; } } ``` @@ -67,17 +67,93 @@ At this point, `d1 > D` or `d2 > D` (but not both -- keep track of how `d1` and `d2` change). So we merely return the desired values: ``` - if(d1 > D) d1 = d2, n1 = n2; + if(d1 > D) { d1 = d2; n1 = n2; } if(!mixed) return [0, n1, d1]; var q = Math.floor(n1/d1); return [q, n1 - q*d1, d1]; }; ``` +## Continued Fraction Method + +The continued fraction technique is employed by various spreadsheet programs. +Note that this technique is inferior to the mediant method (at least, according +to the desired goal of most accurately approximating the floating point number) + +``` +frac.cont = function cont(x, D, mixed) { +``` + +> Record the sign of x, take b0=|x|, p_{-2}=0, p_{-1}=1, q_{-2}=1, q_{-1}=0 + +Note that the variables are implicitly indexed at `k` (so `B` refers to `b_k`): + +``` + var sgn = x < 0 ? -1 : 1; + var B = x * sgn; + var P_2 = 0, P_1 = 1, P = 0; + var Q_2 = 1, Q_1 = 0, Q = 0; + var A = B|0; +``` + +> Iterate + +> ... for k = 0,1,...,K, where K is the first instance of k where +> either q_{k+1} > Q or b_{k+1} is undefined (b_k = a_k). + +``` + while(Q_1 < D) { +``` + +> a_k = [b_k], i.e., the greatest integer <= b_k + +``` + A = B|0; +``` + +> p_k = a_k p_{k-1} + p_{k-2} +> q_k = a_k q_{k-1} + q_{k-2} + +``` + P = A * P_1 + P_2; + Q = A * Q_1 + Q_2; +``` + +> b_{k+1} = (b_{k} - a_{k})^{-1} + +``` + if((B - A) < 0.00000001) break; +``` + +At the end of each iteration, advance `k` by one step: + +``` + B = 1 / (B - A); + P_2 = P_1; P_1 = P; + Q_2 = Q_1; Q_1 = Q; + } +``` + +In case we end up overstepping, walk back an iteration or two: + +``` + if(Q > D) { Q = Q_1; P = P_1; } + if(Q > D) { Q = Q_2; P = P_2; } +``` + +The final result is `r = (sgn x)p_k / q_k`: + +``` + if(!mixed) return [0, sgn * P, Q]; + var q = Math.floor(sgn * P/Q); + return [q, sgn*P - q*Q, Q]; +}; +``` + Finally we put some export jazz: ``` -if(typeof module !== undefined) module.exports = frac; +if(typeof module !== 'undefined') module.exports = frac; ``` # Tests @@ -103,7 +179,7 @@ test: ```json>package.json { "name": "frac", - "version": "0.1.0", + "version": "0.2.0", "author": "SheetJS", "description": "Rational approximation with bounded denominator", "keywords": [ "math", "fraction", "rational", "approximation" ], diff --git a/package.json b/package.json index 064e49c..6106043 100644 --- a/package.json +++ b/package.json @@ -1,6 +1,6 @@ { "name": "frac", - "version": "0.1.0", + "version": "0.2.0", "author": "SheetJS", "description": "Rational approximation with bounded denominator", "keywords": [ "math", "fraction", "rational", "approximation" ],