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
|
/*
* This file is part of tela the Tensor Language.
* Copyright (c) 1994-1999 Pekka Janhunen
*/
#ifdef __GNUC__
# pragma implementation "Tcomplex.H"
#endif
#include "def.H"
#include "Tcomplex.H"
Tcomplex& Tcomplex::operator/=(const Tcomplex& y)
{
// Treat special cases
if (y.imag() == 0.0) {
re /= y.real();
im /= y.real();
return *this;
}
if (y.real() == 0.0) {
re /= y.imag();
im /= -y.imag();
return *this;
}
// Generic division algorithm
if (fabs(y.real()) >= fabs(y.imag())) {
const Treal r = y.imag()/y.real();
const Treal d = y.real() + r*y.imag();
const Treal re_tmp = (re + im*r)/d;
im = (im - re*r)/d;
re = re_tmp;
} else {
const Treal r = y.real()/y.imag();
const Treal d = y.imag() + r*y.real();
const Treal re_tmp = (re*r + im)/d;
im = (im*r - re)/d;
re = re_tmp;
}
return *this;
}
Tcomplex pow(const Tcomplex& x, const Tcomplex& p) {
const Treal h = hypot(x.real(), x.imag());
const Treal a = atan2(x.imag(), x.real());
Treal lr = pow(h, p.real());
Treal li = p.real()*a;
if (p.imag() != 0.0) {
lr/= exp(p.imag() * a);
li+= p.imag() * log(h);
}
return Tcomplex(lr*cos(li), lr*sin(li));
}
Tcomplex pow(const Tcomplex& x, int p)
{
if (p == 0)
return Tcomplex(1.0, 0.0);
else if (x == 0.0)
return Tcomplex(0.0, 0.0);
else if (abs(p) > 20)
return pow(x,Treal(p));
else {
Tcomplex res(1.0, 0.0);
Tcomplex b = x;
if (p < 0) {
p = -p;
b = 1.0 / b;
}
for(;;) {
if (p & 1)
res *= b;
if ((p >>= 1) == 0)
return res;
else
b *= b;
}
return res;
}
}
Tcomplex sqrt(const Tcomplex& x) {
if (x.real() == 0.0 && x.imag() == 0.0)
return Tcomplex(0.0, 0.0);
else {
const Treal s = sqrt((fabs(x.real()) + hypot(x.real(), x.imag()))*0.5);
const Treal d = (x.imag()/s) * 0.5;
if (x.real() > 0.0)
return Tcomplex(s, d);
else if (x.imag() >= 0.0)
return Tcomplex(d, s);
else
return Tcomplex(-d, -s);
}
}
Tcomplex asin(const Tcomplex& z) {
const Tcomplex I(0,1);
return -I*log(sqrt(1-z*z)+I*z);
}
Tcomplex acos(const Tcomplex& z) {
const Tcomplex I(0,1);
return -I*log(z+I*sqrt(1-z*z));
}
Tcomplex atan(const Tcomplex& z) {
const Tcomplex I(0,1);
return 0.5*I*log((1-I*z)/(1+I*z));
}
|
