From: John B. Matthews on 29 Jun 2010 17:47 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 29 Jun 2010 17:52 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 29 Jun 2010 18:22 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 30 Jun 2010 01:01 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 30 Jun 2010 10:29
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 |