Subversion Repositories HelenOS-historic

Rev

Rev 737 | Rev 828 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 737 Rev 804
Line 34... Line 34...
34
 *
34
 *
35
 */
35
 */
36
float32 mulFloat32(float32 a, float32 b)
36
float32 mulFloat32(float32 a, float32 b)
37
{
37
{
38
    float32 result;
38
    float32 result;
39
    __u64 mant1, mant2;
39
    __u64 frac1, frac2;
40
    __s32 exp;
40
    __s32 exp;
41
 
41
 
42
    result.parts.sign = a.parts.sign ^ b.parts.sign;
42
    result.parts.sign = a.parts.sign ^ b.parts.sign;
43
   
43
   
44
    if (isFloat32NaN(a) || isFloat32NaN(b) ) {
44
    if (isFloat32NaN(a) || isFloat32NaN(b) ) {
45
        /* TODO: fix SigNaNs */
45
        /* TODO: fix SigNaNs */
46
        if (isFloat32SigNaN(a)) {
46
        if (isFloat32SigNaN(a)) {
47
            result.parts.mantisa = a.parts.mantisa;
47
            result.parts.fraction = a.parts.fraction;
48
            result.parts.exp = a.parts.exp;
48
            result.parts.exp = a.parts.exp;
49
            return result;
49
            return result;
50
        };
50
        };
51
        if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
51
        if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
52
            result.parts.mantisa = b.parts.mantisa;
52
            result.parts.fraction = b.parts.fraction;
53
            result.parts.exp = b.parts.exp;
53
            result.parts.exp = b.parts.exp;
54
            return result;
54
            return result;
55
        };
55
        };
56
        /* set NaN as result */
56
        /* set NaN as result */
57
        result.binary = FLOAT32_NAN;
57
        result.binary = FLOAT32_NAN;
Line 62... Line 62...
62
        if (isFloat32Zero(b)) {
62
        if (isFloat32Zero(b)) {
63
            /* FIXME: zero * infinity */
63
            /* FIXME: zero * infinity */
64
            result.binary = FLOAT32_NAN;
64
            result.binary = FLOAT32_NAN;
65
            return result;
65
            return result;
66
        }
66
        }
67
        result.parts.mantisa = a.parts.mantisa;
67
        result.parts.fraction = a.parts.fraction;
68
        result.parts.exp = a.parts.exp;
68
        result.parts.exp = a.parts.exp;
69
        return result;
69
        return result;
70
    }
70
    }
71
 
71
 
72
    if (isFloat32Infinity(b)) {
72
    if (isFloat32Infinity(b)) {
73
        if (isFloat32Zero(a)) {
73
        if (isFloat32Zero(a)) {
74
            /* FIXME: zero * infinity */
74
            /* FIXME: zero * infinity */
75
            result.binary = FLOAT32_NAN;
75
            result.binary = FLOAT32_NAN;
76
            return result;
76
            return result;
77
        }
77
        }
78
        result.parts.mantisa = b.parts.mantisa;
78
        result.parts.fraction = b.parts.fraction;
79
        result.parts.exp = b.parts.exp;
79
        result.parts.exp = b.parts.exp;
80
        return result;
80
        return result;
81
    }
81
    }
82
 
82
 
83
    /* exp is signed so we can easy detect underflow */
83
    /* exp is signed so we can easy detect underflow */
Line 93... Line 93...
93
    };
93
    };
94
   
94
   
95
    if (exp < 0) {
95
    if (exp < 0) {
96
        /* FIXME: underflow */
96
        /* FIXME: underflow */
97
        /* return signed zero */
97
        /* return signed zero */
98
        result.parts.mantisa = 0x0;
98
        result.parts.fraction = 0x0;
99
        result.parts.exp = 0x0;
99
        result.parts.exp = 0x0;
100
        return result;
100
        return result;
101
    };
101
    };
102
   
102
   
103
    mant1 = a.parts.mantisa;
103
    frac1 = a.parts.fraction;
104
    if (a.parts.exp > 0) {
104
    if (a.parts.exp > 0) {
105
        mant1 |= FLOAT32_HIDDEN_BIT_MASK;
105
        frac1 |= FLOAT32_HIDDEN_BIT_MASK;
106
    } else {
106
    } else {
107
        ++exp;
107
        ++exp;
108
    };
108
    };
109
   
109
   
110
    mant2 = b.parts.mantisa;
110
    frac2 = b.parts.fraction;
111
 
111
 
112
    if (b.parts.exp > 0) {
112
    if (b.parts.exp > 0) {
113
        mant2 |= FLOAT32_HIDDEN_BIT_MASK;
113
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
114
    } else {
114
    } else {
115
        ++exp;
115
        ++exp;
116
    };
116
    };
117
 
117
 
118
    mant1 <<= 1; /* one bit space for rounding */
118
    frac1 <<= 1; /* one bit space for rounding */
119
 
119
 
120
    mant1 = mant1 * mant2;
120
    frac1 = frac1 * frac2;
121
/* round and return */
121
/* round and return */
122
   
122
   
123
    while ((exp < FLOAT32_MAX_EXPONENT) && (mant1 >= ( 1 << (FLOAT32_MANTISA_SIZE + 2)))) {
123
    while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
124
        /* 23 bits of mantisa + one more for hidden bit (all shifted 1 bit left)*/
124
        /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
125
        ++exp;
125
        ++exp;
126
        mant1 >>= 1;
126
        frac1 >>= 1;
127
    };
127
    };
128
 
128
 
129
    /* rounding */
129
    /* rounding */
130
    //++mant1; /* FIXME: not works - without it is ok */
130
    /* ++frac1; FIXME: not works - without it is ok */
131
    mant1 >>= 1; /* shift off rounding space */
131
    frac1 >>= 1; /* shift off rounding space */
132
   
132
   
133
    if ((exp < FLOAT32_MAX_EXPONENT) && (mant1 >= (1 << (FLOAT32_MANTISA_SIZE + 1)))) {
133
    if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
134
        ++exp;
134
        ++exp;
135
        mant1 >>= 1;
135
        frac1 >>= 1;
136
    };
136
    };
137
 
137
 
138
    if (exp >= FLOAT32_MAX_EXPONENT ) {
138
    if (exp >= FLOAT32_MAX_EXPONENT ) {
139
        /* TODO: fix overflow */
139
        /* TODO: fix overflow */
140
        /* return infinity*/
140
        /* return infinity*/
141
        result.parts.exp = FLOAT32_MAX_EXPONENT;
141
        result.parts.exp = FLOAT32_MAX_EXPONENT;
142
        result.parts.mantisa = 0x0;
142
        result.parts.fraction = 0x0;
143
        return result;
143
        return result;
144
    }
144
    }
145
   
145
   
146
    exp -= FLOAT32_MANTISA_SIZE;
146
    exp -= FLOAT32_FRACTION_SIZE;
147
 
147
 
148
    if (exp <= FLOAT32_MANTISA_SIZE) {
148
    if (exp <= FLOAT32_FRACTION_SIZE) {
149
        /* denormalized number */
149
        /* denormalized number */
150
        mant1 >>= 1; /* denormalize */
150
        frac1 >>= 1; /* denormalize */
151
        while ((mant1 > 0) && (exp < 0)) {
151
        while ((frac1 > 0) && (exp < 0)) {
152
            mant1 >>= 1;
152
            frac1 >>= 1;
153
            ++exp;
153
            ++exp;
154
        };
154
        };
155
        if (mant1 == 0) {
155
        if (frac1 == 0) {
156
            /* FIXME : underflow */
156
            /* FIXME : underflow */
157
        result.parts.exp = 0;
157
        result.parts.exp = 0;
158
        result.parts.mantisa = 0;
158
        result.parts.fraction = 0;
159
        return result;
159
        return result;
160
        };
160
        };
161
    };
161
    };
162
    result.parts.exp = exp;
162
    result.parts.exp = exp;
163
    result.parts.mantisa = mant1 & ( (1 << FLOAT32_MANTISA_SIZE) - 1);
163
    result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
164
   
164
   
165
    return result; 
165
    return result; 
166
   
166
   
167
}
167
}
168
 
168
 
Line 170... Line 170...
170
 *
170
 *
171
 */
171
 */
172
float64 mulFloat64(float64 a, float64 b)
172
float64 mulFloat64(float64 a, float64 b)
173
{
173
{
174
    float64 result;
174
    float64 result;
175
    __u64 mant1, mant2;
175
    __u64 frac1, frac2;
176
    __s32 exp;
176
    __s32 exp;
177
 
177
 
178
    result.parts.sign = a.parts.sign ^ b.parts.sign;
178
    result.parts.sign = a.parts.sign ^ b.parts.sign;
179
   
179
   
180
    if (isFloat64NaN(a) || isFloat64NaN(b) ) {
180
    if (isFloat64NaN(a) || isFloat64NaN(b) ) {
181
        /* TODO: fix SigNaNs */
181
        /* TODO: fix SigNaNs */
182
        if (isFloat64SigNaN(a)) {
182
        if (isFloat64SigNaN(a)) {
183
            result.parts.mantisa = a.parts.mantisa;
183
            result.parts.fraction = a.parts.fraction;
184
            result.parts.exp = a.parts.exp;
184
            result.parts.exp = a.parts.exp;
185
            return result;
185
            return result;
186
        };
186
        };
187
        if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
187
        if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
188
            result.parts.mantisa = b.parts.mantisa;
188
            result.parts.fraction = b.parts.fraction;
189
            result.parts.exp = b.parts.exp;
189
            result.parts.exp = b.parts.exp;
190
            return result;
190
            return result;
191
        };
191
        };
192
        /* set NaN as result */
192
        /* set NaN as result */
193
        result.binary = FLOAT64_NAN;
193
        result.binary = FLOAT64_NAN;
Line 198... Line 198...
198
        if (isFloat64Zero(b)) {
198
        if (isFloat64Zero(b)) {
199
            /* FIXME: zero * infinity */
199
            /* FIXME: zero * infinity */
200
            result.binary = FLOAT64_NAN;
200
            result.binary = FLOAT64_NAN;
201
            return result;
201
            return result;
202
        }
202
        }
203
        result.parts.mantisa = a.parts.mantisa;
203
        result.parts.fraction = a.parts.fraction;
204
        result.parts.exp = a.parts.exp;
204
        result.parts.exp = a.parts.exp;
205
        return result;
205
        return result;
206
    }
206
    }
207
 
207
 
208
    if (isFloat64Infinity(b)) {
208
    if (isFloat64Infinity(b)) {
209
        if (isFloat64Zero(a)) {
209
        if (isFloat64Zero(a)) {
210
            /* FIXME: zero * infinity */
210
            /* FIXME: zero * infinity */
211
            result.binary = FLOAT64_NAN;
211
            result.binary = FLOAT64_NAN;
212
            return result;
212
            return result;
213
        }
213
        }
214
        result.parts.mantisa = b.parts.mantisa;
214
        result.parts.fraction = b.parts.fraction;
215
        result.parts.exp = b.parts.exp;
215
        result.parts.exp = b.parts.exp;
216
        return result;
216
        return result;
217
    }
217
    }
218
 
218
 
219
    /* exp is signed so we can easy detect underflow */
219
    /* exp is signed so we can easy detect underflow */
Line 229... Line 229...
229
    };
229
    };
230
   
230
   
231
    if (exp < 0) {
231
    if (exp < 0) {
232
        /* FIXME: underflow */
232
        /* FIXME: underflow */
233
        /* return signed zero */
233
        /* return signed zero */
234
        result.parts.mantisa = 0x0;
234
        result.parts.fraction = 0x0;
235
        result.parts.exp = 0x0;
235
        result.parts.exp = 0x0;
236
        return result;
236
        return result;
237
    };
237
    };
238
   
238
   
239
    mant1 = a.parts.mantisa;
239
    frac1 = a.parts.fraction;
240
    if (a.parts.exp > 0) {
240
    if (a.parts.exp > 0) {
241
        mant1 |= FLOAT64_HIDDEN_BIT_MASK;
241
        frac1 |= FLOAT64_HIDDEN_BIT_MASK;
242
    } else {
242
    } else {
243
        ++exp;
243
        ++exp;
244
    };
244
    };
245
   
245
   
246
    mant2 = b.parts.mantisa;
246
    frac2 = b.parts.fraction;
247
 
247
 
248
    if (b.parts.exp > 0) {
248
    if (b.parts.exp > 0) {
249
        mant2 |= FLOAT64_HIDDEN_BIT_MASK;
249
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
250
    } else {
250
    } else {
251
        ++exp;
251
        ++exp;
252
    };
252
    };
253
 
253
 
254
    mant1 <<= 1; /* one bit space for rounding */
254
    frac1 <<= 1; /* one bit space for rounding */
255
 
255
 
256
    mul64integers(mant1, mant2, &mant1, &mant2);
256
    mul64integers(frac1, frac2, &frac1, &frac2);
257
 
257
 
258
/* round and return */
258
/* round and return */
259
    /* FIXME: ugly soulution is to shift whole mant2 >> as in 32bit version
259
    /* FIXME: ugly soulution is to shift whole frac2 >> as in 32bit version
260
     * Here is is more slower because we have to shift two numbers with carry
260
     * Here is is more slower because we have to shift two numbers with carry
261
     * Better is find first nonzero bit and make only one shift
261
     * Better is find first nonzero bit and make only one shift
262
     * Third version is to shift both numbers a bit to right and result will be then
262
     * Third version is to shift both numbers a bit to right and result will be then
263
     * placed in higher part of result. Then lower part will be good only for rounding.
263
     * placed in higher part of result. Then lower part will be good only for rounding.
264
     */
264
     */
265
   
265
   
266
    while ((exp < FLOAT64_MAX_EXPONENT) && (mant2 > 0 )) {
266
    while ((exp < FLOAT64_MAX_EXPONENT) && (frac2 > 0 )) {
267
        mant1 >>= 1;
267
        frac1 >>= 1;
268
        mant1 &= ((mant2 & 0x1) << 63);
268
        frac1 &= ((frac2 & 0x1) << 63);
269
        mant2 >>= 1;
269
        frac2 >>= 1;
270
        ++exp;
270
        ++exp;
271
    }
271
    }
272
   
272
   
273
    while ((exp < FLOAT64_MAX_EXPONENT) && (mant1 >= ( (__u64)1 << (FLOAT64_MANTISA_SIZE + 2)))) {
273
    while ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ( (__u64)1 << (FLOAT64_FRACTION_SIZE + 2)))) {
274
        ++exp;
274
        ++exp;
275
        mant1 >>= 1;
275
        frac1 >>= 1;
276
    };
276
    };
277
 
277
 
278
    /* rounding */
278
    /* rounding */
279
    //++mant1; /* FIXME: not works - without it is ok */
279
    /* ++frac1;  FIXME: not works - without it is ok */
280
    mant1 >>= 1; /* shift off rounding space */
280
    frac1 >>= 1; /* shift off rounding space */
281
   
281
   
282
    if ((exp < FLOAT64_MAX_EXPONENT) && (mant1 >= ((__u64)1 << (FLOAT64_MANTISA_SIZE + 1)))) {
282
    if ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ((__u64)1 << (FLOAT64_FRACTION_SIZE + 1)))) {
283
        ++exp;
283
        ++exp;
284
        mant1 >>= 1;
284
        frac1 >>= 1;
285
    };
285
    };
286
 
286
 
287
    if (exp >= FLOAT64_MAX_EXPONENT ) {
287
    if (exp >= FLOAT64_MAX_EXPONENT ) {
288
        /* TODO: fix overflow */
288
        /* TODO: fix overflow */
289
        /* return infinity*/
289
        /* return infinity*/
290
        result.parts.exp = FLOAT64_MAX_EXPONENT;
290
        result.parts.exp = FLOAT64_MAX_EXPONENT;
291
        result.parts.mantisa = 0x0;
291
        result.parts.fraction = 0x0;
292
        return result;
292
        return result;
293
    }
293
    }
294
   
294
   
295
    exp -= FLOAT64_MANTISA_SIZE;
295
    exp -= FLOAT64_FRACTION_SIZE;
296
 
296
 
297
    if (exp <= FLOAT64_MANTISA_SIZE) {
297
    if (exp <= FLOAT64_FRACTION_SIZE) {
298
        /* denormalized number */
298
        /* denormalized number */
299
        mant1 >>= 1; /* denormalize */
299
        frac1 >>= 1; /* denormalize */
300
        while ((mant1 > 0) && (exp < 0)) {
300
        while ((frac1 > 0) && (exp < 0)) {
301
            mant1 >>= 1;
301
            frac1 >>= 1;
302
            ++exp;
302
            ++exp;
303
        };
303
        };
304
        if (mant1 == 0) {
304
        if (frac1 == 0) {
305
            /* FIXME : underflow */
305
            /* FIXME : underflow */
306
        result.parts.exp = 0;
306
        result.parts.exp = 0;
307
        result.parts.mantisa = 0;
307
        result.parts.fraction = 0;
308
        return result;
308
        return result;
309
        };
309
        };
310
    };
310
    };
311
    result.parts.exp = exp;
311
    result.parts.exp = exp;
312
    result.parts.mantisa = mant1 & ( ((__u64)1 << FLOAT64_MANTISA_SIZE) - 1);
312
    result.parts.fraction = frac1 & ( ((__u64)1 << FLOAT64_FRACTION_SIZE) - 1);
313
   
313
   
314
    return result; 
314
    return result; 
315
   
315
   
316
}
316
}
317
 
317
 
Line 336... Line 336...
336
    middle1 = a * blow;
336
    middle1 = a * blow;
337
    middle2 = alow * b;
337
    middle2 = alow * b;
338
    high = a * b;
338
    high = a * b;
339
 
339
 
340
    middle1 += middle2;
340
    middle1 += middle2;
341
    high += ((__u64)(middle1 < middle2) << 32) + middle1>>32;
341
    high += ((__u64)(middle1 < middle2) << 32) + (middle1 >> 32);
342
    middle1 << 32;
342
    middle1 <<= 32;
343
    low += middle1;
343
    low += middle1;
344
    high += (low < middle1);
344
    high += (low < middle1);
345
    *lo = low;
345
    *lo = low;
346
    *hi = high;
346
    *hi = high;
347
    return;
347
    return;