Subversion Repositories HelenOS-historic

Rev

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

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