← Back to branch summary

~john-koepi/ubuntu/trusty/golang/default

« back to all changes in this revision

Viewing changes to src/cmd/gc/mparith3.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondřej Surý
  • Date: 2011-04-20 17:36:48 UTC
  • Revision ID: james.westby@ubuntu.com-20110420173648-ifergoxyrm832trd
Tags: upstream-2011.03.07.1
Import upstream version 2011.03.07.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/cmd/gc/mparith3.c
 
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
}