Subversion Repositories HelenOS-historic

Rev

Rev 1657 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 1657 Rev 1740
1
/*
1
/*
2
 * Copyright (C) 2005 Josef Cejka
2
 * Copyright (C) 2005 Josef Cejka
3
 * All rights reserved.
3
 * All rights reserved.
4
 *
4
 *
5
 * Redistribution and use in source and binary forms, with or without
5
 * Redistribution and use in source and binary forms, with or without
6
 * modification, are permitted provided that the following conditions
6
 * modification, are permitted provided that the following conditions
7
 * are met:
7
 * are met:
8
 *
8
 *
9
 * - Redistributions of source code must retain the above copyright
9
 * - Redistributions of source code must retain the above copyright
10
 *   notice, this list of conditions and the following disclaimer.
10
 *   notice, this list of conditions and the following disclaimer.
11
 * - Redistributions in binary form must reproduce the above copyright
11
 * - Redistributions in binary form must reproduce the above copyright
12
 *   notice, this list of conditions and the following disclaimer in the
12
 *   notice, this list of conditions and the following disclaimer in the
13
 *   documentation and/or other materials provided with the distribution.
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
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.
15
 *   derived from this software without specific prior written permission.
16
 *
16
 *
17
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
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
18
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
19
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
20
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
21
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
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
23
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
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
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.
26
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27
 */
27
 */
28
 
28
 
29
 /** @addtogroup softfloat 
29
/** @addtogroup softfloat  
30
 * @{
30
 * @{
31
 */
31
 */
32
/** @file
32
/** @file
33
 */
33
 */
34
 
34
 
35
#include<sftypes.h>
35
#include<sftypes.h>
36
#include<add.h>
36
#include<add.h>
37
#include<comparison.h>
37
#include<comparison.h>
38
 
38
 
39
/** Add two Float32 numbers with same signs
39
/** Add two Float32 numbers with same signs
40
 */
40
 */
41
float32 addFloat32(float32 a, float32 b)
41
float32 addFloat32(float32 a, float32 b)
42
{
42
{
43
    int expdiff;
43
    int expdiff;
44
    uint32_t exp1, exp2,frac1, frac2;
44
    uint32_t exp1, exp2,frac1, frac2;
45
   
45
   
46
    expdiff = a.parts.exp - b.parts.exp;
46
    expdiff = a.parts.exp - b.parts.exp;
47
    if (expdiff < 0) {
47
    if (expdiff < 0) {
48
        if (isFloat32NaN(b)) {
48
        if (isFloat32NaN(b)) {
49
            /* TODO: fix SigNaN */
49
            /* TODO: fix SigNaN */
50
            if (isFloat32SigNaN(b)) {
50
            if (isFloat32SigNaN(b)) {
51
            };
51
            };
52
 
52
 
53
            return b;
53
            return b;
54
        };
54
        };
55
       
55
       
56
        if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
56
        if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
57
            return b;
57
            return b;
58
        }
58
        }
59
       
59
       
60
        frac1 = b.parts.fraction;
60
        frac1 = b.parts.fraction;
61
        exp1 = b.parts.exp;
61
        exp1 = b.parts.exp;
62
        frac2 = a.parts.fraction;
62
        frac2 = a.parts.fraction;
63
        exp2 = a.parts.exp;
63
        exp2 = a.parts.exp;
64
        expdiff *= -1;
64
        expdiff *= -1;
65
    } else {
65
    } else {
66
        if ((isFloat32NaN(a)) || (isFloat32NaN(b))) {
66
        if ((isFloat32NaN(a)) || (isFloat32NaN(b))) {
67
            /* TODO: fix SigNaN */
67
            /* TODO: fix SigNaN */
68
            if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
68
            if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
69
            };
69
            };
70
            return (isFloat32NaN(a)?a:b);
70
            return (isFloat32NaN(a)?a:b);
71
        };
71
        };
72
       
72
       
73
        if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
73
        if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
74
            return a;
74
            return a;
75
        }
75
        }
76
       
76
       
77
        frac1 = a.parts.fraction;
77
        frac1 = a.parts.fraction;
78
        exp1 = a.parts.exp;
78
        exp1 = a.parts.exp;
79
        frac2 = b.parts.fraction;
79
        frac2 = b.parts.fraction;
80
        exp2 = b.parts.exp;
80
        exp2 = b.parts.exp;
81
    };
81
    };
82
   
82
   
83
    if (exp1 == 0) {
83
    if (exp1 == 0) {
84
        /* both are denormalized */
84
        /* both are denormalized */
85
        frac1 += frac2;
85
        frac1 += frac2;
86
        if (frac1 & FLOAT32_HIDDEN_BIT_MASK ) {
86
        if (frac1 & FLOAT32_HIDDEN_BIT_MASK ) {
87
            /* result is not denormalized */
87
            /* result is not denormalized */
88
            a.parts.exp = 1;
88
            a.parts.exp = 1;
89
        };
89
        };
90
        a.parts.fraction = frac1;
90
        a.parts.fraction = frac1;
91
        return a;
91
        return a;
92
    };
92
    };
93
   
93
   
94
    frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
94
    frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
95
 
95
 
96
    if (exp2 == 0) {
96
    if (exp2 == 0) {
97
        /* second operand is denormalized */
97
        /* second operand is denormalized */
98
        --expdiff;
98
        --expdiff;
99
    } else {
99
    } else {
100
        /* add hidden bit to second operand */
100
        /* add hidden bit to second operand */
101
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
101
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
102
    };
102
    };
103
   
103
   
104
    /* create some space for rounding */
104
    /* create some space for rounding */
105
    frac1 <<= 6;
105
    frac1 <<= 6;
106
    frac2 <<= 6;
106
    frac2 <<= 6;
107
   
107
   
108
    if (expdiff < (FLOAT32_FRACTION_SIZE + 2) ) {
108
    if (expdiff < (FLOAT32_FRACTION_SIZE + 2) ) {
109
        frac2 >>= expdiff;
109
        frac2 >>= expdiff;
110
        frac1 += frac2;
110
        frac1 += frac2;
111
    } else {
111
    } else {
112
        a.parts.exp = exp1;
112
        a.parts.exp = exp1;
113
        a.parts.fraction = (frac1 >> 6) & (~(FLOAT32_HIDDEN_BIT_MASK));
113
        a.parts.fraction = (frac1 >> 6) & (~(FLOAT32_HIDDEN_BIT_MASK));
114
        return a;
114
        return a;
115
    }
115
    }
116
   
116
   
117
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
117
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
118
        ++exp1;
118
        ++exp1;
119
        frac1 >>= 1;
119
        frac1 >>= 1;
120
    };
120
    };
121
   
121
   
122
    /* rounding - if first bit after fraction is set then round up */
122
    /* rounding - if first bit after fraction is set then round up */
123
    frac1 += (0x1 << 5);
123
    frac1 += (0x1 << 5);
124
   
124
   
125
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
125
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
126
        /* rounding overflow */
126
        /* rounding overflow */
127
        ++exp1;
127
        ++exp1;
128
        frac1 >>= 1;
128
        frac1 >>= 1;
129
    };
129
    };
130
   
130
   
131
   
131
   
132
    if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
132
    if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
133
            /* overflow - set infinity as result */
133
            /* overflow - set infinity as result */
134
            a.parts.exp = FLOAT32_MAX_EXPONENT;
134
            a.parts.exp = FLOAT32_MAX_EXPONENT;
135
            a.parts.fraction = 0;
135
            a.parts.fraction = 0;
136
            return a;
136
            return a;
137
            }
137
            }
138
   
138
   
139
    a.parts.exp = exp1;
139
    a.parts.exp = exp1;
140
   
140
   
141
    /*Clear hidden bit and shift */
141
    /*Clear hidden bit and shift */
142
    a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
142
    a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
143
    return a;
143
    return a;
144
}
144
}
145
 
145
 
146
 
146
 
147
/** Add two Float64 numbers with same signs
147
/** Add two Float64 numbers with same signs
148
 */
148
 */
149
float64 addFloat64(float64 a, float64 b)
149
float64 addFloat64(float64 a, float64 b)
150
{
150
{
151
    int expdiff;
151
    int expdiff;
152
    uint32_t exp1, exp2;
152
    uint32_t exp1, exp2;
153
    uint64_t frac1, frac2;
153
    uint64_t frac1, frac2;
154
   
154
   
155
    expdiff = ((int )a.parts.exp) - b.parts.exp;
155
    expdiff = ((int )a.parts.exp) - b.parts.exp;
156
    if (expdiff < 0) {
156
    if (expdiff < 0) {
157
        if (isFloat64NaN(b)) {
157
        if (isFloat64NaN(b)) {
158
            /* TODO: fix SigNaN */
158
            /* TODO: fix SigNaN */
159
            if (isFloat64SigNaN(b)) {
159
            if (isFloat64SigNaN(b)) {
160
            };
160
            };
161
 
161
 
162
            return b;
162
            return b;
163
        };
163
        };
164
       
164
       
165
        /* b is infinity and a not */  
165
        /* b is infinity and a not */  
166
        if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
166
        if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
167
            return b;
167
            return b;
168
        }
168
        }
169
       
169
       
170
        frac1 = b.parts.fraction;
170
        frac1 = b.parts.fraction;
171
        exp1 = b.parts.exp;
171
        exp1 = b.parts.exp;
172
        frac2 = a.parts.fraction;
172
        frac2 = a.parts.fraction;
173
        exp2 = a.parts.exp;
173
        exp2 = a.parts.exp;
174
        expdiff *= -1;
174
        expdiff *= -1;
175
    } else {
175
    } else {
176
        if (isFloat64NaN(a)) {
176
        if (isFloat64NaN(a)) {
177
            /* TODO: fix SigNaN */
177
            /* TODO: fix SigNaN */
178
            if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
178
            if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
179
            };
179
            };
180
            return a;
180
            return a;
181
        };
181
        };
182
       
182
       
183
        /* a is infinity and b not */
183
        /* a is infinity and b not */
184
        if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
184
        if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
185
            return a;
185
            return a;
186
        }
186
        }
187
       
187
       
188
        frac1 = a.parts.fraction;
188
        frac1 = a.parts.fraction;
189
        exp1 = a.parts.exp;
189
        exp1 = a.parts.exp;
190
        frac2 = b.parts.fraction;
190
        frac2 = b.parts.fraction;
191
        exp2 = b.parts.exp;
191
        exp2 = b.parts.exp;
192
    };
192
    };
193
   
193
   
194
    if (exp1 == 0) {
194
    if (exp1 == 0) {
195
        /* both are denormalized */
195
        /* both are denormalized */
196
        frac1 += frac2;
196
        frac1 += frac2;
197
        if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
197
        if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
198
            /* result is not denormalized */
198
            /* result is not denormalized */
199
            a.parts.exp = 1;
199
            a.parts.exp = 1;
200
        };
200
        };
201
        a.parts.fraction = frac1;
201
        a.parts.fraction = frac1;
202
        return a;
202
        return a;
203
    };
203
    };
204
   
204
   
205
    /* add hidden bit - frac1 is sure not denormalized */
205
    /* add hidden bit - frac1 is sure not denormalized */
206
    frac1 |= FLOAT64_HIDDEN_BIT_MASK;
206
    frac1 |= FLOAT64_HIDDEN_BIT_MASK;
207
 
207
 
208
    /* second operand ... */
208
    /* second operand ... */
209
    if (exp2 == 0) {
209
    if (exp2 == 0) {
210
        /* ... is denormalized */
210
        /* ... is denormalized */
211
        --expdiff; 
211
        --expdiff; 
212
    } else {
212
    } else {
213
        /* is not denormalized */
213
        /* is not denormalized */
214
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
214
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
215
    };
215
    };
216
   
216
   
217
    /* create some space for rounding */
217
    /* create some space for rounding */
218
    frac1 <<= 6;
218
    frac1 <<= 6;
219
    frac2 <<= 6;
219
    frac2 <<= 6;
220
   
220
   
221
    if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
221
    if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
222
        frac2 >>= expdiff;
222
        frac2 >>= expdiff;
223
        frac1 += frac2;
223
        frac1 += frac2;
224
    } else {
224
    } else {
225
        a.parts.exp = exp1;
225
        a.parts.exp = exp1;
226
        a.parts.fraction = (frac1 >> 6) & (~(FLOAT64_HIDDEN_BIT_MASK));
226
        a.parts.fraction = (frac1 >> 6) & (~(FLOAT64_HIDDEN_BIT_MASK));
227
        return a;
227
        return a;
228
    }
228
    }
229
   
229
   
230
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
230
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
231
        ++exp1;
231
        ++exp1;
232
        frac1 >>= 1;
232
        frac1 >>= 1;
233
    };
233
    };
234
   
234
   
235
    /* rounding - if first bit after fraction is set then round up */
235
    /* rounding - if first bit after fraction is set then round up */
236
    frac1 += (0x1 << 5);
236
    frac1 += (0x1 << 5);
237
   
237
   
238
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
238
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
239
        /* rounding overflow */
239
        /* rounding overflow */
240
        ++exp1;
240
        ++exp1;
241
        frac1 >>= 1;
241
        frac1 >>= 1;
242
    };
242
    };
243
   
243
   
244
    if ((exp1 == FLOAT64_MAX_EXPONENT ) || (exp2 > exp1)) {
244
    if ((exp1 == FLOAT64_MAX_EXPONENT ) || (exp2 > exp1)) {
245
            /* overflow - set infinity as result */
245
            /* overflow - set infinity as result */
246
            a.parts.exp = FLOAT64_MAX_EXPONENT;
246
            a.parts.exp = FLOAT64_MAX_EXPONENT;
247
            a.parts.fraction = 0;
247
            a.parts.fraction = 0;
248
            return a;
248
            return a;
249
            }
249
            }
250
   
250
   
251
    a.parts.exp = exp1;
251
    a.parts.exp = exp1;
252
    /*Clear hidden bit and shift */
252
    /*Clear hidden bit and shift */
253
    a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
253
    a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
254
   
254
   
255
    return a;
255
    return a;
256
}
256
}
257
 
257
 
258
 
-
 
259
 
-
 
260
 /** @}
258
/** @}
261
 */
259
 */
262
 
-
 
263
 
260