1 /* 2 * This file is part of the MicroPython project, http://micropython.org/ 3 * 4 * These math functions are taken from newlib-nano-2, the newlib/libm/math 5 * directory, available from https://github.com/32bitmicro/newlib-nano-2. 6 * 7 * Appropriate copyright headers are reproduced below. 8 */ 9 10 /* ef_rem_pio2.c -- float version of e_rem_pio2.c 11 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 12 */ 13 14 /* 15 * ==================================================== 16 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 17 * 18 * Developed at SunPro, a Sun Microsystems, Inc. business. 19 * Permission to use, copy, modify, and distribute this 20 * software is freely granted, provided that this notice 21 * is preserved. 22 * ==================================================== 23 * 24 */ 25 26 /* __ieee754_rem_pio2f(x,y) 27 * 28 * return the remainder of x rem pi/2 in y[0]+y[1] 29 * use __kernel_rem_pio2f() 30 */ 31 32 #include "fdlibm.h" 33 34 /* 35 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 36 */ 37 #ifdef __STDC__ 38 static const __uint8_t two_over_pi[] = { 39 #else 40 static __uint8_t two_over_pi[] = { 41 #endif 42 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC, 43 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62, 44 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63, 45 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A, 46 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09, 47 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29, 48 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44, 49 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41, 50 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C, 51 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8, 52 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11, 53 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF, 54 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E, 55 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5, 56 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92, 57 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08, 58 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0, 59 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3, 60 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85, 61 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80, 62 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA, 63 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B, 64 }; 65 66 /* This array is like the one in e_rem_pio2.c, but the numbers are 67 single precision and the last 8 bits are forced to 0. */ 68 #ifdef __STDC__ 69 static const __int32_t npio2_hw[] = { 70 #else 71 static __int32_t npio2_hw[] = { 72 #endif 73 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00, 74 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00, 75 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100, 76 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00, 77 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00, 78 0x4242c700, 0x42490f00 79 }; 80 81 /* 82 * invpio2: 24 bits of 2/pi 83 * pio2_1: first 17 bit of pi/2 84 * pio2_1t: pi/2 - pio2_1 85 * pio2_2: second 17 bit of pi/2 86 * pio2_2t: pi/2 - (pio2_1+pio2_2) 87 * pio2_3: third 17 bit of pi/2 88 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 89 */ 90 91 #ifdef __STDC__ 92 static const float 93 #else 94 static float 95 #endif 96 zero = 0.0000000000e+00f, /* 0x00000000 */ 97 half = 5.0000000000e-01f, /* 0x3f000000 */ 98 two8 = 2.5600000000e+02f, /* 0x43800000 */ 99 invpio2 = 6.3661980629e-01f, /* 0x3f22f984 */ 100 pio2_1 = 1.5707855225e+00f, /* 0x3fc90f80 */ 101 pio2_1t = 1.0804334124e-05f, /* 0x37354443 */ 102 pio2_2 = 1.0804273188e-05f, /* 0x37354400 */ 103 pio2_2t = 6.0770999344e-11f, /* 0x2e85a308 */ 104 pio2_3 = 6.0770943833e-11f, /* 0x2e85a300 */ 105 pio2_3t = 6.1232342629e-17f; /* 0x248d3132 */ 106 107 #ifdef __STDC__ __ieee754_rem_pio2f(float x,float * y)108 __int32_t __ieee754_rem_pio2f(float x, float *y) 109 #else 110 __int32_t __ieee754_rem_pio2f(x,y) 111 float x,y[]; 112 #endif 113 { 114 float z,w,t,r,fn; 115 float tx[3]; 116 __int32_t i,j,n,ix,hx; 117 int e0,nx; 118 119 GET_FLOAT_WORD(hx,x); 120 ix = hx&0x7fffffff; 121 if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */ 122 {y[0] = x; y[1] = 0; return 0;} 123 if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */ 124 if(hx>0) { 125 z = x - pio2_1; 126 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ 127 y[0] = z - pio2_1t; 128 y[1] = (z-y[0])-pio2_1t; 129 } else { /* near pi/2, use 24+24+24 bit pi */ 130 z -= pio2_2; 131 y[0] = z - pio2_2t; 132 y[1] = (z-y[0])-pio2_2t; 133 } 134 return 1; 135 } else { /* negative x */ 136 z = x + pio2_1; 137 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ 138 y[0] = z + pio2_1t; 139 y[1] = (z-y[0])+pio2_1t; 140 } else { /* near pi/2, use 24+24+24 bit pi */ 141 z += pio2_2; 142 y[0] = z + pio2_2t; 143 y[1] = (z-y[0])+pio2_2t; 144 } 145 return -1; 146 } 147 } 148 if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */ 149 t = fabsf(x); 150 n = (__int32_t) (t*invpio2+half); 151 fn = (float)n; 152 r = t-fn*pio2_1; 153 w = fn*pio2_1t; /* 1st round good to 40 bit */ 154 if(n<32&&(ix&0xffffff00)!=npio2_hw[n-1]) { 155 y[0] = r-w; /* quick check no cancellation */ 156 } else { 157 __uint32_t high; 158 j = ix>>23; 159 y[0] = r-w; 160 GET_FLOAT_WORD(high,y[0]); 161 i = j-((high>>23)&0xff); 162 if(i>8) { /* 2nd iteration needed, good to 57 */ 163 t = r; 164 w = fn*pio2_2; 165 r = t-w; 166 w = fn*pio2_2t-((t-r)-w); 167 y[0] = r-w; 168 GET_FLOAT_WORD(high,y[0]); 169 i = j-((high>>23)&0xff); 170 if(i>25) { /* 3rd iteration need, 74 bits acc */ 171 t = r; /* will cover all possible cases */ 172 w = fn*pio2_3; 173 r = t-w; 174 w = fn*pio2_3t-((t-r)-w); 175 y[0] = r-w; 176 } 177 } 178 } 179 y[1] = (r-y[0])-w; 180 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 181 else return n; 182 } 183 /* 184 * all other (large) arguments 185 */ 186 if(!FLT_UWORD_IS_FINITE(ix)) { 187 y[0]=y[1]=x-x; return 0; 188 } 189 /* set z = scalbn(|x|,ilogb(x)-7) */ 190 e0 = (int)((ix>>23)-134); /* e0 = ilogb(z)-7; */ 191 SET_FLOAT_WORD(z, ix - ((__int32_t)e0<<23)); 192 for(i=0;i<2;i++) { 193 tx[i] = (float)((__int32_t)(z)); 194 z = (z-tx[i])*two8; 195 } 196 tx[2] = z; 197 nx = 3; 198 while(tx[nx-1]==zero) nx--; /* skip zero term */ 199 n = __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi); 200 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 201 return n; 202 } 203