Subversion Repositories HelenOS-historic

Rev

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

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