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
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
|
/* See numarray.h for Complex32, Complex64:
typedef struct { Float32 r, i; } Complex32;
typedef struct { Float64 r, i; } Complex64;
*/
typedef struct { Float32 a, theta; } PolarComplex32;
typedef struct { Float64 a, theta; } PolarComplex64;
#define NUM_SQ(x) ((x)*(x))
#define NUM_CABSSQ(p) (NUM_SQ((p).r) + NUM_SQ((p).i))
#define NUM_CABS(p) sqrt(NUM_CABSSQ(p))
#define NUM_C_TO_P(c, p) (p).a = NUM_CABS(c); \
(p).theta = atan2((c).i, (c).r);
#define NUM_P_TO_C(p, c) (c).r = (p).a*cos((p).theta); \
(c).i = (p).a*sin((p).theta);
#define NUM_CASS(p, q) (q).r = (p).r, (q).i = (p).i
#define NUM_CADD(p, q, s) (s).r = (p).r + (q).r, \
(s).i = (p).i + (q).i
#define NUM_CSUB(p, q, s) (s).r = (p).r - (q).r, \
(s).i = (p).i - (q).i
#define NUM_CMUL(p, q, s) \
{ Float64 rp = (p).r; \
Float64 rq = (q).r; \
(s).r = rp*rq - (p).i*(q).i; \
(s).i = rp*(q).i + rq*(p).i; \
}
#define NUM_CDIV(p, q, s) \
{ \
Float64 rp = (p).r; \
Float64 ip = (p).i; \
Float64 rq = (q).r; \
if ((q).i != 0) { \
Float64 temp = NUM_CABSSQ(q); \
(s).r = (rp*rq+(p).i*(q).i)/temp; \
(s).i = (rq*(p).i-(q).i*rp)/temp; \
} else { \
(s).r = rp/rq; \
(s).i = ip/rq; \
} \
}
#define NUM_CREM(p, q, s) \
{ Complex64 r; \
NUM_CDIV(p, q, r); \
r.r = floor(r.r); \
r.i = 0; \
NUM_CMUL(r, q, r); \
NUM_CSUB(p, r, s); \
}
#define NUM_CMINUS(p, s) (s).r = -(p).r; (s).i = -(p).i;
#define NUM_CNEG NUM_CMINUS
#define NUM_CEQ(p, q) (((p).r == (q).r) && ((p).i == (q).i))
#define NUM_CNE(p, q) (((p).r != (q).r) || ((p).i != (q).i))
#define NUM_CLT(p, q) ((p).r < (q).r)
#define NUM_CGT(p, q) ((p).r > (q).r)
#define NUM_CLE(p, q) ((p).r <= (q).r)
#define NUM_CGE(p, q) ((p).r >= (q).r)
/* e**z = e**x * (cos(y)+ i*sin(y)) where z = x + i*y
so e**z = e**x * cos(y) + i * e**x * sin(y)
*/
#define NUM_CEXP(p, s) \
{ Float64 ex = exp((p).r); \
(s).r = ex * cos((p).i); \
(s).i = ex * sin((p).i); \
}
/* e**w = z; w = u + i*v; z = r * e**(i*theta);
e**u * e**(i*v) = r * e**(i*theta);
log(z) = w; log(z) = log(r) + i*theta;
*/
#define NUM_CLOG(p, s) \
{ PolarComplex64 temp; NUM_C_TO_P(p, temp); \
(s).r = num_log(temp.a); \
(s).i = temp.theta; \
}
#define NUM_LOG10_E 0.43429448190325182
#define NUM_CLOG10(p, s) \
{ NUM_CLOG(p, s); \
(s).r *= NUM_LOG10_E; \
(s).i *= NUM_LOG10_E; \
}
/* s = p ** q */
#define NUM_CPOW(p, q, s) { if (NUM_CABSSQ(p) == 0) { \
if ((q).r == 0 && (q).i == 0) { \
(s).r = (s).i = 1; \
} else { \
(s).r = (s).i = 0; \
} \
} else { \
NUM_CLOG(p, s); \
NUM_CMUL(s, q, s); \
NUM_CEXP(s, s); \
} \
}
#define NUM_CSQRT(p, s) { Complex64 temp; temp.r = 0.5; temp.i=0; \
NUM_CPOW(p, temp, s); \
}
#define NUM_CSQR(p, s) { Complex64 temp; temp.r = 2.0; temp.i=0; \
NUM_CPOW(p, temp, s); \
}
#define NUM_CSIN(p, s) { Float64 sp = sin((p).r); \
Float64 cp = cos((p).r); \
(s).r = cosh((p).i) * sp; \
(s).i = sinh((p).i) * cp; \
}
#define NUM_CCOS(p, s) { Float64 sp = sin((p).r); \
Float64 cp = cos((p).r); \
(s).r = cosh((p).i) * cp; \
(s).i = -sinh((p).i) * sp; \
}
#define NUM_CTAN(p, s) { Complex64 ss, cs; \
NUM_CSIN(p, ss); \
NUM_CCOS(p, cs); \
NUM_CDIV(ss, cs, s); \
}
#define NUM_CSINH(p, s) { Float64 sp = sin((p).i); \
Float64 cp = cos((p).i); \
(s).r = sinh((p).r) * cp; \
(s).i = cosh((p).r) * sp; \
}
#define NUM_CCOSH(p, s) { Float64 sp = sin((p).i); \
Float64 cp = cos((p).i); \
(s).r = cosh((p).r) * cp; \
(s).i = sinh((p).r) * sp; \
}
#define NUM_CTANH(p, s) { Complex64 ss, cs; \
NUM_CSINH(p, ss); \
NUM_CCOSH(p, cs); \
NUM_CDIV(ss, cs, s); \
}
#define NUM_CRPOW(p, v, s) { Complex64 cr; cr.r = v; cr.i = 0; \
NUM_CPOW(p,cr,s); \
}
#define NUM_CRMUL(p, v, s) (s).r = (p).r * v; (s).i = (p).i * v;
#define NUM_CIMUL(p, s) { Float64 temp = (s).r; \
(s).r = -(p).i; (s).i = temp; \
}
/* asin(z) = -i * log(i*z + (1 - z**2)**0.5) */
#define NUM_CASIN(p, s) { Complex64 p1; NUM_CASS(p, p1); \
NUM_CIMUL(p, p1); \
NUM_CMUL(p, p, s); \
NUM_CNEG(s, s); \
(s).r += 1; \
NUM_CRPOW(s, 0.5, s); \
NUM_CADD(p1, s, s); \
NUM_CLOG(s, s); \
NUM_CIMUL(s, s); \
NUM_CNEG(s, s); \
}
/* acos(z) = -i * log(z + i*(1 - z**2)**0.5) */
#define NUM_CACOS(p, s) { Complex64 p1; NUM_CASS(p, p1); \
NUM_CMUL(p, p, s); \
NUM_CNEG(s, s); \
(s).r += 1; \
NUM_CRPOW(s, 0.5, s); \
NUM_CIMUL(s, s); \
NUM_CADD(p1, s, s); \
NUM_CLOG(s, s); \
NUM_CIMUL(s, s); \
NUM_CNEG(s, s); \
}
/* atan(z) = i/2 * log( (i+z) / (i - z) ) */
#define NUM_CATAN(p, s) { Complex64 p1, p2; \
NUM_CASS(p, p1); NUM_CNEG(p, p2); \
p1.i += 1; \
p2.i += 1; \
NUM_CDIV(p1, p2, s); \
NUM_CLOG(s, s); \
NUM_CIMUL(s, s); \
NUM_CRMUL(s, 0.5, s); \
}
/* asinh(z) = log( z + (z**2 + 1)**0.5 ) */
#define NUM_CASINH(p, s) { Complex64 p1; NUM_CASS(p, p1); \
NUM_CMUL(p, p, s); \
(s).r += 1; \
NUM_CRPOW(s, 0.5, s); \
NUM_CADD(p1, s, s); \
NUM_CLOG(s, s); \
}
/* acosh(z) = log( z + (z**2 - 1)**0.5 ) */
#define NUM_CACOSH(p, s) { Complex64 p1; NUM_CASS(p, p1); \
NUM_CMUL(p, p, s); \
(s).r -= 1; \
NUM_CRPOW(s, 0.5, s); \
NUM_CADD(p1, s, s); \
NUM_CLOG(s, s); \
}
/* atanh(z) = 1/2 * log( (1+z)/(1-z) ) */
#define NUM_CATANH(p, s) { Complex64 p1, p2; \
NUM_CASS(p, p1); NUM_CNEG(p, p2); \
p1.r += 1; \
p2.r += 1; \
NUM_CDIV(p1, p2, s); \
NUM_CLOG(s, s); \
NUM_CRMUL(s, 0.5, s); \
}
#define NUM_CMIN(p, q) (NUM_CLE(p, q) ? p : q)
#define NUM_CMAX(p, q) (NUM_CGE(p, q) ? p : q)
#define NUM_CNZ(p) (((p).r != 0) || ((p).i != 0))
#define NUM_CLAND(p, q) (NUM_CNZ(p) & NUM_CNZ(q))
#define NUM_CLOR(p, q) (NUM_CNZ(p) | NUM_CNZ(q))
#define NUM_CLXOR(p, q) (NUM_CNZ(p) ^ NUM_CNZ(q))
#define NUM_CLNOT(p) (!NUM_CNZ(p))
#define NUM_CFLOOR(p, s) (s).r = floor((p).r); (s).i = floor((p).i);
#define NUM_CCEIL(p, s) (s).r = ceil((p).r); (s).i = ceil((p).i);
#define NUM_CFABS(p, s) (s).r = fabs((p).r); (s).i = fabs((p).i);
#define NUM_CROUND(p, s) (s).r = num_round((p).r); (s).i = num_round((p).i);
#define NUM_CHYPOT(p, q, s) { Complex64 t; \
NUM_CSQR(p, s); NUM_CSQR(q, t); \
NUM_CADD(s, t, s); \
NUM_CSQRT(s, s); \
}
|
