Atlas - e_fmod.c
Home / ext / SDL / src / libm
Lines: 1 | Size: 3608 bytes
[Download] [Show on GitHub] [Search similar files] [Raw] [Raw (proxy)] 
[FILE BEGIN]
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * __ieee754_fmod(x,y)
15 * Return x mod y in exact arithmetic
16 * Method: shift and subtract
17 */
18
19#include "math_libm.h"
20#include "math_private.h"
21
22static const double one = 1.0, Zero[] = {0.0, -0.0,};
23
24double attribute_hidden __ieee754_fmod(double x, double y)
25{
26	int32_t n,hx,hy,hz,ix,iy,sx,i;
27	u_int32_t lx,ly,lz;
28
29	EXTRACT_WORDS(hx,lx,x);
30	EXTRACT_WORDS(hy,ly,y);
31	sx = hx&0x80000000;		/* sign of x */
32	hx ^=sx;		/* |x| */
33	hy &= 0x7fffffff;	/* |y| */
34
35    /* purge off exception values */
36	if((hy|ly)==0||(hx>=0x7ff00000)||	/* y=0,or x not finite */
37	  ((hy|((ly|-(int32_t)ly)>>31))>0x7ff00000))	/* or y is NaN */
38	    return (x*y)/(x*y);
39	if(hx<=hy) {
40	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
41	    if(lx==ly)
42		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
43	}
44
45    /* determine ix = ilogb(x) */
46	if(hx<0x00100000) {	/* subnormal x */
47	    if(hx==0) {
48		for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
49	    } else {
50		for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
51	    }
52	} else ix = (hx>>20)-1023;
53
54    /* determine iy = ilogb(y) */
55	if(hy<0x00100000) {	/* subnormal y */
56	    if(hy==0) {
57		for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
58	    } else {
59		for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
60	    }
61	} else iy = (hy>>20)-1023;
62
63    /* set up {hx,lx}, {hy,ly} and align y to x */
64	if(ix >= -1022)
65	    hx = 0x00100000|(0x000fffff&hx);
66	else {		/* subnormal x, shift x to normal */
67	    n = -1022-ix;
68	    if(n<=31) {
69	        hx = (hx<<n)|(lx>>(32-n));
70	        lx <<= n;
71	    } else {
72		hx = lx<<(n-32);
73		lx = 0;
74	    }
75	}
76	if(iy >= -1022)
77	    hy = 0x00100000|(0x000fffff&hy);
78	else {		/* subnormal y, shift y to normal */
79	    n = -1022-iy;
80	    if(n<=31) {
81	        hy = (hy<<n)|(ly>>(32-n));
82	        ly <<= n;
83	    } else {
84		hy = ly<<(n-32);
85		ly = 0;
86	    }
87	}
88
89    /* fix point fmod */
90	n = ix - iy;
91	while(n--) {
92	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
93	    if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
94	    else {
95	    	if((hz|lz)==0) 		/* return sign(x)*0 */
96		    return Zero[(u_int32_t)sx>>31];
97	    	hx = hz+hz+(lz>>31); lx = lz+lz;
98	    }
99	}
100	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
101	if(hz>=0) {hx=hz;lx=lz;}
102
103    /* convert back to floating value and restore the sign */
104	if((hx|lx)==0) 			/* return sign(x)*0 */
105	    return Zero[(u_int32_t)sx>>31];
106	while(hx<0x00100000) {		/* normalize x */
107	    hx = hx+hx+(lx>>31); lx = lx+lx;
108	    iy -= 1;
109	}
110	if(iy>= -1022) {	/* normalize output */
111	    hx = ((hx-0x00100000)|((iy+1023)<<20));
112	    INSERT_WORDS(x,hx|sx,lx);
113	} else {		/* subnormal output */
114	    n = -1022 - iy;
115	    if(n<=20) {
116		lx = (lx>>n)|((u_int32_t)hx<<(32-n));
117		hx >>= n;
118	    } else if (n<=31) {
119		lx = (hx<<(32-n))|(lx>>n); hx = sx;
120	    } else {
121		lx = hx>>(n-32); hx = sx;
122	    }
123	    INSERT_WORDS(x,hx|sx,lx);
124	    x *= one;		/* create necessary signal */
125	}
126	return x;		/* exact output */
127}
128
129/*
130 * wrapper fmod(x,y)
131 */
132#ifndef _IEEE_LIBM
133double fmod(double x, double y)
134{
135	double z = __ieee754_fmod(x, y);
136	if (_LIB_VERSION == _IEEE_ || isnan(y) || isnan(x))
137		return z;
138	if (y == 0.0)
139		return __kernel_standard(x, y, 27); /* fmod(x,0) */
140	return z;
141}
142#else
143strong_alias(__ieee754_fmod, fmod)
144#endif
145libm_hidden_def(fmod)
146
[FILE END]
(C) 2025 0x4248
(C) 2025 4248 Media and 4248 Systems, All part of 0x4248
See LICENCE files for more information. Not all files are by 0x4248 always check Licencing.