[prev in list] [next in list] [prev in thread] [next in thread] 

List:       gcc-patches
Subject:    [PATCH, libstdc++, complex] complex multiplication algorithm improved
From:       Wang Feng <wanng.fenng () gmail ! com>
Date:       2010-09-15 8:10:40
Message-ID: AANLkTikfU1+E9K8B8-6XxvxsCLdu01TcznOTURKaE0vF () mail ! gmail ! com
[Download RAW message or body]

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;
     }

["complex.patch" (application/octet-stream)]

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;
     }
     


[prev in list] [next in list] [prev in thread] [next in thread]