Square root

Then of course there's the "babylonian" method:

final int ITERATIONS = 10000; //Any large number - The larger the more acccurate

double num = 5; //Number to take square root of

double r = 1; //Any number actually works for r

for (int i = 0; i < ITERATIONS; i++)

r = ((num / r) + r) / 2.0;

System.out.println("Square root of "+num+" = "+r);

-Sam

colesburya at 2007-7-16 2:34:18 >

If you have a quantum computer this should work for you:

double sqrt( double value) {

Random r = new Random();

double tolerance = 0.0625;

double sqrt = 0;

while (true) {

long randomLong = r.nextLong();

sqrt = Double.longBitsToDouble( randomLong );

if (Math.abs(sqrt*sqrt - value) < tolerance) {

return sqrt;

}

}

}

sethawa at 2007-7-16 2:34:18 >

> I guess the difference is that with Newton's method

> you can get an iterative approximation to any

> function, while the Babylonian method was a good

> algorithm that would be hard to extend without general

> principles.

>

> I'm not a mathematician (just an ex-engineer), nor a

> mathematics historian, so I can't say.

>

> It's still Newton's method to me. - MOD

The Babylonian method is identical to a Newton/Raphson iteration for x^2-p = 0.

Have a look -- x' = x-f(x)/f'(x) where x' is the new approximation and f' denotes the derivative of f w.r.t. x.

This is the NR iteration. Now fill in the blanks -- x' = x-(x^2-p)/(2*x) --> x' = 1/2*(x+p/x), whick is exactly

what the Babylonian method does. Both methods take the same iteration steps. Those old guys

where smart back in those days.

kind regards,

Jos

JosHorsmeiera at 2007-7-16 2:34:18 >

>

> The Babylonian method is identical to a

> Newton/Raphson iteration for x^2-p = 0.

> Have a look -- x' = x-f(x)/f'(x) where x' is the new

> approximation and f' denotes the derivative of f

> w.r.t. x.

> This is the NR iteration. Now fill in the blanks -- x'

> = x-(x^2-p)/(2*x) --> x' = 1/2*(x+p/x), whick is

> exactly

> what the Babylonian method does. Both methods take the

> same iteration steps. Those old guys

> where smart back in those days.

>

> kind regards,

>

> Jos

Absolutely right, Jos.

I'm just wondering if they managed the same trick for cube roots or arbitrary roots. You can do it with Newton's method.Did those Babylonians figure out the general rule?

Pretty darned smart. - MOD

duffymoa at 2007-7-16 2:34:18 >

> That must have been a hell of a calculation for the Babylonians, using a number system with a mixture of

> base 10 and base 60 and using integers only. Just the division part alone would be horribly complicated.

Don't be surprised; the Babylonians used a decimal notation (not positional) for their scientific

calculations; when they needed to perform calculations in base 60, they resorted to reciprocals of

factors of 60. Another form of decimal notation was later known by the Vedics (1000 BC), including

a notion of zero and brought back a bit to the west by the Arabs.

Note that students of Fibonacci were brought to death by the catholic church somewhere around the

12th century because the used decimal notation instead of roman numerals. It's quite a history ...

kind regards,

Jos

JosHorsmeiera at 2007-7-16 2:34:18 >

> I'm just wondering if they managed the same trick for

> cube roots or arbitrary roots.

When I was in high school we learned to calculate square roots with paper and pencil. The method we used looked like division, but you used the decimal places of the number you were square-rooting two at a time instead of the one at a time that ordinary division does, and you made up the number you were "dividing by" as you went along so that the result of the "division" was the same as it. There was a multiplication by 2 in there somewhere as well.

(This was in the days before calculators with square-root buttons.)

Anyway, I thought this was really cool and I figured out how to do cube roots using a similar but more complicated method. That was the day I realized I was a mathematician.

DrClapa at 2007-7-16 2:34:18 >

> was this it?

> http://www.nist.gov/dads/HTML/squareRoot.html

Yes, that's exactly it. Cute isn't it? Anyone interested in finding the decimal expansion of 1/p where p

happens to be prime? It can even be done without pencil and paper after a bit of practicing ...

kind regards,

Jos

JosHorsmeiera at 2007-7-19 18:29:58 >

http://www.nist.gov/dads/HTML/squareRoot.html

A scalar binary implementation.class Sqrt {

public static void main(String[] argv) {

long l;

try {

l = Long.parseLong(argv[0]);

} catch (NumberFormatException ex) {

ex.printStackTrace();

return;

}

System.out.println(Sqrt.sqrt(l));

}

public static int sqrt(long n) {

int bits = 0;

while (n> 1<<bits ) bits += 2;

int scale = 1 ><< (bits>>1);

int s = 0, t;

while (scale>0) {

scale >>>= 1;

t = s+scale;

if (t*(long)t<=n) s = t;

}

return s;

}

}

tschodta at 2007-7-19 18:29:58 >

Math.sqrt() gives you an approximation of the square root that happens to be equals to the exact value in some cases (for instance, in the cases 4 and 5.0625 you've mentioned).

This is an excerpt of j2se\src\share\native\java\lang\fdlibm\src\e_sqrt.c, that is the source code of the native implementation of Math.sqrt. It explains the algorithm that it uses.

/* __ieee754_sqrt(x)

* Return correctly rounded sqrt.

*

*| Use the hardware sqrt if you have one |

*

* Method:

*Bit by bit method using integer arithmetic. (Slow, but portable)

*1. Normalization

*Scale x to y in [1,4) with even powers of 2:

*find an integer k such that 1 <= (y=x*2^(2k)) < 4, then

*sqrt(x) = 2^k * sqrt(y)

*2. Bit by bit computation

*Let q = sqrt(y) truncated to i bit after binary point (q = 1),

*i 0

* i+1 2

*s = 2*q , andy = 2* ( y - q ).(1)

*iiii

*

*To compute qfrom q , one checks whether

*i+1i

*

*-(i+1) 2

*(q + 2) <= y.(2)

*i

*-(i+1)

*If (2) is false, then q= q ; otherwise q= q + 2.

*i+1i i+1i

*

*With some algebric manipulation, it is not difficult to see

*that (2) is equivalent to

*-(i+1)

*s + 2<= y(3)

* ii

*

*The advantage of (3) is that s and y can be computed by

*ii

*the following recurrence formula:

*if (3) is false

*

*s= s ,y= y;(4)

*i+1i i+1i

*

*otherwise,

* -i -(i+1)

*s = s + 2 , y= y - s - 2 (5)

*i+1i i+1ii

*

*One may easily use induction to prove (4) and (5).

*Note. Since the left hand side of (3) contain only i+2 bits,

*it does not necessary to do a full (53-bit) comparison

*in (3).

*3. Final rounding

*After generating the 53 bits result, we compute one more bit.

*Together with the remainder, we can decide whether the

*result is exact, bigger than 1/2ulp, or less than 1/2ulp

*(it will never equal to 1/2ulp).

*The rounding mode can be detected by checking whether

*huge + tiny is equal to huge, and whether huge - tiny is

*equal to huge for some floating point number "huge" and "tiny".

*

* Special cases:

*sqrt(+-0) = +-0 ... exact

*sqrt(inf) = inf

*sqrt(-ve) = NaN... with invalid signal

*sqrt(NaN) = NaN... with invalid signal for signaling NaN

*

* Other methods : see the appended file at the end of the program below.

*

*/

.....

.....

/*

Other methods (use floating-point arithmetic)

-

(This is a copy of a drafted paper by Prof W. Kahan

and K.C. Ng, written in May, 1986)

Two algorithms are given here to implement sqrt(x)

(IEEE double precision arithmetic) in software.

Both supply sqrt(x) correctly rounded. The first algorithm (in

Section A) uses newton iterations and involves four divisions.

The second one uses reciproot iterations to avoid division, but

requires more multiplications. Both algorithms need the ability

to chop results of arithmetic operations instead of round them,

and the INEXACT flag to indicate when an arithmetic operation

is executed exactly with no roundoff error, all part of the

standard (IEEE 754-1985). The ability to perform shift, add,

subtract and logical AND operations upon 32-bit words is needed

too, though not part of the standard.

A. sqrt(x) by Newton Iteration

(1)Initial approximation

Let x0 and x1 be the leading and the trailing 32-bit words of

a floating point number x (in IEEE double format) respectively

11152 ...widths

x: |s| e|f|

msblsb msblsb ...order

x0: |s|e|f1| x1: | f2|

By performing shifts and subtracts on x0 and x1 (both regarded

as integers), we obtain an 8-bit approximation of sqrt(x) as

follows.

k := (x0>>1) + 0x1ff80000;

y0 := k - T1[31&(k>>15)].... y ~ sqrt(x) to 8 bits

Here k is a 32-bit integer and T1[] is an integer array containing

correction terms. Now magically the floating value of y (y's

leading 32-bit word is y0, the value of its trailing word is 0)

approximates sqrt(x) to almost 8-bit.

Value of T1:

static int T1[32]= {

0,1024,3062,5746,9193,13348,18162,23592,

29598,36145,43202,50740,58733,67158,75992,85215,

83599,71378,60428,50647,41945,34246,27478,21581,

16499,12183,8588,5674,3403,1742,661,130,};

(2)Iterative refinement

Apply Heron's rule three times to y, we have y approximates

sqrt(x) to within 1 ulp (Unit in the Last Place):

y := (y+x/y)/2... almost 17 sig. bits

y := (y+x/y)/2... almost 35 sig. bits

y := y-(y-x/y)/2... within 1 ulp

Remark 1.

Another way to improve y to within 1 ulp is:

y := (y+x/y)... almost 17 sig. bits to 2*sqrt(x)

y := y - 0x00100006... almost 18 sig. bits to sqrt(x)

2

(x-y )*y

y := y + 2* -...within 1 ulp

2

3y + x

This formula has one division fewer than the one above; however,

it requires more multiplications and additions. Also x must be

scaled in advance to avoid spurious overflow in evaluating the

expression 3y*y+x. Hence it is not recommended uless division

is slow. If division is very slow, then one should use the

reciproot algorithm given in section B.

(3) Final adjustment

By twiddling y's last bit it is possible to force y to be

correctly rounded according to the prevailing rounding mode

as follows. Let r and i be copies of the rounding mode and

inexact flag before entering the square root program. Also we

use the expression y+-ulp for the next representable floating

numbers (up and down) of y. Note that y+-ulp = either fixed

point y+-1, or multiply y by nextafter(1,+-inf) in chopped

mode.

I := FALSE;... reset INEXACT flag I

R := RZ;... set rounding mode to round-toward-zero

z := x/y;... chopped quotient, possibly inexact

If(not I) then {... if the quotient is exact

if(z=y) {

I := i; ... restore inexact flag

R := r; ... restore rounded mode

return sqrt(x):=y.

} else {

z := z - ulp;... special rounding

}

}

i := TRUE;... sqrt(x) is inexact

If (r=RN) then z=z+ulp... rounded-to-nearest

If (r=RP) then {... round-toward-+inf

y = y+ulp; z=z+ulp;

}

y := y+z;... chopped sum

y0:=y0-0x00100000;... y := y/2 is correctly rounded.

I := i; ... restore inexact flag

R := r; ... restore rounded mode

return sqrt(x):=y.

(4)Special cases

Square root of +inf, +-0, or NaN is itself;

Square root of a negative number is NaN with invalid signal.

B. sqrt(x) by Reciproot Iteration

(1)Initial approximation

Let x0 and x1 be the leading and the trailing 32-bit words of

a floating point number x (in IEEE double format) respectively

(see section A). By performing shifs and subtracts on x0 and y0,

we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.

k := 0x5fe80000 - (x0>>1);

y0:= k - T2[63&(k>>14)].... y ~ 1/sqrt(x) to 7.8 bits

Here k is a 32-bit integer and T2[] is an integer array

containing correction terms. Now magically the floating

value of y (y's leading 32-bit word is y0, the value of

its trailing word y1 is set to zero) approximates 1/sqrt(x)

to almost 7.8-bit.

Value of T2:

static int T2[64]= {

0x1500,0x2ef8,0x4d67,0x6b02,0x87be,0xa395,0xbe7a,0xd866,

0xf14a,0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,

0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,

0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,

0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,

0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,

0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,

0x1527f,0x1334a,0x11051,0xe951,0xbe01,0x8e0d,0x5924,0x1edd,};

(2)Iterative refinement

Apply Reciproot iteration three times to y and multiply the

result by x to get an approximation z that matches sqrt(x)

to about 1 ulp. To be exact, we will have

-1ulp < sqrt(x)-z<1.0625ulp.

... set rounding mode to Round-to-nearest

y := y*(1.5-0.5*x*y*y)... almost 15 sig. bits to 1/sqrt(x)

y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)

... special arrangement for better accuracy

z := x*y... 29 bits to sqrt(x), with z*y<1

z := z + 0.5*z*(1-z*y)... about 1 ulp to sqrt(x)

Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that

(a) the term z*y in the final iteration is always less than 1;

(b) the error in the final result is biased upward so that

-1 ulp < sqrt(x) - z < 1.0625 ulp

instead of |sqrt(x)-z|<1.03125ulp.

(3)Final adjustment

By twiddling y's last bit it is possible to force y to be

correctly rounded according to the prevailing rounding mode

as follows. Let r and i be copies of the rounding mode and

inexact flag before entering the square root program. Also we

use the expression y+-ulp for the next representable floating

numbers (up and down) of y. Note that y+-ulp = either fixed

point y+-1, or multiply y by nextafter(1,+-inf) in chopped

mode.

R := RZ;... set rounding mode to round-toward-zero

switch(r) {

case RN:... round-to-nearest

if(x<= z*(z-ulp)...chopped) z = z - ulp; else

if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;

break;

case RZ:case RM:... round-to-zero or round-to--inf

R:=RP;... reset rounding mod to round-to-+inf

if(x<z*z ... rounded up) z = z - ulp; else

if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;

break;

case RP:... round-to-+inf

if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else

if(x>z*z ...chopped) z = z+ulp;

break;

}

Remark 3. The above comparisons can be done in fixed point. For

example, to compare x and w=z*z chopped, it suffices to compare

x1 and w1 (the trailing parts of x and w), regarding them as

two's complement integers.

...Is z an exact square root?

To determine whether z is an exact square root of x, let z1 be the

trailing part of z, and also let x0 and x1 be the leading and

trailing parts of x.

If ((z1&0x03ffffff)!=0)... not exact if trailing 26 bits of z!=0

I := 1;... Raise Inexact flag: z is not exact

else {

j := 1 - [(x0>>20)&1]... j = logb(x) mod 2

k := z1 >> 26;... get z's 25-th and 26-th

fraction bits

I := i or (k&j) or ((k&(j+j+1))!=(x1&3));

}

R:= r... restore rounded mode

return sqrt(x):=z.

If multiplication is cheaper then the foregoing red tape, the

Inexact flag can be evaluated by

I := i;

I := (z*z!=x) or I.

Note that z*z can overwrite I; this value must be sensed if it is

True.

Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be

zero.

--

z1: |f2|

--

bit 31bit 0

Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd

or even of logb(x) have the following relations:

-

bit 27,26 of z1bit 1,0 of x1logb(x)

-

0000odd and even

0101even

1010odd

1000even

1101even

-

(4)Special cases (see (4) of Section A).

*/

edsonwa at 2007-7-19 18:29:58 >

> You know, sqrt can be computed in constant time.

> Rather, time that is linear in the number of bits in

> your number (32 or 64). This can be done with a simple

> guess-and-check algorithm, if you understand how the

> IEEE 754 float layout works. Hard to beat constant

> time. :-)

Could you provide a link to this algorithm?

I don't mean to sound flip, but I don't see how mere knowledge of the IEEE 754 float layout will help calculate a square root or anything else.

Does it work for cube roots? How about transcendental functions? Complimentary error function? Bessel functions?

Wouldn't time depend on the number and quality of your guesses? Isn't Newton-Raphson iterative method just a way to guide that guess process? What if you guess far from the root and things diverge? (It's been known to happen.)

I might learn something here, but I'm very skeptical. Please post the link. Thanks - %

duffymoa at 2007-7-19 18:29:58 >

But these are integer square roots. That's not the problem I'm thinking of. I thought this discussion was about roots of general doubles, not just integers.

There's no doubt that a table lookup will be pretty darned fast (constant time), but that won't help if you want the square root of an arbitrary double. That's what Newton-Raphson gives you.

Is the IEEE floating point standard algorithm only for integers? If so, it's not applicable to this discussion, IMO.

%

duffymoa at 2007-7-19 18:29:58 >

> I start the first of September but I quit the old job

> the first of July and they pay 'till September, so I

> have

> a luxury two months sabbatical leave now ;-) I'm

> sitting at my picknick table in my back garden;

> it's sunny and I'm playing with my new Java primal

> simplex solver ... I know, I know, it's a rotten job

> but someone's got to do it! ;-)

>

> kind regards,

>

> Jos

Very nice. The idea of two months of is incredible. Europe can be so much more enlightened sometimes.

%

duffymoa at 2007-7-19 18:30:05 >

Very nice !

I spent a bunch of time doing some LP/MIQP solving a while back using ILog's tools - really fun and interesting stuff (yes, I know I need to get out more).

I'd really love to see this - perhaps you could find it a home on codehaus/sourceforge/other (if you need some webspace I'd be delighted to oblige) ? And if you make it available, please do let us know via the forums.

Dave.

dcmintera at 2007-7-19 18:30:05 >

Sounds like great work, as usual. I'd like very much to see how you approached the design.

Your preprocessor and builder approach to solving make me think of the usual finite element approach: local element matricies assembled into a global matrix for solution. Assemble, apply constraints, solve. You can even save yourself the trouble of assembling a global matrix if you use a wavefront solver. That generates element matricies and brings them into the wavefront as needed.

But your sparse mapper might make all that superfluous. Very elegant indeed. Probably easy on memory, too.

Brilliant. Thanks, JosAH.

%

duffymoa at 2007-7-19 18:30:05 >

> Maybe this is a silly question but ... what exactly does this thing solve :) Seems

> interesting but there are far too many acronyms in there for me to understand at the moment :)

Yes, LP (Linear Programming) carries its own bag of jargon. It optimizes (minimizes or maximizes)

a linear problem, i.e. a bunch of variables given a set of linear equations or inequalities.

Let A be a n*m matrix (m columns and n rows), let b be an n vector (also called the 'rhs' or Right

Hand Side vector) and let c be an m vector, the cost vector. The standard problem can be formulated as -- minimize sum(j= 1,m: c_j*x_j) w.r.t.

A*x= b, lo_i <= x_i <= hi_i (i in 1,n)

This is an NP complete problem but the simplex method (the thing I'm implementing now) uses a

very favourable big-Oh number of iterations on average to solve the problem.

Here is where [url=http://www-fp.mcs.anl.gov/otc/Guide/CaseStudies/simplex/]the simplex method[/url] is explained, and [url=http://www-fp.mcs.anl.gov/otc/Guide/TestProblems/LPtest/]here[/url] are some very nasty

test problems; as for now my implementation can just solve some of these problems, due to numerical

inaccuracies, bad problem formulation (the row rank of matrix A has to be equal to n), etc. etc. Note that

matrix A happens to be huge most of the time ( > 1000 rows and columns ), but very sparse.

kind regards,

Jos

JosAHa at 2007-7-19 18:30:05 >

From a much more simple perspective, Linear Programming is a topic I discovered on my AS Level maths course. We had a discrete maths module along with pure and statistics and it covered amongst other things The Simplex Method and Linear Programming. This is a problem I believe is best understood mathematically before attempting to write a program.

I found it very difficult to swallow but slowly came to terms with the fact that it is as usual a very very simple way to solve what is a difficult problem. In fact when I grasped it I found it really satisfying. Often when studying maths I find the material abstract but this problem seemed really rooted in practicallity.

You should have a look at it. The text book I was using at the time was Cambridge Advanced Level Mathematics, Discrete Mathematics 1. However my Maths teacher didn't like this book. He didn't recommend anything else, but for this level maybe a Shoam's Outline book would be good.

stanton_iana at 2007-7-19 18:30:05 >