computing t = pow(-11.5, .333)
John Smith
negative and the exponent is a fraction. In standard C this would
be something like t = pow(-11.5, .333), but with this combination
of arguments there is a domain error and the result is a
percolating NaN.
I worked around the problem by putting the standard function in a
wrapper:
double xpow(double b, double x)
{
double p;
if(b < 0)
{
p = pow(-b, x);
p *= -1.0;
}
else p = pow(b, x);
return p;
}
Is there any particular reason for this limitation in the
standard library function? Has anyone else had a problem with
this and found a better solution?
JS
pete
-11.5 has a third root which is negative,
if that's what you're really trying to find.
However, if you consider that
a negative number raised to an even integer power, is positive,
and that
a negative number raised to an odd integer power, is negative,
then it just seems natural that
the sign of a negative number raised to a non integer,
should be undefined.
--
pete
Robert Gamble
I don't know why the pow function is defined this way and there is
little insight in the rationale. You can obtain the expected result of
pow(-11.5, .333) by using the cube root function from C99 if your
implementation provide it: cbrt(-11.5);.
Robert Gamble
Robert Gamble
> John Smith wrote:
> > In a C program I need to do exponentiation where the base is
> > negative and the exponent is a fraction. In standard C this would
> > be something like t = pow(-11.5, .333), but with this combination
> > of arguments there is a domain error and the result is a
> > percolating NaN.
> >
> > I worked around the problem by putting the standard function in a
> > wrapper:
> >
> > double xpow(double b, double x)
> > {
> > double p;
> >
> > if(b < 0)
> > {
> > p = pow(-b, x);
> > p *= -1.0;
> > }
> > else p = pow(b, x);
> >
> > return p;
> > }
> >
> > Is there any particular reason for this limitation in the
> > standard library function? Has anyone else had a problem with
> > this and found a better solution?
>
> I don't know why the pow function is defined this way and there is
> little insight in the rationale.
I obviously didn't give enough thought to that. Any negative number
raised to a fractional power that is not the reciprocal of an odd
number (1/3, 1/5, 1/7, etc) will have a complex result. C99 provides
support for complex numbers and you can perform such calculations with
the cpow family of functions.
Eric Sosman
The reason for the limitation is mathematical: The proper
result of pow(-1,1.5) is -i, where i is the square root of -1.
That imaginary value is not in the range of double, so the pow
function cannot return it.
You define xpow(-1,1.5) as -pow(1,1.5), or -1. C allows
you to do so, but notice that you've now violated the usual
laws of exponents:
pow(x,a) * pow(x,b) == pow(x, a+b) /* mathematically */
so
pow(-1,1.5) * pow(-1,1.5)
== pow(-1,3.0)
== -1
but
xpow(-1,1.5) * xpow(-1,1.5)
== -pow(1,1.5) * -pow(1,1.5)
== -1 * -1
== 1
If you're happy with this state of affairs, go ahead and
be happy -- but don't ask the rest of the world to join in
your rejoicing.
--
Eric Sosman
eso...@acm-dot-org.invalid
Howard Hinnant
John Smith <JSm...@mail.net> wrote:
You might try something more along the lines of:
double root(double b, int x)
{
if(b < 0 && !(x & 1))
return 0./0.;
if (b < 0)
return -pow(-b, 1./x);
return pow(b, 1./x);
}
This will give you correct results for both:
root(-11.5, 3);
root(-11.5, 2);
root( 11.5, 2);
root( 11.5, 3);
-Howard
CBFalconer
Yes. What you are doing has no meaning. You are actually
returning -(b**x) instead of (-b)**x, where I am using ** to
represent exponentiation. What is the real problem?
Consider: what is the square root of -1/2? There is no answer in
the real plane. You need complex variables for this. Yet this is
just the value of (-0.5)**0.5.
--
"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/>
Dik T. Winter
> In a C program I need to do exponentiation where the base is
> negative and the exponent is a fraction. In standard C this would
> be something like t = pow(-11.5, .333), but with this combination
> of arguments there is a domain error and the result is a
> percolating NaN.
Rightly so.
> I worked around the problem by putting the standard function in a
> wrapper:
> double xpow(double b, double x)
> {
> double p;
>
> if(b < 0)
> {
> p = pow(-b, x);
> p *= -1.0;
> }
> else p = pow(b, x);
>
> return p;
> }
What is xpow(-1, 2) with your function? Does it divert from reality?
Or (talking about fractions) what about xpow(-1, 0.5)?
> Is there any particular reason for this limitation in the
> standard library function?
Yes, there is a pretty good reason. When the 'b' argument is negative
there may not even be a solution within the reals when the 'x' argument
is a fraction.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/
Jordan Abel
Or with such reciprocals - I believe that, traditionally, raising
a negative number to the fractional power yields the root with the
lowest [counterclockwise] angle on the complex plane from the positive
real axis - which is going to be less than pi (the angle for negative
numbers), while the nth-root operator yields the negative real root.
Julian V. Noble
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.
Although I have not looked into <math.h> for the details, I would
not be surprised to find that in many C implementations pow can
tell whether the exponent b is a floating point number or a positive
(or negative) integer smaller in magnitude than some value N, for
which a^b = a*...*a (computed via the minimal-multiplication algorithm)
is faster than computing two transcendental functions.
That is, I would expect a properly written pow function to compute
pow(a,2) = a*a;
pow(a,-2) = 1/(a*a);
whereas
pow(a,2.0) = exp(2.0*ln(a));
--
Julian V. Noble
Professor Emeritus of Physics
University of Virginia
pete
> That is, I would expect a properly written pow function to compute
>
> pow(a,2) = a*a;
> pow(a,-2) = 1/(a*a);
>
> whereas
>
> pow(a,2.0) = exp(2.0*ln(a));
It would take some kind of compiler trick,
as well as proper function writing, for a called function
to be able to distinguish between two arguments
which compare equal when converted to the type of the parameter,
but which have different types.
--
pete
Kenneth Brody
>
> John Smith wrote:
> >
> > In a C program I need to do exponentiation where the base is
> > negative and the exponent is a fraction. In standard C this would
> > be something like t = pow(-11.5, .333), but with this combination
> > of arguments there is a domain error and the result is a
> > percolating NaN.
> > Is there any particular reason for this limitation in the
> > standard library function? Has anyone else had a problem with
> > this and found a better solution?
It's not a "limitation in the standard library function", but rather
a "limitation" of math itself. You cannot express the fractional power
of a negative number as a real number in the general case.
> -11.5 has a third root which is negative,
> if that's what you're really trying to find.
Of course, a power of ".333" is _not_ a third root, as ".333" is only
_close_ to 1/3. (For large enough values of "close".)
> However, if you consider that
> a negative number raised to an even integer power, is positive,
> and that
> a negative number raised to an odd integer power, is negative,
> then it just seems natural that
> the sign of a negative number raised to a non integer,
> should be undefined.
Now, his solution of negating the number, and again negating the result,
does "work" for odd integral roots and odd integral powers, but not in
the general case.
[...]
--
+-------------------------+--------------------+-----------------------------+
| Kenneth J. Brody | www.hvcomputer.com | |
| kenbrody/at\spamcop.net | www.fptech.com | #include <std_disclaimer.h> |
+-------------------------+--------------------+-----------------------------+
Don't e-mail me at: <mailto:ThisIsA...@gmail.com>
John Smith
> John Smith wrote:
>
>>In a C program I need to do exponentiation where the base is
>>negative and the exponent is a fraction. In standard C this would
>>be something like t = pow(-11.5, .333), but with this combination
>>of arguments there is a domain error and the result is a
>>percolating NaN.
>>
>>I worked around the problem by putting the standard function in a
>>wrapper:
>>
>> double xpow(double b, double x)
>> {
>> double p;
>>
>> if(b < 0)
>> {
>> p = pow(-b, x);
>> p *= -1.0;
>> }
>> else p = pow(b, x);
>>
>> return p;
>> }
>>
>>Is there any particular reason for this limitation in the
>>standard library function? Has anyone else had a problem with
>>this and found a better solution?
>
>
> Yes. What you are doing has no meaning.
This is what happens when a non-mathematician attempts
mathematical programming.
You are actually
> returning -(b**x) instead of (-b)**x, where I am using ** to
> represent exponentiation. What is the real problem?
>
> Consider: what is the square root of -1/2? There is no answer in
> the real plane. You need complex variables for this. Yet this is
> just the value of (-0.5)**0.5.
>
My project involves solving cubic and quartic equations which
sometimes have complex roots and I am feeling my way along in the
dark. Occasionally there is a faint glimmer of light in the
distance. thanks to all for your responses. They have been helpful.
JS
Julian V. Noble
Not really, Fortran did it by looking for the decimal point. That is,
if it could interpret the literal exponent as an integer it would. If
it were forced to interpret it as a real (by the presence of a dp, e.g.)
it would do that. How it dealt with symbolic exponents doubtless had
to do with the type of the exponent name. Remember this is all done
at compile-time, not run-time. I once wrote a system where a variable
contained a token identifying its type so it could be done at run-time.
It really isn't that hard.
CBFalconer
> CBFalconer wrote:
>
... snip ...
>>
>> Yes. What you are doing has no meaning. You are actually
>> returning -(b**x) instead of (-b)**x, where I am using ** to
>> represent exponentiation. What is the real problem?
>>
>> Consider: what is the square root of -1/2? There is no answer in
>> the real plane. You need complex variables for this. Yet this is
>> just the value of (-0.5)**0.5.
>
> My project involves solving cubic and quartic equations which
> sometimes have complex roots and I am feeling my way along in the
> dark. Occasionally there is a faint glimmer of light in the
> distance. thanks to all for your responses. They have been helpful.
If you want real roots (which may not exist for a quartic, but must
exist for a cubic) consider a Newton-Raphson solution. Off to the
books with you. This is now no longer a C question, but
algorithmic, and better suited to comp.programming.
pete
>
> pete wrote:
> > Julian V. Noble wrote:
> >
> >> That is, I would expect a properly written pow function to compute
> >>
> >> pow(a,2) = a*a;
> >> pow(a,-2) = 1/(a*a);
> >>
> >> whereas
> >>
> >> pow(a,2.0) = exp(2.0*ln(a));
> >
> > It would take some kind of compiler trick,
> > as well as proper function writing, for a called function
> > to be able to distinguish between two arguments
> > which compare equal when converted to the type of the parameter,
> > but which have different types.
> >
>
> Not really, Fortran did it by looking for the decimal point.
How does C do it?
--
pete
P.J. Plauger
news:446A46...@mindspring.com...
Essentially by calling trunc(x) and seeing if the return value
differs from x. (More precisely, we use the _Dint function that
was discussed at length in these pages a few months ago.)
P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com
pete
There's no function that can be called from pow,
that will distinguish
pow(a, 2);
from
pow(a, 2.0);
--
pete
pete
I'm not very familiar with C99, so for a moment
I thought trunc() was one of your own, like _Dint.
Now that I've read the standard description of trunc,
I can say that there's absolutley no way that calling trunc
can help to distinguish
P.J. Plauger
Right. I was addressing the issue of determining whether you have
an exact integer power, however it was arrived at.
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)
--
pete
pete
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.
--
pete
Julian V. Noble
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.
pete
> >>>>>>> pete wrote:
> >>>>>>>> Julian V. Noble wrote:
> >>>>>>>>
> >>>>>>>>> That is,
> >>>>>>>>> I would expect a properly written pow function to compute
> >>>>>>>>>
> >>>>>>>>> pow(a,2) = a*a;
> >>>>>>>>> pow(a,-2) = 1/(a*a);
> >>>>>>>>>
> >>>>>>>>> whereas
> >>>>>>>>>
> >>>>>>>>> pow(a,2.0) = exp(2.0*ln(a));
> >>>>>>>> It would take some kind of compiler trick,
> >>>>>>>> as well as proper function writing, for a called function
> >>>>>>>> to be able to distinguish between two arguments
> >>>>>>>> which compare equal when converted
> >>>>>>>> to the type of the parameter,
> >>>>>>>> but which have different types.
> >>>>>>>>
> >>>>>>> Not really, Fortran did it by looking for the decimal point.
> >>>>>> How does C do it?
> 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.
--
pete
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.
P.J. Plauger
news:e4fbop$4e9$1...@murdoch.acc.Virginia.EDU...
> 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
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.
pete
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++.
--
pete
Keith Thompson
> pete wrote:
>>
> ... 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.
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 (The_Other_Keith) ks...@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.
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.
Jordan Abel
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
I'm not sure what you mean by that.
--
"Welcome to the wonderful world of undefined behavior, where the demons
are nasal and the DeathStation users are nervous." --Daniel Fox
Keith Thompson
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
news:slrne6nasp.2...@random.yi.org...
> 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.
Jordan Abel
> Jordan Abel <rand...@gmail.com> writes:
>> On 2006-05-17, pete <pfi...@mindspring.com> wrote:
>>> 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.
>
> 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
> "Jordan Abel" <rand...@gmail.com> wrote in message
> news:slrne6nasp.2...@random.yi.org...
>
>> 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.
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
[ 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) !
Charles Richmond
>
> John Smith wrote:
> > CBFalconer wrote:
> >
> ... snip ...
> >>
> >> Yes. What you are doing has no meaning. You are actually
> >> returning -(b**x) instead of (-b)**x, where I am using ** to
> >> represent exponentiation. What is the real problem?
> >>
> >> Consider: what is the square root of -1/2? There is no answer in
> >> the real plane. You need complex variables for this. Yet this is
> >> just the value of (-0.5)**0.5.
> >
> > My project involves solving cubic and quartic equations which
> > sometimes have complex roots and I am feeling my way along in the
> > dark. Occasionally there is a faint glimmer of light in the
> > distance. thanks to all for your responses. They have been helpful.
>
> If you want real roots (which may not exist for a quartic, but must
> exist for a cubic) consider a Newton-Raphson solution. Off to the
> books with you. This is now no longer a C question, but
> algorithmic, and better suited to comp.programming.
>
OT: I would recommend Bairstow's method. Google it or see Wikipedia.
--
+----------------------------------------------------------------+
| Charles and Francis Richmond richmond at plano dot net |
+----------------------------------------------------------------+
pete
> 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.)
My portable freestanding version of pow,
which uses my freestanding versions of exp and log,
has trouble matching the accuracy my implementation's pow.
pow(0.0001, -0.25) - 10 is 0.000000e+000
fs_pow(0.0001, -0.25) - 10 is 3.552714e-015
There's a lot of squaring in my fs_exp function,
which I suspect is the main problem.
Portable C code doesn't have access to any types
which are guaranteed to be more precise than double.
double fs_exp(double x)
{
unsigned n, square;
double b, e;
static double x_max;
if (1 > x_max) {
x_max = fs_log(DBL_MAX);
}
if (x_max >= x && x >= -x_max) {
for (square = 0; x > 1; x /= 2) {
++square;
}
while (-1 > x) {
++square;
x /= 2;
}
e = b = n = 1;
do {
b /= n++;
b *= x;
e += b;
b /= n++;
b *= x;
e += b;
} while (b > DBL_EPSILON / 4);
while (square-- != 0) {
e *= e;
}
} else {
e = x > 0 ? DBL_MAX : 0;
}
return e;
}
--
pete
pete
> > Portable C code doesn't have access to any types
> > which are guaranteed to be more precise than double.
>
> Unless you write a portable set of extra precision functions...
Thank you.
That's like the next level of a hard video game for me.
But I only write these things for my own amusement anyway,
so I guess it's time to go to the next level.
--
pete