![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multprec_floating_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 multprec_floating_numbers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers |
1.1 maekawa 1: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
2:
3: package body Multprec_Floating_Numbers is
4:
5: -- NOTES ON THE IMPLEMENTATION :
6: -- 0) See the notes in the bodies of Multprec_{Natural,Integer}_Numbers.
7: -- 1) 0.0 must have zero fraction and zero exponent.
8: -- 2) The truncate and round operations all depend on the decimal radix.
9: -- 3) The function "Decimal_Places" occupies a central place.
10: -- 4) The normalized representation has a nonzero least significant digit.
11:
12: tol : constant double_float := 10.0**(-8);
13:
14: -- AUXILIARIES :
15:
16: function Truncate ( f : double_float ) return integer is
17:
18: i : integer := integer(f);
19:
20: begin
21: if i >= 0
22: then if double_float(i) > f + tol
23: then i := i-1;
24: end if;
25: else if double_float(i) < f - tol
26: then i := i+1;
27: end if;
28: end if;
29: return i;
30: end Truncate;
31:
32: function Max_Size ( i1,i2 : Integer_Number ) return natural is
33:
34: -- DESCRIPTION :
35: -- Return the maximal size of i1 and i2.
36:
37: sz1 : constant natural := Size(i1);
38: sz2 : constant natural := Size(i2);
39:
40: begin
41: if sz1 > sz2
42: then return sz1;
43: else return sz2;
44: end if;
45: end Max_Size;
46:
47: procedure Shift_Left ( n : in out Natural_Number; ns : out natural ) is
48:
49: -- DESCRIPTION :
50: -- Shifts n until the leading coefficient is larger than Base/10.
51: -- The number ns on return equals the number of shifts.
52:
53: begin
54: if (Empty(n) or else Equal(n,0))
55: then ns := 0;
56: else ns := 0;
57: while (Coefficient(n,Size(n)) < Base/10) loop
58: Mul(n,10);
59: ns := ns+1;
60: end loop;
61: end if;
62: end Shift_Left;
63:
64: procedure Shift_Right ( n : in out Natural_Number; ns : out natural ) is
65:
66: -- DESCRIPTION :
67: -- Shifts n until the last digit in the least significant coefficient
68: -- of n is different from zero.
69: -- The number ns on return equals the number of shifts.
70:
71: cff : natural;
72:
73: begin
74: if (Empty(n) or else Equal(n,0))
75: then ns := 0;
76: else ns := 0;
77: cff := Coefficient(n,0) mod 10;
78: while (cff = 0) loop
79: Div(n,10);
80: ns := ns+1;
81: cff := Coefficient(n,0) mod 10;
82: end loop;
83: end if;
84: end Shift_Right;
85:
86: procedure Shift_Right ( i : in out Integer_Number; ns : out natural ) is
87:
88: -- DESCRIPTION :
89: -- Applies the Shift_Right to the number representation of i.
90:
91: ni : Natural_Number;
92:
93: begin
94: if Empty(i)
95: then ns := 0;
96: else --ni := Unsigned(i);
97: Copy(Unsigned(i),ni); -- avoid sharing
98: Shift_Right(ni,ns);
99: if Sign(i) < 0
100: then Clear(i); i := Create(ni); Min(i);
101: else Clear(i); i := Create(ni);
102: end if;
103: Clear(ni);
104: end if;
105: end Shift_Right;
106:
107: procedure Normalize ( f : in out Floating_Number ) is
108:
109: -- DESCRIPTION :
110: -- Normalizes the representation of f so that its fraction has as
111: -- rightmost digit a nonzero number.
112:
113: ns : natural;
114:
115: begin
116: Shift_Right(f.fraction,ns);
117: Add(f.exponent,ns);
118: end Normalize;
119:
120: procedure Adjust_Expo ( f : in out double_float;
121: expo : in Integer_Number ) is
122:
123: -- DESCRIPTION :
124: -- Multiplies f with 10**expo.
125:
126: cnt : Integer_Number;
127:
128: begin
129: Copy(expo,cnt);
130: while cnt > 0 loop
131: f := f*10.0;
132: Sub(cnt,1);
133: end loop;
134: while cnt < 0 loop
135: f := f/10.0;
136: Add(cnt,1);
137: end loop;
138: Clear(cnt);
139: end Adjust_Expo;
140:
141: -- CONSTRUCTORS :
142:
143: function Create ( i : integer ) return Floating_Number is
144:
145: res : Floating_Number;
146:
147: begin
148: res.fraction := Create(i);
149: res.exponent := Create(0);
150: Normalize(res);
151: return res;
152: end Create;
153:
154: function Create ( i : Integer_Number ) return Floating_Number is
155:
156: res : Floating_Number;
157:
158: begin
159: res.fraction := +i;
160: res.exponent := Create(0);
161: Normalize(res);
162: return res;
163: end Create;
164:
165: function Create ( fraction,exponent : integer ) return Floating_Number is
166:
167: res : Floating_Number;
168:
169: begin
170: res.fraction := Create(fraction);
171: res.exponent := Create(exponent);
172: Normalize(res);
173: return res;
174: end Create;
175:
176: function Create ( fraction : Integer_Number;
177: exponent : integer ) return Floating_Number is
178:
179: res : Floating_Number;
180:
181: begin
182: res.fraction := +fraction;
183: res.exponent := Create(exponent);
184: Normalize(res);
185: return res;
186: end Create;
187:
188: function Create ( fraction : integer;
189: exponent : Integer_Number ) return Floating_Number is
190:
191: res : Floating_Number;
192:
193: begin
194: res.fraction := Create(fraction);
195: res.exponent := +exponent;
196: Normalize(res);
197: return res;
198: end Create;
199:
200: function Create ( fraction,exponent : Integer_Number )
201: return Floating_Number is
202:
203: res : Floating_Number;
204:
205: begin
206: res.fraction := +fraction;
207: res.exponent := +exponent;
208: Normalize(res);
209: return res;
210: end Create;
211:
212: function Create ( f : double_float ) return Floating_Number is
213:
214: -- NOTE :
215: -- To avoid rounding errors, comparison operations are prefered
216: -- above relying on intermediate results of calculations, although
217: -- rounding errors still occur!
218:
219: res : Floating_Number;
220: ten : Integer_Number;
221: wrkflt,fltexp,lowflt,uppflt,incflt : double_float;
222: intexp : integer;
223: fraction,exponent,intinc : Integer_Number;
224:
225: begin
226: if f = 0.0
227: then res := Create(0);
228: return res;
229: end if;
230: if f > 0.0
231: then wrkflt := f;
232: else wrkflt := -f;
233: end if;
234: fltexp := LOG10(wrkflt);
235: intexp := Truncate(fltexp);
236: if intexp < 15
237: then wrkflt := wrkflt*(10.0**(15-intexp));
238: exponent := Create(intexp-15);
239: intexp := 15;
240: else exponent := Create(0);
241: end if;
242: lowflt := 10.0**intexp;
243: uppflt := 10.0**(intexp+1);
244: while lowflt > wrkflt loop
245: intexp := intexp-1;
246: lowflt := 10.0**intexp;
247: uppflt := 10.0**(intexp+1);
248: end loop;
249: ten := Create(10);
250: fraction := ten**intexp;
251: for k in 1..16 loop
252: incflt := 10.0**intexp;
253: intinc := ten**intexp;
254: uppflt := lowflt;
255: for i in 1..10 loop
256: uppflt := uppflt + incflt;
257: exit when (uppflt > wrkflt);
258: lowflt := uppflt;
259: Add(fraction,intinc);
260: end loop;
261: intexp := intexp - 1;
262: Clear(intinc);
263: end loop;
264: Clear(ten);
265: if f < 0.0
266: then Min(fraction);
267: end if;
268: res := Create(fraction,exponent);
269: Normalize(res);
270: return res;
271: end Create;
272:
273: -- SELECTORS :
274:
275: function Fraction ( f : Floating_Number ) return Integer_Number is
276: begin
277: return f.fraction;
278: end Fraction;
279:
280: function Exponent ( f : Floating_Number ) return Integer_Number is
281: begin
282: return f.exponent;
283: end Exponent;
284:
285: function Size_Fraction ( f : Floating_Number ) return natural is
286: begin
287: return Size(Fraction(f));
288: end Size_Fraction;
289:
290: function Size_Exponent ( f : Floating_Number ) return natural is
291: begin
292: return Size(Exponent(f));
293: end Size_Exponent;
294:
295: function Decimal_Places_Fraction ( f : Floating_Number ) return natural is
296: begin
297: return Decimal_Places(Fraction(f));
298: end Decimal_Places_Fraction;
299:
300: function Decimal_Places_Exponent ( f : Floating_Number ) return natural is
301: begin
302: return Decimal_Places(Exponent(f));
303: end Decimal_Places_Exponent;
304:
305: function Decimal_to_Size ( deci : natural ) return natural is
306:
307: res : natural;
308: exp : constant natural := Multprec_Natural_Numbers.Exponent;
309:
310: begin
311: res := deci/exp;
312: if res*exp < deci
313: then res := res+1;
314: end if;
315: return res;
316: end Decimal_to_Size;
317:
318: function AbsVal ( f : Floating_Number ) return Floating_Number is
319:
320: res : Floating_Number;
321:
322: begin
323: if Negative(f.fraction)
324: then res.fraction := -f.fraction;
325: else res.fraction := +f.fraction;
326: end if;
327: res.exponent := +f.exponent;
328: return res;
329: end AbsVal;
330:
331: -- COMPARISON AND COPYING :
332:
333: function Equal ( f1 : Floating_Number; f2 : double_float ) return boolean is
334:
335: ff2 : Floating_Number := Create(f2);
336: res : boolean := Equal(f1,ff2);
337:
338: begin
339: Clear(ff2);
340: return res;
341: end Equal;
342:
343: function Equal ( f1,f2 : Floating_Number ) return boolean is
344:
345: res : boolean;
346: f1deci : constant natural := Decimal_Places(f1.fraction);
347: f2deci : constant natural := Decimal_Places(f2.fraction);
348: f1expo : Integer_Number := f1.exponent + f1deci;
349: f2expo : Integer_Number := f2.exponent + f2deci;
350: mulfra : Integer_Number;
351:
352: begin
353: if Equal(f1expo,f2expo)
354: then if f1deci = f2deci
355: then res := Equal(f1.fraction,f2.fraction);
356: else if f1deci < f2deci
357: then mulfra := f1.fraction*10;
358: for i in 1..(f2deci-f1deci-1) loop
359: Mul(mulfra,10);
360: end loop;
361: res := Equal(mulfra,f2.fraction);
362: else mulfra := f2.fraction*10;
363: for i in 1..(f1deci-f2deci-1) loop
364: Mul(mulfra,10);
365: end loop;
366: res := Equal(f1.fraction,mulfra);
367: end if;
368: Clear(mulfra);
369: end if;
370: else res := false;
371: end if;
372: Clear(f1expo); Clear(f2expo);
373: return res;
374: end Equal;
375:
376: function "<" ( f1 : double_float; f2 : Floating_Number ) return boolean is
377:
378: ff1 : Floating_Number := Create(f1);
379: res : boolean := (ff1 < f2);
380:
381: begin
382: Clear(ff1);
383: return res;
384: end "<";
385:
386: function "<" ( f1 : Floating_Number; f2 : double_float ) return boolean is
387:
388: ff2 : Floating_Number := Create(f2);
389: res : boolean := (f1 < ff2);
390:
391: begin
392: Clear(ff2);
393: return res;
394: end "<";
395:
396: function "<" ( f1,f2 : Floating_Number ) return boolean is
397:
398: res : boolean;
399: f1deci : constant natural := Decimal_Places(f1.fraction);
400: f2deci : constant natural := Decimal_Places(f2.fraction);
401: f1expo : Integer_Number := f1.exponent + f1deci;
402: f2expo : Integer_Number := f2.exponent + f2deci;
403: mulfra : Integer_Number;
404:
405: begin
406: if (f1expo < f2expo)
407: then if Sign(f2.fraction) > 0
408: then res := true;
409: elsif Sign(f2.fraction) < 0
410: then res := false;
411: else res := (Sign(f1.fraction) < 0);
412: end if;
413: elsif (f1expo > f2expo)
414: then if Sign(f1.fraction) > 0
415: then res := false;
416: elsif Sign(f1.fraction) < 0
417: then res := true;
418: else res := (Sign(f2.fraction) > 0);
419: end if;
420: else if f1deci = f2deci
421: then res := (f1.fraction < f2.fraction);
422: else if f1deci < f2deci
423: then mulfra := f1.fraction*10;
424: for i in 1..(f2deci-f1deci-1) loop
425: Mul(mulfra,10);
426: end loop;
427: res := (mulfra < f2.fraction);
428: else mulfra := f2.fraction*10;
429: for i in 1..(f1deci-f2deci-1) loop
430: Mul(mulfra,10);
431: end loop;
432: res := (f1.fraction < mulfra);
433: end if;
434: Clear(mulfra);
435: end if;
436: end if;
437: Clear(f1expo); Clear(f2expo);
438: return res;
439: end "<";
440:
441: function ">" ( f1 : double_float; f2 : Floating_Number ) return boolean is
442:
443: ff1 : Floating_Number := Create(f1);
444: res : boolean := (ff1 > f2);
445:
446: begin
447: Clear(ff1);
448: return res;
449: end ">";
450:
451: function ">" ( f1 : Floating_Number; f2 : double_float ) return boolean is
452:
453: ff2 : Floating_Number := Create(f2);
454: res : boolean := (f1 > ff2);
455:
456: begin
457: Clear(ff2);
458: return res;
459: end ">";
460:
461: function ">" ( f1,f2 : Floating_Number ) return boolean is
462:
463: res : boolean;
464: f1deci : constant natural := Decimal_Places(f1.fraction);
465: f2deci : constant natural := Decimal_Places(f2.fraction);
466: f1expo : Integer_Number := f1.exponent + f1deci;
467: f2expo : Integer_Number := f2.exponent + f2deci;
468: mulfra : Integer_Number;
469:
470: begin
471: if (f1expo > f2expo)
472: then if Sign(f1.fraction) > 0
473: then res := true;
474: elsif Sign(f1.fraction) < 0
475: then res := false;
476: else res := (Sign(f2.fraction) < 0);
477: end if;
478: elsif (f1expo < f2expo)
479: then if Sign(f2.fraction) > 0
480: then res := false;
481: elsif Sign(f2.fraction) < 0
482: then res := true;
483: else res := (Sign(f1.fraction) > 0);
484: end if;
485: else if f1deci = f2deci
486: then res := (f1.fraction > f2.fraction);
487: else if f1deci < f2deci
488: then mulfra := f1.fraction*10;
489: for i in 1..(f2deci-f1deci-1) loop
490: Mul(mulfra,10);
491: end loop;
492: res := (mulfra > f2.fraction);
493: else mulfra := f2.fraction*10;
494: for i in 1..(f1deci-f2deci-1) loop
495: Mul(mulfra,10);
496: end loop;
497: res := (f1.fraction > mulfra);
498: end if;
499: Clear(mulfra);
500: end if;
501: end if;
502: Clear(f1expo); Clear(f2expo);
503: return res;
504: end ">";
505:
506: procedure Copy ( f1 : in Floating_Number; f2 : in out Floating_Number ) is
507: begin
508: Clear(f2);
509: Copy(f1.fraction,f2.fraction);
510: Copy(f1.exponent,f2.exponent);
511: end Copy;
512:
513: -- EXPANDING AND SHORTENING THE MANTISSA :
514:
515: function Expand ( i : Integer_Number; k : natural )
516: return Array_of_Naturals is
517:
518: -- DESCRIPTION :
519: -- Returns the expanded coefficient vector of i.
520:
521: cff : constant Array_of_Naturals := Coefficients(Unsigned(i));
522: res : Array_of_Naturals(cff'first..cff'last+k);
523:
524: begin
525: res(cff'range) := cff(cff'range);
526: res(cff'last+1..res'last) := (cff'last+1..res'last => 0);
527: return res;
528: end Expand;
529:
530: procedure Truncate ( i : in Integer_Number; k : in natural;
531: trc : out Array_of_Naturals; exp : out integer ) is
532:
533: -- DESCRIPTION :
534: -- On return is the truncated coefficient vector of i and the number
535: -- of decimal places that are lost, or gained when exp is positive.
536:
537: n : Natural_Number;
538: ns : natural; -- #shifts = #decimal places gained
539:
540: begin
541: Copy(Unsigned(i),n);
542: Shift_Left(n,ns);
543: exp := k*Multprec_Natural_Numbers.Exponent - ns;
544: declare
545: cff : constant Array_of_Naturals := Coefficients(n);
546: begin
547: for i in trc'range loop
548: trc(i) := cff(i+k);
549: end loop;
550: Clear(n);
551: end;
552: end Truncate;
553:
554: procedure Truncate ( i : in out Integer_Number; k : in natural;
555: exp : out integer ) is
556:
557: n : Natural_Number := Unsigned(i);
558: neg : constant boolean := Negative(i);
559: ns : natural;
560: newint : Integer_Number;
561:
562: begin
563: Shift_Left(n,ns);
564: exp := k*Multprec_Natural_Numbers.Exponent - ns;
565: declare
566: cff : constant Array_of_Naturals := Coefficients(n);
567: trc : Array_of_Naturals(cff'first..cff'last-k);
568: begin
569: for i in trc'range loop
570: trc(i) := cff(i+k);
571: end loop;
572: newint := Create(trc);
573: i := newint; --Copy(newint,i); --Clear(newint);
574: if neg
575: then Min(i);
576: end if;
577: end;
578: Clear(n);
579: end Truncate;
580:
581: procedure Round ( i : in Integer_Number; k : in natural;
582: trc : out Array_of_Naturals; exp : out integer ) is
583:
584: -- DESCRIPTION :
585: -- On return is the rounded coefficient vector of i and the number
586: -- of decimal places that are lost, or gained when exp is positive.
587:
588: n : Natural_Number := Unsigned(i);
589: ns : natural; -- #shifts = #decimal places gained
590:
591: begin
592: Shift_Left(n,ns);
593: exp := k*Multprec_Natural_Numbers.Exponent - ns;
594: declare
595: cff : constant Array_of_Naturals := Coefficients(n);
596: begin
597: for i in trc'range loop
598: trc(i) := cff(i+k);
599: end loop;
600: if k > 0
601: then if cff(k-1) >= Multprec_Natural_Numbers.Base/2
602: then trc(0) := trc(0)+1;
603: end if;
604: end if;
605: end;
606: Clear(n);
607: end Round;
608:
609: function Expand ( f : Floating_Number; k : natural ) return Floating_Number is
610:
611: res : Floating_Number;
612: cff : constant Array_of_Naturals := Expand(f.fraction,k);
613:
614: begin
615: res.fraction := Create(Create(cff));
616: if Negative(f.fraction)
617: then Min(res.fraction);
618: end if;
619: Copy(f.exponent,res.exponent);
620: Normalize(res);
621: return res;
622: end Expand;
623:
624: procedure Expand ( f : in out Floating_Number; k : in natural ) is
625:
626: cff : constant Array_of_Naturals := Expand(f.fraction,k);
627: neg : constant boolean := Negative(f.fraction);
628:
629: begin
630: Clear(f.fraction);
631: f.fraction := Create(Create(cff));
632: if neg
633: then Min(f.fraction);
634: end if;
635: Normalize(f);
636: end Expand;
637:
638: function Round ( f : Floating_Number ) return double_float is
639:
640: res : double_float;
641: szf,ns : natural;
642: n : Natural_Number;
643: expo : Integer_Number;
644:
645: begin
646: if (Empty(f.fraction) or else Equal(f.fraction,0))
647: then res := 0.0;
648: else Copy(Unsigned(f.fraction),n);
649: Shift_Left(n,ns);
650: expo := f.exponent - ns;
651: szf := Size(n);
652: while Coefficient(n,szf) = 0 and szf > 0 loop
653: szf := szf - 1;
654: end loop;
655: res := double_float(Coefficient(n,szf));
656: if szf >= 1
657: then res := res*double_float(Multprec_Natural_Numbers.Base)
658: + double_float(Coefficient(n,szf-1));
659: if szf >= 2
660: then if Coefficient(n,szf-2) > Multprec_Natural_Numbers.Base/2
661: then res := res+1.0;
662: end if;
663: ns := (szf - 1)*Multprec_Natural_Numbers.Exponent;
664: Add(expo,ns);
665: end if;
666: end if;
667: Adjust_Expo(res,expo);
668: Clear(n); Clear(expo);
669: end if;
670: return res;
671: end Round;
672:
673: function Round ( f : Floating_Number; k : natural ) return Floating_Number is
674:
675: res : Floating_Number;
676: cff : Array_of_Naturals(0..(Size(f.fraction)-k));
677: exp : integer;
678:
679: begin
680: Round(f.fraction,k,cff,exp);
681: res.fraction := Create(Create(cff));
682: if Negative(f.fraction)
683: then Min(res.fraction);
684: end if;
685: Copy(f.exponent,res.exponent);
686: Add(res.exponent,exp);
687: Normalize(res);
688: return res;
689: end Round;
690:
691: procedure Round ( f : in out Floating_Number; k : in natural ) is
692:
693: exp : integer;
694: n : Natural_Number;
695: neg : constant boolean := Negative(f.fraction);
696: ns : natural;
697:
698: begin
699: Copy(Unsigned(f.fraction),n);
700: Shift_Left(n,ns);
701: exp := k*Multprec_Natural_Numbers.Exponent - ns;
702: declare
703: cff : constant Array_of_Naturals := Coefficients(n);
704: trc : Array_of_Naturals(cff'first..cff'last-k);
705: begin
706: for i in trc'range loop
707: trc(i) := cff(i+k);
708: end loop;
709: if k > 0
710: then if cff(k-1) >= Multprec_Natural_Numbers.Base/2
711: then trc(0) := trc(0)+1;
712: end if;
713: end if;
714: Clear(f.fraction);
715: f.fraction := Create(trc);
716: if neg
717: then Min(f.fraction);
718: end if;
719: end;
720: Clear(n);
721: Add(f.exponent,exp);
722: Normalize(f);
723: end Round;
724:
725: function Trunc ( f : Floating_Number ) return double_float is
726:
727: res : double_float;
728: szf,ns : natural;
729: n : Natural_Number;
730: expo : Integer_Number;
731:
732: begin
733: if (Empty(f.fraction) or else Equal(f.fraction,0))
734: then res := 0.0;
735: else n := Unsigned(f.fraction);
736: Shift_Left(n,ns);
737: expo := f.exponent - ns;
738: szf := Size(n);
739: while Coefficient(n,szf) = 0 and szf > 0 loop
740: szf := szf - 1;
741: end loop;
742: res := double_float(Coefficient(f.fraction,szf));
743: if szf >= 1
744: then res := res*double_float(Multprec_Natural_Numbers.Base)
745: + double_float(Coefficient(f.fraction,szf-1));
746: if szf >= 2
747: then ns := (szf - 1)*Multprec_Natural_Numbers.Exponent;
748: Add(expo,ns);
749: end if;
750: end if;
751: Adjust_Expo(res,expo);
752: Clear(n); Clear(expo);
753: end if;
754: return res;
755: end Trunc;
756:
757: function Trunc ( f : Floating_Number; k : natural ) return Floating_Number is
758:
759: res : Floating_Number;
760: cff : Array_of_Naturals(0..(Size(f.fraction)-k));
761: exp : integer;
762:
763: begin
764: Truncate(f.fraction,k,cff,exp);
765: res.fraction := Create(Create(cff));
766: if Negative(f.fraction)
767: then Min(res.fraction);
768: end if;
769: Copy(f.exponent,res.exponent);
770: Add(res.exponent,exp);
771: return res;
772: end Trunc;
773:
774: procedure Trunc ( f : in out Floating_Number; k : in natural ) is
775:
776: exp : integer;
777: n : Natural_Number;
778: neg : constant boolean := Negative(f.fraction);
779: ns : natural;
780:
781: begin
782: Copy(Unsigned(f.fraction),n);
783: Shift_Left(n,ns);
784: exp := k*Multprec_Natural_Numbers.Exponent - ns;
785: declare
786: cff : constant Array_of_Naturals := Coefficients(n);
787: trc : Array_of_Naturals(cff'first..cff'last-k);
788: begin
789: for i in trc'range loop
790: trc(i) := cff(i+k);
791: end loop;
792: Clear(f.fraction);
793: f.fraction := Create(trc);
794: if neg
795: then Min(f.fraction);
796: end if;
797: end;
798: Clear(n);
799: Add(f.exponent,exp);
800: end Trunc;
801:
802: procedure Set_Size ( f : in out Floating_Number; k : in natural ) is
803:
804: sz : constant natural := Size_Fraction(f);
805:
806: begin
807: if sz > k
808: then Round(f,sz-k);
809: elsif sz < k
810: then Expand(f,k-sz);
811: end if;
812: end Set_Size;
813:
814: -- ARITHMETIC OPERATIONS as functions (no data sharing) :
815:
816: function "+" ( f1 : Floating_Number; f2 : double_float )
817: return Floating_Number is
818:
819: ff2 : Floating_Number := Create(f2);
820: res : Floating_Number := f1+ff2;
821:
822: begin
823: Clear(ff2);
824: return res;
825: end "+";
826:
827: function "+" ( f1 : double_float; f2 : Floating_Number )
828: return Floating_Number is
829: begin
830: return f2+f1;
831: end "+";
832:
833: function "+" ( f1,f2 : Floating_Number ) return Floating_Number is
834:
835: res : Floating_Number;
836: maxsz,diff_size,max_diff : integer;
837: diff_expo : Integer_Number;
838: cnt : natural;
839:
840: begin
841: if (Empty(f1.fraction) or else Equal(f1.fraction,0))
842: then Copy(f2,res);
843: elsif (Empty(f2.fraction) or else Equal(f2.fraction,0))
844: then Copy(f1,res);
845: else maxsz := Max_Size(f1.fraction,f2.fraction) + 1;
846: if Equal(f1.exponent,f2.exponent)
847: then res.fraction := f1.fraction + f2.fraction;
848: Copy(f1.exponent,res.exponent);
849: else diff_expo := f1.exponent - f2.exponent;
850: max_diff := maxsz*Multprec_Natural_Numbers.Exponent;
851: if diff_expo > 0
852: then Copy(f1.fraction,res.fraction);
853: if diff_expo < 2*max_diff
854: then cnt := Create(diff_expo);
855: for i in 1..cnt loop
856: Mul(res.fraction,10);
857: end loop;
858: Add(res.fraction,f2.fraction);
859: Copy(f2.exponent,res.exponent);
860: else Copy(f1.exponent,res.exponent);
861: end if;
862: else Copy(f2.fraction,res.fraction);
863: Min(diff_expo);
864: if diff_expo < 2*max_diff
865: then cnt := Create(diff_expo);
866: for i in 1..cnt loop
867: Mul(res.fraction,10);
868: end loop;
869: Add(res.fraction,f1.fraction);
870: Copy(f1.exponent,res.exponent);
871: else Copy(f2.exponent,res.exponent);
872: end if;
873: end if;
874: Clear(diff_expo);
875: end if;
876: diff_size := Size(res.fraction) + 1 - maxsz;
877: if diff_size > 0
878: then Round(res,diff_size);
879: else Normalize(res);
880: end if;
881: end if;
882: return res;
883: end "+";
884:
885: function "+" ( f : Floating_Number ) return Floating_Number is
886:
887: res : Floating_Number;
888:
889: begin
890: Copy(f,res);
891: return res;
892: end "+";
893:
894: function "-" ( f : Floating_Number ) return Floating_Number is
895:
896: res : Floating_Number;
897:
898: begin
899: Copy(f,res);
900: Min(res.fraction);
901: return res;
902: end "-";
903:
904: function "-" ( f1 : Floating_Number; f2 : double_float )
905: return Floating_Number is
906:
907: minf2 : double_float := -f2;
908: res : Floating_Number := f1+minf2;
909:
910: begin
911: return res;
912: end "-";
913:
914: function "-" ( f1 : double_float; f2 : Floating_Number )
915: return Floating_Number is
916:
917: ff1 : Floating_Number := Create(f1);
918: res : Floating_Number := ff1 - f2;
919:
920: begin
921: Clear(ff1);
922: return res;
923: end "-";
924:
925: function "-" ( f1,f2 : Floating_Number ) return Floating_Number is
926:
927: res,minf2 : Floating_Number;
928:
929: begin
930: minf2.exponent := f2.exponent;
931: minf2.fraction := -f2.fraction;
932: res := f1 + minf2;
933: Clear(minf2.fraction);
934: return res;
935: end "-";
936:
937: function "*" ( f1 : Floating_Number; f2 : double_float )
938: return Floating_Number is
939:
940: ff2 : Floating_Number := Create(f2);
941: res : Floating_Number := f1*ff2;
942:
943: begin
944: Clear(ff2);
945: return res;
946: end "*";
947:
948: function "*" ( f1 : double_float; f2 : Floating_Number )
949: return Floating_Number is
950:
951: ff1 : Floating_Number := Create(f1);
952: res : Floating_Number := ff1*f2;
953:
954: begin
955: Clear(ff1);
956: return res;
957: end "*";
958:
959: function "*" ( f1,f2 : Floating_Number ) return Floating_Number is
960:
961: res : Floating_Number;
962: diffsize : integer;
963:
964: begin
965: if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
966: then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
967: then res.fraction := f1.fraction*f2.fraction;
968: res.exponent := f1.exponent+f2.exponent;
969: diffsize := Size(res.fraction)
970: - Max_Size(f1.fraction,f2.fraction);
971: if diffsize > 0
972: then Round(res,diffsize);
973: end if;
974: Normalize(res);
975: end if;
976: end if;
977: return res;
978: end "*";
979:
980: function "**" ( f : double_float; n : Natural_Number )
981: return Floating_Number is
982:
983: ff : Floating_Number := Create(f);
984: res : Floating_Number := ff**n;
985:
986: begin
987: Clear(ff);
988: return res;
989: end "**";
990:
991: function "**" ( f : Floating_Number; n : Natural_Number )
992: return Floating_Number is
993:
994: res : Floating_Number;
995: cnt : Natural_Number;
996:
997: begin
998: if (Empty(n) or else Equal(n,0))
999: then res := Create(1);
1000: else Copy(f,res);
1001: cnt := Create(1);
1002: while cnt < n loop
1003: Mul(res,f);
1004: Add(cnt,1);
1005: end loop;
1006: Clear(cnt);
1007: Normalize(res);
1008: end if;
1009: return res;
1010: end "**";
1011:
1012: function "**" ( f : Floating_Number; i : integer ) return Floating_Number is
1013:
1014: res : Floating_Number;
1015:
1016: begin
1017: if i = 0
1018: then res := Create(1.0);
1019: else if i > 0
1020: then Copy(f,res);
1021: for j in 1..(i-1) loop
1022: Mul(res,f);
1023: end loop;
1024: else res := Create(1);
1025: for j in 1..(-i) loop
1026: Div(res,f);
1027: end loop;
1028: end if;
1029: Normalize(res);
1030: end if;
1031: return res;
1032: end "**";
1033:
1034: function "**" ( f : double_float; i : Integer_Number )
1035: return Floating_Number is
1036:
1037: ff : Floating_Number := Create(f);
1038: res : Floating_Number := ff**i;
1039:
1040: begin
1041: Clear(ff);
1042: return res;
1043: end "**";
1044:
1045: function "**" ( f : Floating_Number; i : Integer_Number )
1046: return Floating_Number is
1047:
1048: res : Floating_Number;
1049: cnt,n : Natural_Number;
1050:
1051: begin
1052: if (Empty(i) or else Equal(i,0))
1053: then res := Create(1);
1054: else n := Unsigned(i);
1055: if i > 0
1056: then Copy(f,res);
1057: cnt := Create(1);
1058: while cnt < n loop
1059: Mul(res,f);
1060: Add(cnt,1);
1061: end loop;
1062: else res := Create(1);
1063: cnt := Create(0);
1064: while cnt < n loop
1065: Div(res,f);
1066: Add(cnt,1);
1067: end loop;
1068: end if;
1069: Clear(cnt);
1070: Normalize(res);
1071: end if;
1072: return res;
1073: end "**";
1074:
1075: function "/" ( f1 : Floating_Number; f2 : double_float )
1076: return Floating_Number is
1077:
1078: ff2 : Floating_Number := Create(f2);
1079: res : Floating_Number := f1/ff2;
1080:
1081: begin
1082: Clear(ff2);
1083: return res;
1084: end "/";
1085:
1086: function "/" ( f1 : double_float; f2 : Floating_Number )
1087: return Floating_Number is
1088:
1089: ff1 : Floating_Number := Create(f1);
1090: res : Floating_Number := ff1/f2;
1091:
1092: begin
1093: Clear(ff1);
1094: return res;
1095: end "/";
1096:
1097: function Pos_Div ( f1,f2 : Floating_Number ) return Floating_Number is
1098:
1099: -- DESCRIPTION :
1100: -- This division only applies when f1 and f2 are both positive.
1101:
1102: res : Floating_Number;
1103: rmd_frac,new_rmd_frac : Integer_Number;
1104: diffsize : integer;
1105: maxsz : natural := Max_Size(f1.fraction,f2.fraction);
1106:
1107: begin
1108: Div(f1.fraction,f2.fraction,res.fraction,rmd_frac);
1109: res.exponent := f1.exponent - f2.exponent;
1110: diffsize := Size(res.fraction) - maxsz;
1111: if not Equal(rmd_frac,0)
1112: then while diffsize <= 0 loop
1113: while rmd_frac < f2.fraction loop
1114: Mul(rmd_frac,10);
1115: Mul(res.fraction,10);
1116: Sub(res.exponent,1);
1117: end loop;
1118: Div(rmd_frac,f2.fraction,new_rmd_frac);
1119: Add(res.fraction,rmd_frac);
1120: Clear(rmd_frac); rmd_frac := new_rmd_frac;
1121: diffsize := Size(res.fraction) - maxsz;
1122: exit when Equal(rmd_frac,0);
1123: end loop;
1124: end if;
1125: Clear(rmd_frac);
1126: if diffsize > 0
1127: then Round(res,diffsize);
1128: end if;
1129: Normalize(res);
1130: return res;
1131: end Pos_Div;
1132:
1133: function "/" ( f1,f2 : Floating_Number ) return Floating_Number is
1134:
1135: res,minf1,minf2 : Floating_Number;
1136:
1137: begin
1138: if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
1139: then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
1140: then if Multprec_Integer_Numbers.Positive(f1.fraction)
1141: then if Multprec_Integer_Numbers.Positive(f2.fraction)
1142: then res := Pos_Div(f1,f2);
1143: else minf2.fraction := -f2.fraction;
1144: minf2.exponent := f2.exponent;
1145: res := Pos_Div(f1,minf2);
1146: Clear(minf2.fraction);
1147: Min(res);
1148: end if;
1149: else minf1.fraction := -f1.fraction;
1150: minf1.exponent := f1.exponent;
1151: if Multprec_Integer_Numbers.Positive(f2.fraction)
1152: then res := Pos_Div(minf1,f2);
1153: Min(res);
1154: else minf2.fraction := -f2.fraction;
1155: minf2.exponent := f2.exponent;
1156: res := Pos_Div(minf1,minf2);
1157: Clear(minf2.fraction);
1158: end if;
1159: Clear(minf1.fraction);
1160: end if;
1161: else raise NUMERIC_ERROR;
1162: end if;
1163: end if;
1164: return res;
1165: end "/";
1166:
1167: -- ARITHMETIC OPERATIONS as procedures for memory management :
1168:
1169: procedure Add ( f1 : in out Floating_Number; f2 : in double_float ) is
1170:
1171: ff2 : Floating_Number := Create(f2);
1172:
1173: begin
1174: Add(f1,ff2);
1175: Clear(ff2);
1176: end Add;
1177:
1178: procedure Add ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
1179:
1180: res : Floating_Number;
1181: maxsz,diff_size,max_diff : integer;
1182: diff_expo : Integer_Number;
1183: cnt : natural;
1184:
1185: begin
1186: if (Empty(f1.fraction) or else Equal(f1.fraction,0))
1187: then Copy(f2,f1);
1188: elsif (Empty(f2.fraction) or else Equal(f2.fraction,0))
1189: then null;
1190: else maxsz := Max_Size(f1.fraction,f2.fraction) + 1;
1191: if Equal(f1.exponent,f2.exponent)
1192: then Add(f1.fraction,f2.fraction);
1193: else diff_expo := f1.exponent - f2.exponent;
1194: max_diff := maxsz*Multprec_Natural_Numbers.Exponent;
1195: if diff_expo > 0
1196: then if diff_expo < 2*max_diff
1197: then cnt := Create(diff_expo);
1198: for i in 1..cnt loop
1199: Mul(f1.fraction,10);
1200: end loop;
1201: Add(f1.fraction,f2.fraction);
1202: Copy(f2.exponent,f1.exponent);
1203: else null;
1204: end if;
1205: else Copy(f2.fraction,res.fraction);
1206: Min(diff_expo);
1207: if diff_expo < 2*max_diff
1208: then cnt := Create(diff_expo);
1209: for i in 1..cnt loop
1210: Mul(res.fraction,10);
1211: end loop;
1212: Add(f1.fraction,res.fraction);
1213: Clear(res.fraction);
1214: else Copy(f2.exponent,f1.exponent);
1215: Clear(f1.fraction);
1216: f1.fraction := res.fraction;
1217: end if;
1218: end if;
1219: Clear(diff_expo);
1220: end if;
1221: diff_size := Size(f1.fraction) + 1 - maxsz;
1222: if diff_size > 0
1223: then Round(f1,diff_size);
1224: else Normalize(f1);
1225: end if;
1226: end if;
1227: end Add;
1228:
1229: procedure Sub ( f1 : in out Floating_Number; f2 : in double_float ) is
1230:
1231: minf2 : double_float := -f2;
1232:
1233: begin
1234: Add(f1,minf2);
1235: end Sub;
1236:
1237: procedure Sub ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
1238:
1239: minf2 : Floating_Number;
1240:
1241: begin
1242: minf2.exponent := f2.exponent;
1243: minf2.fraction := -f2.fraction;
1244: Add(f1,minf2);
1245: Clear(minf2.fraction);
1246: end Sub;
1247:
1248: procedure Min ( f : in out Floating_Number ) is
1249: begin
1250: Min(f.fraction);
1251: end Min;
1252:
1253: procedure Mul ( f1 : in out Floating_Number; f2 : in double_float ) is
1254:
1255: ff2 : Floating_Number := Create(f2);
1256:
1257: begin
1258: Mul(f1,ff2);
1259: Clear(ff2);
1260: end Mul;
1261:
1262: procedure Mul ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
1263:
1264: diffsize,maxsize : integer;
1265:
1266: begin
1267: if Empty(f2.fraction)
1268: then Clear(f1);
1269: elsif Equal(f2.fraction,0)
1270: then Clear(f1);
1271: elsif Empty(f1.fraction)
1272: then null;
1273: elsif Equal(f1.fraction,0)
1274: then null;
1275: else maxsize := Max_Size(f1.fraction,f2.fraction);
1276: Mul(f1.fraction,f2.fraction);
1277: Add(f1.exponent,f2.exponent);
1278: diffsize := Size(f1.fraction) - maxsize;
1279: if diffsize > 0
1280: then Round(f1,diffsize);
1281: end if;
1282: Normalize(f1);
1283: end if;
1284: end Mul;
1285:
1286: procedure Div ( f1 : in out Floating_Number; f2 : in double_float ) is
1287:
1288: ff2 : Floating_Number := Create(f2);
1289:
1290: begin
1291: Div(f1,ff2);
1292: Clear(ff2);
1293: end Div;
1294:
1295: procedure Pos_Div ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
1296:
1297: -- DESCRIPTION :
1298: -- Does the division in case f1 and f2 are both positive.
1299:
1300: rmd_frac,new_rmd_frac : Integer_Number;
1301: diffsize : integer;
1302: maxsz : natural := Max_Size(f1.fraction,f2.fraction);
1303:
1304: begin
1305: Div(f1.fraction,f2.fraction,rmd_frac);
1306: Sub(f1.exponent,f2.exponent);
1307: diffsize := Size(f1.fraction) - maxsz;
1308: if not Equal(rmd_frac,0)
1309: then while diffsize <= 0 loop
1310: while rmd_frac < f2.fraction loop
1311: Mul(rmd_frac,10);
1312: Mul(f1.fraction,10);
1313: Sub(f1.exponent,1);
1314: end loop;
1315: Div(rmd_frac,f2.fraction,new_rmd_frac);
1316: Add(f1.fraction,rmd_frac);
1317: Clear(rmd_frac); rmd_frac := new_rmd_frac;
1318: diffsize := Size(f1.fraction) - maxsz;
1319: exit when Equal(rmd_frac,0);
1320: end loop;
1321: end if;
1322: Clear(rmd_frac);
1323: if diffsize > 0
1324: then Round(f1,diffsize);
1325: end if;
1326: Normalize(f1);
1327: end Pos_Div;
1328:
1329: procedure Div ( f1 : in out Floating_Number; f2 : in Floating_Number ) is
1330:
1331: minf1,minf2 : Floating_Number;
1332:
1333: begin
1334: if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
1335: then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
1336: then if Multprec_Integer_Numbers.Positive(f1.fraction)
1337: then if Multprec_Integer_Numbers.Positive(f2.fraction)
1338: then Pos_Div(f1,f2);
1339: else minf2.fraction := -f2.fraction;
1340: minf2.exponent := f2.exponent;
1341: Pos_Div(f1,minf2);
1342: Clear(minf2.fraction); Min(f1);
1343: end if;
1344: else minf1 := -f1;
1345: if Multprec_Integer_Numbers.Positive(f2.fraction)
1346: then Pos_Div(minf1,f2);
1347: Clear(f1); f1 := minf1; Min(f1);
1348: else minf2.fraction := -f2.fraction;
1349: minf2.exponent := f2.exponent;
1350: Pos_Div(minf1,minf2);
1351: Clear(minf2.fraction);
1352: Clear(f1); f1 := minf1;
1353: end if;
1354: end if;
1355: else raise NUMERIC_ERROR;
1356: end if;
1357: end if;
1358: end Div;
1359:
1360: -- DESTRUCTOR :
1361:
1362: procedure Clear ( f : in out Floating_Number ) is
1363: begin
1364: Clear(f.fraction);
1365: Clear(f.exponent);
1366: end Clear;
1367:
1368: end Multprec_Floating_Numbers;