Subversion Repositories HelenOS-historic

Rev

Rev 737 | Rev 829 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
731 cejka 1
/*
2
 * Copyright (C) 2005 Josef Cejka
3
 * All rights reserved.
4
 *
5
 * Redistribution and use in source and binary forms, with or without
6
 * modification, are permitted provided that the following conditions
7
 * are met:
8
 *
9
 * - Redistributions of source code must retain the above copyright
10
 *   notice, this list of conditions and the following disclaimer.
11
 * - Redistributions in binary form must reproduce the above copyright
12
 *   notice, this list of conditions and the following disclaimer in the
13
 *   documentation and/or other materials provided with the distribution.
14
 * - The name of the author may not be used to endorse or promote products
15
 *   derived from this software without specific prior written permission.
16
 *
17
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
18
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27
 */
28
 
29
#include<sftypes.h>
30
#include<mul.h>
31
#include<comparison.h>
32
 
33
/** Multiply two 32 bit float numbers
34
 *
35
 */
36
float32 mulFloat32(float32 a, float32 b)
37
{
38
    float32 result;
804 cejka 39
    __u64 frac1, frac2;
731 cejka 40
    __s32 exp;
41
 
42
    result.parts.sign = a.parts.sign ^ b.parts.sign;
43
 
737 cejka 44
    if (isFloat32NaN(a) || isFloat32NaN(b) ) {
731 cejka 45
        /* TODO: fix SigNaNs */
46
        if (isFloat32SigNaN(a)) {
804 cejka 47
            result.parts.fraction = a.parts.fraction;
731 cejka 48
            result.parts.exp = a.parts.exp;
49
            return result;
50
        };
51
        if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
804 cejka 52
            result.parts.fraction = b.parts.fraction;
731 cejka 53
            result.parts.exp = b.parts.exp;
54
            return result;
55
        };
56
        /* set NaN as result */
737 cejka 57
        result.binary = FLOAT32_NAN;
731 cejka 58
        return result;
59
    };
60
 
61
    if (isFloat32Infinity(a)) {
62
        if (isFloat32Zero(b)) {
63
            /* FIXME: zero * infinity */
737 cejka 64
            result.binary = FLOAT32_NAN;
731 cejka 65
            return result;
66
        }
804 cejka 67
        result.parts.fraction = a.parts.fraction;
731 cejka 68
        result.parts.exp = a.parts.exp;
69
        return result;
70
    }
71
 
72
    if (isFloat32Infinity(b)) {
73
        if (isFloat32Zero(a)) {
74
            /* FIXME: zero * infinity */
737 cejka 75
            result.binary = FLOAT32_NAN;
731 cejka 76
            return result;
77
        }
804 cejka 78
        result.parts.fraction = b.parts.fraction;
731 cejka 79
        result.parts.exp = b.parts.exp;
80
        return result;
81
    }
82
 
83
    /* exp is signed so we can easy detect underflow */
84
    exp = a.parts.exp + b.parts.exp;
85
    exp -= FLOAT32_BIAS;
86
 
737 cejka 87
    if (exp >= FLOAT32_MAX_EXPONENT) {
731 cejka 88
        /* FIXME: overflow */
89
        /* set infinity as result */
737 cejka 90
        result.binary = FLOAT32_INF;
91
        result.parts.sign = a.parts.sign ^ b.parts.sign;
731 cejka 92
        return result;
93
    };
94
 
95
    if (exp < 0) {
96
        /* FIXME: underflow */
97
        /* return signed zero */
804 cejka 98
        result.parts.fraction = 0x0;
731 cejka 99
        result.parts.exp = 0x0;
100
        return result;
101
    };
102
 
804 cejka 103
    frac1 = a.parts.fraction;
737 cejka 104
    if (a.parts.exp > 0) {
804 cejka 105
        frac1 |= FLOAT32_HIDDEN_BIT_MASK;
731 cejka 106
    } else {
107
        ++exp;
108
    };
109
 
804 cejka 110
    frac2 = b.parts.fraction;
737 cejka 111
 
112
    if (b.parts.exp > 0) {
804 cejka 113
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
731 cejka 114
    } else {
115
        ++exp;
116
    };
117
 
804 cejka 118
    frac1 <<= 1; /* one bit space for rounding */
731 cejka 119
 
804 cejka 120
    frac1 = frac1 * frac2;
731 cejka 121
/* round and return */
122
 
804 cejka 123
    while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
124
        /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
731 cejka 125
        ++exp;
804 cejka 126
        frac1 >>= 1;
731 cejka 127
    };
128
 
129
    /* rounding */
804 cejka 130
    /* ++frac1; FIXME: not works - without it is ok */
131
    frac1 >>= 1; /* shift off rounding space */
731 cejka 132
 
804 cejka 133
    if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
731 cejka 134
        ++exp;
804 cejka 135
        frac1 >>= 1;
731 cejka 136
    };
137
 
737 cejka 138
    if (exp >= FLOAT32_MAX_EXPONENT ) {
731 cejka 139
        /* TODO: fix overflow */
140
        /* return infinity*/
737 cejka 141
        result.parts.exp = FLOAT32_MAX_EXPONENT;
804 cejka 142
        result.parts.fraction = 0x0;
731 cejka 143
        return result;
144
    }
145
 
804 cejka 146
    exp -= FLOAT32_FRACTION_SIZE;
731 cejka 147
 
804 cejka 148
    if (exp <= FLOAT32_FRACTION_SIZE) {
731 cejka 149
        /* denormalized number */
804 cejka 150
        frac1 >>= 1; /* denormalize */
151
        while ((frac1 > 0) && (exp < 0)) {
152
            frac1 >>= 1;
731 cejka 153
            ++exp;
154
        };
804 cejka 155
        if (frac1 == 0) {
731 cejka 156
            /* FIXME : underflow */
157
        result.parts.exp = 0;
804 cejka 158
        result.parts.fraction = 0;
731 cejka 159
        return result;
160
        };
161
    };
162
    result.parts.exp = exp;
804 cejka 163
    result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
731 cejka 164
 
165
    return result; 
166
 
167
}
168
 
737 cejka 169
/** Multiply two 64 bit float numbers
170
 *
171
 */
172
float64 mulFloat64(float64 a, float64 b)
173
{
174
    float64 result;
804 cejka 175
    __u64 frac1, frac2;
737 cejka 176
    __s32 exp;
731 cejka 177
 
737 cejka 178
    result.parts.sign = a.parts.sign ^ b.parts.sign;
179
 
180
    if (isFloat64NaN(a) || isFloat64NaN(b) ) {
181
        /* TODO: fix SigNaNs */
182
        if (isFloat64SigNaN(a)) {
804 cejka 183
            result.parts.fraction = a.parts.fraction;
737 cejka 184
            result.parts.exp = a.parts.exp;
185
            return result;
186
        };
187
        if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
804 cejka 188
            result.parts.fraction = b.parts.fraction;
737 cejka 189
            result.parts.exp = b.parts.exp;
190
            return result;
191
        };
192
        /* set NaN as result */
193
        result.binary = FLOAT64_NAN;
194
        return result;
195
    };
196
 
197
    if (isFloat64Infinity(a)) {
198
        if (isFloat64Zero(b)) {
199
            /* FIXME: zero * infinity */
200
            result.binary = FLOAT64_NAN;
201
            return result;
202
        }
804 cejka 203
        result.parts.fraction = a.parts.fraction;
737 cejka 204
        result.parts.exp = a.parts.exp;
205
        return result;
206
    }
731 cejka 207
 
737 cejka 208
    if (isFloat64Infinity(b)) {
209
        if (isFloat64Zero(a)) {
210
            /* FIXME: zero * infinity */
211
            result.binary = FLOAT64_NAN;
212
            return result;
213
        }
804 cejka 214
        result.parts.fraction = b.parts.fraction;
737 cejka 215
        result.parts.exp = b.parts.exp;
216
        return result;
217
    }
218
 
219
    /* exp is signed so we can easy detect underflow */
220
    exp = a.parts.exp + b.parts.exp;
221
    exp -= FLOAT64_BIAS;
222
 
223
    if (exp >= FLOAT64_MAX_EXPONENT) {
224
        /* FIXME: overflow */
225
        /* set infinity as result */
226
        result.binary = FLOAT64_INF;
227
        result.parts.sign = a.parts.sign ^ b.parts.sign;
228
        return result;
229
    };
230
 
231
    if (exp < 0) {
232
        /* FIXME: underflow */
233
        /* return signed zero */
804 cejka 234
        result.parts.fraction = 0x0;
737 cejka 235
        result.parts.exp = 0x0;
236
        return result;
237
    };
238
 
804 cejka 239
    frac1 = a.parts.fraction;
737 cejka 240
    if (a.parts.exp > 0) {
804 cejka 241
        frac1 |= FLOAT64_HIDDEN_BIT_MASK;
737 cejka 242
    } else {
243
        ++exp;
244
    };
245
 
804 cejka 246
    frac2 = b.parts.fraction;
737 cejka 247
 
248
    if (b.parts.exp > 0) {
804 cejka 249
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
737 cejka 250
    } else {
251
        ++exp;
252
    };
253
 
804 cejka 254
    frac1 <<= 1; /* one bit space for rounding */
737 cejka 255
 
804 cejka 256
    mul64integers(frac1, frac2, &frac1, &frac2);
737 cejka 257
 
258
/* round and return */
804 cejka 259
    /* FIXME: ugly soulution is to shift whole frac2 >> as in 32bit version
737 cejka 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
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.
264
     */
265
 
804 cejka 266
    while ((exp < FLOAT64_MAX_EXPONENT) && (frac2 > 0 )) {
267
        frac1 >>= 1;
268
        frac1 &= ((frac2 & 0x1) << 63);
269
        frac2 >>= 1;
737 cejka 270
        ++exp;
271
    }
272
 
804 cejka 273
    while ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ( (__u64)1 << (FLOAT64_FRACTION_SIZE + 2)))) {
737 cejka 274
        ++exp;
804 cejka 275
        frac1 >>= 1;
737 cejka 276
    };
277
 
278
    /* rounding */
804 cejka 279
    /* ++frac1;  FIXME: not works - without it is ok */
280
    frac1 >>= 1; /* shift off rounding space */
737 cejka 281
 
804 cejka 282
    if ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ((__u64)1 << (FLOAT64_FRACTION_SIZE + 1)))) {
737 cejka 283
        ++exp;
804 cejka 284
        frac1 >>= 1;
737 cejka 285
    };
286
 
287
    if (exp >= FLOAT64_MAX_EXPONENT ) {
288
        /* TODO: fix overflow */
289
        /* return infinity*/
290
        result.parts.exp = FLOAT64_MAX_EXPONENT;
804 cejka 291
        result.parts.fraction = 0x0;
737 cejka 292
        return result;
293
    }
294
 
804 cejka 295
    exp -= FLOAT64_FRACTION_SIZE;
737 cejka 296
 
804 cejka 297
    if (exp <= FLOAT64_FRACTION_SIZE) {
737 cejka 298
        /* denormalized number */
804 cejka 299
        frac1 >>= 1; /* denormalize */
300
        while ((frac1 > 0) && (exp < 0)) {
301
            frac1 >>= 1;
737 cejka 302
            ++exp;
303
        };
804 cejka 304
        if (frac1 == 0) {
737 cejka 305
            /* FIXME : underflow */
306
        result.parts.exp = 0;
804 cejka 307
        result.parts.fraction = 0;
737 cejka 308
        return result;
309
        };
310
    };
311
    result.parts.exp = exp;
804 cejka 312
    result.parts.fraction = frac1 & ( ((__u64)1 << FLOAT64_FRACTION_SIZE) - 1);
737 cejka 313
 
314
    return result; 
315
 
316
}
317
 
318
/** Multiply two 64 bit numbers and return result in two parts
319
 * @param a first operand
320
 * @param b second operand
321
 * @param lo lower part from result
322
 * @param hi higher part of result
323
 */
324
void mul64integers(__u64 a,__u64 b, __u64 *lo, __u64 *hi)
325
{
326
    __u64 low, high, middle1, middle2;
327
    __u32 alow, blow;
328
 
329
    alow = a & 0xFFFFFFFF;
330
    blow = b & 0xFFFFFFFF;
331
 
332
    a <<= 32;
333
    b <<= 32;
334
 
335
    low = (__u64)alow * blow;
336
    middle1 = a * blow;
337
    middle2 = alow * b;
338
    high = a * b;
339
 
340
    middle1 += middle2;
804 cejka 341
    high += ((__u64)(middle1 < middle2) << 32) + (middle1 >> 32);
342
    middle1 <<= 32;
737 cejka 343
    low += middle1;
344
    high += (low < middle1);
345
    *lo = low;
346
    *hi = high;
347
    return;
348
}
349
 
350