[PATCH, libstdc++, complex] complex multiplication algorithm improved

[Wang Feng]

Hi,
   Complex number multiplication done with 4 multiplications, but
this one only need 3:

            (a+ib)(c+id) = ac - bd + i( bc + ad ) -- 4 multiplies
                               = ac - bd + i[ ( a + b )( c + d ) - ac
- bd ] -- 3 multiplies





diff -ru gcc/libstdc++-v3/include/std/complex
gcc.working.copy/libstdc++-v3//include/std/complex
--- gcc/libstdc++-v3/include/std/complex        2010-08-13
16:48:07.000000000 +0800
+++ gcc.working.copy/libstdc++-v3//include/std/complex  2010-09-11
13:21:00.000000000 +0800
@@ -291,9 +291,13 @@
    complex<_Tp>&
    complex<_Tp>::operator*=(const complex<_Up>& __z)
    {
-      const _Tp __r = _M_real * __z.real() - _M_imag * __z.imag();
-      _M_imag = _M_real * __z.imag() + _M_imag * __z.real();
-      _M_real = __r;
+      //(a+ib)(c+id) = ac - bd + i( bc + ad ) -- 4 multiplies
+      //             = ac - bd + i[ ( a + b )( c + d ) - ac - bd ] --
3 multiplies
+      const _Tp ac = _M_real * __z.real();
+      const _Tp bd = _M_imag * __z.imag();
+      const _Tp ab_cd = (_M_real+_M_imag) * (__z.real()+__z.imag());
+      _M_real = ac - bd;
+      _M_imag = ab_cd - ac - bd;
      return *this;
    }

@@ -304,10 +308,11 @@
    complex<_Tp>&
    complex<_Tp>::operator/=(const complex<_Up>& __z)
    {
-      const _Tp __r =  _M_real * __z.real() + _M_imag * __z.imag();
-      const _Tp __n = std::norm(__z);
-      _M_imag = (_M_imag * __z.real() - _M_real * __z.imag()) / __n;
-      _M_real = __r / __n;
+      // a + ib    ( a + ib ) ( c - id )
+      //-------- = ---------------------  -- 6 multiplies
+      // c + id          c^2 + d^2
+      *this *= conj(__z);
+      *this /= norm(__z);
      return *this;
    }
-------------- next part --------------
A non-text attachment was scrubbed...
Name: complex.patch
Type: application/octet-stream
Size: 1448 bytes
Desc: not available
URL: <http://gcc.gnu.org/pipermail/libstdc++/attachments/20100915/c6d65c3f/attachment.obj>