|
|
1.1 root 1:
2: #include "complex.h"
3: #include "const.h"
4:
5: complex operator*(complex a1,complex a2)
6: {
7: return complex(a1.re*a2.re-a1.im*a2.im, a1.re*a2.im+a1.im*a2.re);
8: }
9:
10:
11: complex operator/(complex a1, complex a2)
12: {
13: register double r = a2.re; /* (r,i) */
14: register double i = a2.im;
15: register double ti; /* (tr,ti) */
16: register double tr;
17:
18: tr = ABS(r);
19: ti = ABS(i);
20:
21: if (tr <= ti) {
22: ti = r/i;
23: tr = i * (1 + ti*ti);
24: r = a1.re;
25: i = a1.im;
26: }
27: else {
28: ti = -i/r;
29: tr = r * (1 + ti*ti);
30: r = -a1.im;
31: i = a1.re;
32: }
33:
34: return complex( (r*ti + i)/tr, (i*ti - r)/tr );
35: }
36:
37: void complex.operator*=(complex a)
38: {
39: register double r = re*a.re - im*a.im;
40: register double i = re*a.im + im*a.re;
41: re = r;
42: im = i;
43: }
44:
45: void complex.operator/=(complex a)
46: {
47: complex quot, temp1, temp2;
48:
49: if ( (temp2.re = a.re) < 0 ) temp2.re = -temp2.re;
50: if ( (temp2.im = a.im) < 0 ) temp2.im = -temp2.im;
51: if ( temp2.re <= temp2.im) {
52: temp2.im = a.re/a.im;
53: temp2.re = a.im * (1 + temp2.im*temp2.im);
54: temp1 = *this;
55: }
56: else {
57: temp2.im = -a.im/a.re;
58: temp2.re = a.re * (1 + temp2.im*temp2.im);
59: temp1.re = -im;
60: temp1.im = re;
61: }
62: re = (temp1.re * temp2.im + temp1.im) / temp2.re;
63: im = (temp1.im * temp2.im - temp1.re) / temp2.re;
64: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.