On Wednesday, 95/12/06, Mike Bilow wrote to Mike Fisher about "Cube
root" as follows:

MB>  MF>                root:=exp(ln(P)/Q);
MB> 
MB> Am I really missing something here, or does Pascal have no 
MB> support for direct exponentiation?

Hi Mike,

You have just stumbled on one reason why I have felt for 20 years that
Pascal is a little too limited to be a real-world programming
language.

MB> The Qth root of P is 
MB> equal to P raised to the 1/Q power.  Using C, it is a heck 
MB> of a lot more efficient to compute a Qth root as such:
MB> 
MB>    double P, Q, root;
MB> 
MB>    root = pow(P, 1.0 / Q);

This is correct about pow() being more efficient than exp(ln()), but
it is not necessarily optimal.

MB> The direct exponentiation form uses a double precision 
MB> special case division -- taking a reciprocal -- which 
MB> corresponds to a machine instruction on the Intel math 
MB> coprocessor, while the logarithmic form uses a double 
MB> precision division that is not a special case.

The savings of using pow() are much greater than the handfull of CPU
cycles difference between a reciprocal and a full division.

MB> The direct 
MB> exponentiation form also incurs one function call, while 
MB> the logarithmic form incurs two.

This is also true. But there is more ...

Since I had nothing better to do today, I cooked up a couple of
programs that implemented the Newton-Raphson algorithm for solving
cube roots and compared this with raising the number to a one-third
power.

Firstly, I coded it in FORTRAN IV and used MS FORTRAN 5.1 to produce a
16:16 OS/2 executable. This yielded a slightly unexpected result that
the N-R method took about 25% longer than exponentiation. From this, I
gathered that exponentiation is handled by means other than calls to
log() and exp().

I then re-coded it in C and used Borland C/C++ 3.1 to produce a 16:16
DOS program. This gave me the more expected result that the N-R method
took about 40% less time than using pow(). I think that pow() is
implemented in the BCAF 3.1 run-time as exp(log()) without any further
thought on Borland's part.

So, I copied the code to use Borland C/C++ for OS/2 1.5 to produce a
0:32 OS/2 executable. This gave a very surprising result indeed! The
N-R method took 150% more time (two and a half times as long) to
execute than did pow().

Finally, I de-Borlandised the C source and compiled using C Set ++
2.01. This gave the even more extreme result that the N-R method took
more than 3 times as long as using pow(). What is more, the pow()
supplied by IBM ran almost twice as fast as the pow() supplied by
Borland for OS/2.

From the foregoing, I would conclude that implemnters of C/C++ do not
use exp(log()) any longer for coding pow(). Instead, they use some
more specialised algorithm, possibly with special cases for [small]
integer roots.

But, Pascal programmers cannot readily exploit this, since there is no
exponentiation supported directly by the language. So their situation
is like that of BCAF 3.1 programmers. Indeed, I would expect the code
in the BP7 RTL to be nearly identical to that in the BCAF 3.1 RTL, but
with no means of condensing the logarithm and exponentiation into a
single call. In fact, Wirth himself suggests using EXP(y*LN(x)) to
calculate x^y. Hence, the Newton-Raphson algorithm would be more
efficient for obtaining cube roots. [Strangely enough, Wirth
implemented the SQR() function to calculate x*x. When designing Pascal,
he seems to have ignored the significant to focus on the trivial.]

There is still one gotcha. The pow() function of C/C++ and the **
operator of FORTRAN do not like non-integer exponents of negative
numbers. You have to take absolute value and then re-apply the sign,
which is a little clumsy. The N-R method works just peachy keen with
negative arguments.

Therefore, I attach [Pascal] source code for the Newton-Raphson cube
root.

BTW, 'fair dinkum' means 'genuine'. ... :-)

Regards

Dave


========================== CUBRT.PAS =========================

{ A function to calculate cube roots
  using the Newton-Raphson method.}

FUNCTION CUBRT(N : REAL) : REAL;

CONST
     eps = 1E-12; { Adjust for hardware precision }
VAR
     x, xnew : REAL;
BEGIN
     IF N = 0 THEN
          xnew := 0
     ELSE IF N > 0 THEN
          xnew := 1
     ELSE
          xnew := -1;

     IF xnew  0 THEN
     REPEAT
          x := xnew;
          xnew := (x + x + N/SQR(x))*(1E0/3E0)
     UNTIL ABS(x-xnew) >= ABS(x)*eps;

     CUBRT := xnew
END;

==============================================================

 * KWQ/2 1.2i * Program call to load Windows NT- "Here_piggy_piggy_piggy"
--- Maximus/2 3.00