Subversion Repositories HelenOS

Rev

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

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