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

Subversion Repositories HelenOS-historic

(root)/uspace/trunk/softfloat/generic/div.c – Rev 835 → 1031