Subversion Repositories HelenOS-historic

Rev

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

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