WebSVN – HelenOS-historic – Diff – /uspace/trunk/softfloat/generic/mul.c


/*
 * Copyright (C) 2005 Josef Cejka
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * - Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * - Redistributions in binary form must reproduce the above copyright
 *   notice, this list of conditions and the following disclaimer in the
 *   documentation and/or other materials provided with the distribution.
 * - The name of the author may not be used to endorse or promote products
 *   derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
#include<sftypes.h>
#include<mul.h>
#include<comparison.h>
#include<common.h>
 
/** Multiply two 32 bit float numbers
 *
 */
float32 mulFloat32(float32 a, float32 b)
{
    float32 result;
    uint64_t frac1, frac2;
    int32_t exp;
 
    result.parts.sign = a.parts.sign ^ b.parts.sign;
   
    if (isFloat32NaN(a) || isFloat32NaN(b) ) {
        /* TODO: fix SigNaNs */
        if (isFloat32SigNaN(a)) {
            result.parts.fraction = a.parts.fraction;
            result.parts.exp = a.parts.exp;
            return result;
        };
        if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
            result.parts.fraction = b.parts.fraction;
            result.parts.exp = b.parts.exp;
            return result;
        };
        /* set NaN as result */
        result.binary = FLOAT32_NAN;
        return result;
    };
       
    if (isFloat32Infinity(a)) {
        if (isFloat32Zero(b)) {
            /* FIXME: zero * infinity */
            result.binary = FLOAT32_NAN;
            return result;
        }
        result.parts.fraction = a.parts.fraction;
        result.parts.exp = a.parts.exp;
        return result;
    }
 
    if (isFloat32Infinity(b)) {
        if (isFloat32Zero(a)) {
            /* FIXME: zero * infinity */
            result.binary = FLOAT32_NAN;
            return result;
        }
        result.parts.fraction = b.parts.fraction;
        result.parts.exp = b.parts.exp;
        return result;
    }
 
    /* exp is signed so we can easy detect underflow */
    exp = a.parts.exp + b.parts.exp;
    exp -= FLOAT32_BIAS;
   
    if (exp >= FLOAT32_MAX_EXPONENT) {
        /* FIXME: overflow */
        /* set infinity as result */
        result.binary = FLOAT32_INF;
        result.parts.sign = a.parts.sign ^ b.parts.sign;
        return result;
    };
   
    if (exp < 0) {
        /* FIXME: underflow */
        /* return signed zero */
        result.parts.fraction = 0x0;
        result.parts.exp = 0x0;
        return result;
    };
   
    frac1 = a.parts.fraction;
    if (a.parts.exp > 0) {
        frac1 |= FLOAT32_HIDDEN_BIT_MASK;
    } else {
        ++exp;
    };
   
    frac2 = b.parts.fraction;
 
    if (b.parts.exp > 0) {
        frac2 |= FLOAT32_HIDDEN_BIT_MASK;
    } else {
        ++exp;
    };
 
    frac1 <<= 1; /* one bit space for rounding */
 
    frac1 = frac1 * frac2;
/* round and return */
   
    while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
        /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
        ++exp;
        frac1 >>= 1;
    };
 
    /* rounding */
    /* ++frac1; FIXME: not works - without it is ok */
    frac1 >>= 1; /* shift off rounding space */
   
    if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
        ++exp;
        frac1 >>= 1;
    };
 
    if (exp >= FLOAT32_MAX_EXPONENT ) {
        /* TODO: fix overflow */
        /* return infinity*/
        result.parts.exp = FLOAT32_MAX_EXPONENT;
        result.parts.fraction = 0x0;
        return result;
    }
   
    exp -= FLOAT32_FRACTION_SIZE;
 
    if (exp <= FLOAT32_FRACTION_SIZE) {
        /* denormalized number */
        frac1 >>= 1; /* denormalize */
        while ((frac1 > 0) && (exp < 0)) {
            frac1 >>= 1;
            ++exp;
        };
        if (frac1 == 0) {
            /* FIXME : underflow */
        result.parts.exp = 0;
        result.parts.fraction = 0;
        return result;
        };
    };
    result.parts.exp = exp;
    result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
   
    return result; 
   
}
 
/** Multiply two 64 bit float numbers
 *
 */
float64 mulFloat64(float64 a, float64 b)
{
    float64 result;
    uint64_t frac1, frac2;
    int32_t exp;
 
    result.parts.sign = a.parts.sign ^ b.parts.sign;
   
    if (isFloat64NaN(a) || isFloat64NaN(b) ) {
        /* TODO: fix SigNaNs */
        if (isFloat64SigNaN(a)) {
            result.parts.fraction = a.parts.fraction;
            result.parts.exp = a.parts.exp;
            return result;
        };
        if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
            result.parts.fraction = b.parts.fraction;
            result.parts.exp = b.parts.exp;
            return result;
        };
        /* set NaN as result */
        result.binary = FLOAT64_NAN;
        return result;
    };
       
    if (isFloat64Infinity(a)) {
        if (isFloat64Zero(b)) {
            /* FIXME: zero * infinity */
            result.binary = FLOAT64_NAN;
            return result;
        }
        result.parts.fraction = a.parts.fraction;
        result.parts.exp = a.parts.exp;
        return result;
    }
 
    if (isFloat64Infinity(b)) {
        if (isFloat64Zero(a)) {
            /* FIXME: zero * infinity */
            result.binary = FLOAT64_NAN;
            return result;
        }
        result.parts.fraction = b.parts.fraction;
        result.parts.exp = b.parts.exp;
        return result;
    }
 
    /* exp is signed so we can easy detect underflow */
    exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
   
    frac1 = a.parts.fraction;
 
    if (a.parts.exp > 0) {
        frac1 |= FLOAT64_HIDDEN_BIT_MASK;
    } else {
        ++exp;
    };
   
    frac2 = b.parts.fraction;
 
    if (b.parts.exp > 0) {
        frac2 |= FLOAT64_HIDDEN_BIT_MASK;
    } else {
        ++exp;
    };
 
    frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
    frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
 
    mul64integers(frac1, frac2, &frac1, &frac2);
 
    frac2 |= (frac1 != 0);
    if (frac2 & (0x1ll << 62)) {
        frac2 <<= 1;
        exp--;
    }
 
    result = finishFloat64(exp, frac2, result.parts.sign);
    return result;
}
 
/** Multiply two 64 bit numbers and return result in two parts
 * @param a first operand
 * @param b second operand
 * @param lo lower part from result
 * @param hi higher part of result
 */
void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi)
{
    uint64_t low, high, middle1, middle2;
    uint32_t alow, blow;
 
    alow = a & 0xFFFFFFFF;
    blow = b & 0xFFFFFFFF;
   
    a >>= 32;
    b >>= 32;
   
    low = ((uint64_t)alow) * blow;
    middle1 = a * blow;
    middle2 = alow * b;
    high = a * b;
 
    middle1 += middle2;
    high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
    middle1 <<= 32;
    low += middle1;
    high += (low < middle1);
    *lo = low;
    *hi = high;
   
    return;
}
 
 
 

Rev 829	Rev 1031
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	#include<sftypes.h>	29	#include<sftypes.h>
30	#include<mul.h>	30	#include<mul.h>
31	#include<comparison.h>	31	#include<comparison.h>
32	#include<common.h>	32	#include<common.h>
33		33
34	/** Multiply two 32 bit float numbers	34	/** Multiply two 32 bit float numbers
35	*	35	*
36	*/	36	*/
37	float32 mulFloat32(float32 a, float32 b)	37	float32 mulFloat32(float32 a, float32 b)
38	{	38	{
39	float32 result;	39	float32 result;
40	__u64 frac1, frac2;	40	uint64_t frac1, frac2;
41	__s32 exp;	41	int32_t exp;
42		42
43	result.parts.sign = a.parts.sign ^ b.parts.sign;	43	result.parts.sign = a.parts.sign ^ b.parts.sign;
44		44
45	if (isFloat32NaN(a) \|\| isFloat32NaN(b) ) {	45	if (isFloat32NaN(a) \|\| isFloat32NaN(b) ) {
46	/* TODO: fix SigNaNs */	46	/* TODO: fix SigNaNs */
47	if (isFloat32SigNaN(a)) {	47	if (isFloat32SigNaN(a)) {
48	result.parts.fraction = a.parts.fraction;	48	result.parts.fraction = a.parts.fraction;
49	result.parts.exp = a.parts.exp;	49	result.parts.exp = a.parts.exp;
50	return result;	50	return result;
51	};	51	};
52	if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */	52	if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
53	result.parts.fraction = b.parts.fraction;	53	result.parts.fraction = b.parts.fraction;
54	result.parts.exp = b.parts.exp;	54	result.parts.exp = b.parts.exp;
55	return result;	55	return result;
56	};	56	};
57	/* set NaN as result */	57	/* set NaN as result */
58	result.binary = FLOAT32_NAN;	58	result.binary = FLOAT32_NAN;
59	return result;	59	return result;
60	};	60	};
61		61
62	if (isFloat32Infinity(a)) {	62	if (isFloat32Infinity(a)) {
63	if (isFloat32Zero(b)) {	63	if (isFloat32Zero(b)) {
64	/* FIXME: zero * infinity */	64	/* FIXME: zero * infinity */
65	result.binary = FLOAT32_NAN;	65	result.binary = FLOAT32_NAN;
66	return result;	66	return result;
67	}	67	}
68	result.parts.fraction = a.parts.fraction;	68	result.parts.fraction = a.parts.fraction;
69	result.parts.exp = a.parts.exp;	69	result.parts.exp = a.parts.exp;
70	return result;	70	return result;
71	}	71	}
72		72
73	if (isFloat32Infinity(b)) {	73	if (isFloat32Infinity(b)) {
74	if (isFloat32Zero(a)) {	74	if (isFloat32Zero(a)) {
75	/* FIXME: zero * infinity */	75	/* FIXME: zero * infinity */
76	result.binary = FLOAT32_NAN;	76	result.binary = FLOAT32_NAN;
77	return result;	77	return result;
78	}	78	}
79	result.parts.fraction = b.parts.fraction;	79	result.parts.fraction = b.parts.fraction;
80	result.parts.exp = b.parts.exp;	80	result.parts.exp = b.parts.exp;
81	return result;	81	return result;
82	}	82	}
83		83
84	/* exp is signed so we can easy detect underflow */	84	/* exp is signed so we can easy detect underflow */
85	exp = a.parts.exp + b.parts.exp;	85	exp = a.parts.exp + b.parts.exp;
86	exp -= FLOAT32_BIAS;	86	exp -= FLOAT32_BIAS;
87		87
88	if (exp >= FLOAT32_MAX_EXPONENT) {	88	if (exp >= FLOAT32_MAX_EXPONENT) {
89	/* FIXME: overflow */	89	/* FIXME: overflow */
90	/* set infinity as result */	90	/* set infinity as result */
91	result.binary = FLOAT32_INF;	91	result.binary = FLOAT32_INF;
92	result.parts.sign = a.parts.sign ^ b.parts.sign;	92	result.parts.sign = a.parts.sign ^ b.parts.sign;
93	return result;	93	return result;
94	};	94	};
95		95
96	if (exp < 0) {	96	if (exp < 0) {
97	/* FIXME: underflow */	97	/* FIXME: underflow */
98	/* return signed zero */	98	/* return signed zero */
99	result.parts.fraction = 0x0;	99	result.parts.fraction = 0x0;
100	result.parts.exp = 0x0;	100	result.parts.exp = 0x0;
101	return result;	101	return result;
102	};	102	};
103		103
104	frac1 = a.parts.fraction;	104	frac1 = a.parts.fraction;
105	if (a.parts.exp > 0) {	105	if (a.parts.exp > 0) {
106	frac1 \|= FLOAT32_HIDDEN_BIT_MASK;	106	frac1 \|= FLOAT32_HIDDEN_BIT_MASK;
107	} else {	107	} else {
108	++exp;	108	++exp;
109	};	109	};
110		110
111	frac2 = b.parts.fraction;	111	frac2 = b.parts.fraction;
112		112
113	if (b.parts.exp > 0) {	113	if (b.parts.exp > 0) {
114	frac2 \|= FLOAT32_HIDDEN_BIT_MASK;	114	frac2 \|= FLOAT32_HIDDEN_BIT_MASK;
115	} else {	115	} else {
116	++exp;	116	++exp;
117	};	117	};
118		118
119	frac1 <<= 1; /* one bit space for rounding */	119	frac1 <<= 1; /* one bit space for rounding */
120		120
121	frac1 = frac1 * frac2;	121	frac1 = frac1 * frac2;
122	/* round and return */	122	/* round and return */
123		123
124	while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {	124	while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
125	/* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/	125	/* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
126	++exp;	126	++exp;
127	frac1 >>= 1;	127	frac1 >>= 1;
128	};	128	};
129		129
130	/* rounding */	130	/* rounding */
131	/* ++frac1; FIXME: not works - without it is ok */	131	/* ++frac1; FIXME: not works - without it is ok */
132	frac1 >>= 1; /* shift off rounding space */	132	frac1 >>= 1; /* shift off rounding space */
133		133
134	if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {	134	if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
135	++exp;	135	++exp;
136	frac1 >>= 1;	136	frac1 >>= 1;
137	};	137	};
138		138
139	if (exp >= FLOAT32_MAX_EXPONENT ) {	139	if (exp >= FLOAT32_MAX_EXPONENT ) {
140	/* TODO: fix overflow */	140	/* TODO: fix overflow */
141	/* return infinity*/	141	/* return infinity*/
142	result.parts.exp = FLOAT32_MAX_EXPONENT;	142	result.parts.exp = FLOAT32_MAX_EXPONENT;
143	result.parts.fraction = 0x0;	143	result.parts.fraction = 0x0;
144	return result;	144	return result;
145	}	145	}
146		146
147	exp -= FLOAT32_FRACTION_SIZE;	147	exp -= FLOAT32_FRACTION_SIZE;
148		148
149	if (exp <= FLOAT32_FRACTION_SIZE) {	149	if (exp <= FLOAT32_FRACTION_SIZE) {
150	/* denormalized number */	150	/* denormalized number */
151	frac1 >>= 1; /* denormalize */	151	frac1 >>= 1; /* denormalize */
152	while ((frac1 > 0) && (exp < 0)) {	152	while ((frac1 > 0) && (exp < 0)) {
153	frac1 >>= 1;	153	frac1 >>= 1;
154	++exp;	154	++exp;
155	};	155	};
156	if (frac1 == 0) {	156	if (frac1 == 0) {
157	/* FIXME : underflow */	157	/* FIXME : underflow */
158	result.parts.exp = 0;	158	result.parts.exp = 0;
159	result.parts.fraction = 0;	159	result.parts.fraction = 0;
160	return result;	160	return result;
161	};	161	};
162	};	162	};
163	result.parts.exp = exp;	163	result.parts.exp = exp;
164	result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);	164	result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
165		165
166	return result;	166	return result;
167		167
168	}	168	}
169		169
170	/** Multiply two 64 bit float numbers	170	/** Multiply two 64 bit float numbers
171	*	171	*
172	*/	172	*/
173	float64 mulFloat64(float64 a, float64 b)	173	float64 mulFloat64(float64 a, float64 b)
174	{	174	{
175	float64 result;	175	float64 result;
176	__u64 frac1, frac2;	176	uint64_t frac1, frac2;
177	__s32 exp;	177	int32_t exp;
178		178
179	result.parts.sign = a.parts.sign ^ b.parts.sign;	179	result.parts.sign = a.parts.sign ^ b.parts.sign;
180		180
181	if (isFloat64NaN(a) \|\| isFloat64NaN(b) ) {	181	if (isFloat64NaN(a) \|\| isFloat64NaN(b) ) {
182	/* TODO: fix SigNaNs */	182	/* TODO: fix SigNaNs */
183	if (isFloat64SigNaN(a)) {	183	if (isFloat64SigNaN(a)) {
184	result.parts.fraction = a.parts.fraction;	184	result.parts.fraction = a.parts.fraction;
185	result.parts.exp = a.parts.exp;	185	result.parts.exp = a.parts.exp;
186	return result;	186	return result;
187	};	187	};
188	if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */	188	if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
189	result.parts.fraction = b.parts.fraction;	189	result.parts.fraction = b.parts.fraction;
190	result.parts.exp = b.parts.exp;	190	result.parts.exp = b.parts.exp;
191	return result;	191	return result;
192	};	192	};
193	/* set NaN as result */	193	/* set NaN as result */
194	result.binary = FLOAT64_NAN;	194	result.binary = FLOAT64_NAN;
195	return result;	195	return result;
196	};	196	};
197		197
198	if (isFloat64Infinity(a)) {	198	if (isFloat64Infinity(a)) {
199	if (isFloat64Zero(b)) {	199	if (isFloat64Zero(b)) {
200	/* FIXME: zero * infinity */	200	/* FIXME: zero * infinity */
201	result.binary = FLOAT64_NAN;	201	result.binary = FLOAT64_NAN;
202	return result;	202	return result;
203	}	203	}
204	result.parts.fraction = a.parts.fraction;	204	result.parts.fraction = a.parts.fraction;
205	result.parts.exp = a.parts.exp;	205	result.parts.exp = a.parts.exp;
206	return result;	206	return result;
207	}	207	}
208		208
209	if (isFloat64Infinity(b)) {	209	if (isFloat64Infinity(b)) {
210	if (isFloat64Zero(a)) {	210	if (isFloat64Zero(a)) {
211	/* FIXME: zero * infinity */	211	/* FIXME: zero * infinity */
212	result.binary = FLOAT64_NAN;	212	result.binary = FLOAT64_NAN;
213	return result;	213	return result;
214	}	214	}
215	result.parts.fraction = b.parts.fraction;	215	result.parts.fraction = b.parts.fraction;
216	result.parts.exp = b.parts.exp;	216	result.parts.exp = b.parts.exp;
217	return result;	217	return result;
218	}	218	}
219		219
220	/* exp is signed so we can easy detect underflow */	220	/* exp is signed so we can easy detect underflow */
221	exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;	221	exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
222		222
223	frac1 = a.parts.fraction;	223	frac1 = a.parts.fraction;
224		224
225	if (a.parts.exp > 0) {	225	if (a.parts.exp > 0) {
226	frac1 \|= FLOAT64_HIDDEN_BIT_MASK;	226	frac1 \|= FLOAT64_HIDDEN_BIT_MASK;
227	} else {	227	} else {
228	++exp;	228	++exp;
229	};	229	};
230		230
231	frac2 = b.parts.fraction;	231	frac2 = b.parts.fraction;
232		232
233	if (b.parts.exp > 0) {	233	if (b.parts.exp > 0) {
234	frac2 \|= FLOAT64_HIDDEN_BIT_MASK;	234	frac2 \|= FLOAT64_HIDDEN_BIT_MASK;
235	} else {	235	} else {
236	++exp;	236	++exp;
237	};	237	};
238		238
239	frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);	239	frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
240	frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);	240	frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
241		241
242	mul64integers(frac1, frac2, &frac1, &frac2);	242	mul64integers(frac1, frac2, &frac1, &frac2);
243		243
244	frac2 \|= (frac1 != 0);	244	frac2 \|= (frac1 != 0);
245	if (frac2 & (0x1ll << 62)) {	245	if (frac2 & (0x1ll << 62)) {
246	frac2 <<= 1;	246	frac2 <<= 1;
247	exp--;	247	exp--;
248	}	248	}
249		249
250	result = finishFloat64(exp, frac2, result.parts.sign);	250	result = finishFloat64(exp, frac2, result.parts.sign);
251	return result;	251	return result;
252	}	252	}
253		253
254	/** Multiply two 64 bit numbers and return result in two parts	254	/** Multiply two 64 bit numbers and return result in two parts
255	* @param a first operand	255	* @param a first operand
256	* @param b second operand	256	* @param b second operand
257	* @param lo lower part from result	257	* @param lo lower part from result
258	* @param hi higher part of result	258	* @param hi higher part of result
259	*/	259	*/
260	void mul64integers(__u64 a,__u64 b, __u64 lo, __u64 hi)	260	void mul64integers(uint64_t a,uint64_t b, uint64_t lo, uint64_t hi)
261	{	261	{
262	__u64 low, high, middle1, middle2;	262	uint64_t low, high, middle1, middle2;
263	__u32 alow, blow;	263	uint32_t alow, blow;
264		264
265	alow = a & 0xFFFFFFFF;	265	alow = a & 0xFFFFFFFF;
266	blow = b & 0xFFFFFFFF;	266	blow = b & 0xFFFFFFFF;
267		267
268	a >>= 32;	268	a >>= 32;
269	b >>= 32;	269	b >>= 32;
270		270
271	low = ((__u64)alow) * blow;	271	low = ((uint64_t)alow) * blow;
272	middle1 = a * blow;	272	middle1 = a * blow;
273	middle2 = alow * b;	273	middle2 = alow * b;
274	high = a * b;	274	high = a * b;
275		275
276	middle1 += middle2;	276	middle1 += middle2;
277	high += (((__u64)(middle1 < middle2)) << 32) + (middle1 >> 32);	277	high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
278	middle1 <<= 32;	278	middle1 <<= 32;
279	low += middle1;	279	low += middle1;
280	high += (low < middle1);	280	high += (low < middle1);
281	*lo = low;	281	*lo = low;
282	*hi = high;	282	*hi = high;
283		283
284	return;	284	return;
285	}	285	}
286		286
287		287
288		288

Subversion Repositories HelenOS-historic

(root)/uspace/trunk/softfloat/generic/mul.c – Rev 829 → 1031