From: John B. Matthews on
In article <Xns9DA6ACEFDCF8BWarrensBlatherings(a)81.169.183.62>,
Warren <ve3wwg(a)gmail.com> wrote:

> Warren expounded in news:Xns9DA6AC2A2C8BWarrensBlatherings(a)81.169.183.62:
>
> > Adam Beneschan expounded in news:73a0af1d-7213-44af-90fa-ed6de4c64ce8
> > @b4g2000pra.googlegroups.com:
> >
> >> On Jun 29, 12:31 pm, Warren <ve3...(a)gmail.com> wrote:
> >>>
> >>> > So, what is the "missing" function?
> >>>
> >>> > <http://www.adaic.com/standards/05rm/html/RM-G-1-2.html>
> >>>
> >>> 1) The "inverse" of a complex number.
> >>
> >> Isn't that just 1.0 / X?
> >
> > Ok, it seems to be, as the GSL (Gnu Scientific Library)
> > defines it as:
> >
> > gsl_complex
> > gsl_complex_inverse (gsl_complex a)
> > { /* z=1/a */
> > double s = 1.0 / gsl_complex_abs (a);
> >
> > gsl_complex z;
> > GSL_SET_COMPLEX (
> > &z, (GSL_REAL (a) * s) * s, -(GSL_IMAG (a) * s) * s);
> > return z;
> > }
>
> Apologies for following up my own post, but
> looking at the GSL implementation of complex
> division:
> [...]
> Given a=1.0, the inverse function is definitely
> higher performance. It's simpler calculation
> probably implies more accuracy as well.

The more comparable Ada implementation (from GNAT) would be

function "/" (Left : Real'Base; Right : Complex) return Complex is
a : constant R := Left;
c : constant R := Right.Re;
d : constant R := Right.Im;
begin
return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
Im => -((a * d) / (c ** 2 + d ** 2)));
end "/";

--
John B. Matthews
trashgod at gmail dot com
<http://sites.google.com/site/drjohnbmatthews>
From: Damien Carbonne on
Le 29/06/2010 23:00, Warren a �crit :
> Warren expounded in news:Xns9DA6AC2A2C8BWarrensBlatherings(a)81.169.183.62:
>>>>
>>>> 1) The "inverse" of a complex number.
>>>
>>> Isn't that just 1.0 / X?
>>
<snip>
>
> Given a=1.0, the inverse function is definitely
> higher performance. It's simpler calculation
> probably implies more accuracy as well.
>
> Warren

Ada has several overloaded versions of "/".
GNAT seems to implement "/" (Real, Complex) using a simplified version
of "/" (Complex, Complex).
It looks (almost) like the example you gave.

function "/" (Left : Real'Base; Right : Complex) return Complex is
a : constant R := Left;
c : constant R := Right.Re;
d : constant R := Right.Im;
begin
return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
Im => -((a * d) / (c ** 2 + d ** 2)));
end "/";

There are differences that may have an impact on speed.
But the above version is simpler than this one:

function "/" (Left, Right : Complex) return Complex is
a : constant R := Left.Re;
b : constant R := Left.Im;
c : constant R := Right.Re;
d : constant R := Right.Im;

begin
if c = 0.0 and then d = 0.0 then
raise Constraint_Error;
else
return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
end if;
end "/";

Damien
From: Adam Beneschan on
On Jun 29, 2:00 pm, Warren <ve3...(a)gmail.com> wrote:
> Warren expounded innews:Xns9DA6AC2A2C8BWarrensBlatherings(a)81.169.183.62:
>
>
>
>
>
>
>
> > Adam Beneschan expounded in news:73a0af1d-7213-44af-90fa-ed6de4c64ce8
> > @b4g2000pra.googlegroups.com:
>
> >> On Jun 29, 12:31 pm, Warren <ve3...(a)gmail.com> wrote:
>
> >>> > So, what is the "missing" function?
>
> >>> > <http://www.adaic.com/standards/05rm/html/RM-G-1-2.html>
>
> >>> 1) The "inverse" of a complex number.
>
> >> Isn't that just 1.0 / X?  
>
> > Ok, it seems to be, as the GSL (Gnu Scientific Library)
> > defines it as:
>
> > gsl_complex
> > gsl_complex_inverse (gsl_complex a)
> > {                               /* z=1/a */
> >   double s = 1.0 / gsl_complex_abs (a);
>
> >   gsl_complex z;
> >   GSL_SET_COMPLEX (&z, (GSL_REAL (a) * s) * s, -(GSL_IMAG (a) * s) *
> s);
> >   return z;
> > }
>
> Apologies for following up my own post, but
> looking at the GSL implementation of complex
> division:
>
> gsl_complex
> gsl_complex_div (gsl_complex a, gsl_complex b)
> {                               /* z=a/b */
>   double ar = GSL_REAL (a), ai = GSL_IMAG (a);
>   double br = GSL_REAL (b), bi = GSL_IMAG (b);
>
>   double s = 1.0 / gsl_complex_abs (b);
>
>   double sbr = s * br;
>   double sbi = s * bi;
>
>   double zr = (ar * sbr + ai * sbi) * s;
>   double zi = (ai * sbr - ar * sbi) * s;
>
>   gsl_complex z;
>   GSL_SET_COMPLEX (&z, zr, zi);
>   return z;
>
> }
>
> Given a=1.0, the inverse function is definitely
> higher performance.  It's simpler calculation
> probably implies more accuracy as well.

That last makes no sense to me. The code to compute A / Z for real A
has got to be essentially the same as computing 1.0 / Z and then
multiplying the real and imaginary parts of the result by A---how else
would it be done? So while you might gain slightly in efficiency by
not performing an extra step, I don't see how you can gain in
accuracy---I've never heard of accuracy ever getting lost by
multiplying a floating-point number by 1.0. Not according to
everything I know about how floating-point arithmetic works.

Even if you implement A / Z by converting A to a complex number A+0i
and then performing a complex division, I still don't see any way to
do this except by computing 1.0/Z and performing a complex
multiplication; again, while there will be some wasted processor
cycles, you're still going to be multiplying numbers by 1.0 and then
adding 0.0, which again is not going to result in any lost accuracy.

-- Adam

From: Simon Wright on
Adam Beneschan <adam(a)irvine.com> writes:

> By the way, I'm not sure that the error is the main reason for having
> this function. The result of the two-argument arctan function is in
> the range -pi .. +pi, as opposed to -pi/2 .. +pi/2 for the one-
> argument variation; this allows you to get a result in any of the four
> quadrants (instead of just two), giving you the correct result when
> you want to compute the angle corresponding to a point (x, y) on a
> Cartesian graph.

And avoiding problems with the extremes; arctan (y / x) is a Bad Idea!
From: Adam Beneschan on
On Jun 29, 10:01 pm, Simon Wright <si...(a)pushface.org> wrote:
> Adam Beneschan <a...(a)irvine.com> writes:
> > By the way, I'm not sure that the error is the main reason for having
> > this function.  The result of the two-argument arctan function is in
> > the range -pi .. +pi, as opposed to -pi/2 .. +pi/2 for the one-
> > argument variation; this allows you to get a result in any of the four
> > quadrants (instead of just two), giving you the correct result when
> > you want to compute the angle corresponding to a point (x, y) on a
> > Cartesian graph.
>
> And avoiding problems with the extremes; arctan (y / x) is a Bad Idea!

I assume you mean because x could be 0.0 (or close to it)? Yes,
that's another case where the two-argument arctan will give you the
result you probably want (i.e. +pi/2 or -pi/2).

-- Adam