Subversion Repositories HelenOS

Rev

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

Rev Author Line No. Line
731 cejka 1
/*
2071 jermar 2
 * Copyright (c) 2005 Josef Cejka
731 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
 
1740 jermar 29
/** @addtogroup softfloat  
1657 cejka 30
 * @{
31
 */
32
/** @file
33
 */
34
 
731 cejka 35
#include<sftypes.h>
36
#include<mul.h>
37
#include<comparison.h>
829 cejka 38
#include<common.h>
731 cejka 39
 
40
/** Multiply two 32 bit float numbers
41
 *
42
 */
43
float32 mulFloat32(float32 a, float32 b)
44
{
45
    float32 result;
1031 cejka 46
    uint64_t frac1, frac2;
47
    int32_t exp;
731 cejka 48
 
49
    result.parts.sign = a.parts.sign ^ b.parts.sign;
50
 
737 cejka 51
    if (isFloat32NaN(a) || isFloat32NaN(b) ) {
731 cejka 52
        /* TODO: fix SigNaNs */
53
        if (isFloat32SigNaN(a)) {
804 cejka 54
            result.parts.fraction = a.parts.fraction;
731 cejka 55
            result.parts.exp = a.parts.exp;
56
            return result;
57
        };
58
        if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
804 cejka 59
            result.parts.fraction = b.parts.fraction;
731 cejka 60
            result.parts.exp = b.parts.exp;
61
            return result;
62
        };
63
        /* set NaN as result */
737 cejka 64
        result.binary = FLOAT32_NAN;
731 cejka 65
        return result;
66
    };
67
 
68
    if (isFloat32Infinity(a)) {
69
        if (isFloat32Zero(b)) {
70
            /* FIXME: zero * infinity */
737 cejka 71
            result.binary = FLOAT32_NAN;
731 cejka 72
            return result;
73
        }
804 cejka 74
        result.parts.fraction = a.parts.fraction;
731 cejka 75
        result.parts.exp = a.parts.exp;
76
        return result;
77
    }
78
 
79
    if (isFloat32Infinity(b)) {
80
        if (isFloat32Zero(a)) {
81
            /* FIXME: zero * infinity */
737 cejka 82
            result.binary = FLOAT32_NAN;
731 cejka 83
            return result;
84
        }
804 cejka 85
        result.parts.fraction = b.parts.fraction;
731 cejka 86
        result.parts.exp = b.parts.exp;
87
        return result;
88
    }
89
 
90
    /* exp is signed so we can easy detect underflow */
91
    exp = a.parts.exp + b.parts.exp;
92
    exp -= FLOAT32_BIAS;
93
 
737 cejka 94
    if (exp >= FLOAT32_MAX_EXPONENT) {
731 cejka 95
        /* FIXME: overflow */
96
        /* set infinity as result */
737 cejka 97
        result.binary = FLOAT32_INF;
98
        result.parts.sign = a.parts.sign ^ b.parts.sign;
731 cejka 99
        return result;
100
    };
101
 
102
    if (exp < 0) {
103
        /* FIXME: underflow */
104
        /* return signed zero */
804 cejka 105
        result.parts.fraction = 0x0;
731 cejka 106
        result.parts.exp = 0x0;
107
        return result;
108
    };
109
 
804 cejka 110
    frac1 = a.parts.fraction;
737 cejka 111
    if (a.parts.exp > 0) {
804 cejka 112
        frac1 |= FLOAT32_HIDDEN_BIT_MASK;
731 cejka 113
    } else {
114
        ++exp;
115
    };
116
 
804 cejka 117
    frac2 = b.parts.fraction;
737 cejka 118
 
119
    if (b.parts.exp > 0) {
804 cejka 120
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
731 cejka 121
    } else {
122
        ++exp;
123
    };
124
 
804 cejka 125
    frac1 <<= 1; /* one bit space for rounding */
731 cejka 126
 
804 cejka 127
    frac1 = frac1 * frac2;
731 cejka 128
/* round and return */
129
 
804 cejka 130
    while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
131
        /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
731 cejka 132
        ++exp;
804 cejka 133
        frac1 >>= 1;
731 cejka 134
    };
135
 
136
    /* rounding */
804 cejka 137
    /* ++frac1; FIXME: not works - without it is ok */
138
    frac1 >>= 1; /* shift off rounding space */
731 cejka 139
 
804 cejka 140
    if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
731 cejka 141
        ++exp;
804 cejka 142
        frac1 >>= 1;
731 cejka 143
    };
144
 
737 cejka 145
    if (exp >= FLOAT32_MAX_EXPONENT ) {
731 cejka 146
        /* TODO: fix overflow */
147
        /* return infinity*/
737 cejka 148
        result.parts.exp = FLOAT32_MAX_EXPONENT;
804 cejka 149
        result.parts.fraction = 0x0;
731 cejka 150
        return result;
151
    }
152
 
804 cejka 153
    exp -= FLOAT32_FRACTION_SIZE;
731 cejka 154
 
804 cejka 155
    if (exp <= FLOAT32_FRACTION_SIZE) {
731 cejka 156
        /* denormalized number */
804 cejka 157
        frac1 >>= 1; /* denormalize */
158
        while ((frac1 > 0) && (exp < 0)) {
159
            frac1 >>= 1;
731 cejka 160
            ++exp;
161
        };
804 cejka 162
        if (frac1 == 0) {
731 cejka 163
            /* FIXME : underflow */
164
        result.parts.exp = 0;
804 cejka 165
        result.parts.fraction = 0;
731 cejka 166
        return result;
167
        };
168
    };
169
    result.parts.exp = exp;
804 cejka 170
    result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
731 cejka 171
 
172
    return result; 
173
 
174
}
175
 
737 cejka 176
/** Multiply two 64 bit float numbers
177
 *
178
 */
179
float64 mulFloat64(float64 a, float64 b)
180
{
181
    float64 result;
1031 cejka 182
    uint64_t frac1, frac2;
183
    int32_t exp;
731 cejka 184
 
737 cejka 185
    result.parts.sign = a.parts.sign ^ b.parts.sign;
186
 
187
    if (isFloat64NaN(a) || isFloat64NaN(b) ) {
188
        /* TODO: fix SigNaNs */
189
        if (isFloat64SigNaN(a)) {
804 cejka 190
            result.parts.fraction = a.parts.fraction;
737 cejka 191
            result.parts.exp = a.parts.exp;
192
            return result;
193
        };
194
        if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
804 cejka 195
            result.parts.fraction = b.parts.fraction;
737 cejka 196
            result.parts.exp = b.parts.exp;
197
            return result;
198
        };
199
        /* set NaN as result */
200
        result.binary = FLOAT64_NAN;
201
        return result;
202
    };
203
 
204
    if (isFloat64Infinity(a)) {
205
        if (isFloat64Zero(b)) {
206
            /* FIXME: zero * infinity */
207
            result.binary = FLOAT64_NAN;
208
            return result;
209
        }
804 cejka 210
        result.parts.fraction = a.parts.fraction;
737 cejka 211
        result.parts.exp = a.parts.exp;
212
        return result;
213
    }
731 cejka 214
 
737 cejka 215
    if (isFloat64Infinity(b)) {
216
        if (isFloat64Zero(a)) {
217
            /* FIXME: zero * infinity */
218
            result.binary = FLOAT64_NAN;
219
            return result;
220
        }
804 cejka 221
        result.parts.fraction = b.parts.fraction;
737 cejka 222
        result.parts.exp = b.parts.exp;
223
        return result;
224
    }
225
 
226
    /* exp is signed so we can easy detect underflow */
829 cejka 227
    exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
737 cejka 228
 
804 cejka 229
    frac1 = a.parts.fraction;
829 cejka 230
 
737 cejka 231
    if (a.parts.exp > 0) {
804 cejka 232
        frac1 |= FLOAT64_HIDDEN_BIT_MASK;
737 cejka 233
    } else {
234
        ++exp;
235
    };
236
 
804 cejka 237
    frac2 = b.parts.fraction;
737 cejka 238
 
239
    if (b.parts.exp > 0) {
804 cejka 240
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
737 cejka 241
    } else {
242
        ++exp;
243
    };
244
 
829 cejka 245
    frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
246
    frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
737 cejka 247
 
804 cejka 248
    mul64integers(frac1, frac2, &frac1, &frac2);
737 cejka 249
 
829 cejka 250
    frac2 |= (frac1 != 0);
251
    if (frac2 & (0x1ll << 62)) {
252
        frac2 <<= 1;
253
        exp--;
737 cejka 254
    }
255
 
829 cejka 256
    result = finishFloat64(exp, frac2, result.parts.sign);
257
    return result;
737 cejka 258
}
259
 
260
/** Multiply two 64 bit numbers and return result in two parts
261
 * @param a first operand
262
 * @param b second operand
263
 * @param lo lower part from result
264
 * @param hi higher part of result
265
 */
1031 cejka 266
void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi)
737 cejka 267
{
1031 cejka 268
    uint64_t low, high, middle1, middle2;
269
    uint32_t alow, blow;
829 cejka 270
 
737 cejka 271
    alow = a & 0xFFFFFFFF;
272
    blow = b & 0xFFFFFFFF;
273
 
828 cejka 274
    a >>= 32;
275
    b >>= 32;
737 cejka 276
 
1031 cejka 277
    low = ((uint64_t)alow) * blow;
737 cejka 278
    middle1 = a * blow;
279
    middle2 = alow * b;
280
    high = a * b;
281
 
282
    middle1 += middle2;
1031 cejka 283
    high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
804 cejka 284
    middle1 <<= 32;
737 cejka 285
    low += middle1;
286
    high += (low < middle1);
287
    *lo = low;
288
    *hi = high;
828 cejka 289
 
737 cejka 290
    return;
291
}
292
 
1740 jermar 293
/** @}
1657 cejka 294
 */