cubic root subroutine

Do you have a question? Post it now! No Registration Necessary.  Now with pictures!

•  Subject
• Author
• Posted on

Hi.

I needed a cubic root subroutine but didnīt find one at a quick search.
So I wrote a simple one.
It seems pretty fast and accurate enough.
Any ideas  how to improve it,
or point me to better solutions.
/dg

#---------------------------------------------------------------------------
# Windows Perl

for (\$number=2; \$number < 126;\$number++) {
\$resp = cubicroot(\$number);
print "3 root of \$number = \$resp\n";
}

exit;
#------------------------------------------------------------------
sub cubicroot{

my (\$number) = @_;
\$number = abs(\$number);

my \$reply = int(sqrt(sqrt(\$number)));   #Need to start somewhere. Any
better idea ?

my (\$diff,\$decimal,\$fraction);

# first find the nearest higher whole number
while (\$reply**3  < \$number ) {
}

# We start with 1.0 to catch the even cubic ( 8 27 64 125 ...)
# and then count down to the nerast lower by 0.1
\$decimal  = 1.0;
\$fraction = 0.1;
\$diff = 1;

WHILE:
while (\$diff > 0)  {
while ( (\$reply + \$decimal)**3 > \$number ) {
\$decimal -= \$fraction;
#we need a precision check here
if (\$decimal - \$fraction == \$decimal) {last WHILE}
}
# Now we divide fraction by ten and go up again
\$fraction  = \$fraction / 10;
while ( (\$reply + \$decimal)**3 < \$number ) {
\$decimal +=\$fraction;
#we need a precision check here
if (\$decimal + \$fraction == \$decimal) {last WHILE}
}
# Divide fraction by ten for next loop
\$fraction  = \$fraction / 10;
# Calculate diff before next loop
\$diff = abs( (\$reply + \$decimal)**3 -  \$number);
}
}

Re: cubic root subroutine

Dan van Ginhoven wrote:
> Hi.
>
> I needed a cubic root subroutine but didnīt find one at a quick
> search. So I wrote a simple one. It seems pretty fast and accurate
> enough. Any ideas  how to improve it, or point me to better
> solutions.

<code snipped>

If you remember back to your first year algebra class, you will recall
that the third root of a number is the same as that number raised to the
1/3 power.

use strict;
use warnings;

for my \$number (2 .. 125) {
my \$resp = cubicroot(\$number);
print "3 root of \$number = \$resp\n";
}

sub cubicroot {
return \$_[0] ** (1/3);
}

Re: cubic root subroutine

gw@comcast.com:

> If you remember back to your first year algebra class, you will recall
> that the third root of a number is the same as that number raised to the
> 1/3 power.

Oh, God! I was so mesmerized by the mish-mash the OP posted, that I missed
the very very obvious.

I have a feeling I will not be able to sleep tonight :)

Thanks

Sinan.

Re: cubic root subroutine

> I needed a cubic root subroutine but didnīt find one at a quick
> search. So I wrote a simple one.

You mean "a simplistic and naive one"

> It seems pretty fast and accurate enough.

How did you test it?

> Any ideas  how to improve it,

Using strict and warnings will reveal some problems but you basically have
an algorithmic problem.

D:\Home\asu1\UseNet\clpmisc> ttt
Cubic root of -1 = 1

Sinan.

Re: cubic root subroutine

> or point me to better solutions.

I haven't tested it myself, but I would give

http://search.cpan.org/~jonathan/Math-Calculus-NewtonRaphson-
0.1/NewtonRaphson.pm

a shot first.

Sinan

Re: cubic root subroutine

Dan van Ginhoven wrote:
> I needed a cubic root subroutine but didnīt find one at a quick
> search. So I wrote a simple one.
> It seems pretty fast and accurate enough.
> Any ideas  how to improve it,
> or point me to better solutions.

Is there anything wrong with just using the exponentiation operator?

\$cubicroot = \$foo ** (1/3);

jue

Re: cubic root subroutine

> Dan van Ginhoven wrote:
> > I needed a cubic root subroutine but didnīt find one at a quick
> > search. So I wrote a simple one.
> > It seems pretty fast and accurate enough.
> > Any ideas  how to improve it,
> > or point me to better solutions.
>
> Is there anything wrong with just using the exponentiation operator?
>
>     \$cubicroot = \$foo ** (1/3);

You may want to handle the sign separately, since letting \$foo be
negative results in NaN ("Not a number"), at least here, on
Solaris 2.8 x86.

Also, it may be worth making sure that the cube of an integer always
comes back as an integer, and not with a little rounding-error, if
that is important for the application. The IEEE floating-point
specification doesn't guarantee that, or does it?

Re: cubic root subroutine

Arndt Jonasson wrote:

> You may want to handle the sign separately, since letting \$foo be
> negative results in NaN ("Not a number"), at least here, on
> Solaris 2.8 x86.

Same on Win32.

> Also, it may be worth making sure that the cube of an integer always
> comes back as an integer, and not with a little rounding-error,

I assume you mean that the cube root of a cube comes back an integer.
That doesn't appear to be a problem on my machine.

use strict;
use warnings;

# Demonstrate that this works for numbers that are cubes,
# including negative cubes, without rounding error.
for my \$number (-20 .. 20) {
my \$x = \$number ** 3;
my \$resp = cubicroot(\$x);
print "3 root of \$x = \$resp\n";
}

sub cubicroot {

my \$x = shift;

if (\$x < 0) {
return - abs(\$x) ** (1/3);
}
else {
return \$x ** (1/3);
}
}

There are actually 3 cube roots of any number, but finding the other two
is a more complex problem.

Re: cubic root subroutine

> Arndt Jonasson wrote:
>
> > You may want to handle the sign separately, since letting \$foo be
> > negative results in NaN ("Not a number"), at least here, on
> > Solaris 2.8 x86.
>
> Same on Win32.
>
> > Also, it may be worth making sure that the cube of an integer always
> > comes back as an integer, and not with a little rounding-error,
>
> I assume you mean that the cube root of a cube comes back an
> integer. That doesn't appear to be a problem on my machine.

Yes, I meant "comes back from the cube root subroutine". But it is a
problem on my machine:

perl -e 'printf "%d\n", int(64**(1/3))'
3

Re: cubic root subroutine

>
>> Arndt Jonasson wrote:
>>
>> > You may want to handle the sign separately, since letting \$foo be
>> > negative results in NaN ("Not a number"), at least here, on
>> > Solaris 2.8 x86.
>>
>> Same on Win32.
>>
>> > Also, it may be worth making sure that the cube of an integer
>> > always comes back as an integer, and not with a little
>> > rounding-error,
>>
>> I assume you mean that the cube root of a cube comes back an
>> integer. That doesn't appear to be a problem on my machine.
>
> Yes, I meant "comes back from the cube root subroutine". But it is a
> problem on my machine:
>
>     perl -e 'printf "%d\n", int(64**(1/3))'
>     3

int EXPR
int     Returns the integer portion of EXPR. If EXPR is omitted,
uses \$_. You should not use this function for rounding: one
because it truncates towards 0, and two because machine
representations of floating point numbers can sometimes
produce counterintuitive results. For example,
"int(-6.725/0.025)" produces -268 rather
than the correct -269; that's because it's really more like
"printf", or the "POSIX::floor" and "POSIX::ceil" functions
will serve you better than will int().

These calculations necessary to calculate the power of a number are
performed in floating point by Perl. If you are going to use the number
in other calculations, you should leave it as it is instead of jumping
through a bunch of hoops.

On the other hand, if you just wanted to report results in integers:

D:\> perl -e "printf q, ((64)**(1/3))"
4

Sinan.

Re: cubic root subroutine

>
> >
> >> Arndt Jonasson wrote:
> >>
> >> > You may want to handle the sign separately, since letting \$foo be
> >> > negative results in NaN ("Not a number"), at least here, on
> >> > Solaris 2.8 x86.
> >>
> >> Same on Win32.
> >>
> >> > Also, it may be worth making sure that the cube of an integer
> >> > always comes back as an integer, and not with a little
> >> > rounding-error,
> >>
> >> I assume you mean that the cube root of a cube comes back an
> >> integer. That doesn't appear to be a problem on my machine.
> >
> > Yes, I meant "comes back from the cube root subroutine". But it is a
> > problem on my machine:
> >
> >     perl -e 'printf "%d\n", int(64**(1/3))'
> >     3
>
>     int EXPR
>     int     Returns the integer portion of EXPR. If EXPR is omitted,
>             uses \$_. You should not use this function for rounding: one
>             because it truncates towards 0, and two because machine
>             representations of floating point numbers can sometimes
>             produce counterintuitive results. For example,
>             "int(-6.725/0.025)" produces -268 rather
>             than the correct -269; that's because it's really more like
>             -268.99999999999994315658 instead. Usually, the "sprintf",
>             "printf", or the "POSIX::floor" and "POSIX::ceil" functions
>             will serve you better than will int().
>
> These calculations necessary to calculate the power of a number are
> performed in floating point by Perl. If you are going to use the number
> in other calculations, you should leave it as it is instead of jumping
> through a bunch of hoops.
>
> On the other hand, if you just wanted to report results in integers:
>
> D:\> perl -e "printf q, ((64)**(1/3))"
> 4

Call me dense, but I'm not sure if you're taking issue with my
statement, or just providing additional information.

Re: cubic root subroutine

>

>> >> Arndt Jonasson wrote:
>> > Yes, I meant "comes back from the cube root subroutine". But it is
>> > a problem on my machine:
>> >
>> >     perl -e 'printf "%d\n", int(64**(1/3))'
>> >     3
>>
>>     int EXPR
>>     int     Returns the integer portion of EXPR. If EXPR is omitted,
>>             uses \$_. You should not use this function for rounding:

....

>> These calculations necessary to calculate the power of a number are
>> performed in floating point by Perl. If you are going to use the
>>  number in other calculations, you should leave it as it is instead
>> of jumping through a bunch of hoops.
>>
>> On the other hand, if you just wanted to report results in integers:
>>
>> D:\> perl -e "printf q, ((64)**(1/3))"
>> 4
>
> Call me dense, but I'm not sure if you're taking issue with my
> statement, or just providing additional information.

I am taking issue with your statement "it is a problem on my machine".
You seem to think there is a problem when in actuality, the behavior you
are observing is only due to the fact that you used the '%d' format for
printf instead of the more appropriate '%.0f' format.

There is nothing you can do about the fact that:

perl -e "printf q, ((64)**(1/3))"

will print something like

3.9999999999999996

i.e. it is impossible to:

<quoted from another message of yours>
.... it may be worth making sure that the cube of an integer always
comes back as an integer, and not with a little rounding-error,
</quoted from another message of yours>

For values used in interim calculations, you should not bother "fixing"
anything. For values that are printed, you should use the appropriate
printf format.

Sinan.

Re: cubic root subroutine

>
> >
>
> >> >> Arndt Jonasson wrote:
> >> > Yes, I meant "comes back from the cube root subroutine". But it is
> >> > a problem on my machine:
> >> >
> >> >     perl -e 'printf "%d\n", int(64**(1/3))'
> >> >     3
> >>
> >>     int EXPR
> >>     int     Returns the integer portion of EXPR. If EXPR is omitted,
> >>             uses \$_. You should not use this function for rounding:
>
> ...
>
> >> These calculations necessary to calculate the power of a number are
> >> performed in floating point by Perl. If you are going to use the
> >>  number in other calculations, you should leave it as it is instead
> >> of jumping through a bunch of hoops.
> >>
> >> On the other hand, if you just wanted to report results in integers:
> >>
> >> D:\> perl -e "printf q, ((64)**(1/3))"
> >> 4
> >
> > Call me dense, but I'm not sure if you're taking issue with my
> > statement, or just providing additional information.
>
> I am taking issue with your statement "it is a problem on my machine".
> You seem to think there is a problem when in actuality, the behavior you
> are observing is only due to the fact that you used the '%d' format for
> printf instead of the more appropriate '%.0f' format.
>
> There is nothing you can do about the fact that:
>
> perl -e "printf q, ((64)**(1/3))"
>
> will print something like
>
> 3.9999999999999996
>
> i.e. it is impossible to:
>
> <quoted from another message of yours>
> ... it may be worth making sure that the cube of an integer always
> comes back as an integer, and not with a little rounding-error,
> </quoted from another message of yours>
>
> For values used in interim calculations, you should not bother "fixing"
> anything. For values that are printed, you should use the appropriate
> printf format.

I think you are misunderstanding what I perceive as the problem.  I
stand by my original warning (and did not use the word "fixing", so
your quotation marks are not warranted). "It may be worth" was to be
read in the sense that you are fully entitled to consider it not worth
it, if the results are good enough for you - not that the results are
wrong.  My print statement only showed that the suggested cube root
calculation is not what you want if you want the cube root of the cube
of an integer to be the original integer. If printing out was all I
wanted, there is no problem. But you cut out the end of my warning,
which was something like "depending on the application". If I remember
correctly (which I maybe don't), the original poster didn't say what
he wanted to use the cube root for.

(Seeing it from another angle: if I had just come into the group with
appropriate. Thank you.)

A "fix" would be an actual integer cube root routine (one that returns
\$x when given \$x**3 as input, \$x being an integer). We haven't seen
one yet, but it's left as an exercise for the interested reader.

Re: cubic root subroutine

Arndt Jonasson wrote:

> A "fix" would be an actual integer cube root routine (one that returns
> \$x when given \$x**3 as input, \$x being an integer). We haven't seen
> one yet, but it's left as an exercise for the interested reader.

Did I miss something? I thought I posted one.

Re: cubic root subroutine

>

....

>> For values used in interim calculations, you should not bother
>> "fixing" anything. For values that are printed, you should use the
>> appropriate printf format.

....

> (and did not use the word "fixing", so your quotation marks are not
> warranted).

In this instance, I did not use the quotation marks to indicate I was
quoting you, but to qualify the word.

Sinan

Re: cubic root subroutine

>
>On the other hand, if you just wanted to report results in integers:
>
>D:\> perl -e "printf q, ((64)**(1/3))"
>4
>

I'm still at a loss as to why Math::Complex's cbrt() isn't
considered sufficient.

c:\users\jgamble>perl -MMath::Complex -le "print cbrt(64);"
4

--
-john

February 28 1997: Last day libraries could order catalogue cards
from the Library of Congress.

Re: cubic root subroutine

jgamble@ripco.com (John M. Gamble) writes:
> >
> >On the other hand, if you just wanted to report results in integers:
> >
> >D:\> perl -e "printf q, ((64)**(1/3))"
> >4
> >
>
> I'm still at a loss as to why Math::Complex's cbrt() isn't
> considered sufficient.
>
> c:\users\jgamble>perl -MMath::Complex -le "print cbrt(64);"
> 4

Somehow I missed your suggestion when I continued to critizise the
integerness (whether that is guaranteed I couldn't see from a quick
glance in the manual):

arndt ~/perl 4359> perl -MMath::Complex -le "print (cbrt(64)==4);"
1
arndt ~/perl 4360> perl -le "print (((64)**(1/3))==4)"

arndt ~/perl 4361>

Re: cubic root subroutine

: I needed a cubic root subroutine but didnīt find one at a quick search.
: So I wrote a simple one.
: It seems pretty fast and accurate enough.
: Any ideas  how to improve it,
: or point me to better solutions.

A good stab at the problem, but as numeric solutions go, the algorithm
sucks.

The classic algorithm, introduced very early in numeric method courses, is
application of Newton's method.  Much easier to code and
orders-of-magnitude faster.

Re: cubic root subroutine

Dan van Ginhoven wrote:

> Hi.
>
> I needed a cubic root subroutine but didnīt find one at a quick search.
> So I wrote a simple one.

Theres a trivial analytic solution
http://mathworld.wolfram.com/RootofUnity.html

gtoomey

Re: cubic root subroutine

>Hi.
>
>I needed a cubic root subroutine but didnīt find one at a quick search.
>So I wrote a simple one.
>It seems pretty fast and accurate enough.
>Any ideas  how to improve it,
>or point me to better solutions.
>

perldoc Math::Complex

Page down until you come to "The following extra operations are
supported on both real and complex numbers:"

Look in particular for cbrt() and root().

Very useful functions.

--
-john

February 28 1997: Last day libraries could order catalogue cards
from the Library of Congress.