Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #80075

Re: Challenge: optimizing isqrt

Show all headers | View raw

Hi Steven,

let me start by answering from reverse:
 > Q3: What is the largest value of n beyond which you can never use the 
float
 > optimization?
 >

A3: There is no such value, besides the upper limit of floats (DBL_MAX~ 
10^308)

P3: If you feed a perfect square into the floating point square root 
algorithm, with a mantissa of the root of length smaller than the 
bitwidth of your float, it will always come out perfectly. I.e., 
computing sqrt(25) in FP math is no different from sqrt(25*2**200):

 >>> 25*2**200
40173451106474756888549052308529065063055074844569820882534400L
 >>> x=int(math.sqrt(25*2**200))
 >>> x
6338253001141147007483516026880L
 >>> x*x
40173451106474756888549052308529065063055074844569820882534400L
 >>>



Am 01.11.14 02:29, schrieb Steven D'Aprano:
> There is an algorithm for calculating the integer square root of any
> positive integer using only integer operations:
>
> def isqrt(n):
>      if n < 0: raise ValueError
>      if n == 0:
>          return 0
>      bits = n.bit_length()
>      a, b = divmod(bits, 2)
>      x = 2**(a+b)
>      while True:
>          y = (x + n//x)//2
>          if y >= x:
>              return x
>          x = y

> Q2: For values above M, is there a way of identifying which values of n are
> okay to use the optimized version?

A2: Do it in a different way.

Your above algorithm is obviously doing Heron- or Newton-Raphson 
iterations, so the same as with floating point math. The first line 
before the while loop computes some approximation to sqrt(n). Instead of 
doing bit shuffling, you could compute this by FP math and get closer to 
the desired result, unless the integer is too large to be represented by 
FP. Now, the terminating condition seems to rely on the fact that the 
initial estimate x>=sqrt(n), but I don't really understand it. My guess 
is that if you do x=int(sqrt(n)), then do the first iteration, then swap 
x and y such that x>y, then enter the loop, you would simply start with 
a better estimate in case that the significant bits can be represented 
by the float.

So this is my try, but not thoroughly tested:

def isqrt(n):
	if n < 0: raise ValueError
	if n == 0:
		return 0
	bits = n.bit_length()
	# the highest exponent in 64bit IEEE is 1023
	if n>2**1022:
		a, b = divmod(bits, 2)
		x = 2**(a+b)
	else:
		x=int(math.sqrt(n))
		y=n//x
		if x<y:
			x,y = (y,x)

	while True:
		y = (x + n//x)//2
		if y >= x:
			return x
		x = y


Christian

Back to comp.lang.python | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread

Thread

Challenge: optimizing isqrt Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2014-11-01 12:29 +1100
  Re: Challenge: optimizing isqrt Chris Angelico <rosuav@gmail.com> - 2014-11-01 13:12 +1100
  Re: Challenge: optimizing isqrt Christian Gollwitzer <auriocus@gmx.de> - 2014-11-01 09:02 +0100
    Re: Challenge: optimizing isqrt Christian Gollwitzer <auriocus@gmx.de> - 2014-11-01 09:09 +0100
    Re: Challenge: optimizing isqrt Chris Angelico <rosuav@gmail.com> - 2014-11-01 19:13 +1100
      Re: Challenge: optimizing isqrt Christian Gollwitzer <auriocus@gmx.de> - 2014-11-01 09:25 +0100
        Re: Challenge: optimizing isqrt Chris Angelico <rosuav@gmail.com> - 2014-11-01 19:33 +1100
          Re: Challenge: optimizing isqrt Christian Gollwitzer <auriocus@gmx.de> - 2014-11-01 09:38 +0100
            Re: Challenge: optimizing isqrt Chris Angelico <rosuav@gmail.com> - 2014-11-01 19:45 +1100
  Re: Challenge: optimizing isqrt Akira Li <4kir4.1i@gmail.com> - 2014-11-01 12:12 +0300
  Re: Challenge: optimizing isqrt Serhiy Storchaka <storchaka@gmail.com> - 2014-11-20 20:00 +0200
  Re: Challenge: optimizing isqrt Stephen Tucker <stephen_tucker@sil.org> - 2014-11-21 08:15 +0000
  Re: Challenge: optimizing isqrt Dave Angel <davea@davea.name> - 2014-11-23 10:43 -0500
  Re: Challenge: optimizing isqrt Serhiy Storchaka <storchaka@gmail.com> - 2014-11-25 21:31 +0200
  Re: Challenge: optimizing isqrt Dave Angel <davea@davea.name> - 2014-11-25 19:46 -0500
  Re: Challenge: optimizing isqrt Stephen Tucker <stephen_tucker@sil.org> - 2014-11-26 07:57 +0000
    Re: Challenge: optimizing isqrt Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2014-11-26 21:25 +1100
  Re: Challenge: optimizing isqrt Serhiy Storchaka <storchaka@gmail.com> - 2014-11-26 18:24 +0200

csiph-web