← Back to branch summary

~martin-decky/helenos/rcu

« back to all changes in this revision

Viewing changes to uspace/lib/softfloat/generic/mul.c

  • Committer: Martin Decky
  • Date: 2009-08-04 11:19:19 UTC
  • Revision ID: martin@uranus.dsrg.hide.ms.mff.cuni.cz-20090804111919-evyclddlr3v5lhmp
Initial import

Show diffs side-by-side

added added

removed removed

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