]>
Commit | Line | Data |
---|---|---|
dd82cadc | 1 | //This file is part of CRAB |
2 | //http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html | |
3 | ||
4 | /*-------------------------------------------------- | |
5 | complex.cxx | |
6 | Version 1.0.1 | |
7 | ||
8 | -- Some modification for speed | |
9 | -- Some protections from overflow and underflow | |
10 | -- Solving of equations of order 2,3,4: | |
11 | Solve2, Solve3, Solve4; | |
12 | --------------------------------------------------*/ | |
13 | ||
14 | #ifndef __COMPLEX_CXX__ | |
15 | #define __COMPLEX_CXX__ | |
16 | ||
17 | #ifndef small_epsilon | |
18 | #define small_epsilon 1e-12 | |
19 | #endif | |
20 | ||
21 | #ifndef PI | |
22 | #define PI 3.141592653589793238462643 | |
23 | #endif | |
24 | ||
f93c53d8 | 25 | |
26 | #include <TMath.h> | |
27 | #include <Riostream.h> | |
28 | //#include <math.h> | |
3294e9ff | 29 | #include <stdlib.h> |
dd82cadc | 30 | |
31 | /////////////////////////// Helper Error Function | |
32 | ||
8e8eae84 | 33 | void ComplexError (const char *msg) |
dd82cadc | 34 | { |
35 | cout << endl << msg << endl; | |
36 | exit (1); | |
37 | } | |
38 | ||
8e8eae84 | 39 | void ComplexWarning (const char *msg) |
dd82cadc | 40 | { |
41 | cout << endl << msg << endl; | |
42 | } | |
43 | ||
44 | class Complex | |
45 | { | |
46 | public: | |
47 | double Re,Im; | |
48 | ||
49 | Complex (double a=0, double b=0): Re(a),Im(b) {} | |
50 | ||
51 | double GetRe () { return Re; } | |
52 | double GetIm () { return Im; } | |
53 | void SetRe (double a) { Re = a; } | |
54 | void SetIm (double a) { Im = a; } | |
55 | void Set (double a, double b) { Re = a; Im = b; } | |
56 | void SetPolar (double, double); | |
57 | void SetAbs (double); | |
58 | void SetArg (double); | |
59 | Complex operator ! () const; | |
60 | Complex operator ! (); | |
61 | Complex operator + () const; | |
62 | Complex operator + (); | |
63 | Complex operator - () const; | |
64 | Complex operator - (); | |
65 | Complex operator += (const Complex &); | |
66 | Complex operator += (double); | |
67 | Complex operator -= (const Complex &); | |
68 | Complex operator -= (double); | |
69 | Complex operator *= (const Complex &); | |
70 | Complex operator *= (double); | |
71 | Complex operator /= (const Complex &); | |
72 | Complex operator /= (double); | |
73 | ||
74 | friend ostream& operator << (ostream&, const Complex &); | |
75 | friend istream& operator >> (istream&, Complex &); | |
76 | }; | |
77 | ||
78 | const Complex ImUnit (0,1); | |
79 | ||
80 | void Zero (Complex &); | |
81 | Complex Polar (double, double); | |
82 | double real (const Complex &); | |
83 | double imag (const Complex &); | |
84 | double Arg (const Complex &); | |
85 | double phase(const Complex &); | |
86 | double abs (const Complex &); | |
87 | void Solve2 (Complex*, const Complex&, const Complex&); | |
88 | void Solve3 (Complex*, const Complex&, const Complex&, const Complex&); | |
89 | void Solve4 (Complex*, const Complex&, const Complex&, const Complex&, const Complex&); | |
90 | Complex Solve2 (const Complex&, const Complex&, int RootNumber); | |
91 | Complex Solve3 (const Complex&, const Complex&, const Complex&, int RootNumber); | |
92 | Complex Solve4 (const Complex&, const Complex&, const Complex&, const Complex&, int RootNumber); | |
93 | Complex conj (const Complex &); | |
94 | Complex sqr (const Complex &); | |
95 | Complex sqrt (const Complex &, int=0); | |
96 | Complex pow (Complex, int); | |
97 | Complex pow (Complex, double); | |
98 | Complex pow (Complex, Complex); | |
99 | Complex root (const Complex &, int, int=0); | |
100 | Complex exp (const Complex &); | |
101 | Complex log (const Complex &); | |
102 | Complex sin (const Complex &); | |
103 | Complex cos (const Complex &); | |
104 | Complex tan (const Complex &); | |
105 | Complex cot (const Complex &); | |
106 | Complex sec (const Complex &); | |
107 | Complex csc (const Complex &); | |
108 | Complex sinh (const Complex &); | |
109 | Complex cosh (const Complex &); | |
110 | Complex tanh (const Complex &); | |
111 | Complex coth (const Complex &); | |
112 | Complex sech (const Complex &); | |
113 | Complex csch (const Complex &); | |
114 | Complex asin (const Complex &, int=0); | |
115 | Complex acos (const Complex &, int=0); | |
116 | Complex atan (const Complex &); | |
117 | Complex acot (const Complex &); | |
118 | Complex asec (const Complex &, int=0); | |
119 | Complex acsc (const Complex &, int=0); | |
120 | Complex asinh(const Complex &, int=0); | |
121 | Complex acosh(const Complex &, int=0); | |
122 | Complex atanh(const Complex &); | |
123 | Complex acoth(const Complex &); | |
124 | Complex asech(const Complex &, int=0); | |
125 | Complex acsch(const Complex &, int=0); | |
126 | ||
127 | Complex operator ^ (const Complex &, int); | |
128 | Complex operator ^ (const Complex &, double); | |
129 | Complex operator ^ (const Complex &, const Complex &); | |
130 | ||
131 | double operator += (double&, const Complex &); | |
132 | double operator -= (double&, const Complex &); | |
133 | double operator *= (double&, const Complex &); | |
134 | double operator /= (double&, const Complex &); | |
135 | ||
136 | Complex operator + (const Complex &, const Complex &); | |
137 | Complex operator + (const Complex &, double); | |
138 | Complex operator + (double, const Complex &); | |
139 | Complex operator - (const Complex &, const Complex &); | |
140 | Complex operator - (const Complex &, double); | |
141 | Complex operator - (double, const Complex &); | |
142 | Complex operator * (const Complex &, const Complex &); | |
143 | Complex operator * (const Complex &, double); | |
144 | Complex operator * (double, const Complex &); | |
145 | Complex operator / (const Complex &, const Complex &); | |
146 | Complex operator / (const Complex &, double); | |
147 | Complex operator / (double, const Complex &); | |
148 | ||
149 | int operator == (const Complex &, double); | |
150 | int operator == (double, const Complex &); | |
151 | int operator == (const Complex &, const Complex &); | |
152 | int operator != (const Complex &, double); | |
153 | int operator != (double, const Complex &); | |
154 | int operator != (const Complex &, const Complex &); | |
155 | int operator >= (const Complex &, double); | |
156 | int operator >= (double, const Complex &); | |
157 | int operator >= (const Complex &, const Complex &); | |
158 | int operator <= (const Complex &, double); | |
159 | int operator <= (double, const Complex &); | |
160 | int operator <= (const Complex &, const Complex &); | |
161 | int operator > (const Complex &, double); | |
162 | int operator > (double, const Complex &); | |
163 | int operator > (const Complex &, const Complex &); | |
164 | int operator < (const Complex &, double); | |
165 | int operator < (double, const Complex &); | |
166 | int operator < (const Complex &, const Complex &); | |
167 | ||
168 | //////////////////////////////////////////////// | |
169 | ||
dd82cadc | 170 | |
171 | inline Complex Complex::operator + () const | |
172 | { | |
173 | return *this; | |
174 | } | |
175 | ||
176 | inline Complex Complex::operator + () | |
177 | { | |
178 | return *this; | |
179 | } | |
180 | ||
181 | inline Complex Complex::operator - () const | |
182 | { | |
183 | return Complex (-Re, -Im); | |
184 | } | |
185 | ||
186 | inline Complex Complex::operator - () | |
187 | { | |
188 | return Complex (-Re, -Im); | |
189 | } | |
190 | ||
191 | inline Complex Complex::operator ! () const // complex conjugate // | |
192 | { | |
193 | return Complex (Re, -Im); | |
194 | } | |
195 | ||
196 | inline Complex Complex::operator ! () // complex conjugate // | |
197 | { | |
198 | return Complex (Re, -Im); | |
199 | } | |
200 | ||
201 | /////////////////////// +=,-=,*=,/= | |
202 | ||
203 | double operator += (double &a, const Complex &b) | |
204 | { | |
3294e9ff | 205 | if (TMath::Abs(b.Im) < small_epsilon) a+=b.Re; |
dd82cadc | 206 | else |
207 | ComplexError ("Error in double+=Complex: Complex is not double number"); | |
208 | return a; | |
209 | } | |
210 | ||
211 | double operator -= (double &a, const Complex &b) | |
212 | { | |
3294e9ff | 213 | if (TMath::Abs(b.Im) < small_epsilon) a-=b.Re; |
dd82cadc | 214 | else |
215 | ComplexError ("Error in double -=Complex: Complex is not double number"); | |
216 | return a; | |
217 | } | |
218 | ||
219 | double operator *= (double &a, const Complex &b) | |
220 | { | |
3294e9ff | 221 | if (TMath::Abs(b.Im) < small_epsilon) a*=b.Re; |
dd82cadc | 222 | else |
223 | ComplexError ("Error in double *=Complex: Complex is not double number"); | |
224 | return a; | |
225 | } | |
226 | ||
227 | double operator /= (double &a, const Complex &b) | |
228 | { | |
3294e9ff | 229 | if (TMath::Abs(b.Im) < small_epsilon) a/=b.Re; |
dd82cadc | 230 | else |
231 | ComplexError ("Error in double /=Complex: Complex is not double number"); | |
232 | return a; | |
233 | } | |
234 | ||
235 | inline Complex Complex::operator += (const Complex &a) | |
236 | { | |
237 | Re+=a.Re; | |
238 | Im+=a.Im; | |
239 | return *this; | |
240 | } | |
241 | ||
242 | inline Complex Complex::operator += (double a) | |
243 | { | |
244 | Re+=a; | |
245 | return *this; | |
246 | } | |
247 | ||
248 | inline Complex Complex::operator -= (const Complex &a) | |
249 | { | |
250 | Re-=a.Re; | |
251 | Im-=a.Im; | |
252 | return *this; | |
253 | } | |
254 | ||
255 | inline Complex Complex::operator -= (double a) | |
256 | { | |
257 | Re-=a; | |
258 | return *this; | |
259 | } | |
260 | ||
261 | Complex Complex::operator *= (const Complex &a) | |
262 | { | |
263 | double t1=Re*a.Re, t2=Im*a.Im; | |
264 | Im = (Re+Im) * (a.Re+a.Im) - t1 - t2; | |
265 | Re = t1 - t2; | |
266 | return *this; | |
267 | } | |
268 | ||
269 | inline Complex Complex::operator *= (double a) | |
270 | { | |
271 | Re*=a; | |
272 | Im*=a; | |
273 | return *this; | |
274 | } | |
275 | ||
276 | Complex Complex::operator /= (const Complex &a) | |
277 | { | |
278 | double t1, t2, temp; | |
3294e9ff | 279 | if (TMath::Abs(a.Re) >= TMath::Abs(a.Im)) |
dd82cadc | 280 | { |
281 | t1 = a.Im / a.Re; | |
282 | t2 = a.Re + a.Im * t1; | |
283 | temp = (Re + Im * t1) / t2; | |
284 | Im = (Im - Re * t1) / t2; | |
285 | Re = temp; | |
286 | } | |
287 | else | |
288 | { | |
289 | t1 = a.Re / a.Im; | |
290 | t2 = a.Re * t1 + a.Im; | |
291 | temp = (Re * t1 + Im) / t2; | |
292 | Im = (Im * t1 - Re) / t2; | |
293 | Re = temp; | |
294 | } | |
295 | return *this; | |
296 | } | |
297 | ||
298 | inline Complex Complex::operator /= (double a) | |
299 | { | |
300 | Re/=a; | |
301 | Im/=a; | |
302 | return *this; | |
303 | } | |
304 | ||
305 | /////////////////////////// +,-,*,/ | |
306 | ||
307 | inline Complex operator + (const Complex &a, const Complex &b) | |
308 | { | |
309 | return Complex (a.Re+b.Re, a.Im+b.Im); | |
310 | } | |
311 | ||
312 | inline Complex operator + (const Complex &a, double b) | |
313 | { | |
314 | return Complex (a.Re+b, a.Im); | |
315 | } | |
316 | ||
317 | inline Complex operator + (double b, const Complex &a) | |
318 | { | |
319 | return Complex (a.Re+b, a.Im); | |
320 | } | |
321 | ||
322 | inline Complex operator - (const Complex &a, const Complex &b) | |
323 | { | |
324 | return Complex (a.Re-b.Re, a.Im-b.Im); | |
325 | } | |
326 | ||
327 | inline Complex operator - (const Complex &a, double b) | |
328 | { | |
329 | return Complex (a.Re-b, a.Im); | |
330 | } | |
331 | ||
332 | inline Complex operator - (double b, const Complex &a) | |
333 | { | |
334 | return Complex (b-a.Re, -a.Im); | |
335 | } | |
336 | ||
337 | Complex operator * (const Complex &a, const Complex &b) | |
338 | { | |
339 | double t1 = a.Re * b.Re; | |
340 | double t2 = a.Im * b.Im; | |
341 | return Complex (t1 - t2, (a.Re+a.Im) * (b.Re+b.Im) - t1 - t2); | |
342 | } | |
343 | ||
344 | inline Complex operator * (const Complex &a, double b) | |
345 | { | |
346 | return Complex (a.Re*b, a.Im*b); | |
347 | } | |
348 | ||
349 | inline Complex operator * (double b, const Complex &a) | |
350 | { | |
351 | return Complex (a.Re*b, a.Im*b); | |
352 | } | |
353 | ||
354 | inline Complex operator / (const Complex &a, double b) | |
355 | { | |
356 | return Complex (a.Re/b, a.Im/b); | |
357 | } | |
358 | ||
359 | inline Complex operator / (const Complex &a, const Complex &b) | |
360 | { | |
361 | double t1, t2; | |
3294e9ff | 362 | if (TMath::Abs(b.Re) >= TMath::Abs(b.Im)) |
dd82cadc | 363 | { |
364 | t1 = b.Im / b.Re; | |
365 | t2 = b.Re + b.Im * t1; | |
366 | return Complex ((a.Re + a.Im * t1) / t2, (a.Im - a.Re * t1) / t2); | |
367 | } | |
368 | else | |
369 | { | |
370 | t1 = b.Re / b.Im; | |
371 | t2 = b.Re * t1 + b.Im; | |
372 | return Complex ((a.Re * t1 + a.Im) / t2, (a.Im * t1 - a.Re) / t2); | |
373 | } | |
374 | } | |
375 | ||
376 | inline Complex operator / (double b, const Complex &a) | |
377 | { | |
378 | return (Complex (b,0)) / a; | |
379 | } | |
380 | ||
381 | //////////////////////// <<,>> | |
382 | ||
383 | ostream& operator << (ostream &stream, const Complex &a) | |
384 | { | |
385 | stream<<" "<<a.Re<<" "<<a.Im<<" "; | |
386 | return stream; | |
387 | } | |
388 | ||
389 | istream& operator >> (istream &stream, Complex &a) | |
390 | { | |
391 | stream>>a.Re>>a.Im; | |
392 | return stream; | |
393 | } | |
394 | ||
395 | /////////////////////// ==,!= | |
396 | ||
397 | inline int operator == (const Complex &a, double b) | |
398 | { | |
3294e9ff | 399 | return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon); |
dd82cadc | 400 | } |
401 | ||
402 | inline int operator == (double b, const Complex &a) | |
403 | { | |
3294e9ff | 404 | return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon); |
dd82cadc | 405 | } |
406 | ||
407 | inline int operator == (const Complex &a, const Complex &b) | |
408 | { | |
409 | return (a.Re==b.Re) && (a.Im==b.Im); | |
410 | } | |
411 | ||
412 | inline int operator != (const Complex &a, double b) | |
413 | { | |
3294e9ff | 414 | return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon); |
dd82cadc | 415 | } |
416 | ||
417 | inline int operator != (double b, const Complex &a) | |
418 | { | |
3294e9ff | 419 | return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon); |
dd82cadc | 420 | } |
421 | ||
422 | inline int operator != (const Complex &a, const Complex &b) | |
423 | { | |
424 | return (a.Re!=b.Re) || (a.Im!=b.Im); | |
425 | } | |
426 | ||
427 | /////////////////////// >=,<= | |
428 | ||
429 | inline int operator >= (const Complex &a, double b) | |
430 | { | |
431 | return abs(a) >= b; | |
432 | } | |
433 | ||
434 | inline int operator >= (double b, const Complex &a) | |
435 | { | |
436 | return b >= abs(a); | |
437 | } | |
438 | ||
439 | inline int operator >= (const Complex &a, const Complex &b) | |
440 | { | |
441 | return abs(a) >= abs(b); | |
442 | } | |
443 | ||
444 | inline int operator <= (const Complex &a, double b) | |
445 | { | |
446 | return abs(a) <= b; | |
447 | } | |
448 | ||
449 | inline int operator <= (double b, const Complex &a) | |
450 | { | |
451 | return b <= abs(a); | |
452 | } | |
453 | ||
454 | inline int operator <= (const Complex &a, const Complex &b) | |
455 | { | |
456 | return abs(a) <= abs(b); | |
457 | } | |
458 | ||
459 | /////////////////////// >,< | |
460 | ||
461 | inline int operator > (const Complex &a, double b) | |
462 | { | |
463 | return abs(a) > b; | |
464 | } | |
465 | ||
466 | inline int operator > (double b, const Complex &a) | |
467 | { | |
468 | return b > abs(a); | |
469 | } | |
470 | ||
471 | inline int operator > (const Complex &a, const Complex &b) | |
472 | { | |
473 | return abs(a) > abs(b); | |
474 | } | |
475 | ||
476 | inline int operator < (const Complex &a, double b) | |
477 | { | |
478 | return abs(a) < b; | |
479 | } | |
480 | ||
481 | inline int operator < (double b, const Complex &a) | |
482 | { | |
483 | return b < abs(a); | |
484 | } | |
485 | ||
486 | inline int operator < (const Complex &a, const Complex &b) | |
487 | { | |
488 | return abs(a) < abs(b); | |
489 | } | |
490 | ||
491 | ///////////////////////// Functions | |
492 | double real (const Complex &a) | |
493 | { | |
494 | return a.Re; | |
495 | } | |
496 | double imag (const Complex &a) | |
497 | { | |
498 | return a.Im; | |
499 | } | |
500 | double abs (const Complex &a) | |
501 | { | |
3294e9ff | 502 | if (a.Im == 0) return TMath::Abs(a.Re); |
503 | if (a.Re == 0) return TMath::Abs(a.Im); | |
504 | double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im); | |
dd82cadc | 505 | return (R >= I) ? |
ec094451 | 506 | (R * sqrt (1 + a.Im*a.Im/a.Re/a.Re)) : |
507 | (I * sqrt (1 + a.Re*a.Re/a.Im/a.Im)) ; | |
dd82cadc | 508 | } |
509 | ||
510 | double Arg (const Complex &a) | |
511 | { | |
512 | return atan2(a.Im,a.Re); | |
513 | } | |
514 | ||
515 | double phase (const Complex &a) | |
516 | { | |
517 | return atan2(a.Im,a.Re); | |
518 | } | |
519 | ||
520 | inline void Zero (Complex &a) { a.Re=a.Im=0; } | |
521 | ||
522 | Complex conj (const Complex &a) { return !a; } | |
523 | Complex pow (Complex a, int n) { return a^n; } | |
524 | Complex pow (Complex a, double n) { return a^n; } | |
525 | Complex pow (Complex a, Complex b) { return a^b; } | |
526 | ||
527 | inline Complex sqr (const Complex &a) { return a*a; } | |
528 | ||
529 | Complex sqrt (const Complex &a, int flag) | |
530 | { | |
3294e9ff | 531 | if ((a.Re>=0) && (TMath::Abs(a.Im)<small_epsilon)) |
dd82cadc | 532 | return flag ? -sqrt(a.Re) : sqrt(a.Re); |
3294e9ff | 533 | double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im); |
dd82cadc | 534 | double w = (R >= I) ? |
ec094451 | 535 | sqrt (R/2 * ( 1 + sqrt (1 + a.Im*a.Im/a.Re/a.Re))): |
536 | sqrt (I/2 * (R/I + sqrt (1 + a.Re*a.Re/a.Im/a.Im))); | |
dd82cadc | 537 | Complex c; |
538 | if (a.Re >= 0) | |
539 | { | |
540 | c.Re = w; | |
541 | c.Im = a.Im / (2*w); | |
542 | } | |
543 | else | |
544 | { | |
545 | c.Re = I / (2*w); | |
546 | c.Im = (a.Im >= 0) ? w : -w; | |
547 | } | |
548 | return ((flag && (c.Re<0)) || (!flag && (c.Re>=0))) ? c : -c; | |
549 | } | |
550 | ||
551 | Complex operator ^ (const Complex &a, int n) | |
552 | { | |
553 | Complex c(1,0); | |
554 | if (n==0) return c; | |
555 | if (n>0) for (int i=0;i<n;i++) c*=a; | |
556 | else for (int j=0;j>n;j--) c*=a; | |
557 | if (n>0) return c; | |
558 | else return 1/c; | |
559 | } | |
560 | ||
561 | Complex operator ^ (const Complex &a, double n) | |
562 | { | |
563 | return exp(n*log(a)); | |
564 | } | |
565 | ||
566 | Complex operator ^ (const Complex &a, const Complex &b) | |
567 | { | |
568 | return exp(b*log(a)); | |
569 | } | |
570 | ||
571 | Complex root (const Complex &z, int n, int k) | |
572 | { | |
573 | double c=exp(log(abs(z))/n); | |
574 | double t=(Arg(z)+2*PI*k)/n; | |
575 | return Complex (c*cos(t), c*sin(t)); | |
576 | } | |
577 | ||
578 | Complex exp (const Complex &a) | |
579 | { | |
580 | double t=exp(a.Re); | |
581 | return Complex (t*cos(a.Im), t*sin(a.Im)); | |
582 | } | |
583 | ||
584 | Complex log (const Complex &a) | |
585 | { | |
586 | if (a==0) | |
587 | ComplexError("Error in function log(Complex): argument is 0"); | |
588 | return Complex (log(abs(a)), Arg(a)); | |
589 | } | |
590 | ||
591 | Complex sin (const Complex &a) | |
592 | { | |
593 | return Complex (sin(a.Re)*cosh(a.Im), cos(a.Re)*sinh(a.Im)); | |
594 | } | |
595 | ||
596 | Complex cos (const Complex &a) | |
597 | { | |
598 | return Complex (cos(a.Re)*cosh(a.Im), -sin(a.Re)*sinh(a.Im)); | |
599 | } | |
600 | ||
601 | Complex tan (const Complex &a) | |
602 | { | |
603 | return sin(a)/cos(a); | |
604 | } | |
605 | ||
606 | Complex cot (const Complex &a) | |
607 | { | |
608 | return cos(a)/sin(a); | |
609 | } | |
610 | ||
611 | Complex sec (const Complex &a) | |
612 | { | |
613 | return 1/cos(a); | |
614 | } | |
615 | ||
616 | Complex csc (const Complex &a) | |
617 | { | |
618 | return 1/sin(a); | |
619 | } | |
620 | ||
621 | Complex sinh (const Complex &a) | |
622 | { | |
623 | return Complex (sinh(a.Re)*cos(a.Im), cosh(a.Re)*sin(a.Im)); | |
624 | } | |
625 | ||
626 | Complex cosh (const Complex &a) | |
627 | { | |
628 | return Complex (cosh(a.Re)*cos(a.Im), sinh(a.Re)*sin(a.Im)); | |
629 | } | |
630 | ||
631 | Complex tanh (const Complex &a) | |
632 | { | |
633 | return sinh(a)/cosh(a); | |
634 | } | |
635 | ||
636 | Complex coth (const Complex &a) | |
637 | { | |
638 | return cosh(a)/sinh(a); | |
639 | } | |
640 | ||
641 | Complex sech (const Complex &a) | |
642 | { | |
643 | return 1/cosh(a); | |
644 | } | |
645 | ||
646 | Complex csch (const Complex &a) | |
647 | { | |
648 | return 1/sinh(a); | |
649 | } | |
650 | //////////////////////// Inverce trigonometric functions | |
651 | ||
652 | Complex asin (const Complex &a, int flag) | |
653 | { | |
654 | return -ImUnit * log(ImUnit*a + sqrt(1-sqr(a), flag)); | |
655 | } | |
656 | ||
657 | Complex acos (const Complex &a, int flag) | |
658 | { | |
659 | return -ImUnit * log(a + ImUnit*sqrt(1-sqr(a), flag)); | |
660 | } | |
661 | ||
662 | Complex atan (const Complex &a) | |
663 | { | |
664 | return ImUnit/2 * log((ImUnit+a)/(ImUnit-a)); | |
665 | } | |
666 | ||
667 | Complex acot (const Complex &a) | |
668 | { | |
669 | return ImUnit/2 * log((a-ImUnit)/(a+ImUnit)); | |
670 | } | |
671 | ||
672 | Complex asec (const Complex &a, int flag) | |
673 | { | |
674 | return acos(1/a, flag); | |
675 | } | |
676 | ||
677 | Complex acsc (const Complex &a, int flag) | |
678 | { | |
679 | return asin(1/a, flag); | |
680 | } | |
681 | ||
682 | Complex asinh (const Complex &a, int flag) | |
683 | { | |
684 | return log(a + sqrt(sqr(a)+1, flag)); | |
685 | } | |
686 | ||
687 | Complex acosh (const Complex &a, int flag) | |
688 | { | |
689 | return log(a + sqrt(sqr(a)-1, flag)); | |
690 | } | |
691 | ||
692 | Complex atanh (const Complex &a) | |
693 | { | |
694 | return log((1+a)/(1-a)) / 2; | |
695 | } | |
696 | ||
697 | Complex acoth (const Complex &a) | |
698 | { | |
699 | return log((a+1)/(a-1)) / 2; | |
700 | } | |
701 | ||
702 | Complex asech (const Complex &a, int flag) | |
703 | { | |
704 | return acosh(1/a, flag); | |
705 | } | |
706 | ||
707 | Complex acsch (const Complex &a, int flag) | |
708 | { | |
709 | return asinh(1/a, flag); | |
710 | } | |
711 | ||
712 | Complex Polar (double a, double b) | |
713 | { | |
714 | return Complex (a*cos(b), a*sin(b)); | |
715 | } | |
716 | ||
717 | void Complex::SetPolar (double a, double b) | |
718 | { | |
719 | Re=a*cos(b); | |
720 | Im=a*sin(b); | |
721 | } | |
722 | ||
723 | void Complex::SetAbs (double a) | |
724 | { | |
725 | double b=Arg(*this); | |
726 | Re=a*cos(b); | |
727 | Im=a*sin(b); | |
728 | } | |
729 | ||
730 | void Complex::SetArg (double b) | |
731 | { | |
732 | double a=abs(*this); | |
733 | Re=a*cos(b); | |
734 | Im=a*sin(b); | |
735 | } | |
736 | ||
737 | void Solve2 (Complex* z, const Complex &b, const Complex &c) | |
738 | // finding z of equation: z^2 + b*z + c = 0 | |
739 | { | |
740 | Complex t = sqrt(sqr(b)-4*c); | |
741 | Complex q = ((!b * t).Re >= 0) ? (-(b + t) / 2) : (-(b - t) / 2); | |
742 | z[0] = q; | |
743 | z[1] = c/q; | |
744 | } | |
745 | ||
746 | Complex Solve2 (const Complex &b, const Complex &c, int RootNumber) | |
747 | { | |
748 | if ((RootNumber < 0) || (RootNumber > 1)) | |
749 | ComplexError ("Error in Solve2: wrong root number"); | |
750 | Complex z[2]; | |
751 | Solve2 (z,b,c); | |
752 | return z[RootNumber]; | |
753 | } | |
754 | ||
755 | void Solve3 (Complex* z, const Complex &a2, const Complex &a1, const Complex &a0) | |
756 | // finding z of equation: z^3 + a2*z^2 + a1*z + a0 = 0 | |
757 | { | |
758 | Complex q, r, t, a, b, zero(0,0); | |
759 | q = (sqr(a2) - 3*a1) / 9; | |
760 | r = (2*(a2^3) - 9*a2*a1 + 27*a0) / 54; | |
761 | t = sqrt((r^2) - (q^3)); | |
762 | a = ((!r * t) >=0.0) ? -((r+t)^(1.0/3)) : -((r-t)^(1.0/3)); | |
763 | b = ((a == zero) ? zero : (q/a)); | |
f93c53d8 | 764 | z[0] = -(a+b)/2 - a2/3 + ImUnit*sqrt(3.0)*(a-b)/2; |
765 | z[1] = -(a+b)/2 - a2/3 - ImUnit*sqrt(3.0)*(a-b)/2; | |
dd82cadc | 766 | z[2] = a + b - a2/3; |
767 | } | |
768 | ||
769 | Complex Solve3 (const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber) | |
770 | { | |
771 | if ((RootNumber < 0) || (RootNumber > 2)) | |
772 | ComplexError ("Error in Solve3: wrong root number"); | |
773 | Complex z[3]; | |
774 | Solve3 (z,a2,a1,a0); | |
775 | return z[RootNumber]; | |
776 | } | |
777 | ||
778 | void Solve4 (Complex* z, const Complex &a3, const Complex &a2, const Complex &a1, const Complex &a0) | |
779 | // finding z of equation: z^4 + a3*z^3 + a2*z^2 + a1*z + a0 = 0 | |
780 | { | |
781 | Complex u[3], t1, t2, t; | |
782 | Solve3 (u, -a2, (a1*a3 - 4*a0), -((a1^2) + a0*(a3^2) - 4*a0*a2)); | |
783 | t = *u; | |
784 | t1 = sqrt((a3^2)/4 - a2 + t); | |
785 | t2 = sqrt((t^2)/4 - a0); | |
786 | Solve2 ( z, (a3/2 - t1), (t/2 + t2)); | |
787 | Solve2 (&(z[2]), (a3/2 + t1), (t/2 - t2)); | |
788 | } | |
789 | ||
790 | Complex Solve4 (const Complex &a3, const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber) | |
791 | { | |
792 | if ((RootNumber < 0) || (RootNumber > 3)) | |
793 | ComplexError ("Error in Solve4: wrong root number"); | |
794 | Complex z[4]; | |
795 | Solve4 (z,a3,a2,a1,a0); | |
796 | return z[RootNumber]; | |
797 | } | |
798 | ||
799 | ||
800 | #endif // __COMPLEX_CXX__ // |