1 // Copyright 2009 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
4
5 #include "go.h"
6
7 /*
8 * returns the leading non-zero
9 * word of the number
10 */
11 int
12 sigfig(Mpflt *a)
13 {
14 int i;
15
16 for(i=Mpprec-1; i>=0; i--)
17 if(a->val.a[i] != 0)
18 break;
19 //print("sigfig %d %d\n", i-z+1, z);
20 return i+1;
21 }
22
23 /*
24 * shifts the leading non-zero
25 * word of the number to Mpnorm
26 */
27 void
28 mpnorm(Mpflt *a)
29 {
30 int s, os;
31 long x;
32
33 os = sigfig(a);
34 if(os == 0) {
35 // zero
36 a->exp = 0;
37 a->val.neg = 0;
38 return;
39 }
40
41 // this will normalize to the nearest word
42 x = a->val.a[os-1];
43 s = (Mpnorm-os) * Mpscale;
44
45 // further normalize to the nearest bit
46 for(;;) {
47 x <<= 1;
48 if(x & Mpbase)
49 break;
50 s++;
51 if(x == 0) {
52 // this error comes from trying to
53 // convert an Inf or something
54 // where the initial x=0x80000000
55 s = (Mpnorm-os) * Mpscale;
56 break;
57 }
58 }
59
60 mpshiftfix(&a->val, s);
61 a->exp -= s;
62 }
63
64 /// implements float arihmetic
65
66 void
67 mpaddfltflt(Mpflt *a, Mpflt *b)
68 {
69 int sa, sb, s;
70 Mpflt c;
71
72 if(Mpdebug)
73 print("\n%F + %F", a, b);
74
75 sa = sigfig(a);
76 if(sa == 0) {
77 mpmovefltflt(a, b);
78 goto out;
79 }
80
81 sb = sigfig(b);
82 if(sb == 0)
83 goto out;
84
85 s = a->exp - b->exp;
86 if(s > 0) {
87 // a is larger, shift b right
88 mpmovefltflt(&c, b);
89 mpshiftfix(&c.val, -s);
90 mpaddfixfix(&a->val, &c.val);
91 goto out;
92 }
93 if(s < 0) {
94 // b is larger, shift a right
95 mpshiftfix(&a->val, s);
96 a->exp -= s;
97 mpaddfixfix(&a->val, &b->val);
98 goto out;
99 }
100 mpaddfixfix(&a->val, &b->val);
101
102 out:
103 mpnorm(a);
104 if(Mpdebug)
105 print(" = %F\n\n", a);
106 }
107
108 void
109 mpmulfltflt(Mpflt *a, Mpflt *b)
110 {
111 int sa, sb;
112
113 if(Mpdebug)
114 print("%F\n * %F\n", a, b);
115
116 sa = sigfig(a);
117 if(sa == 0) {
118 // zero
119 a->exp = 0;
120 a->val.neg = 0;
121 return;
122 }
123
124 sb = sigfig(b);
125 if(sb == 0) {
126 // zero
127 mpmovefltflt(a, b);
128 return;
129 }
130
131 mpmulfract(&a->val, &b->val);
132 a->exp = (a->exp + b->exp) + Mpscale*Mpprec - Mpscale - 1;
133
134 mpnorm(a);
135 if(Mpdebug)
136 print(" = %F\n\n", a);
137 }
138
139 void
140 mpdivfltflt(Mpflt *a, Mpflt *b)
141 {
142 int sa, sb;
143 Mpflt c;
144
145 if(Mpdebug)
146 print("%F\n / %F\n", a, b);
147
148 sb = sigfig(b);
149 if(sb == 0) {
150 // zero and ovfl
151 a->exp = 0;
152 a->val.neg = 0;
153 a->val.ovf = 1;
154 yyerror("mpdivfltflt divide by zero");
155 return;
156 }
157
158 sa = sigfig(a);
159 if(sa == 0) {
160 // zero
161 a->exp = 0;
162 a->val.neg = 0;
163 return;
164 }
165
166 // adjust b to top
167 mpmovefltflt(&c, b);
168 mpshiftfix(&c.val, Mpscale);
169
170 // divide
171 mpdivfract(&a->val, &c.val);
172 a->exp = (a->exp-c.exp) - Mpscale*(Mpprec-1) + 1;
173
174 mpnorm(a);
175 if(Mpdebug)
176 print(" = %F\n\n", a);
177 }
178
179 double
180 mpgetflt(Mpflt *a)
181 {
182 int s, i, e;
183 uvlong v, vm;
184 double f;
185
186 if(a->val.ovf)
187 yyerror("mpgetflt ovf");
188
189 s = sigfig(a);
190 if(s == 0)
191 return 0;
192
193 if(s != Mpnorm) {
194 yyerror("mpgetflt norm");
195 mpnorm(a);
196 }
197
198 while((a->val.a[Mpnorm-1] & Mpsign) == 0) {
199 mpshiftfix(&a->val, 1);
200 a->exp -= 1;
201 }
202
203 // the magic numbers (64, 63, 53, 10, -1074) are
204 // IEEE specific. this should be done machine
205 // independently or in the 6g half of the compiler
206
207 // pick up the mantissa and a rounding bit in a uvlong
208 s = 53+1;
209 v = 0;
210 for(i=Mpnorm-1; s>=Mpscale; i--) {
211 v = (v<<Mpscale) | a->val.a[i];
212 s -= Mpscale;
213 }
214 vm = v;
215 if(s > 0)
216 vm = (vm<<s) | (a->val.a[i]>>(Mpscale-s));
217
218 // continue with 64 more bits
219 s += 64;
220 for(; s>=Mpscale; i--) {
221 v = (v<<Mpscale) | a->val.a[i];
222 s -= Mpscale;
223 }
224 if(s > 0)
225 v = (v<<s) | (a->val.a[i]>>(Mpscale-s));
226
227 // gradual underflow
228 e = Mpnorm*Mpscale + a->exp - 53;
229 if(e < -1074) {
230 s = -e - 1074;
231 if(s > 54)
232 s = 54;
233 v |= vm & ((1ULL<<s) - 1);
234 vm >>= s;
235 e = -1074;
236 }
237
238 //print("vm=%.16llux v=%.16llux\n", vm, v);
239 // round toward even
240 if(v != 0 || (vm&2ULL) != 0)
241 vm = (vm>>1) + (vm&1ULL);
242 else
243 vm >>= 1;
244
245 f = (double)(vm);
246 f = ldexp(f, e);
247
248 if(a->val.neg)
249 f = -f;
250 return f;
251 }
252
253 void
254 mpmovecflt(Mpflt *a, double c)
255 {
256 int i;
257 double f;
258 long l;
259
260 if(Mpdebug)
261 print("\nconst %g", c);
262 mpmovecfix(&a->val, 0);
263 a->exp = 0;
264 if(c == 0)
265 goto out;
266 if(c < 0) {
267 a->val.neg = 1;
268 c = -c;
269 }
270
271 f = frexp(c, &i);
272 a->exp = i;
273
274 for(i=0; i<10; i++) {
275 f = f*Mpbase;
276 l = floor(f);
277 f = f - l;
278 a->exp -= Mpscale;
279 a->val.a[0] = l;
280 if(f == 0)
281 break;
282 mpshiftfix(&a->val, Mpscale);
283 }
284
285 out:
286 mpnorm(a);
287 if(Mpdebug)
288 print(" = %F\n", a);
289 }
290
291 void
292 mpnegflt(Mpflt *a)
293 {
294 a->val.neg ^= 1;
295 }
296
297 int
298 mptestflt(Mpflt *a)
299 {
300 int s;
301
302 if(Mpdebug)
303 print("\n%F?", a);
304 s = sigfig(a);
305 if(s != 0) {
306 s = +1;
307 if(a->val.neg)
308 s = -1;
309 }
310 if(Mpdebug)
311 print(" = %d\n", s);
312 return s;
313 }