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
|
/*
* This file is part of tela the Tensor Language.
* Copyright (c) 1994-1996 Pekka Janhunen
*/
#ifdef __GNUC__
# pragma implementation "Tcomplex.H"
#endif
#include "def.H"
#include "Tcomplex.H"
Tcomplex& Tcomplex::operator/=(const Tcomplex& y)
{
Treal invden = 1.0/(fabs(y.real()) + fabs(y.imag()));
Treal xrden = re*invden;
Treal xiden = im*invden;
Treal yrden = y.real()*invden;
Treal yiden = y.imag()*invden;
Treal invnrm = 1.0/(yrden*yrden + yiden*yiden);
re = (xrden*yrden + xiden*yiden)*invnrm;
im = (xiden*yrden - xrden*yiden)*invnrm;
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));
}
|
