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
30extern "C++" {
31template <class FLOAT> complex<FLOAT>
32cos (const complex<FLOAT>& x)
33{
34  return complex<FLOAT> (cos (real (x)) * cosh (imag (x)),
35			   - sin (real (x)) * sinh (imag (x)));
36}
37
38template <class FLOAT> complex<FLOAT>
39cosh (const complex<FLOAT>& x)
40{
41  return complex<FLOAT> (cosh (real (x)) * cos (imag (x)),
42			   sinh (real (x)) * sin (imag (x)));
43}
44
45template <class FLOAT> complex<FLOAT>
46exp (const complex<FLOAT>& x)
47{
48  return polar (FLOAT (exp (real (x))), imag (x));
49}
50
51template <class FLOAT> complex<FLOAT>
52log (const complex<FLOAT>& x)
53{
54  return complex<FLOAT> (log (abs (x)), arg (x));
55}
56
57template <class FLOAT> complex<FLOAT>
58pow (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
67template <class FLOAT> complex<FLOAT>
68pow (const complex<FLOAT>& x, FLOAT y)
69{
70  return exp (FLOAT (y) * log (x));
71}
72
73template <class FLOAT> complex<FLOAT>
74pow (FLOAT x, const complex<FLOAT>& y)
75{
76  return exp (y * FLOAT (log (x)));
77}
78
79template <class FLOAT> complex<FLOAT>
80sin (const complex<FLOAT>& x)
81{
82  return complex<FLOAT> (sin (real (x)) * cosh (imag (x)),
83			   cos (real (x)) * sinh (imag (x)));
84}
85
86template <class FLOAT> complex<FLOAT>
87sinh (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
95template <class FLOAT> istream&
96operator >> (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
123template <class FLOAT> ostream&
124operator << (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/****************************************************************
133Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
134
135Permission to use, copy, modify, and distribute this software
136and its documentation for any purpose and without fee is hereby
137granted, provided that the above copyright notice appear in all
138copies and that both that the copyright notice and this
139permission notice and warranty disclaimer appear in supporting
140documentation, and that the names of AT&T Bell Laboratories or
141Bellcore or any of their entities not be used in advertising or
142publicity pertaining to distribution of the software without
143specific, written prior permission.
144
145AT&T and Bellcore disclaim all warranties with regard to this
146software, including all implied warranties of merchantability
147and fitness.  In no event shall AT&T or Bellcore be liable for
148any special, indirect or consequential damages or any damages
149whatsoever resulting from loss of use, data or profits, whether
150in an action of contract, negligence or other tortious action,
151arising out of or in connection with the use or performance of
152this software.
153****************************************************************/
154
155template <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
181template <class FLOAT> complex<FLOAT>
182operator / (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
205template <class FLOAT> complex<FLOAT>
206operator / (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
229template <class FLOAT> complex<FLOAT>
230pow (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
252template <class FLOAT> complex<FLOAT>
253sqrt (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++"
274