cubic root subroutine

[Dan van Ginhoven]

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) = @_;
my ($reply);
$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 ) {
$reply +=1;
}
$reply -=1;

# 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);
}
return $reply+$decimal;
}

 

[A. Sinan Unur]

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.

 

[Scott Bryce]

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);
}

 

[A. Sinan Unur]

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

 

[A. Sinan Unur]

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.

 

[Gregory Toomey]

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

 

[Jürgen Exner]

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

 

[Jay Tilton]

: 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.

 

[John M. Gamble]

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.

 

[Arndt Jonasson]

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?

 

[Scott Bryce]

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.

 

[Arndt Jonasson]

Same on Win32.


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

 

[A. Sinan Unur]

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{%.0f}, ((64)**(1/3))"
4

Sinan.

 

[Richard Gration]

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

<SNIP>

Here's one using Newton-Raphson:
(http://www.sosmath.com/calculus/diff/der07/der07.html)

#!/usr/bin/perl

use strict;
use warnings;

# Program to find the cubic root of a number using the Newton-Raphson
# method for approximating the roots of polynomials
#
# Here, f(x) = x**3 - n where n is the number to be rooted
# thus f'(x) = 3x**2
#
# [ f'(x) is the first differential of f(x) wrt x ]

my $findcubicrootof = $ARGV[0] || 27;
my $first_guess = $findcubicrootof / 3;
my $accuracy = 0.000000000000001;

# Arbitrarily set to what works on my machine. YMMV
die "Accuracy exceeds machine limits\n" if ($accuracy < 1e-15);

print "The cubic root of $findcubicrootof is ",
nr($first_guess,$findcubicrootof,$accuracy),
" +/- $accuracy\n";

sub nr {
my ($guess,$n,$accuracy) = @_;

# Find delta
my $delta = fofx($guess,$n) / fdashofx($guess,$n);

if (abs($delta) < $accuracy) {
return $guess;
} else {
return nr($guess - $delta,$n,$accuracy);
}
}

# Returns f(x)
sub fofx {
my ($x,$n) = @_;
return $x ** 3 - $n;
}

# Returns f'(x)
sub fdashofx {
my ($x,$n) = @_;
return 3 * $x ** 2;
}

 

[Arndt Jonasson]

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{%.0f}, ((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.

 

[A. Sinan Unur]

....


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{%.16f}, ((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.

 

[Arndt Jonasson]

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{%.16f}, ((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
my "problem" as the original posting, your reply would be 100%
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.

 

[Scott Bryce]

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.

 

[A. Sinan Unur]

....

....

(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

 

[John M. Gamble]

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

D:\> perl -e "printf q{%.0f}, ((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