1: // public domain
2:
3: #include <math.h>
4:
5: #include "bcomplex.h"
6:
7: /*********************************************************/
8:
9: Complex_ld Complex_ld::div(Complex_ld &x, Complex_ld &y)
10: {
11: Complex_ld q;
12: long double r;
13: long double den;
14:
15: if (fabs(y.re) < fabs(y.im))
16: {
17: r = y.re / y.im;
18: den = y.im + r * y.re;
19: q.re = (x.re * r + x.im) / den;
20: q.im = (x.im * r - x.re) / den;
21: }
22: else
23: {
24: r = y.im / y.re;
25: den = y.re + r * y.im;
26: q.re = (x.re + r * x.im) / den;
27: q.im = (x.im - r * x.re) / den;
28: }
29: return q;
30: }
31:
32: Complex_ld Complex_ld::mul(Complex_ld &x, Complex_ld &y)
33: {
34: Complex_ld p;
35:
36: p.re = x.re * y.re - x.im * y.im;
37: p.im = x.im * y.re + x.re * y.im;
38: return p;
39: }
40:
41: long double Complex_ld::abs(Complex_ld &z)
42: {
43: long double x,y,ans,temp;
44:
45: x = fabs(z.re);
46: y = fabs(z.im);
47: if (x == 0)
48: ans = y;
49: else if (y == 0)
50: ans = x;
51: else if (x > y)
52: {
53: temp = y / x;
54: ans = x * sqrt(1 + temp * temp);
55: }
56: else
57: {
58: temp = x / y;
59: ans = y * sqrt(1 + temp * temp);
60: }
61: return ans;
62: }
63:
64: Complex_ld Complex_ld::sqrtc(Complex_ld &z)
65: {
66: Complex_ld c;
67: long double x,y,w,r;
68:
69: if (z.re == 0 && z.im == 0)
70: {
71: c.re = 0;
72: c.im = 0;
73: }
74: else
75: {
76: x = fabs(z.re);
77: y = fabs(z.im);
78: if (x >= y)
79: {
80: r = y / x;
81: w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
82: }
83: else
84: {
85: r = x / y;
86: w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
87: }
88: if (z.re >= 0)
89: {
90: c.re = w;
91: c.im = z.im / (w + w);
92: }
93: else
94: {
95: c.im = (z.im >= 0) ? w : -w;
96: c.re = z.im / (c.im + c.im);
97: }
98: }
99: return c;
100: }
101:
102: /*********************************************************/
103:
104: Complex_d Complex_d::div(Complex_d &x, Complex_d &y)
105: {
106: Complex_d q;
107: long double r;
108: long double den;
109:
110: if (fabs(y.re) < fabs(y.im))
111: {
112: r = y.re / y.im;
113: den = y.im + r * y.re;
114: q.re = (x.re * r + x.im) / den;
115: q.im = (x.im * r - x.re) / den;
116: }
117: else
118: {
119: r = y.im / y.re;
120: den = y.re + r * y.im;
121: q.re = (x.re + r * x.im) / den;
122: q.im = (x.im - r * x.re) / den;
123: }
124: return q;
125: }
126:
127: Complex_d Complex_d::mul(Complex_d &x, Complex_d &y)
128: {
129: Complex_d p;
130:
131: p.re = x.re * y.re - x.im * y.im;
132: p.im = x.im * y.re + x.re * y.im;
133: return p;
134: }
135:
136: long double Complex_d::abs(Complex_d &z)
137: {
138: long double x,y,ans,temp;
139:
140: x = fabs(z.re);
141: y = fabs(z.im);
142: if (x == 0)
143: ans = y;
144: else if (y == 0)
145: ans = x;
146: else if (x > y)
147: {
148: temp = y / x;
149: ans = x * sqrt(1 + temp * temp);
150: }
151: else
152: {
153: temp = x / y;
154: ans = y * sqrt(1 + temp * temp);
155: }
156: return ans;
157: }
158:
159: Complex_d Complex_d::sqrtc(Complex_d &z)
160: {
161: Complex_d c;
162: long double x,y,w,r;
163:
164: if (z.re == 0 && z.im == 0)
165: {
166: c.re = 0;
167: c.im = 0;
168: }
169: else
170: {
171: x = fabs(z.re);
172: y = fabs(z.im);
173: if (x >= y)
174: {
175: r = y / x;
176: w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
177: }
178: else
179: {
180: r = x / y;
181: w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
182: }
183: if (z.re >= 0)
184: {
185: c.re = w;
186: c.im = z.im / (w + w);
187: }
188: else
189: {
190: c.im = (z.im >= 0) ? w : -w;
191: c.re = z.im / (c.im + c.im);
192: }
193: }
194: return c;
195: }
196:
197: /*********************************************************/
198:
199: Complex_f Complex_f::div(Complex_f &x, Complex_f &y)
200: {
201: Complex_f q;
202: long double r;
203: long double den;
204:
205: if (fabs(y.re) < fabs(y.im))
206: {
207: r = y.re / y.im;
208: den = y.im + r * y.re;
209: q.re = (x.re * r + x.im) / den;
warning C4244: '=' : conversion from 'long double' to 'float', possible loss of data
210: q.im = (x.im * r - x.re) / den;
warning C4244: '=' : conversion from 'long double' to 'float', possible loss of data
211: }
212: else
213: {
214: r = y.im / y.re;
215: den = y.re + r * y.im;
216: q.re = (x.re + r * x.im) / den;
warning C4244: '=' : conversion from 'long double' to 'float', possible loss of data
217: q.im = (x.im - r * x.re) / den;
warning C4244: '=' : conversion from 'long double' to 'float', possible loss of data
218: }
219: return q;
220: }
221:
222: Complex_f Complex_f::mul(Complex_f &x, Complex_f &y)
223: {
224: Complex_f p;
225:
226: p.re = x.re * y.re - x.im * y.im;
227: p.im = x.im * y.re + x.re * y.im;
228: return p;
229: }
230:
231: long double Complex_f::abs(Complex_f &z)
232: {
233: long double x,y,ans,temp;
234:
235: x = fabs(z.re);
236: y = fabs(z.im);
237: if (x == 0)
238: ans = y;
239: else if (y == 0)
240: ans = x;
241: else if (x > y)
242: {
243: temp = y / x;
244: ans = x * sqrt(1 + temp * temp);
245: }
246: else
247: {
248: temp = x / y;
249: ans = y * sqrt(1 + temp * temp);
250: }
251: return ans;
252: }
253:
254: Complex_f Complex_f::sqrtc(Complex_f &z)
255: {
256: Complex_f c;
257: long double x,y,w,r;
258:
259: if (z.re == 0 && z.im == 0)
260: {
261: c.re = 0;
262: c.im = 0;
263: }
264: else
265: {
266: x = fabs(z.re);
267: y = fabs(z.im);
268: if (x >= y)
269: {
270: r = y / x;
271: w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
272: }
273: else
274: {
275: r = x / y;
276: w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
277: }
278: if (z.re >= 0)
279: {
280: c.re = w;
warning C4244: '=' : conversion from 'long double' to 'float', possible loss of data
281: c.im = z.im / (w + w);
warning C4244: '=' : conversion from 'long double' to 'float', possible loss of data
282: }
283: else
284: {
285: c.im = (z.im >= 0) ? w : -w;
warning C4244: '=' : conversion from 'long double' to 'float', possible loss of data
286: c.re = z.im / (c.im + c.im);
287: }
288: }
289: return c;
290: }
291:
292:
293: