contrib/libstdc++/std/complext.cc

   1 // Member templates for the -*- C++ -*- complex number classes.
   2 // Copyright (C) 1994 Free Software Foundation
   3
   4 // This file is part of the GNU ANSI C++ Library.  This library is free
   5 // software; you can redistribute it and/or modify it under the
   6 // terms of the GNU General Public License as published by the
   7 // Free Software Foundation; either version 2, or (at your option)
   8 // any later version.
   9
  10 // This library is distributed in the hope that it will be useful,
  11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 // GNU General Public License for more details.
  14
  15 // You should have received a copy of the GNU General Public License
  16 // along with this library; see the file COPYING.  If not, write to the Free
  17 // Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  18
  19 // As a special exception, if you link this library with files
  20 // compiled with a GNU compiler to produce an executable, this does not cause
  21 // the resulting executable to be covered by the GNU General Public License.
  22 // This exception does not however invalidate any other reasons why
  23 // the executable file might be covered by the GNU General Public License.
  24
  25 // Written by Jason Merrill based upon the specification in the 27 May 1994
  26 // C++ working paper, ANSI document X3J16/94-0098.
  27
  28 #include <complex>
  29
  30 extern "C++" {
  31 template <class FLOAT> complex<FLOAT>
  32 cos (const complex<FLOAT>& x)
  33 {
  34   return complex<FLOAT> (cos (real (x)) * cosh (imag (x)),
  35                            - sin (real (x)) * sinh (imag (x)));
  36 }
  37
  38 template <class FLOAT> complex<FLOAT>
  39 cosh (const complex<FLOAT>& x)
  40 {
  41   return complex<FLOAT> (cosh (real (x)) * cos (imag (x)),
  42                            sinh (real (x)) * sin (imag (x)));
  43 }
  44
  45 template <class FLOAT> complex<FLOAT>
  46 exp (const complex<FLOAT>& x)
  47 {
  48   return polar (FLOAT (exp (real (x))), imag (x));
  49 }
  50
  51 template <class FLOAT> complex<FLOAT>
  52 log (const complex<FLOAT>& x)
  53 {
  54   return complex<FLOAT> (log (abs (x)), arg (x));
  55 }
  56
  57 template <class FLOAT> complex<FLOAT>
  58 pow (const complex<FLOAT>& x, const complex<FLOAT>& y)
  59 {
  60   FLOAT logr = log (abs (x));
  61   FLOAT t = arg (x);
  62
  63   return polar (FLOAT (exp (logr * real (y) - imag (y) * t)),
  64                 FLOAT (imag (y) * logr + real (y) * t));
  65 }
  66
  67 template <class FLOAT> complex<FLOAT>
  68 pow (const complex<FLOAT>& x, FLOAT y)
  69 {
  70   return exp (FLOAT (y) * log (x));
  71 }
  72
  73 template <class FLOAT> complex<FLOAT>
  74 pow (FLOAT x, const complex<FLOAT>& y)
  75 {
  76   return exp (y * FLOAT (log (x)));
  77 }
  78
  79 template <class FLOAT> complex<FLOAT>
  80 sin (const complex<FLOAT>& x)
  81 {
  82   return complex<FLOAT> (sin (real (x)) * cosh (imag (x)),
  83                            cos (real (x)) * sinh (imag (x)));
  84 }
  85
  86 template <class FLOAT> complex<FLOAT>
  87 sinh (const complex<FLOAT>& x)
  88 {
  89   return complex<FLOAT> (sinh (real (x)) * cos (imag (x)),
  90                            cosh (real (x)) * sin (imag (x)));
  91 }
  92
  93 #include <iostream.h>
  94
  95 template <class FLOAT> istream&
  96 operator >> (istream& is, complex<FLOAT>& x)
  97 {
  98   FLOAT re, im = 0;
  99   char ch = 0;
 100
 101   if (is.ipfx0 ())
 102     {
 103       if (is.peek () == '(')
 104         is >> ch;
 105       is >> re;
 106       if (ch == '(')
 107         {
 108           is >> ch;
 109           if (ch == ',')
 110             is >> im >> ch;
 111         }
 112     }
 113   is.isfx ();
 114
 115   if (ch != 0 && ch != ')')
 116     is.setstate (ios::failbit);
 117   else if (is.good ())
 118     x = complex<FLOAT> (re, im);
 119
 120   return is;
 121 }
 122
 123 template <class FLOAT> ostream&
 124 operator << (ostream& os, const complex<FLOAT>& x)
 125 {
 126   return os << '(' << real (x) << ',' << imag (x) << ')';
 127 }
 128
 129 // The code below is adapted from f2c's libF77, and is subject to this
 130 // copyright:
 131
 132 /****************************************************************
 133 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
 134
 135 Permission to use, copy, modify, and distribute this software
 136 and its documentation for any purpose and without fee is hereby
 137 granted, provided that the above copyright notice appear in all
 138 copies and that both that the copyright notice and this
 139 permission notice and warranty disclaimer appear in supporting
 140 documentation, and that the names of AT&T Bell Laboratories or
 141 Bellcore or any of their entities not be used in advertising or
 142 publicity pertaining to distribution of the software without
 143 specific, written prior permission.
 144
 145 AT&T and Bellcore disclaim all warranties with regard to this
 146 software, including all implied warranties of merchantability
 147 and fitness.  In no event shall AT&T or Bellcore be liable for
 148 any special, indirect or consequential damages or any damages
 149 whatsoever resulting from loss of use, data or profits, whether
 150 in an action of contract, negligence or other tortious action,
 151 arising out of or in connection with the use or performance of
 152 this software.
 153 ****************************************************************/
 154
 155 template <class FLOAT> complex<FLOAT>&
 156 __doadv (complex<FLOAT>* ths, const complex<FLOAT>& y)
 157 {
 158   FLOAT ar = abs (y.re);
 159   FLOAT ai = abs (y.im);
 160   FLOAT nr, ni;
 161   FLOAT t, d;
 162   if (ar <= ai)
 163     {
 164       t = y.re / y.im;
 165       d = y.im * (1 + t*t);
 166       nr = (ths->re * t + ths->im) / d;
 167       ni = (ths->im * t - ths->re) / d;
 168     }
 169   else
 170     {
 171       t = y.im / y.re;
 172       d = y.re * (1 + t*t);
 173       nr = (ths->re + ths->im * t) / d;
 174       ni = (ths->im - ths->re * t) / d;
 175     }
 176   ths->re = nr;
 177   ths->im = ni;
 178   return *ths;
 179 }
 180
 181 template <class FLOAT> complex<FLOAT>
 182 operator / (const complex<FLOAT>& x, const complex<FLOAT>& y)
 183 {
 184   FLOAT ar = abs (real (y));
 185   FLOAT ai = abs (imag (y));
 186   FLOAT nr, ni;
 187   FLOAT t, d;
 188   if (ar <= ai)
 189     {
 190       t = real (y) / imag (y);
 191       d = imag (y) * (1 + t*t);
 192       nr = (real (x) * t + imag (x)) / d;
 193       ni = (imag (x) * t - real (x)) / d;
 194     }
 195   else
 196     {
 197       t = imag (y) / real (y);
 198       d = real (y) * (1 + t*t);
 199       nr = (real (x) + imag (x) * t) / d;
 200       ni = (imag (x) - real (x) * t) / d;
 201     }
 202   return complex<FLOAT> (nr, ni);
 203 }
 204
 205 template <class FLOAT> complex<FLOAT>
 206 operator / (FLOAT x, const complex<FLOAT>& y)
 207 {
 208   FLOAT ar = abs (real (y));
 209   FLOAT ai = abs (imag (y));
 210   FLOAT nr, ni;
 211   FLOAT t, d;
 212   if (ar <= ai)
 213     {
 214       t = real (y) / imag (y);
 215       d = imag (y) * (1 + t*t);
 216       nr = x * t / d;
 217       ni = -x / d;
 218     }
 219   else
 220     {
 221       t = imag (y) / real (y);
 222       d = real (y) * (1 + t*t);
 223       nr = x / d;
 224       ni = -x * t / d;
 225     }
 226   return complex<FLOAT> (nr, ni);
 227 }
 228
 229 template <class FLOAT> complex<FLOAT>
 230 pow (const complex<FLOAT>& xin, int y)
 231 {
 232   if (y == 0)
 233     return complex<FLOAT> (1.0);
 234   complex<FLOAT> r (1.0);
 235   complex<FLOAT> x (xin);
 236   if (y < 0)
 237     {
 238       y = -y;
 239       x = FLOAT(1)/x;
 240     }
 241   for (;;)
 242     {
 243       if (y & 1)
 244         r *= x;
 245       if (y >>= 1)
 246         x *= x;
 247       else
 248         return r;
 249     }
 250 }
 251
 252 template <class FLOAT> complex<FLOAT>
 253 sqrt (const complex<FLOAT>& x)
 254 {
 255   FLOAT r = abs (x);
 256   FLOAT nr, ni;
 257   if (r == 0.0)
 258     nr = ni = r;
 259   else if (real (x) > 0)
 260     {
 261       nr = sqrt (0.5 * (r + real (x)));
 262       ni = imag (x) / nr / 2;
 263     }
 264   else
 265     {
 266       ni = sqrt (0.5 * (r - real (x)));
 267       if (imag (x) < 0)
 268         ni = - ni;
 269       nr = imag (x) / ni / 2;
 270     }
 271   return complex<FLOAT> (nr, ni);
 272 }
 273 } // extern "C++"