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

Subversion Repositories HelenOS-historic

(root)/uspace/trunk/softfloat/generic/div.c – Rev 828 → 829