Subversion Repositories HelenOS

Rev

Rev 2479 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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