Subversion Repositories HelenOS-historic

Rev

Rev 828 | Rev 1031 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 828 Rev 834
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
#include<sftypes.h>
29
#include<sftypes.h>
30
#include<add.h>
30
#include<add.h>
31
#include<comparison.h>
31
#include<comparison.h>
32
 
32
 
33
/** Add two Float32 numbers with same signs
33
/** Add two Float32 numbers with same signs
34
 */
34
 */
35
float32 addFloat32(float32 a, float32 b)
35
float32 addFloat32(float32 a, float32 b)
36
{
36
{
37
    int expdiff;
37
    int expdiff;
38
    __u32 exp1, exp2,frac1, frac2;
38
    __u32 exp1, exp2,frac1, frac2;
39
   
39
   
40
    expdiff = a.parts.exp - b.parts.exp;
40
    expdiff = a.parts.exp - b.parts.exp;
41
    if (expdiff < 0) {
41
    if (expdiff < 0) {
42
        if (isFloat32NaN(b)) {
42
        if (isFloat32NaN(b)) {
43
            /* TODO: fix SigNaN */
43
            /* TODO: fix SigNaN */
44
            if (isFloat32SigNaN(b)) {
44
            if (isFloat32SigNaN(b)) {
45
            };
45
            };
46
 
46
 
47
            return b;
47
            return b;
48
        };
48
        };
49
       
49
       
50
        if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
50
        if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
51
            return b;
51
            return b;
52
        }
52
        }
53
       
53
       
54
        frac1 = b.parts.fraction;
54
        frac1 = b.parts.fraction;
55
        exp1 = b.parts.exp;
55
        exp1 = b.parts.exp;
56
        frac2 = a.parts.fraction;
56
        frac2 = a.parts.fraction;
57
        exp2 = a.parts.exp;
57
        exp2 = a.parts.exp;
58
        expdiff *= -1;
58
        expdiff *= -1;
59
    } else {
59
    } else {
60
        if ((isFloat32NaN(a)) || (isFloat32NaN(b))) {
60
        if ((isFloat32NaN(a)) || (isFloat32NaN(b))) {
61
            /* TODO: fix SigNaN */
61
            /* TODO: fix SigNaN */
62
            if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
62
            if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
63
            };
63
            };
64
            return (isFloat32NaN(a)?a:b);
64
            return (isFloat32NaN(a)?a:b);
65
        };
65
        };
66
       
66
       
67
        if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
67
        if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
68
            return a;
68
            return a;
69
        }
69
        }
70
       
70
       
71
        frac1 = a.parts.fraction;
71
        frac1 = a.parts.fraction;
72
        exp1 = a.parts.exp;
72
        exp1 = a.parts.exp;
73
        frac2 = b.parts.fraction;
73
        frac2 = b.parts.fraction;
74
        exp2 = b.parts.exp;
74
        exp2 = b.parts.exp;
75
    };
75
    };
76
   
76
   
77
    if (exp1 == 0) {
77
    if (exp1 == 0) {
78
        /* both are denormalized */
78
        /* both are denormalized */
79
        frac1 += frac2;
79
        frac1 += frac2;
80
        if (frac1 & FLOAT32_HIDDEN_BIT_MASK ) {
80
        if (frac1 & FLOAT32_HIDDEN_BIT_MASK ) {
81
            /* result is not denormalized */
81
            /* result is not denormalized */
82
            a.parts.exp = 1;
82
            a.parts.exp = 1;
83
        };
83
        };
84
        a.parts.fraction = frac1;
84
        a.parts.fraction = frac1;
85
        return a;
85
        return a;
86
    };
86
    };
87
   
87
   
88
    frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
88
    frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
89
 
89
 
90
    if (exp2 == 0) {
90
    if (exp2 == 0) {
91
        /* second operand is denormalized */
91
        /* second operand is denormalized */
92
        --expdiff;
92
        --expdiff;
93
    } else {
93
    } else {
94
        /* add hidden bit to second operand */
94
        /* add hidden bit to second operand */
95
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
95
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
96
    };
96
    };
97
   
97
   
98
    /* create some space for rounding */
98
    /* create some space for rounding */
99
    frac1 <<= 6;
99
    frac1 <<= 6;
100
    frac2 <<= 6;
100
    frac2 <<= 6;
101
   
101
   
102
    if (expdiff < (FLOAT32_FRACTION_SIZE + 2) ) {
102
    if (expdiff < (FLOAT32_FRACTION_SIZE + 2) ) {
103
        frac2 >>= expdiff;
103
        frac2 >>= expdiff;
104
        frac1 += frac2;
104
        frac1 += frac2;
-
 
105
    } else {
-
 
106
        a.parts.exp = exp1;
-
 
107
        a.parts.fraction = (frac1 >> 6) & (~(FLOAT32_HIDDEN_BIT_MASK));
-
 
108
        return a;
105
        };
109
    }
106
   
110
   
107
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
111
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
108
        ++exp1;
112
        ++exp1;
109
        frac1 >>= 1;
113
        frac1 >>= 1;
110
    };
114
    };
111
   
115
   
112
    /* rounding - if first bit after fraction is set then round up */
116
    /* rounding - if first bit after fraction is set then round up */
113
    frac1 += (0x1 << 5);
117
    frac1 += (0x1 << 5);
114
   
118
   
115
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
119
    if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
116
        /* rounding overflow */
120
        /* rounding overflow */
117
        ++exp1;
121
        ++exp1;
118
        frac1 >>= 1;
122
        frac1 >>= 1;
119
    };
123
    };
120
   
124
   
121
   
125
   
122
    if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
126
    if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
123
            /* overflow - set infinity as result */
127
            /* overflow - set infinity as result */
124
            a.parts.exp = FLOAT32_MAX_EXPONENT;
128
            a.parts.exp = FLOAT32_MAX_EXPONENT;
125
            a.parts.fraction = 0;
129
            a.parts.fraction = 0;
126
            return a;
130
            return a;
127
            }
131
            }
128
   
132
   
129
    a.parts.exp = exp1;
133
    a.parts.exp = exp1;
130
   
134
   
131
    /*Clear hidden bit and shift */
135
    /*Clear hidden bit and shift */
132
    a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
136
    a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
133
    return a;
137
    return a;
134
}
138
}
135
 
139
 
136
 
140
 
137
/** Add two Float64 numbers with same signs
141
/** Add two Float64 numbers with same signs
138
 */
142
 */
139
float64 addFloat64(float64 a, float64 b)
143
float64 addFloat64(float64 a, float64 b)
140
{
144
{
141
    int expdiff;
145
    int expdiff;
142
    __u32 exp1, exp2;
146
    __u32 exp1, exp2;
143
    __u64 frac1, frac2;
147
    __u64 frac1, frac2;
144
   
148
   
145
    expdiff = a.parts.exp - b.parts.exp;
149
    expdiff = ((int )a.parts.exp) - b.parts.exp;
146
    if (expdiff < 0) {
150
    if (expdiff < 0) {
147
        if (isFloat64NaN(b)) {
151
        if (isFloat64NaN(b)) {
148
            /* TODO: fix SigNaN */
152
            /* TODO: fix SigNaN */
149
            if (isFloat64SigNaN(b)) {
153
            if (isFloat64SigNaN(b)) {
150
            };
154
            };
151
 
155
 
152
            return b;
156
            return b;
153
        };
157
        };
154
       
158
       
155
        /* b is infinity and a not */  
159
        /* b is infinity and a not */  
156
        if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
160
        if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
157
            return b;
161
            return b;
158
        }
162
        }
159
       
163
       
160
        frac1 = b.parts.fraction;
164
        frac1 = b.parts.fraction;
161
        exp1 = b.parts.exp;
165
        exp1 = b.parts.exp;
162
        frac2 = a.parts.fraction;
166
        frac2 = a.parts.fraction;
163
        exp2 = a.parts.exp;
167
        exp2 = a.parts.exp;
164
        expdiff *= -1;
168
        expdiff *= -1;
165
    } else {
169
    } else {
166
        if (isFloat64NaN(a)) {
170
        if (isFloat64NaN(a)) {
167
            /* TODO: fix SigNaN */
171
            /* TODO: fix SigNaN */
168
            if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
172
            if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
169
            };
173
            };
170
            return a;
174
            return a;
171
        };
175
        };
172
       
176
       
173
        /* a is infinity and b not */
177
        /* a is infinity and b not */
174
        if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
178
        if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
175
            return a;
179
            return a;
176
        }
180
        }
177
       
181
       
178
        frac1 = a.parts.fraction;
182
        frac1 = a.parts.fraction;
179
        exp1 = a.parts.exp;
183
        exp1 = a.parts.exp;
180
        frac2 = b.parts.fraction;
184
        frac2 = b.parts.fraction;
181
        exp2 = b.parts.exp;
185
        exp2 = b.parts.exp;
182
    };
186
    };
183
   
187
   
184
    if (exp1 == 0) {
188
    if (exp1 == 0) {
185
        /* both are denormalized */
189
        /* both are denormalized */
186
        frac1 += frac2;
190
        frac1 += frac2;
187
        if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
191
        if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
188
            /* result is not denormalized */
192
            /* result is not denormalized */
189
            a.parts.exp = 1;
193
            a.parts.exp = 1;
190
        };
194
        };
191
        a.parts.fraction = frac1;
195
        a.parts.fraction = frac1;
192
        return a;
196
        return a;
193
    };
197
    };
194
   
198
   
195
    /* add hidden bit - frac1 is sure not denormalized */
199
    /* add hidden bit - frac1 is sure not denormalized */
196
    frac1 |= FLOAT64_HIDDEN_BIT_MASK;
200
    frac1 |= FLOAT64_HIDDEN_BIT_MASK;
197
 
201
 
198
    /* second operand ... */
202
    /* second operand ... */
199
    if (exp2 == 0) {
203
    if (exp2 == 0) {
200
        /* ... is denormalized */
204
        /* ... is denormalized */
201
        --expdiff; 
205
        --expdiff; 
202
    } else {
206
    } else {
203
        /* is not denormalized */
207
        /* is not denormalized */
204
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
208
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
205
    };
209
    };
206
   
210
   
207
    /* create some space for rounding */
211
    /* create some space for rounding */
208
    frac1 <<= 6;
212
    frac1 <<= 6;
209
    frac2 <<= 6;
213
    frac2 <<= 6;
210
   
214
   
211
    if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
215
    if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
212
        frac2 >>= expdiff;
216
        frac2 >>= expdiff;
213
        frac1 += frac2;
217
        frac1 += frac2;
-
 
218
    } else {
-
 
219
        a.parts.exp = exp1;
-
 
220
        a.parts.fraction = (frac1 >> 6) & (~(FLOAT64_HIDDEN_BIT_MASK));
-
 
221
        return a;
214
        };
222
    }
215
   
223
   
216
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
224
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
217
        ++exp1;
225
        ++exp1;
218
        frac1 >>= 1;
226
        frac1 >>= 1;
219
    };
227
    };
220
   
228
   
221
    /* rounding - if first bit after fraction is set then round up */
229
    /* rounding - if first bit after fraction is set then round up */
222
    frac1 += (0x1 << 5);
230
    frac1 += (0x1 << 5);
223
   
231
   
224
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
232
    if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
225
        /* rounding overflow */
233
        /* rounding overflow */
226
        ++exp1;
234
        ++exp1;
227
        frac1 >>= 1;
235
        frac1 >>= 1;
228
    };
236
    };
229
   
237
   
230
    if ((a.parts.exp == FLOAT64_MAX_EXPONENT ) || (a.parts.exp < exp1)) {
238
    if ((exp1 == FLOAT64_MAX_EXPONENT ) || (exp2 > exp1)) {
231
            /* overflow - set infinity as result */
239
            /* overflow - set infinity as result */
232
            a.parts.exp = FLOAT64_MAX_EXPONENT;
240
            a.parts.exp = FLOAT64_MAX_EXPONENT;
233
            a.parts.fraction = 0;
241
            a.parts.fraction = 0;
234
            return a;
242
            return a;
235
            }
243
            }
236
   
244
   
237
    a.parts.exp = exp1;
245
    a.parts.exp = exp1;
238
    /*Clear hidden bit and shift */
246
    /*Clear hidden bit and shift */
239
    a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
247
    a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
-
 
248
   
240
    return a;
249
    return a;
241
}
250
}
242
 
251
 
243
 
252
 
244
 
253