comp.lang.ada
 help / color / mirror / Atom feed
* Math error in Annex G ???
@ 1996-07-24  0:00 Joel VanLaven
  1996-07-25  0:00 ` Robert A Duff
  1996-07-25  0:00 ` Karl Cooper {46901}
  0 siblings, 2 replies; 4+ messages in thread
From: Joel VanLaven @ 1996-07-24  0:00 UTC (permalink / raw)



  I was implementing complex types from Annex G in the LRM as a way to learn
some Ada95 when I ran into something that I am sure is incorrect.  It has
to do with exponentiation of a complex number by a negative integer by going
through polar.  Essentially it says that exponentiation can be significantly
improved by converting to polar, using a polar exponentiation trick, and
then converting back.  That is fine, and I knew about that even without the
LRM mentioning it.  The problem is that the method that they give for
exponentiation with the polar of a complex number is wrong for negative
exponents.

  The LRM says (liberty taken with spacing):
"... exponentiating the modulus by the given exponent;
 multiplying the argument by the given exponent,
   when the exponent is positive,
 or dividing the argument by the absolute value of the given exponent,
   when the exponent is negative; ..."

Translating to mathese:
if n > 0 then
   (r,theta)**n=(r**n,n*theta)
if n < 0 then
   (r,theta)**n=(r**n,theta/abs(n))

but i**(-1) = 1/i = -i
and i**(-1) = (1,pi/2)**(-1) ?=? (1**(-1),(pi/2)/abs(-1)) = (1.pi/2) = i
obviously something is wrong there.  

I am almost positive that the desired math is:
for all n
   (r,theta)**n=(r**n,n*theta)

this comes about from how * works for polar complex numbers, namely:

(r1,th1)*(r2,th2)=(r1*r2,th1+th2)

this will prove it by induction for positive n as follows:

-- base case
(r,th)**1 = (r,th) = (r**1,1*th)

--inductive case
if (r,th)**n=(r**n,n*th) then
  (r,th)**(n+1)=((r,th)**n)*(r,th) = (r**n,n*th)*(r,th) =
  ((r**n)*r,(n*th)+th) = (r**(n+1),(n+1)*th)
Proven.

for negative n, simply use the polar definition of * and the part
for positive n as follows:
(r',th') = (r,th)**(-n) = 1/((r,th)**n) = 1/(r**n,n*th)
(r',th')*(r**n,n*th) = 1 = (1,0)
(r'*(r**n),th'+(n*th)) = (1,0)
r'*(r**n) = 1  and  th'+(n*th) = 0
r' = 1/(r**n) = r**(-n)  and  th' = -(n*th) = (-n)*th
(r,th)**(-n) = (r',th') = (r**(-n),(-n)*th)
Proven.

I have purposely not dealt with 0 exponents as they were not included in
this part of the LRM but mathematically it works the same way.

I submitted this to ada-comment a week ago, but I have a feeling that
comp.lang.ada might be a little more active...

If anyone can show me that the LRM is correct I'd love to see it.
I am practically positive that the LRM is just plain wrong math-wise,
but all proofs can have flaws, so if you want to give your complex
arithmetic skills a whirl, have at it.

-- Joel




^ permalink raw reply	[flat|nested] 4+ messages in thread

* Re: Math error in Annex G ???
  1996-07-24  0:00 Math error in Annex G ??? Joel VanLaven
  1996-07-25  0:00 ` Robert A Duff
@ 1996-07-25  0:00 ` Karl Cooper {46901}
  1 sibling, 0 replies; 4+ messages in thread
From: Karl Cooper {46901} @ 1996-07-25  0:00 UTC (permalink / raw)



Joel VanLaven wrote:

>   The LRM says (liberty taken with spacing):
> "... exponentiating the modulus by the given exponent;
>  multiplying the argument by the given exponent,
>    when the exponent is positive,
>  or dividing the argument by the absolute value of the given exponent,
>    when the exponent is negative; ..."
<snip...>
> I am almost positive that the desired math is:
> for all n
>    (r,theta)**n=(r**n,n*theta)
> 

Yes, your mathematics is correct, and the permission given in the LRM
would result in errors if actually adopted by a language implementor.

The LRM also permits the implementor to use repeated complex multiplication,
with a possible final complex reciprocation (when the exponent is negative),
which would result in correct (if potentially inefficient) calculation.

And, of course, no implementor is bound by the permissions.

Does anyone know if this slip in the LRM was actually adopted in any
particular Ada 95 implementation?  Shall we all go out and test our compiler?

I have sometimes wondered why the exponent is limited to integral values,
but perhaps we should be grateful that the LRM didn't try to explain raising
a complex number to a complex power.




^ permalink raw reply	[flat|nested] 4+ messages in thread

* Re: Math error in Annex G ???
  1996-07-24  0:00 Math error in Annex G ??? Joel VanLaven
@ 1996-07-25  0:00 ` Robert A Duff
  1996-07-26  0:00   ` Robert I. Eachus
  1996-07-25  0:00 ` Karl Cooper {46901}
  1 sibling, 1 reply; 4+ messages in thread
From: Robert A Duff @ 1996-07-25  0:00 UTC (permalink / raw)



In article <1996Jul24.214556.24621@ocsystems.com>,
Joel VanLaven <jvl@ocsystems.com> wrote:
>I submitted this to ada-comment a week ago, but I have a feeling that
>comp.lang.ada might be a little more active...

I saw your message to ada-comment, but I was hoping somebody with
expertise in numerics would answer your question (e.g. Robert Eachus
and/or Robert Dewar).  Ada-comment is active, but often somewhat slow
(e.g. months).

- Bob




^ permalink raw reply	[flat|nested] 4+ messages in thread

* Re: Math error in Annex G ???
  1996-07-25  0:00 ` Robert A Duff
@ 1996-07-26  0:00   ` Robert I. Eachus
  0 siblings, 0 replies; 4+ messages in thread
From: Robert I. Eachus @ 1996-07-26  0:00 UTC (permalink / raw)



In article <Dv45up.36o@world.std.com> bobduff@world.std.com (Robert A Duff) writes:

  > I saw your message to ada-comment, but I was hoping somebody with
  > expertise in numerics would answer your question (e.g. Robert Eachus
  > and/or Robert Dewar).  Ada-comment is active, but often somewhat slow
  > (e.g. months).

   There is an error, but not a serious error, and we should take the
time to correct it right.  (Or throw it over the fence to the NRG
people.)  But I wasn't going to respond until I was sure what change
to recommend.  So thanks for finding the problem, but fixing it will
procede with due deliberation.


--

					Robert I. Eachus

with Standard_Disclaimer;
use  Standard_Disclaimer;
function Message (Text: in Clever_Ideas) return Better_Ideas is...




^ permalink raw reply	[flat|nested] 4+ messages in thread

end of thread, other threads:[~1996-07-26  0:00 UTC | newest]

Thread overview: 4+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
1996-07-24  0:00 Math error in Annex G ??? Joel VanLaven
1996-07-25  0:00 ` Robert A Duff
1996-07-26  0:00   ` Robert I. Eachus
1996-07-25  0:00 ` Karl Cooper {46901}

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox