1
2 /*
3 * Complex - Implements a complex number
4 * Copyright (c) 2003-2006 by Mattias Hultgren <mattias_hultgren@tele2.se>
5 *
6 * See complex.h
7 */
8
9 /*
10 News
11 ----
12
13 v3 2006-07-30
14 --
15
16 Replaced all throw_... calls with THROW_ERROR
17 Replaced all std::string usage with utf8_string
18
19 v2 2006-05-07 & 2006-05-21
20 --
21
22 Fixed so that abs returns exact results when the argument is exact with imaginary set to zero
23 Small changes in operator<= and operator>=
24 Small clean up in Complex::append_to_string
25
26 v1 2005-09-25 - 2005-10-10
27 --
28
29 Extracted class Complex from the math2.* system
30 ump_float -> floatx
31 operator < & > doesn't allow imaginary part to be non-zero anymore
32
33 */
34
35 #include "complex.h"
36 #include "keyfile.h"
37
38
append_to_string(utf8_string & str,const Format & fmt) const39 void Complex::append_to_string( utf8_string &str, const Format &fmt ) const throw(error_obj)
40 {
41 try
42 {
43 if( floatx(imaginary) == 0.0 )
44 real.append_to_string( str, fmt );
45 else
46 {
47 if( floatx(real) != 0.0 )
48 {
49 real.append_to_string( str, fmt );
50
51 if( floatx(imaginary) < 0.0 )
52 {
53 if( floatx(imaginary) == -1.0 )
54 str.append( " - i" );
55 else
56 {
57
58 if( !imaginary.isInteger() && imaginary.isExact() )
59 str.append( " - (" );
60 else
61 str.append( " - " );
62 (-imaginary).append_to_string( str, fmt );
63 if( !imaginary.isInteger() && imaginary.isExact() )
64 str.append( ")i" );
65 else
66 str.append( "i" );
67 }
68 }
69 else
70 {
71 if( floatx(imaginary) == 1.0 )
72 str.append( " + i" );
73 else
74 {
75
76 if( !imaginary.isInteger() && imaginary.isExact() )
77 str.append( " + (" );
78 else
79 str.append( " + " );
80 imaginary.append_to_string( str, fmt );
81 if( !imaginary.isInteger() && imaginary.isExact() )
82 str.append( ")i" );
83 else
84 str.append( "i" );
85 }
86 }
87 }
88 else
89 {
90 if( floatx(imaginary) < 0.0 )
91 {
92 if( floatx(imaginary) == -1.0 )
93 str.append( "-i" );
94 else
95 {
96 if( !imaginary.isInteger() && imaginary.isExact() )
97 str.append( "-(" );
98 else
99 str.append( "-" );
100 (-imaginary).append_to_string( str, fmt );
101 if( !imaginary.isInteger() && imaginary.isExact() )
102 str.append( ")i" );
103 else
104 str.append( "i" );
105 }
106 }
107 else
108 {
109 if( floatx(imaginary) == 1.0 )
110 str.append( "i" );
111 else
112 {
113 if( !imaginary.isInteger() && imaginary.isExact() )
114 str.append( "(" );
115 imaginary.append_to_string( str, fmt );
116 if( !imaginary.isInteger() && imaginary.isExact() )
117 str.append( ")i" );
118 else
119 str.append( "i" );
120 }
121 }
122 }
123 }
124 }
125 catch(error_obj error) { throw error; }
126 catch(...) { THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") ); }
127 }
128
operator *(const Complex & var) const129 Complex Complex::operator*(const Complex &var) const
130 {
131 if( floatx(imaginary) == 0.0 && floatx(var.imaginary) == 0.0 )
132 return Complex( real*var.real );
133 else
134 {
135 Complex res( real, real );
136 Real tmp( imaginary );
137 res.real.mul( var.real );
138 tmp.mul( var.imaginary );
139 res.real.sub( tmp );
140
141 res.imaginary.mul( var.imaginary );
142 tmp = imaginary;
143 tmp.mul( var.real );
144 res.imaginary.add( tmp );
145
146 return res;
147 }
148 }
149
operator /(const Complex & var) const150 Complex Complex::operator/(const Complex &var) const throw(error_obj)
151 {
152 if( floatx(var.real) == 0.0 && floatx(var.imaginary) == 0.0 )
153 THROW_ERROR( ErrorType_Divide_by_zero, _("Divide by zero") );
154
155 if( floatx(var.imaginary) == 0.0 )
156 return Complex( real/var.real, imaginary/var.real);
157 else
158 {
159 Real tmp = var.real*var.real + var.imaginary*var.imaginary;
160 return Complex( (real*var.real + imaginary*var.imaginary) / tmp ,
161 (imaginary*var.real - real*var.imaginary) / tmp );
162 }
163 }
164
operator /(const Real & val,const Complex & com)165 Complex operator/(const Real &val,const Complex &com)
166 {
167 if( floatx(com.real) == 0.0 && floatx(com.imaginary) == 0.0 )
168 THROW_ERROR( ErrorType_Divide_by_zero, _("Divide by zero") );
169
170 if(floatx(com.imaginary) == 0.0)
171 return Complex( val/com.real );
172 else
173 return Complex(val*com.real / (com.real*com.real + com.imaginary*com.imaginary),
174 -val*com.imaginary/ (com.real*com.real + com.imaginary*com.imaginary));
175 }
176
operator <(const Complex & var) const177 bool Complex::operator<(const Complex &var) const throw(error_obj)
178 {
179 if( floatx(imaginary) == 0.0 && floatx(var.imaginary) == 0.0 )
180 return (real < var.real);
181 else
182 THROW_ERROR( ErrorType_Domain, _("Domain error: Imaginary part must be zero.") );
183 }
184
operator >(const Complex & var) const185 bool Complex::operator>(const Complex &var) const throw(error_obj)
186 {
187 if( floatx(imaginary) == 0.0 && floatx(var.imaginary) == 0.0 )
188 return (real > var.real);
189 else
190 THROW_ERROR( ErrorType_Domain, _("Domain error: Imaginary part must be zero.") );
191 }
192
sqroot(const Complex & var)193 Complex sqroot(const Complex &var) throw(error_obj)
194 {
195 if( floatx(var.imaginary) != 0.0 )
196 return power( var, Complex( 0.5 ) );
197
198 if( floatx(var.real) >= 0.0 )
199 return Complex( sqrtx(floatx(var.real)) );
200 else
201 return Complex( Real(0), sqrtx(-floatx(var.real)) );
202 }
203
complex_almost_equal(const Complex & left,const Complex & right)204 bool complex_almost_equal(const Complex &left,const Complex &right)
205 {
206 Complex l,r;
207
208 r = right;
209 l = left;
210
211 if( floatx(l.real) <= 1e-10 && floatx(l.real) >= -1e-10 )
212 l.real = Real(1e-10);
213 if( floatx(l.imaginary) <= 1e-10 && floatx(l.imaginary) >= -1e-10 )
214 l.imaginary = Real(1e-10);
215 if( floatx(r.real) <= 1e-10 && floatx(r.real) >= -1e-10 )
216 r.real = Real(1e-10);
217 if( floatx(r.imaginary) <= 1e-10 && floatx(r.imaginary) >= -1e-10 )
218 r.imaginary = Real(1e-10);
219
220 if( (floatx(l.real) * floatx(r.real)) < 0.0 ) // not same sign
221 return false;
222 else if( (floatx(l.imaginary) * floatx(r.imaginary)) < 0.0 )
223 return false;
224 else if( fabsx(log10x(floatx(l.real)) - log10x(floatx(r.real))) > 1e-10 )
225 return false;
226 else if( fabsx(log10x(floatx(l.imaginary)) - log10x(floatx(r.imaginary))) > 1e-10 )
227 return false;
228
229 return true;
230 }
231
power(const Complex & left,const Complex & right)232 Complex power(const Complex &left,const Complex &right) throw(error_obj)
233 {
234 Complex tmp;
235
236 tmp.imaginary = Real(0);
237
238 if( floatx(right.imaginary) != 0.0 )
239 {
240 // let's use the fact that a^b = ( e^(ln a) )^b = e^(b*ln a)
241 tmp = expc( right * logc(left) );
242 }
243 else
244 {
245 if( floatx(left.imaginary) != 0.0 ) // we need to use de Moivres formula
246 {
247 de_moivres:
248 Complex tmp_abs,tmp_arg;
249
250 tmp_abs = abs(left);
251 tmp_arg = arg(left);
252 tmp = powx( floatx(tmp_abs.real), floatx(right.real) ) *
253 (Complex( cosx(floatx(right.real)*floatx(tmp_arg.real)),
254 sinx(floatx(right.real)*floatx(tmp_arg.real)) ) );
255 // tmp = r^n * ( cos(n*phi) + i*sin(n*phi) )
256 }
257 else
258 {
259 if( left.real.isExact() && right.real.isInteger() )
260 {
261 Integer nr;
262 right.real.get_top( &nr );
263 tmp.real = power( left.real, nr );
264 }
265 else if( floatx(left.real) < 0.0 )
266 {
267 if( right.real.isInteger() ) // checks if right.real == N
268 tmp.real = Real( powx(floatx(left.real), floatx(right.real)) );
269 else if((1.0/floatx(right.real)) == (truncx(1.0/floatx(right.real)))) // checks if right.real == 1/N
270 {
271 if(((1.0/floatx(right.real))/4.0) == (truncx((1.0/floatx(right.real))/4.0))) // N is dividable with 4
272 goto de_moivres;
273 else if(((1.0/floatx(right.real))/2.0) == (truncx((1.0/floatx(right.real))/2.0))) // N is dividable with 2
274 {
275 tmp.real = Real(0);
276 tmp.imaginary = Real( powx(-floatx(left.real),floatx(right.real)) );
277 }
278 else
279 tmp.real = Real( -powx(-floatx(left.real),floatx(right.real)) );
280 }
281 else
282 goto de_moivres;
283 }
284 else
285 tmp.real = Real( powx(floatx(left.real),floatx(right.real)) );
286 }
287 }
288 return tmp;
289 }
290
abs(const Complex & var)291 Complex abs(const Complex &var) throw(error_obj)
292 {
293 Real tmp(0);
294
295 if( var.imaginary == tmp )
296 {
297 if( var.real < tmp )
298 return Complex( -var.real );
299 return var;
300 }
301 else
302 return Complex( hypotx( floatx(var.real), floatx(var.imaginary) ) );
303 }
304
305
expc(const Complex & val)306 Complex expc(const Complex &val)
307 {
308 Complex tmp;
309
310 tmp.real = expx( floatx(val.real) );
311 tmp.imaginary = tmp.real * Real( sinx( floatx(val.imaginary) ) );
312 tmp.real = tmp.real * Real( cosx( floatx(val.imaginary) ) );
313
314 return tmp;
315 }
316
logc(const Complex & val)317 Complex logc(const Complex &val)
318 {
319 //ln(z) = ln( abs(z) ) + i*arg(z)
320 return Complex( logx( floatx( abs( val ).real ) ), arg( val ).real );
321 }
322
log2c(const Complex & val)323 Complex log2c(const Complex &val)
324 {
325 //log2(z) = ln(z)/ln(2)
326 return logc( val ) / Complex( logx( 2.0 ) );
327 }
328
log10c(const Complex & val)329 Complex log10c(const Complex &val)
330 {
331 //log10(z) = ln(z)/ln(10)
332 return logc( val ) / Complex( logx( 10.0 ) );
333 }
334
sinc(const Complex & val)335 Complex sinc(const Complex &val)
336 {
337 if( floatx(val.imaginary) == 0.0 )
338 return Complex( sinx( floatx(val.real) ) );
339 else
340 {
341 // sin(z) = ( e^(iz) - e^(-iz) ) / (2i)
342 return ( expc(Complex( 0.0, 1.0 ) * val) - expc( Complex( 0.0, -1.0 ) * val ) ) /
343 Complex( 0.0, 2.0 );
344 }
345 }
346
cosc(const Complex & val)347 Complex cosc(const Complex &val)
348 {
349 if( floatx(val.imaginary) == 0.0 )
350 return Complex( cosx( floatx(val.real) ) );
351 else
352 {
353 // cos(z) = ( e^(iz) + e^(-iz) ) / 2
354 return ( expc(Complex( 0.0, 1.0 ) * val) + expc( Complex( 0.0, -1.0 ) * val ) ) /
355 Complex( 2.0 );
356 }
357 }
358
tanc(const Complex & val)359 Complex tanc(const Complex &val)
360 {
361 if( floatx(val.imaginary) == 0.0 )
362 return Complex( tanx( floatx(val.real) ) );
363 else
364 {
365 // tan(z) = ( e^(2iz) -1 ) / ( i(e^(2iz) +1) )
366 Complex tmp = expc( Complex( 0.0, 2.0 ) * val );
367 return (tmp - Complex( 1.0 )) / ( Complex( 0.0, 1.0 ) * (tmp + Complex( 1.0 )) );
368 }
369 }
370
asinc(const Complex & val)371 Complex asinc(const Complex &val)
372 {
373 floatx tmp = floatx(val.real);
374 if( floatx(val.imaginary) == 0.0 && tmp <= 1.0 && tmp >= -1.0 )
375 return Complex( asinx( tmp ) );
376 else
377 {
378 // asin(z) = -i*ln( iz + sqrt(1-z^2) )
379 return Complex( 0.0, -1.0 ) * logc( Complex( 0.0, 1.0 ) * val +
380 sqroot( Complex( 1.0 ) - val*val ) );
381 }
382 }
383
acosc(const Complex & val)384 Complex acosc(const Complex &val)
385 {
386 floatx tmp = floatx(val.real);
387 if( floatx(val.imaginary) == 0.0 && tmp <= 1.0 && tmp >= -1.0 )
388 return Complex( acosx( tmp ) );
389 else
390 {
391 // acos(z) = pi/2 + i*ln( iz + sqrt(1-z^2) )
392 return Complex( PI/2.0 ) + Complex( 0.0, 1.0 ) * logc( Complex( 0.0, 1.0 ) * val +
393 sqroot( Complex( 1.0 ) - val*val ) );
394 }
395 }
396
atanc(const Complex & val)397 Complex atanc(const Complex &val)
398 {
399 if( floatx(val.imaginary) == 0.0 )
400 return Complex( atanx( floatx(val.real) ) );
401 else
402 {
403 // atan(z) = i/2 * ( ln(1-iz) - ln(1+iz) )
404 return Complex( 0.0, 0.5 ) * ( logc( Complex( 1.0 ) - Complex( 0.0, 1.0 ) * val ) -
405 logc( Complex( 1.0 ) + Complex( 0.0, 1.0 ) * val ) );
406 }
407 }
408
sinhc(const Complex & val)409 Complex sinhc(const Complex &val)
410 {
411 //sinh(z) = (exp(z) - exp(-z))/2
412 return ( expc( val ) - expc( -val ) ) / Complex( 2.0 );
413 }
414
coshc(const Complex & val)415 Complex coshc(const Complex &val)
416 {
417 //cosh(z) = (exp(z) + exp(-z))/2.
418 return ( expc( val ) + expc( -val ) ) / Complex( 2.0 );
419 }
420
tanhc(const Complex & val)421 Complex tanhc(const Complex &val)
422 {
423 //tanh(z) = sinh(z)/cosh(z)
424 return sinhc( val ) / coshc( val );
425 }
426
asinhc(const Complex & val)427 Complex asinhc(const Complex &val)
428 {
429 //asinh(z) = ln(z + sqrt(z^2+1))
430 return logc( val + sqroot( val*val + Complex(1.0) ) );
431 }
432
acoshc(const Complex & val)433 Complex acoshc(const Complex &val)
434 {
435 //acosh(z)=ln(z + sqrt(z^2-1))
436 return logc( val + sqroot( val*val - Complex(1.0) ) );
437 }
438
atanhc(const Complex & val)439 Complex atanhc(const Complex &val)
440 {
441 //atanh(z)=ln( (1+z)/(1-z) ) / 2
442 return logc( (Complex(1.0)+val)/(Complex(1.0)-val) ) / Complex(2.0);
443 }
444