G
Gregory Ewing
Chris said:
It's not constant space, because the nines counter
can grow infinitely large.
It's not constant space, because the nines counter
can grow infinitely large.
I
Ian Kelly
The code for that looks like this:
def cf_sqrt(n):
"""Yield the terms of the square root of n as a continued fraction."""
m = 0
d = 1
a = a0 = floor_sqrt(n)
while True:
yield a
next_m = d * a - m
next_d = (n - next_m * next_m) // d
if next_d == 0:
break
next_a = (a0 + next_m) // next_d
m, d, a = next_m, next_d, next_a
def floor_sqrt(n):
"""Return the integer part of the square root of n."""
n = int(n)
if n == 0: return 0
lower = 2 ** int(math.log(n, 2) // 2)
upper = lower * 2
while upper - lower > 1:
mid = (upper + lower) // 2
if n < mid * mid:
upper = mid
else:
lower = mid
return lower
The floor_sqrt function is merely doing a simple binary search and
could probably be optimized, but then it's only called once during
initialization anyway. The meat of the loop, as you can see, is just
a constant amount of integer arithmetic. If it were desired to halt
once the continued fraction starts to repeat, that would just be a
matter of checking whether the triple (m, d, a) has been seen already.
Going back to your example of adding generated digits though, I don't
know how to add two continued fractions together without evaluating
them.
I
Ian Kelly
Sorry, all that "next" business is totally unnecessary. More simply:
def cf_sqrt(n):
"""Yield the terms of the square root of n as a continued fraction."""
m = 0
d = 1
a = a0 = floor_sqrt(n)
while True:
yield a
m = d * a - m
d = (n - m * m) // d
if d == 0:
break
a = (a0 + m) // d
C
Chris Angelico
Okay, okay, technically yes. But the counter can go a long way up
before it takes up any additional space, so all that's really saying
is that this has a major flaw with anything that produces a long
stream of nines. It can't tell the difference between
..99999999999999998 and .999999999999999999999999[11] where the 11
suddenly carries and it has to flip it all back.
Anyway, that was like ten minutes' knocking-together work, you can't
expect it to be perfect. I'm amazed it even worked.
ChrisA
M
Marko Rauhamaa
Chris Angelico said:
I don't know why the constant space/time requirement is crucial. Anyway,
producing more digits simple: <URL: http://nrich.maths.org/5955>.
I believe producing the nth digit is O(n) in time and space.
Still, there's more to arithmetics than that. For example, if you have
two generated decimal expansions, you don't have an effective algorithm
to generate the decimal expansion of their ratio. That's because there's
no effective algorithm to decide if a < b.
Marko
C
Chris Angelico
The reason for striving for constant space/time is because the obvious
method (cut-and-try) is already O(n) for the nth digit, which means
it's quadratic on the number of digits desired. That gets pretty
nasty. So what I was asking was: By representing values as continued
fractions rather than as decimal digits, are you able to perform a
straight-forward transformation that produces the square root, in
constant time (which means linear in the length)? And I guess the
answer's no. CF representation doesn't have the advantage I was
wondering about.
ChrisA
O
Oscar Benjamin
I don't quite follow your reasoning here. By "cut-and-try" do you mean
bisection? If so it gives the first N decimal digits in N*log2(10)
iterations. However each iteration requires a multiply and when the
number of digits N becomes large the multiplication is worse than
linear. So the result is something like N**2 log(N)log(log(N)),
To me the obvious method is Newton iteration which takes O(sqrt(N))
iterations to obtain N digits of precision. This brings the above
complexity below quadratic:
#!/usr/bin/env python
from decimal import Decimal as D, localcontext
def sqrt(y, prec=1000):
'''Solve x**2 = y'''
assert y > 0
eps = D(10) ** -(prec + 5)
x = D(y)
with localcontext() as ctx:
ctx.prec = prec + 10
while x ** 2 - y > x * eps:
x = (x + y/x) / 2
return x
print(sqrt(2))
Some modification would be required to handle a situation where it
ends in a run of nines or zeros if you really care about the exact
digits rather than having a bounded error.
Oscar
M
Marko Rauhamaa
Oscar Benjamin said:
At a quick glance, I believe x ** 2 is O(N²) and so the total complexity
should be O(N ** 2.5).
Marko
C
Chris Angelico
By "cut and try" I'm talking about the really REALLY simple algorithm
for calculating square roots. It's basically brute force.
epsilon = 0.0001
def sqrt(n):
guess1, guess2 = 1, n
while abs(guess1-guess2) > epsilon:
guess1 = n/guess2
guess2 = (guess1 + guess2)/2
return guess1
It's generally going to take roughly O(n*n) time to generate n digits,
give or take. That's the baseline against which anything else can be
compared. There are plenty of better ways to calculate them.
ChrisA
O
Oscar Benjamin
That's the exact same algorithm I showed! How on earth would you call
that brute force?
It does not take O(n*n) time. This is Newton iteration and for
well-behaved problems such as this it generates more than n digits
after n iterations. I modified my code to show the error (x**2 - y) at
each iteration:
$ python3.3 root.py
2
0.2
0.007
0.000006
5E-12
3E-24
8E-49
8E-98
8E-196
9E-392
1E-783
The number of correct digits doubles at each iteration so after n
iterations you have 2**n digits (I misstated this as n**2 before).
This means that it takes log(N) iterations to get N digits. See here
for more:
http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
See also the section below that:
http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation
Such as?
Oscar
O
Oscar Benjamin
x**2 is just a multiplication which can be done in better than O(N**2):
http://en.wikipedia.org/wiki/Multip...st_multiplication_algorithms_for_large_inputs
Oscar
C
Chris Angelico
It uses a lot of division. There are various refinements that can be
done that replace some of that division with multiplication, but I'd
have to go do some research to figure it out.
This is the purest form of attempted-division algorithm. If you're
describing it on a blackboard, you would write it pretty much like
this. At each iteration, you have to divide by a number that's n
digits long, and then do some additional arithmetic.
It seems I'm partly mistaken, though not entirely. Let's compare two
versions. In the first, you set the precision (I'm talking in terms of
REXX's "NUMERIC DIGITS" statement - anything beyond this many digits
will be rounded (and represented exponentially, if necessary); I'm not
sure if decimal.Decimal precision works this way) such that you get 10
digits. Each iteration requires division by a 10-digit number, which
is an operation that takes a certain amount of time; and it's going to
take some number of iterations to get to the final answer.
Second version, you set the precision so you get 20 digits. Now, it's
going to take you approximately one more iteration to get to the final
answer. (This bit I was mistaken on. I thought it would take something
like 25% more or 50% more iterations.) But each iteration will take
longer. The complexity of division depends on the algorithm - grade
school long division would be O(n) with a fixed-length dividend, I
think, but you could probably pull that down a bit.
So that's why I said it'd be very roughly O(n*n) - because the
division in each step is O(n), and I thought it'd take O(n) steps.
Turns out it's O(n*log(n)), which is a lot better.
Improved versions of the above, and I was under the impression that
there were some completely different techniques that converged a lot
more quickly. But it may be that I was mistaken, as I'd been expecting
this to converge in O(n) steps. Reducing the amount of division would
speed things up significantly, although it probably won't change the
algorithmic complexity.
So, put it down to misremembering and thinking the simple algorithm
was worse than it actually is
ChrisA
O
Oscar Benjamin
There's a description of such a method here:
http://en.wikipedia.org/wiki/Method...Iterative_methods_for_reciprocal_square_roots
I don't know whether that would work out faster (when using decimal -
for float it would probably be slower).
I have no idea what that is.
With the decimal module if you set the precision to 5 digits then it
basically represents the number in "standard form" with 5 digits .e.g:
1.2345 x 10**21.
If we're talking about 10-20 digits then the decimal module is
overkill: just use float. The speed up from hardware arithmetic will
massively out-weigh any other performance considerations. My version
was intended to produce large numbers of digits which is when the
big-O comes in:
$ python3.3 -m timeit -s 'from root import sqrt' 'sqrt(2, 10)'
10000 loops, best of 3: 22.4 usec per loop
$ python3.3 -m timeit -s 'from root import sqrt' 'sqrt(2, 100)'
10000 loops, best of 3: 59.1 usec per loop
$ python3.3 -m timeit -s 'from root import sqrt' 'sqrt(2, 1000)'
1000 loops, best of 3: 1.15 msec per loop
$ python3.3 -m timeit -s 'from root import sqrt' 'sqrt(2, 10000)'
10 loops, best of 3: 85.9 msec per loop
$ python3.3 -m timeit -s 'from root import sqrt' 'sqrt(2, 100000)'
10 loops, best of 3: 1.59 sec per loop
Oscar
C
Chris Angelico
That's how NUMERIC DIGITS works, so we're on the same page. I'm not
familiar enough with decimal.Decimal and how precision is configured,
but it seems to function the same way.
Yeah, I'm just digging into the algorithm. The same concept applies
when going from 100 to 200 digits, or 1000 to 2000, and in each case,
the division will get way slower, but the number of iterations won't
go up as fast as I thought it would.
In theory, it should be possible to do the first few divisions at
lower precision, and scale up as you have need. In practice, would the
churning of precisions mean that you lose all the benefit?
ChrisA
D
Dave Angel
Oscar Benjamin said:
One problem with complexity claims is that it's easy to miss some
contributing time eaters. I haven't done any measuring on modern
machines nor in python, but I'd assume that multiplies take
*much* longer for large integers, and that divides are much
worse. So counting iterations isn't the whole story.
On the assumption that division by 2 is very fast, and that a
general multiply isn't too bad, you could improve on Newton by
observing that the slope is 2.
err = n - guess * guess
guess += err/2
Some 37 years ago I microcoded a math package which included
square root. All the math was decimal, and there was no hardware
multiply or divide. The algorithm I came up with generated the
answer one digit at a time, with no subsequent rounding needed.
And it took just a little less time than two divides. For that
architecture, Newton's method would've been too
slow.
Incidentally, the algorithm did no divides, not even by 2. No
multiplies either. Just repeated subtraction, sorta like divide
was done.
If anyone is curious, I'll be glad to describe the algorithm;
I've never seen it published, before or since. I got my
inspiration from a method used in mechanical, non-motorized,
adding machines. My father had shown me that approach in the
50's.
A
Albert van der Horst
Adding cf's adds all computable numbers in infinite precision.
However that is not even a drop in the ocean, as the computable
numbers have measure zero.
A cf object yielding its coefficients amounts to a program that generates
an infinite amount of data (in infinite time), so it is not
very surprising it can represent any computable number.
Pretty humbling really.
Groetjes Albert
A
Albert van der Horst
That is highly non-trivial indeed. See the gosper.txt reference
I gave in another post.
Groetjes Albert
A
Albert van der Horst
You're very much off the track here. A Turing machine is an abstraction
for a computer were the limitations of size are gone.
The most obvious feature of a Turing machine is an infinite tape.
A Turing machine happily calculates Ackerman functions long after
a real machine runs out of memory to represent it, with as a result
a number of ones on that tape.
But it only happens in the mathematicians mind.
S
Steven D'Aprano
On the other hand, it's not really clear that the non-computable numbers
are useful or necessary for anything. They exist as mathematical
abstractions, but they'll never be the result of any calculation or
measurement that anyone might do.
R
Rustom Mody
There are even more extreme versions of this amounting to roughly this view:
"Any infinity supposedly 'larger' than the natural numbers is a nonsensical notion."
See eg
http://en.wikipedia.org/wiki/Controversy_over_Cantor's_theory
and Weyl/Polya bet (pg 10 of http://research.microsoft.com/en-us/um/people/gurevich/Opera/123.pdf )
I cannot find the exact quote so from memory Weyl says something to this effect:
Cantor's diagonalization PROOF is not in question.
Its CONCLUSION very much is.
The classical/platonic mathematician (subject to wooly thinking) concludes that
the real numbers are a superset of the integers
The constructvist mathematician (who supposedly thinks clearly) only concludes
the obvious, viz that real numbers cannot be enumerated
To go from 'cannot be enumerated' to 'is a proper superset of' requires the
assumption of 'completed infinities' and that is not math but theology