computing t = pow(-11.5, .333)

[pete]

Right. I was addressing the issue of determining whether you have
an exact integer power, however it was arrived at.

That wasn't the issue under discussion.

This issue was a properly written pow
computing a*a, for pow(a,2) and
computing exp(2.0*ln(a)), for pow(a,2.0)

 

[pete]

That wasn't the issue under discussion.

This issue was a properly written pow
computing a*a, for pow(a,2) and
computing exp(2.0*ln(a)), for pow(a,2.0)

Unless I misunderstood,
and computing exp(2.0*ln(a)), for pow(a,2.0)
was supposed to be an example of an improperly written pow.

 

[Julian V. Noble]

Unless I misunderstood,
and computing exp(2.0*ln(a)), for pow(a,2.0)
was supposed to be an example of an improperly written pow.

If you recall the original post (and I am having difficulty doing
so by this time ;-) ) I was pointing out that it is not hard to
make the compiler choose which version of pow to use in a given
situation. I noted that Fortran allowed the programmer to pre-empt
the choice by specifying a real, rather than an integer literal.

For exponents that are variable names, the typing features of the
language are used to specify which version of the function to
compile.

The question of "proper" or "improper" forms of pow boils down to
which is faster, successive multiplication (for integer exponent)
or exponentiation of a logarithm. The first wins (at least on
many Intel cpu's) when the integer power is smaller than about
50-60. The latter wins, even for integers, when the integer is
large enough. Of course if the exponent is a variable integer,
the compiler has no way to tell in advance how big it is. So
either a definite form must be compiled, or a version with a
runtime decision. The latter need not be time-costly, but will
definitely increase the code size and complexity.

 

[CBFalconer]

.... snip ...

That wasn't the issue under discussion.

This issue was a properly written pow
computing a*a, for pow(a,2) and
computing exp(2.0*ln(a)), for pow(a,2.0)

It is impossible to discriminate on those calls within pow, since
the arguments are of type double. Therefore the first call is
automatically converted into the second long before the pow routine
receives control.

--
"If you want to post a followup via groups.google.com, don't use
the broken "Reply" link at the bottom of the article. Click on
"show options" at the top of the article, then click on the
"Reply" at the bottom of the article headers." - Keith Thompson
More details at: <http://cfaj.freeshell.org/google/>
Also see <http://www.safalra.com/special/googlegroupsreply/>

 

[pete]

Julian V. Noble wrote:

If you recall the original post (and I am having difficulty doing
so by this time ;-) ) I was pointing out that it is not hard to
make the compiler choose which version of pow to use in a given
situation.

That would be what I was refering to by
"some kind of compiler trick"

I noted that Fortran allowed the programmer to pre-empt
the choice by specifying a real, rather than an integer literal.

I don't think there's anything like that available in C.
I think the best that pow can do,
concerning the evaluation of it's arguments,
is to determine whether or not a value of type double
has a fractional part.

 

[P.J. Plauger]

The question of "proper" or "improper" forms of pow boils down to
which is faster, successive multiplication (for integer exponent)
or exponentiation of a logarithm. The first wins (at least on
many Intel cpu's) when the integer power is smaller than about
50-60. The latter wins, even for integers, when the integer is
large enough. Of course if the exponent is a variable integer,
the compiler has no way to tell in advance how big it is. So
either a definite form must be compiled, or a version with a
runtime decision. The latter need not be time-costly, but will
definitely increase the code size and complexity.

We used to do successive multiplication, until we discovered that
accuracy degrades quickly with that approach. The only really
safe way to evaluate pow(x, y) accurately is to compute
exp(y * ln(x)) to extra internal precision. (Took us a lot of
rewrites, and a lot of careful testing, to find that out.)

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[P.J. Plauger]

That wasn't the issue under discussion.

This issue was a properly written pow
computing a*a, for pow(a,2) and
computing exp(2.0*ln(a)), for pow(a,2.0)

Really? Why would you want to compute a*a only if the second argument
historically began as an integer 2? By the time you get to the body
of pow, 2 and 2.0 look the same. (And that's what I thought you were
fussing about, that you don't know the early history of the argument.)
So IMO I *was* speaking to the issue -- however you get the second
argument, if it's an exact integer you may well want to compute the
function value differently than if it's not. You certainly need to
know when the exponent is an exact integer for certain error checks.

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[pete]

Really?

Yes.
Julian V. Noble clarified the case elsewhere on this thread,
in a post that you already responded to,
when I said that I might be misunderstanding.

"I was pointing out that it is not hard to
make the compiler choose which version of pow to use in a given
situation." -- Julian V. Noble

I believe he's describing a situation similar
to the "function overloading" feature of C++.

 

[Keith Thompson]

It is impossible to discriminate on those calls within pow, since
the arguments are of type double. Therefore the first call is
automatically converted into the second long before the pow routine
receives control.

A language with function overloading could easily distinguish between
pow(a, 2) and pow(a, 2.0), calling two different functions. C, of
course, doesn't have operator overloading.

A language with an exponentiation operator, say "**", could do
something similar; it could distinguish between a**2 and a**2.0 just
as C distinguishes between a*2 and a*2.0 (a limited form of operator
overloading, available only for built-in operators). C, of course,
doesn't have an exponentiation operator.

C99's <tgmath.h> provides a type-generic macro pow() that invokes one
of powf(), pow(), or powl(), depending on the types of the arguments.
It doesn't distinguish between pow(a, 2) and pow(a, 2.0), but I
suppose it could have been defined to do so.

 

[Keith Thompson]

To understand this you must look under the hood. The power function
with an exponent specified as a decimal fraction almost certainly
will do this:

pow(a,b) = exp(b*ln(a));

since the ordinary (real) logarithm of a negative number is undefined
in many implementations it will return NaN. I myself would make it
abort immediately with an error message such as

"Can't take logarithm of negative number!"

(and I have done just that in floating point libraries I have written)
but de gustibus non est disputandum, and others doubtless had good
reasons for going the NaN route.

exp(b*ln(a)) is the obvious and mathematically correct way to
implement pow(a, b), but I've been told that it's not the best way to
implement it, probably for reasons having to do with numerical
stability. Sorry, I don't have any details. (Elsethread,
P.J. Plauger mentions using this formula with extra precision.)

 

[CBFalconer]

.... snip ...

We used to do successive multiplication, until we discovered that
accuracy degrades quickly with that approach. The only really
safe way to evaluate pow(x, y) accurately is to compute
exp(y * ln(x)) to extra internal precision. (Took us a lot of
rewrites, and a lot of careful testing, to find that out.)

I had the same problem years ago with a mechanism for converting a
floating value to a suitable herd of characters for printing. My
approach was to first convert the float, with a binary exponent, to
a 'decimalfloat' with a power of 10 exponent (always with a binary
significand). This was done by multiplication or division by 10.
The accuracy was not good. So I converted to doing initial
multiplication/division by 1e5 or so (it was the largest exact
power of 10 expressible in my floating point, which had a 16 bit
significand). This got the errors under control. Since the
systems accuracy was only 4.7 decimal digits, I couldn't really
afford to lose one. Everything was reversible, so the same
mechanism was used for input.

--
"If you want to post a followup via groups.google.com, don't use
the broken "Reply" link at the bottom of the article. Click on
"show options" at the top of the article, then click on the
"Reply" at the bottom of the article headers." - Keith Thompson
More details at: <http://cfaj.freeshell.org/google/>
Also see <http://www.safalra.com/special/googlegroupsreply/>

 

[Jordan Abel]

If you recall the original post (and I am having difficulty doing
so by this time ;-) ) I was pointing out that it is not hard to
make the compiler choose which version of pow to use in a given
situation. I noted that Fortran allowed the programmer to pre-empt
the choice by specifying a real, rather than an integer literal.

For exponents that are variable names, the typing features of the
language are used to specify which version of the function to
compile.

The question of "proper" or "improper" forms of pow boils down to
which is faster, successive multiplication (for integer exponent)
or exponentiation of a logarithm. The first wins (at least on
many Intel cpu's) when the integer power is smaller than about
50-60. The latter wins, even for integers, when the integer is
large enough. Of course if the exponent is a variable integer,
the compiler has no way to tell in advance how big it is. So
either a definite form must be compiled, or a version with a
runtime decision. The latter need not be time-costly, but will
definitely increase the code size and complexity.

Aren't there tricks you can do with successive multiplication?

x^23 = x*x^2*x^4*x^16

a=x*x
b=a*a
c=b*b
d=c*c
e=a*b
f=c*d
return e*f

I bet there's a clever way to do this in a loop that generalizes to all
cases.

 

[Jordan Abel]

I believe he's describing a situation similar
to the "function overloading" feature of C++.

C has function overloading. It's just not available to mere mortal
programmers.

 

[Ben Pfaff]

C has function overloading. It's just not available to mere mortal
programmers.

I'm not sure what you mean by that.

 

[Keith Thompson]

C has function overloading. It's just not available to mere mortal
programmers.

C certainly has operator overloading (in a form not available to mere
mortal programmers). I'm not aware of anything in C that looks like
function overloading, unless you count C99's <tgmath.h>.

 

[P.J. Plauger]

Aren't there tricks you can do with successive multiplication?

x^23 = x*x^2*x^4*x^16

a=x*x
b=a*a
c=b*b
d=c*c
e=a*b
f=c*d
return e*f

I bet there's a clever way to do this in a loop that generalizes to all
cases.

Yes, it's the classic squaring loop, where you fold in x^(2^n) if the
n-th bit is nonzero. But you still lose too much precision that way.

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[Jordan Abel]

C certainly has operator overloading (in a form not available to mere
mortal programmers). I'm not aware of anything in C that looks like
function overloading, unless you count C99's <tgmath.h>.

That was precisely what I was talking about.

 

[Jordan Abel]

Yes, it's the classic squaring loop, where you fold in x^(2^n) if the
n-th bit is nonzero. But you still lose too much precision that way.

Could there be more clever ways?

x^3*(((x^3*x^2)^2)^2)
a=x*x //2
b=x*a //3
c=a*b //5
d=c*c //10
e=d*d //20
return e*b

That takes off a step.

 

[Julian V. Noble]

Jordan Abel wrote:

[ snip ]

I bet there's a clever way to do this in a loop that generalizes to all
cases.

Yes. I don't know how clever ;-)

In pseudocode,

pwr(x,n) \\ raise x to integer power -- recursive
\\ algorithm: x^n = (x^2)^[n/2] * x^[n mod 2]
if (n=0) then return 1e0 exit
if (n=1) then return x exit
return pwr(x, n % 2) * pwr(x*x, n/2) exit


pospow(x,y) \\ raise x to y power, x >=0, y >=0
if ((y = integer) & (y < 50))
then return pwr(x,y) exit
else return exp(y*ln(x) exit


pow(x,y) \\ raise x >=0 to power y
if (x < 0) then error_message exit
if (y < 0)
then return 1/pospow(x,abs(y)) exit
else return pospow(x,y) exit

 

[Julian V. Noble]

Jordan Abel wrote:

[ snip ]

I bet there's a clever way to do this in a loop that generalizes to
all cases.

Yes. I don't know how clever ;-)

In pseudocode,

pwr(x,n) \\ raise x to integer power -- recursive
\\ algorithm: x^n = (x^2)^[n/2] * x^[n mod 2]
if (n=0) then return 1e0 exit
if (n=1) then return x exit
return pwr(x, n % 2) * pwr(x*x, n/2) exit


pospow(x,y) \\ raise x to y power, x >=0, y >=0
if ((y = integer) & (y < 50))
then return pwr(x,y) exit
else return exp(y*ln(x) exit


pow(x,y) \\ raise x >=0 to power y
if (x < 0) then error_message exit ^^^^^^^
if (y < 0)
then return 1/pospow(x,abs(y)) exit
else return pospow(x,y) exit

Oops! It should be if (x <= 0) !