WebSVN – HelenOS-historic – Diff – /uspace/trunk/softfloat/generic/div.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.
 */
 
/** @addtogroup softfloat  
 * @{
 */
/** @file
 */
 
#include<sftypes.h>
#include<add.h>
#include<div.h>
#include<comparison.h>
#include<mul.h>
#include<common.h>
 
 
float32 divFloat32(float32 a, float32 b)
{
    float32 result;
    int32_t aexp, bexp, cexp;
    uint64_t afrac, bfrac, cfrac;
   
    result.parts.sign = a.parts.sign ^ b.parts.sign;
   
    if (isFloat32NaN(a)) {
        if (isFloat32SigNaN(a)) {
            /*FIXME: SigNaN*/
        }
        /*NaN*/
        return a;
    }
   
    if (isFloat32NaN(b)) {
        if (isFloat32SigNaN(b)) {
            /*FIXME: SigNaN*/
        }
        /*NaN*/
        return b;
    }
   
    if (isFloat32Infinity(a)) {
        if (isFloat32Infinity(b)) {
            /*FIXME: inf / inf */
            result.binary = FLOAT32_NAN;
            return result;
        }
        /* inf / num */
        result.parts.exp = a.parts.exp;
        result.parts.fraction = a.parts.fraction;
        return result;
    }
 
    if (isFloat32Infinity(b)) {
        if (isFloat32Zero(a)) {
            /* FIXME 0 / inf */
            result.parts.exp = 0;
            result.parts.fraction = 0;
            return result;
        }
        /* FIXME: num / inf*/
        result.parts.exp = 0;
        result.parts.fraction = 0;
        return result;
    }
   
    if (isFloat32Zero(b)) {
        if (isFloat32Zero(a)) {
            /*FIXME: 0 / 0*/
            result.binary = FLOAT32_NAN;
            return result;
        }
        /* FIXME: division by zero */
        result.parts.exp = 0;
        result.parts.fraction = 0;
        return result;
    }
 
   
    afrac = a.parts.fraction;
    aexp = a.parts.exp;
    bfrac = b.parts.fraction;
    bexp = b.parts.exp;
   
    /* denormalized numbers */
    if (aexp == 0) {
        if (afrac == 0) {
        result.parts.exp = 0;
        result.parts.fraction = 0;
        return result;
        }
        /* normalize it*/
       
        afrac <<= 1;
            /* afrac is nonzero => it must stop */ 
        while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
            afrac <<= 1;
            aexp--;
        }
    }
 
    if (bexp == 0) {
        bfrac <<= 1;
            /* bfrac is nonzero => it must stop */ 
        while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
            bfrac <<= 1;
            bexp--;
        }
    }
 
    afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
    bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
 
    if ( bfrac <= (afrac << 1) ) {
        afrac >>= 1;
        aexp++;
    }
   
    cexp = aexp - bexp + FLOAT32_BIAS - 2;
   
    cfrac = (afrac << 32) / bfrac;
    if ((  cfrac & 0x3F ) == 0) {
        cfrac |= ( bfrac * cfrac != afrac << 32 );
    }
   
    /* pack and round */
   
    /* find first nonzero digit and shift result and detect possibly underflow */
    while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
        cexp--;
        cfrac <<= 1;
            /* TODO: fix underflow */
    };
   
    cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
   
    if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
        ++cexp;
        cfrac >>= 1;
        }  
 
    /* check overflow */
    if (cexp >= FLOAT32_MAX_EXPONENT ) {
        /* FIXME: overflow, return infinity */
        result.parts.exp = FLOAT32_MAX_EXPONENT;
        result.parts.fraction = 0;
        return result;
    }
 
    if (cexp < 0) {
        /* FIXME: underflow */
        result.parts.exp = 0;
        if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
            result.parts.fraction = 0;
            return result;
        }
        cfrac >>= 1;
        while (cexp < 0) {
            cexp ++;
            cfrac >>= 1;
        }
       
    } else {
        result.parts.exp = (uint32_t)cexp;
    }
   
    result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
   
    return result; 
}
 
float64 divFloat64(float64 a, float64 b)
{
    float64 result;
    int64_t aexp, bexp, cexp;
    uint64_t afrac, bfrac, cfrac;
    uint64_t remlo, remhi;
   
    result.parts.sign = a.parts.sign ^ b.parts.sign;
   
    if (isFloat64NaN(a)) {
       
        if (isFloat64SigNaN(b)) {
            /*FIXME: SigNaN*/
            return b;
        }
       
        if (isFloat64SigNaN(a)) {
            /*FIXME: SigNaN*/
        }
        /*NaN*/
        return a;
    }
   
    if (isFloat64NaN(b)) {
        if (isFloat64SigNaN(b)) {
            /*FIXME: SigNaN*/
        }
        /*NaN*/
        return b;
    }
   
    if (isFloat64Infinity(a)) {
        if (isFloat64Infinity(b) || isFloat64Zero(b)) {
            /*FIXME: inf / inf */
            result.binary = FLOAT64_NAN;
            return result;
        }
        /* inf / num */
        result.parts.exp = a.parts.exp;
        result.parts.fraction = a.parts.fraction;
        return result;
    }
 
    if (isFloat64Infinity(b)) {
        if (isFloat64Zero(a)) {
            /* FIXME 0 / inf */
            result.parts.exp = 0;
            result.parts.fraction = 0;
            return result;
        }
        /* FIXME: num / inf*/
        result.parts.exp = 0;
        result.parts.fraction = 0;
        return result;
    }
   
    if (isFloat64Zero(b)) {
        if (isFloat64Zero(a)) {
            /*FIXME: 0 / 0*/
            result.binary = FLOAT64_NAN;
            return result;
        }
        /* FIXME: division by zero */
        result.parts.exp = 0;
        result.parts.fraction = 0;
        return result;
    }
 
   
    afrac = a.parts.fraction;
    aexp = a.parts.exp;
    bfrac = b.parts.fraction;
    bexp = b.parts.exp;
   
    /* denormalized numbers */
    if (aexp == 0) {
        if (afrac == 0) {
            result.parts.exp = 0;
            result.parts.fraction = 0;
            return result;
        }
        /* normalize it*/
       
        aexp++;
            /* afrac is nonzero => it must stop */ 
        while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
            afrac <<= 1;
            aexp--;
        }
    }
 
    if (bexp == 0) {
        bexp++;
            /* bfrac is nonzero => it must stop */ 
        while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
            bfrac <<= 1;
            bexp--;
        }
    }
 
    afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
    bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
 
    if ( bfrac <= (afrac << 1) ) {
        afrac >>= 1;
        aexp++;
    }
   
    cexp = aexp - bexp + FLOAT64_BIAS - 2;
   
    cfrac = divFloat64estim(afrac, bfrac);
   
    if ((  cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
        mul64integers( bfrac, cfrac, &remlo, &remhi);
        /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ 
        remhi = afrac - remhi - ( remlo > 0);
        remlo = - remlo;
       
        while ((int64_t) remhi < 0) {
            cfrac--;
            remlo += bfrac;
            remhi += ( remlo < bfrac );
        }
        cfrac |= ( remlo != 0 );
    }
   
    /* round and shift */
    result = finishFloat64(cexp, cfrac, result.parts.sign);
    return result;
 
}
 
uint64_t divFloat64estim(uint64_t a, uint64_t b)
{
    uint64_t bhi;
    uint64_t remhi, remlo;
    uint64_t result;
   
    if ( b <= a ) {
        return 0xFFFFFFFFFFFFFFFFull;
    }
   
    bhi = b >> 32;
    result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
    mul64integers(b, result, &remlo, &remhi);
   
    remhi = a - remhi - (remlo > 0);
    remlo = - remlo;
 
    b <<= 32;
    while ( (int64_t) remhi < 0 ) {
            result -= 0x1ll << 32; 
            remlo += b;
            remhi += bhi + ( remlo < b );
        }
    remhi = (remhi << 32) | (remlo >> 32);
    if (( bhi << 32) <= remhi) {
        result |= 0xFFFFFFFF;
    } else {
        result |= remhi / bhi;
    }
   
   
    return result;
}
 
 
/** @}
 */
 
 

Rev 1657	Rev 1740
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	/** @addtogroup softfloat	29	/** @addtogroup softfloat
30	* @{	30	* @{
31	*/	31	*/
32	/** @file	32	/** @file
33	*/	33	*/
34		34
35	#include<sftypes.h>	35	#include<sftypes.h>
36	#include<add.h>	36	#include<add.h>
37	#include<div.h>	37	#include<div.h>
38	#include<comparison.h>	38	#include<comparison.h>
39	#include<mul.h>	39	#include<mul.h>
40	#include<common.h>	40	#include<common.h>
41		41
42		42
43	float32 divFloat32(float32 a, float32 b)	43	float32 divFloat32(float32 a, float32 b)
44	{	44	{
45	float32 result;	45	float32 result;
46	int32_t aexp, bexp, cexp;	46	int32_t aexp, bexp, cexp;
47	uint64_t afrac, bfrac, cfrac;	47	uint64_t afrac, bfrac, cfrac;
48		48
49	result.parts.sign = a.parts.sign ^ b.parts.sign;	49	result.parts.sign = a.parts.sign ^ b.parts.sign;
50		50
51	if (isFloat32NaN(a)) {	51	if (isFloat32NaN(a)) {
52	if (isFloat32SigNaN(a)) {	52	if (isFloat32SigNaN(a)) {
53	/FIXME: SigNaN/	53	/FIXME: SigNaN/
54	}	54	}
55	/NaN/	55	/NaN/
56	return a;	56	return a;
57	}	57	}
58		58
59	if (isFloat32NaN(b)) {	59	if (isFloat32NaN(b)) {
60	if (isFloat32SigNaN(b)) {	60	if (isFloat32SigNaN(b)) {
61	/FIXME: SigNaN/	61	/FIXME: SigNaN/
62	}	62	}
63	/NaN/	63	/NaN/
64	return b;	64	return b;
65	}	65	}
66		66
67	if (isFloat32Infinity(a)) {	67	if (isFloat32Infinity(a)) {
68	if (isFloat32Infinity(b)) {	68	if (isFloat32Infinity(b)) {
69	/FIXME: inf / inf /	69	/FIXME: inf / inf /
70	result.binary = FLOAT32_NAN;	70	result.binary = FLOAT32_NAN;
71	return result;	71	return result;
72	}	72	}
73	/* inf / num */	73	/* inf / num */
74	result.parts.exp = a.parts.exp;	74	result.parts.exp = a.parts.exp;
75	result.parts.fraction = a.parts.fraction;	75	result.parts.fraction = a.parts.fraction;
76	return result;	76	return result;
77	}	77	}
78		78
79	if (isFloat32Infinity(b)) {	79	if (isFloat32Infinity(b)) {
80	if (isFloat32Zero(a)) {	80	if (isFloat32Zero(a)) {
81	/* FIXME 0 / inf */	81	/* FIXME 0 / inf */
82	result.parts.exp = 0;	82	result.parts.exp = 0;
83	result.parts.fraction = 0;	83	result.parts.fraction = 0;
84	return result;	84	return result;
85	}	85	}
86	/* FIXME: num / inf*/	86	/* FIXME: num / inf*/
87	result.parts.exp = 0;	87	result.parts.exp = 0;
88	result.parts.fraction = 0;	88	result.parts.fraction = 0;
89	return result;	89	return result;
90	}	90	}
91		91
92	if (isFloat32Zero(b)) {	92	if (isFloat32Zero(b)) {
93	if (isFloat32Zero(a)) {	93	if (isFloat32Zero(a)) {
94	/FIXME: 0 / 0/	94	/FIXME: 0 / 0/
95	result.binary = FLOAT32_NAN;	95	result.binary = FLOAT32_NAN;
96	return result;	96	return result;
97	}	97	}
98	/* FIXME: division by zero */	98	/* FIXME: division by zero */
99	result.parts.exp = 0;	99	result.parts.exp = 0;
100	result.parts.fraction = 0;	100	result.parts.fraction = 0;
101	return result;	101	return result;
102	}	102	}
103		103
104		104
105	afrac = a.parts.fraction;	105	afrac = a.parts.fraction;
106	aexp = a.parts.exp;	106	aexp = a.parts.exp;
107	bfrac = b.parts.fraction;	107	bfrac = b.parts.fraction;
108	bexp = b.parts.exp;	108	bexp = b.parts.exp;
109		109
110	/* denormalized numbers */	110	/* denormalized numbers */
111	if (aexp == 0) {	111	if (aexp == 0) {
112	if (afrac == 0) {	112	if (afrac == 0) {
113	result.parts.exp = 0;	113	result.parts.exp = 0;
114	result.parts.fraction = 0;	114	result.parts.fraction = 0;
115	return result;	115	return result;
116	}	116	}
117	/* normalize it*/	117	/* normalize it*/
118		118
119	afrac <<= 1;	119	afrac <<= 1;
120	/* afrac is nonzero => it must stop */	120	/* afrac is nonzero => it must stop */
121	while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {	121	while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
122	afrac <<= 1;	122	afrac <<= 1;
123	aexp--;	123	aexp--;
124	}	124	}
125	}	125	}
126		126
127	if (bexp == 0) {	127	if (bexp == 0) {
128	bfrac <<= 1;	128	bfrac <<= 1;
129	/* bfrac is nonzero => it must stop */	129	/* bfrac is nonzero => it must stop */
130	while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {	130	while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
131	bfrac <<= 1;	131	bfrac <<= 1;
132	bexp--;	132	bexp--;
133	}	133	}
134	}	134	}
135		135
136	afrac = (afrac \| FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );	136	afrac = (afrac \| FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
137	bfrac = (bfrac \| FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );	137	bfrac = (bfrac \| FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
138		138
139	if ( bfrac <= (afrac << 1) ) {	139	if ( bfrac <= (afrac << 1) ) {
140	afrac >>= 1;	140	afrac >>= 1;
141	aexp++;	141	aexp++;
142	}	142	}
143		143
144	cexp = aexp - bexp + FLOAT32_BIAS - 2;	144	cexp = aexp - bexp + FLOAT32_BIAS - 2;
145		145
146	cfrac = (afrac << 32) / bfrac;	146	cfrac = (afrac << 32) / bfrac;
147	if (( cfrac & 0x3F ) == 0) {	147	if (( cfrac & 0x3F ) == 0) {
148	cfrac \|= ( bfrac * cfrac != afrac << 32 );	148	cfrac \|= ( bfrac * cfrac != afrac << 32 );
149	}	149	}
150		150
151	/* pack and round */	151	/* pack and round */
152		152
153	/* find first nonzero digit and shift result and detect possibly underflow */	153	/* find first nonzero digit and shift result and detect possibly underflow */
154	while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {	154	while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
155	cexp--;	155	cexp--;
156	cfrac <<= 1;	156	cfrac <<= 1;
157	/* TODO: fix underflow */	157	/* TODO: fix underflow */
158	};	158	};
159		159
160	cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/	160	cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
161		161
162	if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {	162	if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
163	++cexp;	163	++cexp;
164	cfrac >>= 1;	164	cfrac >>= 1;
165	}	165	}
166		166
167	/* check overflow */	167	/* check overflow */
168	if (cexp >= FLOAT32_MAX_EXPONENT ) {	168	if (cexp >= FLOAT32_MAX_EXPONENT ) {
169	/* FIXME: overflow, return infinity */	169	/* FIXME: overflow, return infinity */
170	result.parts.exp = FLOAT32_MAX_EXPONENT;	170	result.parts.exp = FLOAT32_MAX_EXPONENT;
171	result.parts.fraction = 0;	171	result.parts.fraction = 0;
172	return result;	172	return result;
173	}	173	}
174		174
175	if (cexp < 0) {	175	if (cexp < 0) {
176	/* FIXME: underflow */	176	/* FIXME: underflow */
177	result.parts.exp = 0;	177	result.parts.exp = 0;
178	if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {	178	if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
179	result.parts.fraction = 0;	179	result.parts.fraction = 0;
180	return result;	180	return result;
181	}	181	}
182	cfrac >>= 1;	182	cfrac >>= 1;
183	while (cexp < 0) {	183	while (cexp < 0) {
184	cexp ++;	184	cexp ++;
185	cfrac >>= 1;	185	cfrac >>= 1;
186	}	186	}
187		187
188	} else {	188	} else {
189	result.parts.exp = (uint32_t)cexp;	189	result.parts.exp = (uint32_t)cexp;
190	}	190	}
191		191
192	result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));	192	result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
193		193
194	return result;	194	return result;
195	}	195	}
196		196
197	float64 divFloat64(float64 a, float64 b)	197	float64 divFloat64(float64 a, float64 b)
198	{	198	{
199	float64 result;	199	float64 result;
200	int64_t aexp, bexp, cexp;	200	int64_t aexp, bexp, cexp;
201	uint64_t afrac, bfrac, cfrac;	201	uint64_t afrac, bfrac, cfrac;
202	uint64_t remlo, remhi;	202	uint64_t remlo, remhi;
203		203
204	result.parts.sign = a.parts.sign ^ b.parts.sign;	204	result.parts.sign = a.parts.sign ^ b.parts.sign;
205		205
206	if (isFloat64NaN(a)) {	206	if (isFloat64NaN(a)) {
207		207
208	if (isFloat64SigNaN(b)) {	208	if (isFloat64SigNaN(b)) {
209	/FIXME: SigNaN/	209	/FIXME: SigNaN/
210	return b;	210	return b;
211	}	211	}
212		212
213	if (isFloat64SigNaN(a)) {	213	if (isFloat64SigNaN(a)) {
214	/FIXME: SigNaN/	214	/FIXME: SigNaN/
215	}	215	}
216	/NaN/	216	/NaN/
217	return a;	217	return a;
218	}	218	}
219		219
220	if (isFloat64NaN(b)) {	220	if (isFloat64NaN(b)) {
221	if (isFloat64SigNaN(b)) {	221	if (isFloat64SigNaN(b)) {
222	/FIXME: SigNaN/	222	/FIXME: SigNaN/
223	}	223	}
224	/NaN/	224	/NaN/
225	return b;	225	return b;
226	}	226	}
227		227
228	if (isFloat64Infinity(a)) {	228	if (isFloat64Infinity(a)) {
229	if (isFloat64Infinity(b) \|\| isFloat64Zero(b)) {	229	if (isFloat64Infinity(b) \|\| isFloat64Zero(b)) {
230	/FIXME: inf / inf /	230	/FIXME: inf / inf /
231	result.binary = FLOAT64_NAN;	231	result.binary = FLOAT64_NAN;
232	return result;	232	return result;
233	}	233	}
234	/* inf / num */	234	/* inf / num */
235	result.parts.exp = a.parts.exp;	235	result.parts.exp = a.parts.exp;
236	result.parts.fraction = a.parts.fraction;	236	result.parts.fraction = a.parts.fraction;
237	return result;	237	return result;
238	}	238	}
239		239
240	if (isFloat64Infinity(b)) {	240	if (isFloat64Infinity(b)) {
241	if (isFloat64Zero(a)) {	241	if (isFloat64Zero(a)) {
242	/* FIXME 0 / inf */	242	/* FIXME 0 / inf */
243	result.parts.exp = 0;	243	result.parts.exp = 0;
244	result.parts.fraction = 0;	244	result.parts.fraction = 0;
245	return result;	245	return result;
246	}	246	}
247	/* FIXME: num / inf*/	247	/* FIXME: num / inf*/
248	result.parts.exp = 0;	248	result.parts.exp = 0;
249	result.parts.fraction = 0;	249	result.parts.fraction = 0;
250	return result;	250	return result;
251	}	251	}
252		252
253	if (isFloat64Zero(b)) {	253	if (isFloat64Zero(b)) {
254	if (isFloat64Zero(a)) {	254	if (isFloat64Zero(a)) {
255	/FIXME: 0 / 0/	255	/FIXME: 0 / 0/
256	result.binary = FLOAT64_NAN;	256	result.binary = FLOAT64_NAN;
257	return result;	257	return result;
258	}	258	}
259	/* FIXME: division by zero */	259	/* FIXME: division by zero */
260	result.parts.exp = 0;	260	result.parts.exp = 0;
261	result.parts.fraction = 0;	261	result.parts.fraction = 0;
262	return result;	262	return result;
263	}	263	}
264		264
265		265
266	afrac = a.parts.fraction;	266	afrac = a.parts.fraction;
267	aexp = a.parts.exp;	267	aexp = a.parts.exp;
268	bfrac = b.parts.fraction;	268	bfrac = b.parts.fraction;
269	bexp = b.parts.exp;	269	bexp = b.parts.exp;
270		270
271	/* denormalized numbers */	271	/* denormalized numbers */
272	if (aexp == 0) {	272	if (aexp == 0) {
273	if (afrac == 0) {	273	if (afrac == 0) {
274	result.parts.exp = 0;	274	result.parts.exp = 0;
275	result.parts.fraction = 0;	275	result.parts.fraction = 0;
276	return result;	276	return result;
277	}	277	}
278	/* normalize it*/	278	/* normalize it*/
279		279
280	aexp++;	280	aexp++;
281	/* afrac is nonzero => it must stop */	281	/* afrac is nonzero => it must stop */
282	while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {	282	while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
283	afrac <<= 1;	283	afrac <<= 1;
284	aexp--;	284	aexp--;
285	}	285	}
286	}	286	}
287		287
288	if (bexp == 0) {	288	if (bexp == 0) {
289	bexp++;	289	bexp++;
290	/* bfrac is nonzero => it must stop */	290	/* bfrac is nonzero => it must stop */
291	while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {	291	while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
292	bfrac <<= 1;	292	bfrac <<= 1;
293	bexp--;	293	bexp--;
294	}	294	}
295	}	295	}
296		296
297	afrac = (afrac \| FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );	297	afrac = (afrac \| FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
298	bfrac = (bfrac \| FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);	298	bfrac = (bfrac \| FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
299		299
300	if ( bfrac <= (afrac << 1) ) {	300	if ( bfrac <= (afrac << 1) ) {
301	afrac >>= 1;	301	afrac >>= 1;
302	aexp++;	302	aexp++;
303	}	303	}
304		304
305	cexp = aexp - bexp + FLOAT64_BIAS - 2;	305	cexp = aexp - bexp + FLOAT64_BIAS - 2;
306		306
307	cfrac = divFloat64estim(afrac, bfrac);	307	cfrac = divFloat64estim(afrac, bfrac);
308		308
309	if (( cfrac & 0x1FF ) <= 2) { /FIXME:?? /	309	if (( cfrac & 0x1FF ) <= 2) { /FIXME:?? /
310	mul64integers( bfrac, cfrac, &remlo, &remhi);	310	mul64integers( bfrac, cfrac, &remlo, &remhi);
311	/* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/	311	/* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/
312	remhi = afrac - remhi - ( remlo > 0);	312	remhi = afrac - remhi - ( remlo > 0);
313	remlo = - remlo;	313	remlo = - remlo;
314		314
315	while ((int64_t) remhi < 0) {	315	while ((int64_t) remhi < 0) {
316	cfrac--;	316	cfrac--;
317	remlo += bfrac;	317	remlo += bfrac;
318	remhi += ( remlo < bfrac );	318	remhi += ( remlo < bfrac );
319	}	319	}
320	cfrac \|= ( remlo != 0 );	320	cfrac \|= ( remlo != 0 );
321	}	321	}
322		322
323	/* round and shift */	323	/* round and shift */
324	result = finishFloat64(cexp, cfrac, result.parts.sign);	324	result = finishFloat64(cexp, cfrac, result.parts.sign);
325	return result;	325	return result;
326		326
327	}	327	}
328		328
329	uint64_t divFloat64estim(uint64_t a, uint64_t b)	329	uint64_t divFloat64estim(uint64_t a, uint64_t b)
330	{	330	{
331	uint64_t bhi;	331	uint64_t bhi;
332	uint64_t remhi, remlo;	332	uint64_t remhi, remlo;
333	uint64_t result;	333	uint64_t result;
334		334
335	if ( b <= a ) {	335	if ( b <= a ) {
336	return 0xFFFFFFFFFFFFFFFFull;	336	return 0xFFFFFFFFFFFFFFFFull;
337	}	337	}
338		338
339	bhi = b >> 32;	339	bhi = b >> 32;
340	result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;	340	result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
341	mul64integers(b, result, &remlo, &remhi);	341	mul64integers(b, result, &remlo, &remhi);
342		342
343	remhi = a - remhi - (remlo > 0);	343	remhi = a - remhi - (remlo > 0);
344	remlo = - remlo;	344	remlo = - remlo;
345		345
346	b <<= 32;	346	b <<= 32;
347	while ( (int64_t) remhi < 0 ) {	347	while ( (int64_t) remhi < 0 ) {
348	result -= 0x1ll << 32;	348	result -= 0x1ll << 32;
349	remlo += b;	349	remlo += b;
350	remhi += bhi + ( remlo < b );	350	remhi += bhi + ( remlo < b );
351	}	351	}
352	remhi = (remhi << 32) \| (remlo >> 32);	352	remhi = (remhi << 32) \| (remlo >> 32);
353	if (( bhi << 32) <= remhi) {	353	if (( bhi << 32) <= remhi) {
354	result \|= 0xFFFFFFFF;	354	result \|= 0xFFFFFFFF;
355	} else {	355	} else {
356	result \|= remhi / bhi;	356	result \|= remhi / bhi;
357	}	357	}
358		358
359		359
360	return result;	360	return result;
361	}	361	}
362		362
363		-
364	/** @}	363	/** @}
365	*/	364	*/
366		-
367		365

Subversion Repositories HelenOS-historic

(root)/uspace/trunk/softfloat/generic/div.c – Rev 1657 → 1740