Go to most recent revision | Details | 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; |
||
734 | cejka | 38 | __u32 exp1, exp2, mant1, mant2; |
731 | cejka | 39 | float32 result; |
40 | |||
41 | result.f = 0; |
||
42 | |||
734 | cejka | 43 | expdiff = a.parts.exp - b.parts.exp; |
44 | if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.mantisa < b.parts.mantisa))) { |
||
731 | cejka | 45 | if (isFloat32NaN(b)) { |
46 | //TODO: fix SigNaN |
||
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 | |||
734 | cejka | 59 | mant1 = b.parts.mantisa; |
60 | exp1 = b.parts.exp; |
||
61 | mant2 = a.parts.mantisa; |
||
62 | exp2 = a.parts.exp; |
||
63 | expdiff *= -1; |
||
731 | cejka | 64 | } else { |
65 | if (isFloat32NaN(a)) { |
||
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 */ |
75 | //TODO: fix exception |
||
76 | result.binary = FLOAT32_NAN; |
||
77 | return result; |
||
78 | }; |
||
79 | return a; |
||
80 | } |
||
81 | |||
82 | result.parts.sign = a.parts.sign; |
||
83 | |||
734 | cejka | 84 | mant1 = a.parts.mantisa; |
85 | exp1 = a.parts.exp; |
||
86 | mant2 = b.parts.mantisa; |
||
87 | exp2 = b.parts.exp; |
||
731 | cejka | 88 | }; |
89 | |||
734 | cejka | 90 | if (exp1 == 0) { |
731 | cejka | 91 | //both are denormalized |
734 | cejka | 92 | result.parts.mantisa = mant1-mant2; |
93 | if (result.parts.mantisa > mant1) { |
||
731 | cejka | 94 | //TODO: underflow exception |
95 | return result; |
||
96 | }; |
||
734 | cejka | 97 | result.parts.exp = 0; |
731 | cejka | 98 | return result; |
99 | }; |
||
734 | cejka | 100 | |
101 | /* add hidden bit */ |
||
102 | mant1 |= FLOAT32_HIDDEN_BIT_MASK; |
||
731 | cejka | 103 | |
734 | cejka | 104 | if (exp2 == 0) { |
105 | /* denormalized */ |
||
106 | --expdiff; |
||
107 | } else { |
||
108 | /* normalized */ |
||
109 | mant2 |= FLOAT32_HIDDEN_BIT_MASK; |
||
110 | }; |
||
731 | cejka | 111 | |
734 | cejka | 112 | /* create some space for rounding */ |
113 | mant1 <<= 6; |
||
114 | mant2 <<= 6; |
||
731 | cejka | 115 | |
734 | cejka | 116 | if (expdiff > FLOAT32_MANTISA_SIZE + 1) { |
117 | goto done; |
||
118 | }; |
||
731 | cejka | 119 | |
734 | cejka | 120 | mant1 = mant1 - (mant2 >> expdiff); |
121 | done: |
||
122 | //TODO: find first nonzero digit and shift result and detect possibly underflow |
||
123 | while ((exp1 > 0) && (!(mant1 & (FLOAT32_HIDDEN_BIT_MASK << 6 )))) { |
||
124 | --exp1; |
||
125 | mant1 <<= 1; |
||
126 | /* TODO: fix underflow - mant1 == 0 does not necessary means underflow... */ |
||
127 | }; |
||
128 | |||
129 | /* rounding - if first bit after mantisa is set then round up */ |
||
130 | mant1 += 0x20; |
||
131 | |||
132 | if (mant1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) { |
||
133 | ++exp1; |
||
134 | mant1 >>= 1; |
||
135 | }; |
||
136 | |||
137 | /*Clear hidden bit and shift */ |
||
138 | result.parts.mantisa = ((mant1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); |
||
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; |
||
150 | __u64 mant1, mant2; |
||
151 | float64 result; |
||
152 | |||
153 | result.d = 0; |
||
154 | |||
155 | expdiff = a.parts.exp - b.parts.exp; |
||
156 | if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.mantisa < b.parts.mantisa))) { |
||
157 | if (isFloat64NaN(b)) { |
||
158 | //TODO: fix SigNaN |
||
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 | |||
171 | mant1 = b.parts.mantisa; |
||
172 | exp1 = b.parts.exp; |
||
173 | mant2 = a.parts.mantisa; |
||
174 | exp2 = a.parts.exp; |
||
175 | expdiff *= -1; |
||
176 | } else { |
||
177 | if (isFloat64NaN(a)) { |
||
178 | //TODO: fix SigNaN |
||
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 */ |
||
187 | //TODO: fix exception |
||
188 | result.binary = FLOAT64_NAN; |
||
189 | return result; |
||
190 | }; |
||
191 | return a; |
||
192 | } |
||
193 | |||
194 | result.parts.sign = a.parts.sign; |
||
195 | |||
196 | mant1 = a.parts.mantisa; |
||
197 | exp1 = a.parts.exp; |
||
198 | mant2 = b.parts.mantisa; |
||
199 | exp2 = b.parts.exp; |
||
200 | }; |
||
201 | |||
202 | if (exp1 == 0) { |
||
203 | //both are denormalized |
||
204 | result.parts.mantisa = mant1 - mant2; |
||
205 | if (result.parts.mantisa > mant1) { |
||
206 | //TODO: underflow exception |
||
207 | return result; |
||
208 | }; |
||
209 | result.parts.exp = 0; |
||
210 | return result; |
||
211 | }; |
||
212 | |||
213 | /* add hidden bit */ |
||
214 | mant1 |= FLOAT64_HIDDEN_BIT_MASK; |
||
215 | |||
216 | if (exp2 == 0) { |
||
217 | /* denormalized */ |
||
731 | cejka | 218 | --expdiff; |
219 | } else { |
||
734 | cejka | 220 | /* normalized */ |
221 | mant2 |= FLOAT64_HIDDEN_BIT_MASK; |
||
731 | cejka | 222 | }; |
223 | |||
734 | cejka | 224 | /* create some space for rounding */ |
225 | mant1 <<= 6; |
||
226 | mant2 <<= 6; |
||
227 | |||
228 | if (expdiff > FLOAT64_MANTISA_SIZE + 1) { |
||
731 | cejka | 229 | goto done; |
230 | }; |
||
231 | |||
734 | cejka | 232 | mant1 = mant1 - (mant2 >> expdiff); |
731 | cejka | 233 | done: |
234 | //TODO: find first nonzero digit and shift result and detect possibly underflow |
||
734 | cejka | 235 | while ((exp1 > 0) && (!(mant1 & (FLOAT64_HIDDEN_BIT_MASK << 6 )))) { |
236 | --exp1; |
||
731 | cejka | 237 | mant1 <<= 1; |
238 | /* TODO: fix underflow - mant1 == 0 does not necessary means underflow... */ |
||
239 | }; |
||
240 | |||
734 | cejka | 241 | /* rounding - if first bit after mantisa is set then round up */ |
731 | cejka | 242 | mant1 += 0x20; |
243 | |||
734 | cejka | 244 | if (mant1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) { |
731 | cejka | 245 | ++exp1; |
734 | cejka | 246 | mant1 >>= 1; |
731 | cejka | 247 | }; |
248 | |||
734 | cejka | 249 | /*Clear hidden bit and shift */ |
250 | result.parts.mantisa = ((mant1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK)); |
||
731 | cejka | 251 | result.parts.exp = exp1; |
252 | |||
253 | return result; |
||
734 | cejka | 254 | } |
731 | cejka | 255 |