![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to generic_complex_numbers.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers|
Return to generic_complex_numbers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers |
1.1 maekawa 1: package body Generic_Complex_Numbers is
2:
3: TWO : constant number := one+one;
4:
5: -- CREATORS :
6:
7: function Create ( i : integer ) return Complex_Number is
8:
9: res : Complex_Number;
10:
11: begin
12: if i = 0
13: then res.re := +zero;
14: elsif i = 1
15: then res.re := +one;
16: else res.re := Create(i);
17: end if;
18: res.im := +zero;
19: return res;
20: end Create;
21:
22: function Create ( f : number ) return Complex_Number is
23:
24: res : Complex_Number;
25:
26: begin
27: res.re := +f;
28: res.im := +zero;
29: return res;
30: end Create;
31:
32: function Create ( re,im : number ) return Complex_Number is
33:
34: res : Complex_Number;
35:
36: begin
37: res.RE := +re;
38: res.IM := +im;
39: return res;
40: end Create;
41:
42: function Conjugate ( c : Complex_Number ) return Complex_Number is
43:
44: res : Complex_Number;
45:
46: begin
47: res.RE := +c.RE;
48: res.IM := -c.IM;
49: return res;
50: end Conjugate;
51:
52: -- COMPARISON and COPYING :
53:
54: function Equal ( x,y : Complex_Number ) return boolean is
55: begin
56: return (Equal(x.RE,y.RE) and Equal(x.IM,y.IM));
57: end Equal;
58:
59: procedure Copy ( x : in Complex_Number; y : in out Complex_Number ) is
60: begin
61: Copy(x.RE,y.RE);
62: Copy(x.IM,y.IM);
63: end Copy;
64:
65: function "<" ( x,y : Complex_Number ) return boolean is
66:
67: avx : number := AbsVal(x);
68: avy : number := AbsVal(y);
69: res : boolean := (avx < avy);
70:
71: begin
72: Clear(avx);
73: Clear(avy);
74: return res;
75: end "<";
76:
77: function ">" ( x,y : Complex_Number ) return boolean is
78:
79: avx : number := AbsVal(x);
80: avy : number := AbsVal(y);
81: res : boolean := (avx > avy);
82:
83: begin
84: Clear(avx);
85: Clear(avy);
86: return res;
87: end ">";
88:
89: -- SELECTORS :
90:
91: function REAL_PART ( x : Complex_Number ) return number is
92: begin
93: return +x.RE;
94: end REAL_PART;
95:
96: function IMAG_PART ( x : Complex_Number ) return number is
97: begin
98: return +x.IM;
99: end IMAG_PART;
100:
101: function AbsVal ( x : Complex_Number ) return number is
102:
103: res : number := AbsVal(x.RE);
104: acc : number := AbsVal(x.IM);
105:
106: begin
107: Add(res,acc);
108: Clear(acc);
109: return res;
110: end AbsVal;
111:
112: function AbsVal ( x : Complex_Number ) return Complex_Number is
113:
114: abx : number := AbsVal(x);
115: res : Complex_Number := Create(abx);
116:
117: begin
118: Clear(abx);
119: return res;
120: end AbsVal;
121:
122: -- ARITHMETIC OPERATIONS AS FUNCTIONS :
123:
124: function "+" ( x : Complex_Number; y : number ) return Complex_Number is
125: begin
126: return (x.RE+y,+x.IM);
127: end "+";
128:
129: function "-" ( x : Complex_Number; y : number ) return Complex_Number is
130: begin
131: return (x.RE-y,+x.IM);
132: end "-";
133:
134: function "*" ( x : Complex_Number; y : number ) return Complex_Number is
135: begin
136: return (x.RE*y,x.IM*y);
137: end "*";
138:
139: function "/" ( x : Complex_Number; y : number ) return Complex_Number is
140: begin
141: return (x.RE/y,x.IM/y);
142: end "/";
143:
144: function "+" ( x : number; y : Complex_Number ) return Complex_Number is
145: begin
146: return (x+y.RE,+y.IM);
147: end "+";
148:
149: function "-" ( x : number; y : Complex_Number ) return Complex_Number is
150: begin
151: return (x-y.RE,-y.IM);
152: end "-";
153:
154: function "*" ( x : number; y : Complex_Number ) return Complex_Number is
155: begin
156: return (x*y.RE,x*y.IM);
157: end "*";
158:
159: function "/" ( x : number; y : Complex_Number ) return Complex_Number is
160:
161: res : Complex_Number;
162: acc,avyim,avyre : Number;
163:
164: begin
165: if Equal(y.IM,zero)
166: then res.RE := x/y.RE;
167: res.IM := +zero;
168: elsif Equal(y.RE,zero)
169: then res.RE := +zero;
170: res.IM := x/y.IM; Min(res.IM);
171: else avyim := AbsVal(y.IM); avyre := AbsVal(y.RE);
172: if avyim < avyre
173: then acc := y.IM/y.RE;
174: res.RE := x*acc; Mul(acc,y.IM); Add(acc,y.RE); Div(res.RE,acc);
175: res.IM := x/acc; Min(res.IM);
176: Clear(acc);
177: elsif avyim > avyre
178: then acc := y.RE/y.IM;
179: res.IM := x*acc; Min(res.IM); Mul(acc,y.RE); Add(acc,y.IM);
180: res.RE := x/acc;
181: Clear(acc);
182: elsif Equal(y.IM,y.RE)
183: then acc := TWO*y.IM;
184: res.RE := x/acc;
185: res.IM := -res.RE;
186: Clear(acc);
187: else -- y.IM = -y.RE then
188: acc := TWO*y.IM;
189: res.RE := x/acc; Min(res.RE);
190: res.IM := x/acc; Min(res.IM);
191: Clear(acc);
192: end if;
193: Clear(avyim); Clear(avyre);
194: end if;
195: return res;
196: end "/";
197:
198: function "+" ( x : Complex_Number ) return Complex_Number is
199: begin
200: return (+x.RE,+x.IM);
201: end "+";
202:
203: function "-" ( x : Complex_Number ) return Complex_Number is
204: begin
205: return (-x.RE,-x.IM);
206: end "-";
207:
208: function "+" ( x,y : Complex_Number ) return Complex_Number is
209: begin
210: return (x.RE+y.RE,x.IM+y.IM);
211: end "+";
212:
213: function "-" ( x,y : Complex_Number ) return Complex_Number is
214: begin
215: return (x.RE-y.RE,x.IM-y.IM);
216: end "-";
217:
218: function "*" ( x,y : Complex_Number ) return Complex_Number is
219:
220: res : Complex_Number;
221: acc,avyim,avyre : number;
222:
223: begin
224: if Equal(y.IM,zero)
225: then res.RE := x.RE*y.RE;
226: res.IM := x.IM*y.RE;
227: elsif Equal(y.RE,zero)
228: then res.RE := x.IM*y.IM; Min(res.RE);
229: res.IM := x.RE*y.IM;
230: else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
231: if avyre < avyim
232: then acc := y.RE/y.IM;
233: res.RE := x.RE*acc; Sub(res.RE,x.IM); Mul(res.RE,y.IM);
234: res.IM := x.IM*acc; Add(res.IM,x.RE); Mul(res.IM,y.IM);
235: Clear(acc);
236: elsif avyre > avyim
237: then acc := y.IM/y.RE;
238: res.RE := x.IM*acc; Sub(res.RE,x.RE); Min(res.RE);
239: Mul(res.RE,y.RE);
240: res.IM := x.RE*acc; Add(res.IM,x.IM); Mul(res.IM,y.RE);
241: Clear(acc);
242: elsif Equal(y.RE,y.IM)
243: then res.RE := x.RE - x.IM; Mul(res.RE,y.RE);
244: res.IM := x.IM + x.RE; Mul(res.IM,y.RE);
245: else -- y.RE = -y.IM then
246: res.RE := x.RE + x.IM; Mul(res.RE,y.RE);
247: res.IM := x.IM - x.RE; Mul(res.IM,y.RE);
248: end if;
249: Clear(avyre); Clear(avyim);
250: end if;
251: return res;
252: end "*";
253:
254: function "/" ( x,y : Complex_Number ) return Complex_Number is
255:
256: res : Complex_Number;
257: acc,avyre,avyim : number;
258:
259: begin
260: if Equal(y.IM,zero)
261: then res.RE := x.RE/y.RE;
262: res.IM := x.IM/y.RE;
263: elsif Equal(y.RE,zero)
264: then res.RE := x.IM/y.IM;
265: res.IM := x.RE/y.IM; Min(res.IM);
266: else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
267: if avyre < avyim
268: then acc := y.RE/y.IM;
269: res.RE := x.RE*acc; Add(res.RE,x.IM);
270: res.IM := x.IM*acc; Sub(res.IM,x.RE);
271: Mul(acc,y.RE); Add(acc,y.IM);
272: Div(res.RE,acc);
273: Div(res.IM,acc);
274: Clear(acc);
275: elsif avyre > avyim
276: then acc := y.IM/y.RE;
277: res.RE := x.IM*acc; Add(res.RE,x.RE);
278: res.IM := x.RE*acc; Sub(res.IM,x.IM); Min(res.IM);
279: Mul(acc,y.IM); Add(acc,y.RE);
280: Div(res.RE,acc);
281: Div(res.IM,acc);
282: Clear(acc);
283: elsif Equal(y.RE,y.IM)
284: then acc := TWO*y.RE;
285: res.RE := x.RE + x.IM; Div(res.RE,acc);
286: res.IM := x.IM - x.RE; Div(res.IM,acc);
287: Clear(acc);
288: else -- y.RE = -y.IM then
289: acc := TWO*y.RE;
290: res.RE := x.RE - x.IM; Div(res.RE,acc);
291: res.IM := x.IM + x.RE; Div(res.IM,acc);
292: Clear(acc);
293: end if;
294: Clear(avyre); Clear(avyim);
295: end if;
296: return res;
297: end "/";
298:
299: function "**" ( x : Complex_Number; m : integer ) return Complex_Number is
300:
301: res : Complex_Number;
302:
303: begin
304: if m = 0
305: then res := Create(1);
306: elsif m > 0
307: then res := +x; --Copy(x,res); <- Copy CAUSED BUS ERROR!!
308: for j in 2..m loop
309: Mul(res,x);
310: end loop;
311: else res := Create(1);
312: for j in 1..(-m) loop
313: Div(res,x);
314: end loop;
315: end if;
316: return res;
317: end "**";
318:
319: -- ARITHMETIC OPERATIONS AS PROCEDURES :
320:
321: procedure Add ( x : in out Complex_Number; y : in number ) is
322: begin
323: Add(x.RE,y);
324: end Add;
325:
326: procedure Sub ( x : in out Complex_Number; y : in number ) is
327: begin
328: Sub(x.RE,y);
329: end Sub;
330:
331: procedure Mul ( x : in out Complex_Number; y : in number ) is
332: begin
333: Mul(x.RE,y);
334: Mul(x.IM,y);
335: end Mul;
336:
337: procedure Div ( x : in out Complex_Number; y : in number ) is
338: begin
339: Div(x.RE,y);
340: Div(x.IM,y);
341: end Div;
342:
343: procedure Add ( x : in out Complex_Number; y : in Complex_Number ) is
344: begin
345: Add(x.RE,y.RE);
346: Add(x.IM,y.IM);
347: end Add;
348:
349: procedure Sub ( x : in out Complex_Number; y : in Complex_Number ) is
350: begin
351: Sub(x.RE,y.RE);
352: Sub(x.IM,y.IM);
353: end Sub;
354:
355: procedure Min ( x : in out Complex_Number ) is
356: begin
357: Min(x.RE);
358: Min(x.IM);
359: end Min;
360:
361: procedure Mul ( x : in out Complex_Number; y : in Complex_Number ) is
362:
363: res : Complex_Number;
364: acc,avyim,avyre : number;
365:
366: begin
367: if Equal(y.IM,zero)
368: then Mul(x.RE,y.RE);
369: Mul(x.IM,y.RE);
370: elsif Equal(y.RE,zero)
371: then res.RE := x.IM*y.IM; Min(res.RE);
372: res.IM := x.RE*y.IM;
373: Clear(x); x := res;
374: else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
375: if avyre < avyim
376: then acc := y.RE/y.IM;
377: res.RE := x.RE*acc; Sub(res.RE,x.IM); Mul(res.RE,y.IM);
378: res.IM := x.IM*acc; Add(res.IM,x.RE); Mul(res.IM,y.IM);
379: Clear(acc);
380: elsif avyre > avyim
381: then acc := y.IM/y.RE;
382: res.RE := x.IM*acc; Sub(res.RE,x.RE); Min(res.RE);
383: Mul(res.RE,y.RE);
384: res.IM := x.RE*acc; Add(res.IM,x.IM); Mul(res.IM,y.RE);
385: Clear(acc);
386: elsif Equal(y.RE,y.IM)
387: then res.RE := x.RE - x.IM; Mul(res.RE,y.RE);
388: res.IM := x.IM + x.RE; Mul(res.IM,y.RE);
389: else -- y.RE = -y.IM then
390: res.RE := x.RE + x.IM; Mul(res.RE,y.RE);
391: res.IM := x.IM - x.RE; Mul(res.IM,y.RE);
392: end if;
393: Clear(avyre); Clear(avyim);
394: Clear(x); x := res;
395: end if;
396: end Mul;
397:
398: procedure Div ( x : in out Complex_Number; y : in Complex_Number ) is
399:
400: res : Complex_Number;
401: acc,avyre,avyim : number;
402:
403: begin
404: if Equal(y.IM,zero)
405: then Div(x.RE,y.RE);
406: Div(x.IM,y.RE);
407: elsif Equal(y.IM,zero)
408: then res.RE := x.IM/y.IM;
409: res.IM := x.RE/y.IM; Min(res.IM);
410: Clear(x); x := res;
411: else avyre := AbsVal(y.RE); avyim := AbsVal(y.IM);
412: if avyre < avyim
413: then acc := y.RE/y.IM;
414: res.RE := x.RE*acc; Add(res.RE,x.IM);
415: res.IM := x.IM*acc; Sub(res.IM,x.RE);
416: Mul(acc,y.RE); Add(acc,y.IM);
417: Div(res.RE,acc);
418: Div(res.IM,acc);
419: Clear(acc);
420: elsif avyre > avyim
421: then acc := y.IM/y.RE;
422: res.RE := x.IM*acc; Add(res.RE,x.RE);
423: res.IM := x.RE*acc; Sub(res.IM,x.IM); Min(res.IM);
424: Mul(acc,y.IM); Add(acc,y.RE);
425: Div(res.RE,acc);
426: Div(res.IM,acc);
427: Clear(acc);
428: elsif Equal(y.RE,y.IM)
429: then acc := TWO*y.RE;
430: res.RE := x.RE + x.IM; Div(res.RE,acc);
431: res.IM := x.IM - x.RE; Div(res.IM,acc);
432: Clear(acc);
433: else -- y.RE = -y.IM then
434: acc := TWO*y.RE;
435: res.RE := x.RE - x.IM; Div(res.RE,acc);
436: res.IM := x.IM + x.RE; Div(res.IM,acc);
437: Clear(acc);
438: end if;
439: Clear(avyre); Clear(avyim);
440: Clear(x); x := res;
441: end if;
442: end Div;
443:
444: -- DESTRUCTOR :
445:
446: procedure Clear ( x : in out Complex_Number ) is
447: begin
448: Clear(x.RE);
449: Clear(x.IM);
450: end Clear;
451:
452: end Generic_Complex_Numbers;