1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
|
/*
* This file is part of tela the Tensor Language.
* Copyright (c) 1994-1996 Pekka Janhunen
*/
#ifdef __GNUC__
# pragma interface
#endif
// Treal must be defined before including this
#ifndef TCOMPLEX_H
class Tcomplex {
private:
Treal re;
Treal im;
public:
Treal real() const {return re;}
Treal imag() const {return im;}
Tcomplex() {}
Tcomplex(const Tcomplex& y) :re(y.real()), im(y.imag()) {}
Tcomplex(Treal r, Treal i=0) :re(r), im(i) {}
Tcomplex& operator=(const Tcomplex& y) {re = y.real(); im = y.imag(); return *this;}
Tcomplex& operator+=(const Tcomplex& y) {re+= y.real(); im+= y.imag(); return *this;}
Tcomplex& operator+=(Treal y) {re+= y; return *this;}
Tcomplex& operator-=(const Tcomplex& y) {re-=y.real(); im-=y.imag(); return *this;}
Tcomplex& operator-=(Treal y) {re-= y; return *this;}
Tcomplex& operator*=(const Tcomplex& y) {
Treal r = re*y.real() - im*y.imag();
im = re*y.imag() + im*y.real();
re = r;
return *this;
}
Tcomplex& operator*=(Treal y) {re*= y; im*= y; return *this;}
Tcomplex& operator/=(const Tcomplex& y);
Tcomplex& operator/=(Treal y) {re/= y; im/= y; return *this;}
};
inline int operator==(const Tcomplex& x, const Tcomplex& y) {
return x.real() == y.real() && x.imag() == y.imag();
}
inline int operator==(const Tcomplex& x, Treal y) {
return x.imag() == 0.0 && x.real() == y;
}
inline int operator!=(const Tcomplex& x, const Tcomplex& y) {
return x.real() != y.real() || x.imag() != y.imag();
}
inline int operator!=(const Tcomplex& x, Treal y) {
return x.imag() != 0.0 || x.real() != y;
}
inline Tcomplex operator-(const Tcomplex& x) {
return Tcomplex(-x.real(), -x.imag());
}
inline Tcomplex conj(const Tcomplex& x) {
return Tcomplex(x.real(), -x.imag());
}
inline Tcomplex operator+(const Tcomplex& x, const Tcomplex& y) {
return Tcomplex(x.real() + y.real(), x.imag() + y.imag());
}
inline Tcomplex operator+(const Tcomplex& x, Treal y) {
return Tcomplex(x.real() + y, x.imag());
}
inline Tcomplex operator+(Treal x, const Tcomplex& y) {
return Tcomplex(x + y.real(), y.imag());
}
inline Tcomplex operator-(const Tcomplex& x, const Tcomplex& y) {
return Tcomplex(x.real() - y.real(), x.imag() - y.imag());
}
inline Tcomplex operator-(const Tcomplex& x, Treal y) {
return Tcomplex(x.real() - y, x.imag());
}
inline Tcomplex operator-(Treal x, const Tcomplex& y) {
return Tcomplex(x - y.real(), -y.imag());
}
inline Tcomplex operator*(const Tcomplex& x, const Tcomplex& y) {
return Tcomplex(x.real()*y.real() - x.imag()*y.imag(),
x.real()*y.imag() + x.imag()*y.real());
}
inline Tcomplex operator*(const Tcomplex& x, Treal y) {
return Tcomplex(x.real()*y, x.imag()*y);
}
inline Tcomplex operator*(Treal x, const Tcomplex& y) {
return Tcomplex(x*y.real(), x*y.imag());
}
inline Tcomplex operator/(const Tcomplex& x, const Tcomplex& y) {
const Treal invden = 1.0/(fabs(y.real()) + fabs(y.imag()));
const Treal xrden = x.real()*invden;
const Treal xiden = x.imag()*invden;
const Treal yrden = y.real()*invden;
const Treal yiden = y.imag()*invden;
const Treal invnrm = 1.0/(yrden*yrden + yiden*yiden);
return Tcomplex((xrden*yrden + xiden*yiden)*invnrm,
(xiden*yrden - xrden*yiden)*invnrm);
}
#if 0
inline Tcomplex operator/(const Tcomplex& x, const Tcomplex& y) {
const Treal coeff = 1.0/(y.real()*y.real() + y.imag()*y.imag());
return Tcomplex(
(x.real()*y.real() + x.imag()*y.imag())*coeff,
(x.imag()*y.real() - x.real()*y.imag())*coeff );
}
#endif
#if 0
inline TComplex operator/(const Treal x, const Tcomplex& y) {
const Treal invden = 1.0/(fabs(y.real()) + fabs(y.imag()));
const Treal xrden = x*invden;
const Treal yrden = y.real()*invden;
const Treal yiden = y.imag()*invden;
const Treal coeff = xrden/(yrden*yrden + yiden*yiden);
return Tcomplex(yrden*coeff,-yiden*coeff);
}
#endif
inline Tcomplex operator/(Treal x, const Tcomplex& y) {
const Treal coeff = x/(y.real()*y.real() + y.imag()*y.imag());
return Tcomplex(y.real()*coeff,-y.imag()*coeff);
}
inline Tcomplex operator/(const Tcomplex& x, Treal y) {
const Treal invy = 1.0/y;
return Tcomplex(x.real()*invy, x.imag()*invy);
}
inline Treal real(const Tcomplex& x) {return x.real();}
inline Treal imag(const Tcomplex& x) {return x.imag();}
inline Treal abs(const Tcomplex& x) {
return hypot(x.real(), x.imag());
}
inline Treal arg(const Tcomplex& x) {
return atan2(x.imag(), x.real());
}
inline Tcomplex exp(const Tcomplex& x) {
const Treal r = exp(x.real());
return Tcomplex(r*cos(x.imag()),
r*sin(x.imag()));
}
inline Tcomplex log(const Tcomplex& x) {
Treal h = hypot(x.real(), x.imag());
return Tcomplex(log(h), atan2(x.imag(), x.real()));
}
inline Tcomplex cosh(const Tcomplex& x) {
return Tcomplex(cos(x.imag()) * cosh(x.real()),
sin(x.imag()) * sinh(x.real()));
}
inline Tcomplex sinh(const Tcomplex& x) {
return Tcomplex(cos(x.imag()) * sinh(x.real()),
sin(x.imag()) * cosh(x.real()));
}
inline Tcomplex tanh(const Tcomplex& x) {
// tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x))
// = (1 - exp(-2*x))/(1 + exp(-2*x))
const Tcomplex e = exp(-2*x);
return (1-e)/(1+e);
}
inline Tcomplex cos(const Tcomplex& x) {
return Tcomplex(cos(x.real()) * cosh(x.imag()),
-sin(x.real()) * sinh(x.imag()));
}
inline Tcomplex sin(const Tcomplex& x) {
return Tcomplex(sin(x.real()) * cosh(x.imag()),
cos(x.real()) * sinh(x.imag()));
}
inline Tcomplex pow(const Tcomplex& x, Treal p) {
const Treal h = hypot(x.real(), x.imag());
const Treal lr = pow(h, p);
const Treal a = atan2(x.imag(), x.real());
const Treal li = p * a;
return Tcomplex(lr*cos(li), lr*sin(li));
}
// Non-inline functions
Tcomplex pow(const Tcomplex& x, int p);
Tcomplex pow(const Tcomplex& x, const Tcomplex& p);
Tcomplex sqrt(const Tcomplex& x);
Tcomplex asin(const Tcomplex& x);
Tcomplex acos(const Tcomplex& x);
Tcomplex atan(const Tcomplex& x);
#define TCOMPLEX_H
#endif
|
