Your code isn't calculating the square root of n
. It attempts to calculate the continued fraction of an already calculated √n
. I mean it's OK but, had it been correct, your approach is more suitable for general decimal to rational convertion. However for a sqrt function for the regular (simple) continued fractions (where all numerators are ones) the algorithm is slightly different.
However the problem isn't finished. Yes as a rule, the CF coefficients for √n
is in the form of a repeating palindrome which ends with the double of the first nonzero coefficient. Such as √31 =[5;1,1,3,5,3,1,1,10,1,1,3,5,3,1,1,10..]
. Now there is no simple way to quess the length of the palindrome for every given n
. There are some known patterns however they are far from defining a generalized pattern for all n
. So stoping the iteration at the end of the first palindrome is a very undeterministic approach. Imagine
__
√226 =[15;30]
while
____________________________________________________
√244 =[15;1,1,1,1,1,2,1,5,1,1,9,1,6,1,9,1,1,5,1,2,1,1,1,1,1,30]
If you decide to stop the iteration at 2*f[0]
most of the time you get either a bad approximation like in √226
or an over-calculated one like in √244
case. Besides, once n
grows, chasing the end of the palindrome becomes more meaningless since you will never need such a precision.
___________________________________________________________________________________________________________________________________________________________________________
√7114 = [84;2,1,9,3,1,10,2,23,1,1,1,1,1,2,1,27,2,1,1,3,1,2,1,1,1,16,4,3,1,3,2,1,6,18,1,1,2,6,11,11,6,2,1,1,18,6,1,2,3,1,3,4,16,1,1,1,2,1,3,1,1,2,27,1,2,1,1,1,1,1,23,2,10,1,3,9,1,2,168]
In this case it would be reasonable to stop the iteration once the necessary precision is obtained. As I have mentioned at the beginning, there are two approaches.
- The general Decimal to Rational algorithm obtains simple continued fractions from any decimal number. This will give a CF resolving to exactly to that decimal without any floating point errors. For a detailed explanation of this algorithm you may check a previous answer of mine.
- There happens to be a more direct algorithm for
√n
which is basically the same as 1 but tuned for square roots. In this case you don't provide √n
but n
.
The idea is as follows. We have to define a general form for the input which involves a square root value at the numerator. Then we try to reach to the same expression in the continued fraction part to be able to iterate.
Let our input be of the form
q + √n
______
p
for simple square root operation we can assume q
is 0
and p
is 1
. If we can establish this form at the next stage then we can easily iterate on.
Starting with the initial stage where q = 0
, p = 1
, m
is the integer part of √n
and 1/x
is the floating part, our target is to bring x
into (q + √n) / p
form;
1 1 1 (√n + m) √n + m
√n = m + ___ ⇒ x = _______ ⇒ x = ________ . ________ ⇒ x = ________
x √n - m (√n - m) (√n + m) n - m^2
Now √n
is at the numerator and we have the form of;
√n + q
x = ______
p
where q = m
and p = n - m^2
. Right at this point you can calculate x
and the next m
by flooring x
. The generalized form of the algorithm becomes;
√n + q 1 p p(√n - (q - pm)) p(√n + (pm - q))
x = ______ = m + ___ ⇒ x' = ______________ = _________________ = ________________
p x' √n + (q - pm) n - (q - pm)^2 n - (q - pm)^2
at this point p
is divisible by n - (q - pm)^2
. This is now stable and we can extend it as long as we want. Let's make new assignments for q
and p
as;
q' = pm-q;
p' = (n - q'^2)/p;
√n + q'
x' = ______
p'
m' = Math.floor(x')
Note that when p'
becomes 1 (n - q'^2 = p
) we are at the end of the palindrome. However in order to decide where to stop I am using the same mechanism as described in my toRational algorithm as linked in the alternative 1 above. It basically stops once the JS floating point resolution is reached. The JavaScript code is as follows;
function rationalSqrt(n){
var nr = Math.sqrt(n),
m = Math.floor(nr),
p = n-m**2,
q = m,
cs = [m],
n0 = 1,
d0 = 0,
n1 = m,
d1 = 1,
n2 = 0,
d2 = 1;
if (nr === m) return {n:m,d:1,cs};
while (Math.abs(nr-n2/d2) > Number.EPSILON){
m = Math.floor((nr+q)/p);
q = m*p-q;
p = (n-q**2)/p;
cs.push(m);
n2 = m*n1+n0;
d2 = m*d1+d0;
n0 = n1;
d0 = d1;
n1 = n2;
d1 = d2;
}
return {n:n2,d:d2,cs};
}
These two algorithms differ slightly.
- ~60% of the time they result same fractions.
- ~27% of the time
toRational
gives smaller numerators and denominators within the JS floating point resolution. However surely it's not a better approximation had we have one more digit in JS.
- ~13% of the time
rationalSqrt
(this one) gives smaller numerators and denominators within the JS floating point resolution.
rationalSqrt
will result exact coefficients as one would expect from a square root but will be truncated once the resolution is enough.
toRational
gives expected coefficients of couse but the last one could be totally unrelated to what you would expect from the square root series.
One such example is;
rationalSqrt(511); //returns
{ n : 882184734
, d : 39025555
, cs: [22,1,1,1,1,6,1,14,4,1,21,1,4,14,1,6,1]
}
while
toRational(Math.sqrt(511));
{ n : 1215746799
, d : 53781472
, cs: [22,1,1,1,1,6,1,14,4,1,21,1,4,14,1,10]
}
Further thoughts:
- Consider that we are given the RCF coefficients. Can we reverse the
rationalSqrt
alogrithm to obtain it's(q + √n) / p
form? This
might come as an interesting task.
- Practiacally limiting the CF with the JS decimal resolution would yield a similar precision in CF arithmetics and you would be wasting your time. In order to achieve higher square root resolution I advise to replace
while (Math.abs(nr-n2/d2) > Number.EPSILON)
condition with while (cs.length < 30)
if 30 CF coefficients of resolution would be sufficient for your application. That figure is up to you.
double
. – Veridical{3, 4, 12, 3, 1}
for3.245
. So it looks like it's working correctly. – Veridical