[gcc.git] / libstdc++ / std / complext.cc

// Member templates for the -*- C++ -*- complex number classes.
// Copyright (C) 1994 Free Software Foundation

// This file is part of the GNU ANSI C++ Library.  This library is free
// software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with this library; see the file COPYING.  If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

// As a special exception, if you link this library with files
// compiled with a GNU compiler to produce an executable, this does not cause
// the resulting executable to be covered by the GNU General Public License.
// This exception does not however invalidate any other reasons why
// the executable file might be covered by the GNU General Public License.

// Written by Jason Merrill based upon the specification in the 27 May 1994
// C++ working paper, ANSI document X3J16/94-0098.

#include <complex>

extern "C++" {
template <class FLOAT> complex<FLOAT>
cos (const complex<FLOAT>& x)
{
  return complex<FLOAT> (cos (real (x)) * cosh (imag (x)),
			   - sin (real (x)) * sinh (imag (x)));
}

template <class FLOAT> complex<FLOAT>
cosh (const complex<FLOAT>& x)
{
  return complex<FLOAT> (cosh (real (x)) * cos (imag (x)),
			   sinh (real (x)) * sin (imag (x)));
}

template <class FLOAT> complex<FLOAT>
exp (const complex<FLOAT>& x)
{
  return polar (FLOAT (exp (real (x))), imag (x));
}

template <class FLOAT> complex<FLOAT>
log (const complex<FLOAT>& x)
{
  return complex<FLOAT> (log (abs (x)), arg (x));
}

template <class FLOAT> complex<FLOAT>
pow (const complex<FLOAT>& x, const complex<FLOAT>& y)
{
  FLOAT logr = log (abs (x));
  FLOAT t = arg (x);

  return polar (FLOAT (exp (logr * real (y) - imag (y) * t)),
		FLOAT (imag (y) * logr + real (y) * t));
}

template <class FLOAT> complex<FLOAT>
pow (const complex<FLOAT>& x, FLOAT y)
{
  return exp (FLOAT (y) * log (x));
}

template <class FLOAT> complex<FLOAT>
pow (FLOAT x, const complex<FLOAT>& y)
{
  return exp (y * FLOAT (log (x)));
}

template <class FLOAT> complex<FLOAT>
sin (const complex<FLOAT>& x)
{
  return complex<FLOAT> (sin (real (x)) * cosh (imag (x)),
			   cos (real (x)) * sinh (imag (x)));
}

template <class FLOAT> complex<FLOAT>
sinh (const complex<FLOAT>& x)
{
  return complex<FLOAT> (sinh (real (x)) * cos (imag (x)),
			   cosh (real (x)) * sin (imag (x)));
}

#include <iostream.h>

template <class FLOAT> istream&
operator >> (istream& is, complex<FLOAT>& x)
{
  FLOAT re, im = 0;
  char ch = 0;

  if (is.ipfx0 ())
    {
      if (is.peek () == '(')
	is >> ch;
      is >> re;
      if (ch == '(')
	{
	  is >> ch;
	  if (ch == ',')
	    is >> im >> ch;
	}
    }
  is.isfx ();

  if (ch != 0 && ch != ')')
    is.setstate (ios::failbit);
  else if (is.good ())
    x = complex<FLOAT> (re, im);

  return is;
}

template <class FLOAT> ostream&
operator << (ostream& os, const complex<FLOAT>& x)
{
  return os << '(' << real (x) << ',' << imag (x) << ')';
}

// The code below is adapted from f2c's libF77, and is subject to this
// copyright:

/****************************************************************
Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.

Permission to use, copy, modify, and distribute this software
and its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the names of AT&T Bell Laboratories or
Bellcore or any of their entities not be used in advertising or
publicity pertaining to distribution of the software without
specific, written prior permission.

AT&T and Bellcore disclaim all warranties with regard to this
software, including all implied warranties of merchantability
and fitness.  In no event shall AT&T or Bellcore be liable for
any special, indirect or consequential damages or any damages
whatsoever resulting from loss of use, data or profits, whether
in an action of contract, negligence or other tortious action,
arising out of or in connection with the use or performance of
this software.
****************************************************************/

template <class FLOAT> complex<FLOAT>&
__doadv (complex<FLOAT>* ths, const complex<FLOAT>& y)
{
  FLOAT ar = abs (y.re);
  FLOAT ai = abs (y.im);
  FLOAT nr, ni;
  FLOAT t, d;
  if (ar <= ai)
    {
      t = y.re / y.im;
      d = y.im * (1 + t*t);
      nr = (ths->re * t + ths->im) / d;
      ni = (ths->im * t - ths->re) / d;
    }
  else
    {
      t = y.im / y.re;
      d = y.re * (1 + t*t);
      nr = (ths->re + ths->im * t) / d;
      ni = (ths->im - ths->re * t) / d;
    }
  ths->re = nr;
  ths->im = ni;
  return *ths;
}

template <class FLOAT> complex<FLOAT>
operator / (const complex<FLOAT>& x, const complex<FLOAT>& y)
{
  FLOAT ar = abs (real (y));
  FLOAT ai = abs (imag (y));
  FLOAT nr, ni;
  FLOAT t, d;
  if (ar <= ai)
    {
      t = real (y) / imag (y);
      d = imag (y) * (1 + t*t);
      nr = (real (x) * t + imag (x)) / d;
      ni = (imag (x) * t - real (x)) / d;
    }
  else
    {
      t = imag (y) / real (y);
      d = real (y) * (1 + t*t);
      nr = (real (x) + imag (x) * t) / d;
      ni = (imag (x) - real (x) * t) / d;
    }
  return complex<FLOAT> (nr, ni);
}

template <class FLOAT> complex<FLOAT>
operator / (FLOAT x, const complex<FLOAT>& y)
{
  FLOAT ar = abs (real (y));
  FLOAT ai = abs (imag (y));
  FLOAT nr, ni;
  FLOAT t, d;
  if (ar <= ai)
    {
      t = real (y) / imag (y);
      d = imag (y) * (1 + t*t);
      nr = x * t / d;
      ni = -x / d;
    }
  else
    {
      t = imag (y) / real (y);
      d = real (y) * (1 + t*t);
      nr = x / d;
      ni = -x * t / d;
    }
  return complex<FLOAT> (nr, ni);
}

template <class FLOAT> complex<FLOAT>
pow (const complex<FLOAT>& xin, int y)
{
  if (y == 0)
    return complex<FLOAT> (1.0);
  complex<FLOAT> r (1.0);
  complex<FLOAT> x (xin);
  if (y < 0)
    {
      y = -y;
      x = 1/x;
    }
  for (;;)
    {
      if (y & 1)
	r *= x;
      if (y >>= 1)
	x *= x;
      else
	return r;
    }
}

template <class FLOAT> complex<FLOAT>
sqrt (const complex<FLOAT>& x)
{
  FLOAT r = abs (x);
  FLOAT nr, ni;
  if (r == 0.0)
    nr = ni = r;
  else if (real (x) > 0)
    {
      nr = sqrt (0.5 * (r + real (x)));
      ni = imag (x) / nr / 2;
    }
  else
    {
      ni = sqrt (0.5 * (r - real (x)));
      if (imag (x) < 0)
	ni = - ni;
      nr = imag (x) / ni / 2;
    }
  return complex<FLOAT> (nr, ni); 
}
} // extern "C++"
Commit	Line	Data
6599da04 JM	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;
da8c445d	239	x = 1/x;
6599da04 JM	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++"