Subversion Repositories HelenOS-historic

Rev

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