edw@IUS1.CS.CMU.EDU (Eddie Wyatt) (01/12/88)
In article <521@cresswell.quintus.UUCP>, ok@quintus.UUCP (Richard A. O'Keefe) writes: > In article <622@PT.CS.CMU.EDU>, edw@IUS1.CS.CMU.EDU (Eddie Wyatt) writes: > 1> Actually any operation in a languge need not generate functions calls - > 1> have we heard of inline expansion? > ... > 2> implement. I don't think anyone will question the reasonableness > 2> of an infixed notion for exponentiation. > ... > 3> properties become harder to detect when using a functional format > 3> (integers must be typed coerced into doubles). On machines with > 3> no hardware support for exponetiation, the answer maybe one could > > His points (1) and (3) contradict each other. Points 1 and 3 do not contradict each other. A function call in C differs from an operator in that operators may be OVERLOADED or have IMPLICIT coercion of operands while functions calls must have EXPLICITLY type coercion its parameters. This means that 3.0 + 4 is valid while add(3.0,4) would not be valid (assuming add is a function taking two doubles). And hence determining that the exponent in the function call pow(base, (double) int_exp) could be equivalent to an integer is not as easy as having an integer value given. > In fact his point (3) is the direct opposite of C's defaults: > at present, whenever you do > X <operator> Y > "the usual conversions" are done, which means amongst other things > that if one operand is double, the other *will* be converted (which > Wyatt says *won't* happen), whereas if you call a function integers > will only be promoted to 'int'. In C, functional *syntax* gives us > the freedom to have operands of different types: We are talking about optimizing code. You can generate special code depending on the types of expressions present. There is nothing wrong with that provided that the code generated yields the same results in both the optimized and unoptimized versions. > extern double pow(double base, double expt); > -vs- extern double rtp(double base, int power); > Note that these are quite different functions: rtp(-3.0, 2) is well > defined, but pow(-3.0, 2.0) is illegal and should signal a matherr. > (In Common Lisp, (EXPT -3.0 2.0) --> #C(9.0 0.0), a complex number.) Sorry but there are not two separate functions. exp is just a partial function. It should yield real values for all integer exponents or bases >= 0. If Common Lisp choises to say expt(-3.0,2.0) is 9.0+0i that its problem. In the true math since (2.0 == 2). The behavior of the function should be the same for both 2.0 and 2. There was one good point in all that garbage you posted. That is noting that exp is partial over the domain of real values. Do we want to add an operator to the language that will be so ill-defined. -- Eddie Wyatt e-mail: edw@ius1.cs.cmu.edu
quiroz@cs.rochester.edu (Cesar Quiroz) (01/13/88)
Abstract--
There is no such thing as a unique, universally accepted
definition for the exponentiation. The case of
expt0: RealxInteger->Real
is different (in a mathematical way and from optimization
perspectives) from the one of
expt1: RealxReal->Complex
There is no such thing as a satisfactory
expt-: RealxReal->Real
C may either have two functions (pow for expt1 and the
proposed rtp for expt0) or have an operator that lets the
compiler choose statically between the two. In either case,
pow has to be partially defined, because of the lack of complex
arithmetic in the base language.
As C doesn't support complex arithmetic directly, we currently
cannot do things like making sqrt return complex values for the
negatives and real values for the positives. Instead, we just make
it set EDOM and be happy with it for a negative argument.
Actually, even if we had complex arithmetic, we wouldn't want it to
return differently typed results in different domains (a matter of
taste) so we would be burdened with either returning always complex
numbers (complicating the common case) or keeping with the EDOM error.
Why did I bring sqrt in here? Because the error domain of the real
valued sqrt is well accepted: all the negatives. With pow we don't
have that benefit. O'Keefe explained this recently, but Wyatt
missed the point completely. In a sentence, *there is more than one
possible exponentiation function*.
We cannot just have an operator that fails to compute pow on all the
negative bases, because the basic definition of exponentiation is
that of `repeated multiplication'. So, there are some people who
would rightly desire (pow(-2.0,3) == -8.0) to be true. Regrettably,
they cannot provide any meaningful value for pow(-2.0, 3.1), so they
might want to settle for EDOM or a complex result.
The most accepted (and perhaps most acceptable) extension of the
repeated multiplication definition, that provides meaning for all
combinations of bases and exponents, is the one embedded in Common
Lisp's EXPT function. Wyatt, and other people interested in how it
works, are referred to Steele's "Common Lisp: The Language". Wyatt
seems to believe the choice of such definition was whimsical and
that other languages could `choose' otherwise. Actually, he is
almost right, in that doing the right thing in C requires many more
changes to the language than is worth doing in the coming Standard.
Common Lisp `did it right'. There is a rule of complex
canonicalization, that lets you identify #C(12 0) with 12
(rationals, nor floats, all over). The natural extension to this
rule is to permit their sqrt to return sensible things in the
positive and negative domains. The same happens with expt, it is
defined in such a way that it is permissible to do intelligent
things with rational exponents and this does not break any semantic
properties of the rest of the language. But this is possible
due to the polymorphism of Lisp functions and their ability to query
the types of their arguments (remember that typeof discussion?)
The EXPT solution cannot be implemented in C, however, because C
functions cannot compute the types of their arguments. This out of
the way, there is still the lesser problem of not having complex
arithmetic. All in all, I don't think C as now exists can host a
decent implementation of exponentiation as a single function. It
might still be convenient to optimize the common case of an integer
exponent (to use the repeated multiplication rule) and then the rtp
function proposed by O'Keefe comes handy. No need for an
independent operator, at least not until a time comes to make more
radical changes to the language.
Some one (sorry, I didn't keep the article) already said that a new
operator is not an isolated change, but would require other
compatible changes, to avoid kludginess. (For instance, a
corresponding assignment). I hope people who want to introduce
exponentiation ops in C are willing to do it right and therefore
will refrain from giving the ANSI committee a reason to invent a new
feature without prior art. People interested in the problems of
embedding exponentiation in a language without complex arithmetic
are invited to follow up to comp.lang.misc or correspond by Email, I
will be happy to participate in such discussion. (Past experience
makes me add this: If you don't know why it is true that RealxReal
and Complex are *not* the same structure, please figure that out
first.)
--
Cesar Augusto Quiroz Gonzalez
Department of Computer Science ...allegra!rochester!quiroz
University of Rochester or
Rochester, NY 14627 quiroz@cs.rochester.edugwyn@brl-smoke.ARPA (Doug Gwyn ) (01/13/88)
In article <5825@sol.ARPA>, quiroz@cs.rochester.edu (Cesar Quiroz) writes: > There is no such thing as a unique, universally accepted > definition for the exponentiation. I don't buy this. I was under the impression that in the complex number field, which is the natural extension of the reals and therefore also of the integers, there was a single, well-defined (not necessarily well-behaved!), universally agreed-upon mathematical meaning for "w raised to the power z". If this isn't so, please elaborate. By the way, C does have what you labeled as expt-, the pow() function, even though its domain was restricted so that the result wouldn't be a complex number. I think the error domain of pow() is clear: anything that would have a result not representable as a double, basically, plus some arguments that are representable but wouldn't be correctly computed by the use of log() followed by exp() that is the usual implementation of pow(). The latter exclusion is unfortunate. > All in all, I don't think C as now exists can host a > decent implementation of exponentiation as a single function. Right! But all we have is pow(), unless we provide our own RPow() (real-to-integer, more generally useful than the LPow() I posted, which can be turned into RPow() easily enough) or write more involved code than should have been necessary.
edw@IUS1.CS.CMU.EDU (Eddie Wyatt) (01/14/88)
In article <5825@sol.ARPA>, quiroz@cs.rochester.edu (Cesar Quiroz) writes: > > Abstract-- > There is no such thing as a unique, universally accepted > definition for the exponentiation. The case of > expt0: RealxInteger->Real > is different (in a mathematical way and from optimization > perspectives) from the one of > expt1: RealxReal->Complex > There is no such thing as a satisfactory > expt-: RealxReal->Real > C may either have two functions (pow for expt1 and the > proposed rtp for expt0) or have an operator that lets the > compiler choose statically between the two. In either case, > pow has to be partially defined, because of the lack of complex > arithmetic in the base language. You are wrong. exp: (realXreal) -> real is WELL DEFINED but partial. What does well defined mean in this case : intuitively forever pair (x,y) x,y real, exp(x,y) evaluates a a unique real value or no real value at all. Lets consider the function exp2(complex X complex) -> complex z^x is defined to be e^(x log z) (the logrithm is open to choose). BUT for x integer, e^(x log z) = z**x forever branch of the log function. Where x**n is taken to be defined as : x**n : 1/(x**(-n)) : if (n < 0) x*(x**(n-1)) : if (n > 0) 1 : if (n = 0) What this says is that the definition of z^x is compadable with our definition of z**x when x is an integer. Lets consider the function exp(real X real) -> real y^x is defined to be e^(x log y) for y > 0. A real valued function. exp(y,x) is well defined for both y > 0 and x integer. However in the complex plane, exp(z,x) does not have an expectable definition because of the choose of logs. -- Eddie Wyatt e-mail: edw@ius1.cs.cmu.edu
quiroz@cs.rochester.edu (Cesar Quiroz) (01/14/88)
From article <7071@brl-smoke.ARPA> (gwyn@brl-smoke.ARPA (Doug Gwyn )):
:In article <5825@sol.ARPA>, quiroz@cs.rochester.edu (Cesar Quiroz) writes:
:> There is no such thing as a unique, universally accepted
:> definition for the exponentiation.
:
:I don't buy this.
:
Neither do I! I found myself arguing against my own position (this
bad habit of posting when late and lacking sleep).
The confusion came about this way: Somebody (I think E. Wyatt)
argued that expt(-2,3) could be chosen to be either -8 real or
-8 + 0i complex (which happens to project into a real, of
course). My article tried to defend the point that there is no such
choice, the universally accepted *definition* for exponentiation is
the complex one.
The problem is that the *implementation* of some special cases
(those with rational exponents) can be made very efficient, by
ignoring that it doesn't generalize into a continuous function in
any sensible way. Same goes for sqrt, where the NonNegative->Real
definition can be made very efficient (by a Newton-Raphson
iteration, for instance) while exponentiation to the 0.5 (or
negative arguments for the common sqrt) might involve more expensive
stuff. Catering for efficiency in such special cases I find
defensible (so the power functions you posted are desirable in the
library).
The noise about sqrt reduces to this: There are no negative numbers
for which we might want to define a real square root, so setting
EDOM or signalling a matherr for all the negatives is good enough.
With pow, some negative integer exponents produce complex numbers
that can be identified with reals, and so there are people who find
useful to believe that -2^3 = -8, although they don't know how to
define things like -2^3.1. Catering to this sloppiness, even if it
happens to argue for the same library functions, is what I found
indefensible.
--
Cesar Augusto Quiroz Gonzalez
Department of Computer Science ...allegra!rochester!quiroz
University of Rochester or
Rochester, NY 14627 quiroz@cs.rochester.eduok@quintus.UUCP (Richard A. O'Keefe) (01/14/88)
In article <7071@brl-smoke.ARPA>, gwyn@brl-smoke.ARPA (Doug Gwyn ) writes: > In article <5825@sol.ARPA>, quiroz@cs.rochester.edu (Cesar Quiroz) writes: > > There is no such thing as a unique, universally accepted > > definition for the exponentiation. > > I don't buy this. > > I was under the impression that in the complex number field, which is the > natural extension of the reals and therefore also of the integers, there > was a single, well-defined (not necessarily well-behaved!), universally > agreed-upon mathematical meaning for "w raised to the power z". If this > isn't so, please elaborate. > Sorry, but it *isn't* so. Any good book on complex analysis will explain why. I'll try to sketch the basic problem. We can represent any complex number (X + i*Y) in polar co-ordinates as (R, Theta) for some real R >= 0 and real Theta. Any value 2*pi*k+Theta will do, for any integer k. Now the nth root of a complex number is most easily taken in polar form: (R, Theta) ** (1/n) = (R**(1/n), Theta/n + (2*pi*k)/n) which means that every complex number other than 0 has n distinct nth roots. If we want a single-valued function, we have to pick which of these roots we want. It turns out that this cannot be done without introducing a discontinuity (look up "branch cut" in any good text). In the case of nth root, what we get is a ray (R,Theta) for fixed Theta and R > 0 where the nth root function is continuous from one side only; we can place that ray at any Theta we like. Different choices of branch cut are useful in different cases. Common Lisp spells out where the branch cuts are for all the built-in "elementary functions", which is one of the reasons why I mentioned Common Lisp in an earlier message. But if you want an nth root function on complex numbers, not only is there not a single well-defined meaning, there are uncountably many of them! When we restrict ourselves to the reals, a real number may have 0, 1, or 2 real nth roots. Only when we restrict ourselves to non-negative reals does "nth root" become a function. (R, Theta) ** n = (R ** n, Theta*n) for integer n _is_ well-defined. So is the exponential function exp, which is defined thus: exp(X + i*Y) = (e ** X)*(cos Y + i*sin Y) As you might expect, the break-down of pow() can be traced to the break-down of log(): if you try to find an inverse of exp() you get log((R, Theta)) = log(R) + i*(Theta + 2*pi*k) which has infinitely many values. The usual choice of branch cut for logarithm (yielding the so-called "principal value" of the logarithm) restricts Theta to -pi < Theta < +pi, so that the logarithm of negative numbers is not defined. There is no reason at all why we couldn't put the branch cut somewhere else, say -pi/2 < Theta < 3pi/2, which *would* give (-1 + 0i) a logarithm, but not (0 + 1i). Thus if one tried to define pow(z1, z2) as exp("log"(z1)*z2), "log" being the principal value, pow(-27.0, 1.0/3.0) would not be defined, which would be rather odd, as 4.3BSD and SunOS have a cbrt() function, and cbrt(-27.0) is perfectly well defined (-3.0). {This is why cbrt() exists, of course.} double rtp(x, n) /* this is the Raise To Power */ double x; /* function. The binary method */ int n; /* is used. Comments welcome. */ { double a; if (n < 0) x = 1.0/x, n = -n; if (n == 0) return 1.0; while ((n&1) == 0) x *= x, n >>= 1; for (a = x; --n > 0; a *= x) while ((n&1) == 0) x *= x, n >>= 1; return a; } I tried to measure the accuracy of this, but it turned out that the errors (~4e-16) were in my test data. You might lose about one decimal digit in relative error. It is possible in principle to get 1/2ulp accuracy for rtp(), but I don't know how to do that in C. C programmers on UNIX systems which also support F77 may be able to get at Fortran's (double precision)**(integer) function by doing extern double pow_di(/* double, int */); and then linking with -lF77 -lm Having rtp() {perhaps called pow_di, I wouldn't mind} in the C standard would be much cleaner.