john w. kennedy

nth_root: procedure; arg x, p; x = trunc(x); p = trunc(p)
y = 1; y1 = -1; p1 = p - 1
do while (abs(y - y1) > 1)
y1 = y
y = (p1 * y + x % y ** p1) % p
end
if y ** p < x then
return y1
else
return y

The above is not completely stable against overflow, which will not
happen in REXX, but may happen in PL/I.

(I have to thank you for asking this; if you hadn't, I wouldn't have
discovered Open Object REXX; I've missed having REXX since I had to
abandon OS/2.)

--
John W. Kennedy
"But now is a new thing which is very old--
that the rich make themselves richer and not poorer,
which is the true Gospel, for the poor's sake."
-- Charles Williams. "Judgement at Chelmsford"

robin

Don't recall it being posted,
but it may be that by Heinrich, published in 1996,
in Dr Dobbs' Journal, #4, April.

n. shamsundar

Please correct Line-3 to
do while (abs(y - y1) > 0)

To see why, run the original and corrected versions with x = 1027 and p
= 10.

The code is simply an integer version of Newton's method.

N. Shamsundar
University of Houston

gerard schildberger

|> Gerard Schildberger wrote:
|>> Below is a subroutine to find the integer squart root (I've
|>> forgetten who posted this via this newsgroup):
|>> __________________________________________________ ____
|>>
|>> /*find integer sqrt of a non-negative #.*/
|>>
|>> isqrt: procedure; parse arg x; x=trunc(x); r=0; q=1
|>>
|>> do while q<=x
|>> q=q*4
|>> end
|>>
|>> do while q>1
|>> q=q%4
|>> _=x-r-q
|>> r=r%2
|>>
|>> if _>=0 then do
|>> x=_
|>> r=r+q
|>> end
|>>
|>> end
|>>
|>> return r
|>> __________________________________________________ _____
|>>
|>>
|>> Does anyone know of (or where can I get it) an integer
|>> cube root routine? ______________________________Gerard S.
|> nth_root: procedure; arg x, p; x = trunc(x); p = trunc(p)
|> y = 1; y1 = -1; p1 = p - 1
|> do while (abs(y - y1) > 1)
|> y1 = y
|> y = (p1 * y + x % y ** p1) % p
|> end
|> if y ** p < x then
|> return y1
|> else
|> return y
|>
|> The above is not completely stable against overflow, which will not
|> happen in REXX, but may happen in PL/I.
|>
|> (I have to thank you for asking this; if you hadn't, I wouldn't have
|> discovered Open Object REXX; I've missed having REXX since I had to
|> abandon OS/2.)

| Please correct Line-3 to |
| do while (abs(y - y1) > 0)
|
| To see why, run the original and corrected versions with x = 1027 and p
| = 10.
|
| The code is simply an integer version of Newton's method.

What I'm looking for is a specific routine just for integer cube roots,
similar to the integer square root (above), which isn't (as far as I
know, an integer version of Newton's method). I know there is one, and
I was hoping that somebody whould know of it (or where it is on the www).
The general routine (above) for any root seems like overkill for what I
wanted to do. __________________________________________________ Gerard S.

gerard schildberger

| > Gerard Schildberger wrote:
| >> Below is a subroutine to find the integer squart root (I've
| >> forgetten who posted this via this newsgroup):
| >> __________________________________________________ ____
| >>
| >> /*find integer sqrt of a non-negative #.*/
| >>
| >> isqrt: procedure; parse arg x; x=trunc(x); r=0; q=1
| >>
| >> do while q<=x
| >> q=q*4
| >> end
| >>
| >> do while q>1
| >> q=q%4
| >> _=x-r-q
| >> r=r%2
| >>
| >> if _>=0 then do
| >> x=_
| >> r=r+q
| >> end
| >>
| >> end
| >>
| >> return r
| >> __________________________________________________ _____
| >>
| >>
| >> Does anyone know of (or where can I get it) an integer
| >> cube root routine? ______________________________Gerard S.
| >
| > nth_root: procedure; arg x, p; x = trunc(x); p = trunc(p)
| > y = 1; y1 = -1; p1 = p - 1
| > do while (abs(y - y1) > 1)
| > y1 = y
| > y = (p1 * y + x % y ** p1) % p
| > end
| > if y ** p < x then
| > return y1
| > else
| > return y
| >
| > The above is not completely stable against overflow, which will not
| > happen in REXX, but may happen in PL/I.
| >
| > (I have to thank you for asking this; if you hadn't, I wouldn't have
| > discovered Open Object REXX; I've missed having REXX since I had to
| > abandon OS/2.)
| >
|
| Please correct Line-3 to |
| do while (abs(y - y1) > 0)
|
| To see why, run the original and corrected versions with x = 1027 and p
| = 10.
|
| The code is simply an integer version of Newton's method.

I tried x=99, p=2 and it loops. ______________________________Gerard S.

glen herrmannsfeldt

Find a book an how to use an abacus that includes square
and cube roots. I once bought a cheap abacus in a chinatown
store that included a little book.

For the square root algorithm, similar to the integer square
root algorithm you mention, successive odd integers are subtracted
from the high digit, and the number of such is counted.
(The sum of the first n odd integers is n**2.)
(Though with a shift of the decimal point it will work for
fixed point non-integer numbers, too.)

Then there is a correction when you shift to the next lower digits
(one in the accumulating square root, and two in the ever decreasing
starting number.) The square root algorithm comes from the identity
(10A+B)**2=100A+20AB+B**2

The cube root algorithm takes three sets of abacus columns,
accumulating x and x**2 while reducing the original number by x**3.
(For more digits you put two abaci side by side.)

See: http://hem.passagen.se/ceem/china.htm

The cube root algorithm comes from:
(10A+B)**3=1000A**3+300*A**2*B+30*A*B**2+B**3


-- glen

gerard schildberger

| > What I'm looking for is a specific routine just for integer cube roots,
| > similar to the integer square root (above), which isn't (as far as I
| > know, an integer version of Newton's method). I know there is one, and
| > I was hoping that somebody whould know of it (or where it is on the www).
| > The general routine (above) for any root seems like overkill for what I
| > wanted to do. __________________________________________________ Gerard S.
|
|
| Find a book an how to use an abacus that includes square
| and cube roots. I once bought a cheap abacus in a chinatown
| store that included a little book.

Oh, great. I ask what time it is, and you tell me to buy a book
on how to build a watch. Not everyone has access to a good book
store. _____Gerard S.


| For the square root algorithm, similar to the integer square
| root algorithm you mention, successive odd integers are subtracted
| from the high digit, and the number of such is counted.
| (The sum of the first n odd integers is n**2.)
| (Though with a shift of the decimal point it will work for
| fixed point non-integer numbers, too.)

I don't see that in the code below, a subroutine to find the
integer squart root:
__________________________________________________ ____

/*find integer sqrt of a non-negative #.*/

isqrt: procedure; parse arg x; x=trunc(x); r=0; q=1

do while q<=x
q=q*4 end
do while q>1
q=q%4
_=x-r-q r=r%2
if _>=0 then do
x=_
r=r+q
end

end

return r
__________________________________________________ _____

(% is integer divide in Rexx). __________________Gerard S.


| Then there is a correction when you shift to the next lower digits
| (one in the accumulating square root, and two in the ever decreasing
| starting number.) The square root algorithm comes from the identity
| (10A+B)**2=100A+20AB+B**2
|
| The cube root algorithm takes three sets of abacus columns,
| accumulating x and x**2 while reducing the original number by x**3.
| (For more digits you put two abaci side by side.)
|
| See: http://hem.passagen.se/ceem/china.htm
|
| The cube root algorithm comes from:
| (10A+B)**3=1000A**3+300*A**2*B+30*A*B**2+B**3
| -- glen

john w. kennedy

With your correction, it can fall into a loop. I should have gone with
my first instinct; doing it with integers is a dumb idea.

How to get a cube root:

cube_root: procedure; arg x; x = trunc(x)
return trunc(x ** 1/3) Ta-da!


--
John W. Kennedy
"But now is a new thing which is very old--
that the rich make themselves richer and not poorer,
which is the true Gospel, for the poor's sake."
-- Charles Williams. "Judgement at Chelmsford"

gerard schildberger

|> John W. Kennedy wrote:
|>> Gerard Schildberger wrote:
|>>> Below is a subroutine to find the integer squart root (I've
|>>> forgetten who posted this via this newsgroup):
|>>> __________________________________________________ ____
|>>> /*find integer sqrt of a non-negative #.*/
|>>> isqrt: procedure; parse arg x; x=trunc(x); r=0; q=1
|>>> do while q<=x
|>>> q=q*4
|>>> end
|>>> do while q>1
|>>> q=q%4
|>>> _=x-r-q
|>>> r=r%2
|>>> if _>=0 then do
|>>> x=_
|>>> r=r+q
|>>> end
|>>> end
|>>> return r
|>>> __________________________________________________ _____
|>>> Does anyone know of (or where can I get it) an integer
|>>> cube root routine? ______________________________Gerard S. |
| >> nth_root: procedure; arg x, p; x = trunc(x); p = trunc(p)
| >> y = 1; y1 = -1; p1 = p - 1
| >> do while (abs(y - y1) > 1)
| >> y1 = y
| >> y = (p1 * y + x % y ** p1) % p
| >> end
| >> if y ** p < x then
| >> return y1
| >> else
| >> return y
| >>
| >> The above is not completely stable against overflow, which will not
| >> happen in REXX, but may happen in PL/I.
| >>
| >> (I have to thank you for asking this; if you hadn't, I wouldn't have
| >> discovered Open Object REXX; I've missed having REXX since I had to
| >> abandon OS/2.) |
| > Please correct Line-3 to
| > do while (abs(y - y1) > 0)
| > To see why, run the original and corrected versions with x = 1027 and p
| > = 10.

| With your correction, it can fall into a loop. I should have gone with
| my first instinct; doing it with integers is a dumb idea.
|
| How to get a cube root:
|
| cube_root: procedure; arg x; x = trunc(x)
| return trunc(x ** 1/3)
|
| Ta-da!

Nice joke (if you intended it to be). It simply doesn't work in Rexx,
which is where I originally posted from and wanted a REXX solution.

By the way, Robin Vowels gave me a working version (in PL/1) which is
copyrighted, so I won't post it here.
__________________________________________________ _____________Gerard S.

n. shamsundar

If you are willing to use RxMath, the solution is:

/* cube root using float/truncate */

call RxFuncAdd"MathLoadFuncs","rxmath","MathLoadFuncs"
call MathLoadFuncs

yo=0
do x=1 to 100000 y=trunc(rxCalcPower(x+0.1,1/3.0))
if y<>yo then say 'Cube root of ' || x || ' is ' || y
yo=y
end
return

If you do not, here is the fixed-up Newton-Raphson integer code (no
copyright claimed) :

/* integer cube root by Newton-Raphson */
yo=0
do x=1 to 100000 y=cubrt(x)
if y<>yo then say 'Cube root of ' || x || ' is ' || y
yo=y
end
return

cubrt: procedure; arg x; x = trunc(x)
y = 1; v=1; v1=x;
if x < 8 then return 1
do while (abs(y*(v-v1)) > 2*v)
v1 = v
y = (2 * y + x % v) % 3
v = y * y
end
if v*y <= x then
return y
else
return y-1

-- N. Shamsundar
-- University of Houston