# Working with the set of real numbers

Dave Angel davea at davea.name
Wed Mar 5 13:57:49 CET 2014

``` Oscar Benjamin <oscar.j.benjamin at gmail.com> Wrote in message:
> On 4 March 2014 23:20, Dave Angel <davea at davea.name> wrote:
>>
>> 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.
>
> Agreed but there's a big difference between log(N) iterations and N iterations!
>
>> 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
>
> I gues you mean like this:
>
> def sqrt(n):
>     err = guess = 1
>     while err > 1e-10:
>         err = n - guess * guess
>         guess += err/2
>     return guess
>
> This requires log2(10)*N iterations to get N digits.

No idea how you came up with that,  but I see an error in my
stated algorithm,  which does surely penalize it. The slope isn't
2, but 2x.  So the line should have been
guess += err/(2*guess)

Now if you stop the loop after 3 iterations (or use some other
approach to get a low-precision estimate,  then you can calculate
scale = 1/(2*estimate)

and then for remaining iterations,
guess += err *scale

> So the penalty
> for using division would have to be extreme in order for this to
> better. Using Decimal to get many digits we can write that as:
>
> def sqrt2(n, prec=1000):
>     '''Solve x**2 = y'''
>     eps = D(10) ** -(prec + 5)
>     err = guess = D(1)
>     with localcontext() as ctx:
>         ctx.prec = prec + 10
>         while abs(err) > eps:
>             err = n - guess*guess
>             guess += err/2
>     return guess
>
> This method works out much slower than Newton with division at 10000
> digits: 40s (based on a single trial) vs 80ms (timeit result).

Well considering you did not special-case the divide by 2, I'm not
surprised it's slower.

>
>> 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.
>
> If you're working with a fixed small precision then it might be.
>
>> 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.
>
> I'm curious.
>

A later message,  I guess.  I can't write that much on the tablet.

--
DaveA

```