12 #include "cl_LF_impl.h"
15 #include "cl_xmacros.h"
19 const cl_LF LF_LF_plus_LF (const cl_LF& arg1, const cl_LF& arg2)
21 // Methode (nach [Knuth, II, Seminumerical Algorithms, Abschnitt 4.2.1., S.200]):
22 // Falls e1<e2, vertausche x1 und x2.
24 // Falls e2=0, also x2=0.0, Ergebnis x1.
25 // Falls e1 - e2 >= 16n+2, Ergebnis x1.
26 // Erweitere die Mantissen rechts um 3 Bits (Bit -1 als Schutzbit, Bits -2,-3
27 // als Rundungsbits: 00 exakt, 01 1.Hälfte, 10 exakte Mitte, 11 2.Hälfte.)
28 // Schiebe die Mantisse von x2 um e0-e1 Bits nach rechts. (Dabei die Rundung
29 // ausführen: Bit -3 ist das logische Oder der Bits -3,-4,-5,...)
30 // Falls x1,x2 selbes Vorzeichen haben: Addiere dieses zur Mantisse von x1.
31 // Falls x1,x2 verschiedenes Vorzeichen haben: Subtrahiere dieses von der
32 // Mantisse von x1. <0 -> (Es war e1=e2) Vertausche die Vorzeichen, negiere.
35 // Normalisiere, fertig.
38 var uintL uexp1 = TheLfloat(arg1)->expo;
39 var uintL uexp2 = TheLfloat(arg2)->expo;
41 // x1 und x2 vertauschen
42 { x1 = arg2; x2 = arg1; swap(uintL, uexp1,uexp2); }
44 if (uexp2==0) { return x1; } // x2=0.0 -> x1 als Ergebnis
45 var uintC len = TheLfloat(x1)->len; // Länge n von x1 und x2
46 var uintL expdiff = uexp1-uexp2; // e1-e2
47 if ((expdiff == 0) && (TheLfloat(x1)->sign != TheLfloat(x2)->sign))
48 // verschiedene Vorzeichen, aber gleicher Exponent
49 { // Vorzeichen des Ergebnisses festlegen:
50 var cl_signean erg = // Mantissen (je len Digits) vergleichen
51 compare_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len),arrayMSDptr(TheLfloat(x2)->data,len),len);
52 if (erg==0) // Mantissen gleich
53 { return encode_LF0(len); } // Ergebnis 0.0
54 if (erg<0) // |x1| < |x2|
55 // x1 und x2 vertauschen, expdiff bleibt =0
56 { x1.pointer = arg2.pointer; x2.pointer = arg1.pointer;
57 swap(uintL, uexp1,uexp2);
60 if (expdiff >= intDsize * (uintL)len + 2) // e1-e2 >= 16n+2 ?
61 { return x1; } // ja -> x1 als Ergebnis
62 // neues Long-Float allozieren:
63 var Lfloat y = allocate_lfloat(len,uexp1,TheLfloat(x1)->sign);
64 var uintL i = floor(expdiff,intDsize); // e1-e2 div 16 (>=0, <=n)
65 var uintL j = expdiff % intDsize; // e1-e2 mod 16 (>=0, <16)
66 // Mantisse von x2 muß um intDsize*i+j Bits nach rechts geschoben werden.
67 var uintC x2_len = len - i; // n-i Digits von x2 gebraucht
68 // x2_len Digits um j Bits nach rechts schieben und dabei kopieren:
72 var uintD rounding_bits;
73 num_stack_alloc(x2_len, x2_MSDptr=,x2_LSDptr=); // x2_len Digits Platz
75 { copy_loop_msp(arrayMSDptr(TheLfloat(x2)->data,len),x2_MSDptr,x2_len); rounding_bits = 0; }
77 { rounding_bits = shiftrightcopy_loop_msp(arrayMSDptr(TheLfloat(x2)->data,len),x2_MSDptr,x2_len,j,0); }
78 // x2_MSDptr/x2_len/x2_LSDptr sind die essentiellen Digits von x2.
79 // rounding_bits enthält die letzten j herausgeschobenen Bits.
80 // Aus rounding_bits und den nächsten i Digits die 3 Rundungsbits
81 // (als Bits intDsize-1..intDsize-3 von rounding_bits) aufbauen:
83 // j>=2 -> Bits -1,-2 sind OK, Bit -3 bestimmen:
84 { if ((rounding_bits & (bit(intDsize-3)-1)) ==0)
85 { if (test_loop_msp(arrayMSDptr(TheLfloat(x2)->data,len) mspop x2_len,i))
86 { rounding_bits |= bit(intDsize-3); } // Rundungsbit -3 setzen
89 { rounding_bits |= bit(intDsize-3); // Rundungsbit -3 setzen
90 rounding_bits &= bitm(intDsize)-bit(intDsize-3); // andere Bits löschen
93 // j<=3 -> Bits intDsize-4..0 von rounding_bits sind bereits Null.
94 // nächstes und weitere i-1 Digits heranziehen:
95 { if (i > 0) // i=0 -> Bits -1,-2,-3 sind OK.
96 { var uintD* ptr = arrayMSDptr(TheLfloat(x2)->data,len) mspop x2_len;
97 rounding_bits |= (mspref(ptr,0) >> j); // weitere relevante Bits des nächsten Digit dazu
98 if ((rounding_bits & (bit(intDsize-3)-1)) ==0) // Alle Bits -3,-4,... =0 ?
99 { if ( (!((mspref(ptr,0) & (bit(3)-1)) ==0)) // j (<=3) untere Bits von ptr[0] alle =0 ?
100 || test_loop_msp(ptr mspop 1,i-1)
102 { rounding_bits |= bit(intDsize-3); } // Rundungsbit -3 setzen
105 { rounding_bits |= bit(intDsize-3); // Rundungsbit -3 setzen
106 rounding_bits &= bitm(intDsize)-bit(intDsize-3); // andere Bits löschen
108 // x2 liegt in verschobener Form in der UDS x2_MSDptr/x2_len/x2_LSDptr
109 // vor, mit Rundungsbits in Bit intDsize-1..intDsize-3 von rounding_bits.
110 {var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
111 var uintD* y_mantLSDptr = arrayLSDptr(TheLfloat(y)->data,len);
112 if (TheLfloat(x1)->sign == TheLfloat(x2)->sign)
113 // gleiche Vorzeichen -> Mantissen addieren
114 { // erst rechten Mantissenteil (x2_len Digits) durch Addition:
116 add_loop_lsp(arrayLSDptr(TheLfloat(x1)->data,len),x2_LSDptr,
119 // dann linken Mantissenteil (i Digits) direkt kopieren:
121 copy_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len),y_mantMSDptr,i);
122 // dann Übertrag vom rechten zum linken Mantissenteil addieren:
124 { if ( inc_loop_lsp(ptr,i) )
125 // Übertrag über das erste Digit hinaus
126 { // Exponent von y incrementieren:
127 if ( ++(TheLfloat(y)->expo) == LF_exp_high+1 ) { cl_error_floating_point_overflow(); }
128 // normalisiere durch Schieben um 1 Bit nach rechts:
129 {var uintD carry_rechts =
130 shift1right_loop_msp(y_mantMSDptr,len,~(uintD)0);
131 rounding_bits = rounding_bits>>1; // Rundungsbits mitschieben
132 if (!(carry_rechts==0)) { rounding_bits |= bit(intDsize-1); }
136 // verschiedene Vorzeichen -> Mantissen subtrahieren
137 { // erst rechten Mantissenteil (x2_len Digits) durch Subtraktion:
138 rounding_bits = -rounding_bits;
140 subx_loop_lsp(arrayLSDptr(TheLfloat(x1)->data,len),x2_LSDptr,
141 y_mantLSDptr, x2_len,
142 (rounding_bits==0 ? 0 : ~(uintD)0)
144 // dann linken Mantissenteil (i Digits) direkt kopieren:
146 copy_loop_msp(arrayMSDptr(TheLfloat(x1)->data,len),y_mantMSDptr,i);
147 // dann Übertrag des rechten vom linken Mantissenteil subtrahieren:
149 { if ( dec_loop_lsp(ptr,i) )
150 // Übertrag über das erste Digit hinaus, also e1=e2
151 { NOTREACHED } // diesen Fall haben wir schon behandelt
154 // UDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits normalisieren:
155 {var uintD* ptr = y_mantMSDptr;
159 { if (!(mspref(ptr,0)==0)) goto nonzero_found;
160 ptr = ptr mspop 1; k++;
162 if (!(rounding_bits==0)) goto nonzero_found;
163 // Die UDS ist ganz Null. Also war e1=e2, keine Rundungsbits.
164 { NOTREACHED } // diesen Fall haben wir schon behandelt
165 nonzero_found: // Digit /=0 gefunden
166 // UDS von ptr nach y_mantMSDptr um k Digits nach unten kopieren:
168 // mindestens ein führendes Nulldigit. Also war e1-e2 = 0 oder 1.
169 { ptr = copy_loop_msp(ptr,y_mantMSDptr,len-k); // len-k Digits verschieben
170 msprefnext(ptr) = rounding_bits; // Rundungsbits als weiteres Digit
171 clear_loop_msp(ptr,k-1); // dann k-1 Nulldigits
172 rounding_bits = 0; // und keine weiteren Rundungsbits
173 // Exponenten um intDsize*k erniedrigen:
175 {var uintL uexp = TheLfloat(y)->expo;
177 if (uexp < k+LF_exp_low)
181 { if (underflow_allowed())
182 { cl_error_floating_point_underflow(); }
184 { return encode_LF0(len); } // Ergebnis 0.0
186 TheLfloat(y)->expo = uexp - k;
189 // NUDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits normalisieren:
191 integerlengthD(mspref(y_mantMSDptr,0), s = intDsize - );
192 // s = Anzahl der führenden Nullbits im ersten Word (>=0, <intDsize)
194 { // Muß die NUDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits
195 // um s Bits nach links schieben.
196 // (Bei e1-e2>1 ist dabei zwangsläufig s=1.)
198 { shift1left_loop_lsp(y_mantLSDptr,len);
199 if (rounding_bits & bit(intDsize-1))
200 { lspref(y_mantLSDptr,0) |= bit(0); }
201 rounding_bits = rounding_bits << 1;
203 else // s>1, also e1-e2 <= 1 <= s.
204 { shiftleft_loop_lsp(y_mantLSDptr,len,s,rounding_bits>>(intDsize-s));
205 rounding_bits = 0; // = rounding_bits << s;
207 // Exponenten um s erniedrigen:
208 {var uintL uexp = TheLfloat(y)->expo;
210 if (uexp < s+LF_exp_low)
214 { if (underflow_allowed())
215 { cl_error_floating_point_underflow(); }
217 { return encode_LF0(len); } // Ergebnis 0.0
219 TheLfloat(y)->expo = uexp - s;
222 // Hier enthält rounding_bits Bit -1 als Bit intDsize-1, Bit -2 als
223 // Bit intDsize-2, Bit -3 als Oder(Bits intDsize-3..0) !
224 // Runden. Dazu rounding_bits inspizieren:
225 if ((rounding_bits & bit(intDsize-1)) ==0) goto ab; // Bit -1 gelöscht -> abrunden
226 rounding_bits = rounding_bits<<1; // Bits -2,-3
227 if (!(rounding_bits==0)) goto auf; // Bit -2 oder Bit -3 gesetzt -> aufrunden
229 if ((lspref(y_mantLSDptr,0) & bit(0)) ==0) goto ab;
231 if ( inc_loop_lsp(y_mantLSDptr,len) )
232 { // Übertrag durchs Aufrunden
233 mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
235 if (++(TheLfloat(y)->expo) == LF_exp_high+1) { cl_error_floating_point_overflow(); }