1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 24 * Use is subject to license terms. 25 */ 26 27 /* 28 * Copyright 2011, Richard Lowe 29 */ 30 31 #ifndef _LIBM_INLINES_H 32 #define _LIBM_INLINES_H 33 34 #ifdef __GNUC__ 35 36 #ifdef __cplusplus 37 extern "C" { 38 #endif 39 40 #include <sys/types.h> 41 #include <sys/ieeefp.h> 42 43 #define _LO_WORD(x) ((uint32_t *)&x)[0] 44 #define _HI_WORD(x) ((uint32_t *)&x)[1] 45 #define _HIER_WORD(x) ((uint32_t *)&x)[2] 46 47 extern __GNU_INLINE double 48 __inline_sqrt(double a) 49 { 50 double ret; 51 52 __asm__ __volatile__("fsqrt\n\t" : "=t" (ret) : "0" (a) : "cc"); 53 return (ret); 54 } 55 56 extern __GNU_INLINE double 57 __ieee754_sqrt(double a) 58 { 59 return (__inline_sqrt(a)); 60 } 61 62 extern __GNU_INLINE float 63 __inline_sqrtf(float a) 64 { 65 float ret; 66 67 __asm__ __volatile__("fsqrt\n\t" : "=t" (ret) : "0" (a) : "cc"); 68 return (ret); 69 } 70 71 extern __GNU_INLINE double 72 __inline_rint(double a) 73 { 74 __asm__ __volatile__( 75 "andl $0x7fffffff,%1\n\t" 76 "cmpl $0x43300000,%1\n\t" 77 "jae 1f\n\t" 78 "frndint\n\t" 79 "1: fwait\n\t" 80 : "+t" (a), "+&r" (_HI_WORD(a)) 81 : 82 : "cc"); 83 84 return (a); 85 } 86 87 /* 88 * 00 - 24 bits 89 * 01 - reserved 90 * 10 - 53 bits 91 * 11 - 64 bits 92 */ 93 extern __GNU_INLINE int 94 __swapRP(int i) 95 { 96 int ret; 97 uint16_t cw; 98 99 __asm__ __volatile__("fstcw %0\n\t" : "=m" (cw)); 100 101 ret = (cw >> 8) & 0x3; 102 cw = (cw & 0xfcff) | ((i & 0x3) << 8); 103 104 __asm__ __volatile__("fldcw %0\n\t" : : "m" (cw)); 105 106 return (ret); 107 } 108 109 /* 110 * 00 - Round to nearest, with even preferred 111 * 01 - Round down 112 * 10 - Round up 113 * 11 - Chop 114 */ 115 extern __GNU_INLINE enum fp_direction_type 116 __swap87RD(enum fp_direction_type i) 117 { 118 int ret; 119 uint16_t cw; 120 121 __asm__ __volatile__("fstcw %0\n\t" : "=m" (cw)); 122 123 ret = (cw >> 10) & 0x3; 124 cw = (cw & 0xf3ff) | ((i & 0x3) << 10); 125 126 __asm__ __volatile__("fldcw %0\n\t" : : "m" (cw)); 127 128 return (ret); 129 } 130 131 extern __GNU_INLINE double 132 ceil(double d) 133 { 134 /* 135 * Let's set a Rounding Control (RC) bits from x87 FPU Control Word 136 * to fp_positive and save old bits in rd. 137 */ 138 short rd = __swap87RD(fp_positive); 139 140 /* 141 * The FRNDINT instruction returns a floating-point value that is the 142 * integral value closest to the source value in the direction of the 143 * rounding mode specified in the RC field of the x87 FPU control word. 144 * 145 * Rounds the source value in the ST(0) register to the nearest 146 * integral value, depending on the current rounding mode 147 * (setting of the RC field of the FPU control word), 148 * and stores the result in ST(0). 149 */ 150 __asm__ __volatile__("frndint" : "+t" (d) : : "cc"); 151 152 /* restore old RC bits */ 153 __swap87RD(rd); 154 155 return (d); 156 } 157 158 extern __GNU_INLINE double 159 copysign(double d1, double d2) 160 { 161 __asm__ __volatile__( 162 "andl $0x7fffffff,%0\n\t" /* %0 <-- hi_32(abs(d)) */ 163 "andl $0x80000000,%1\n\t" /* %1[31] <-- sign_bit(d2) */ 164 "orl %1,%0\n\t" /* %0 <-- hi_32(copysign(x,y)) */ 165 : "+&r" (_HI_WORD(d1)), "+r" (_HI_WORD(d2)) 166 : 167 : "cc"); 168 169 return (d1); 170 } 171 172 extern __GNU_INLINE double 173 fabs(double d) 174 { 175 __asm__ __volatile__("fabs\n\t" : "+t" (d) : : "cc"); 176 return (d); 177 } 178 179 extern __GNU_INLINE float 180 fabsf(float d) 181 { 182 __asm__ __volatile__("fabs\n\t" : "+t" (d) : : "cc"); 183 return (d); 184 } 185 186 extern __GNU_INLINE long double 187 fabsl(long double d) 188 { 189 __asm__ __volatile__("fabs\n\t" : "+t" (d) : : "cc"); 190 return (d); 191 } 192 193 extern __GNU_INLINE int 194 finite(double d) 195 { 196 int ret = _HI_WORD(d); 197 198 __asm__ __volatile__( 199 "notl %0\n\t" 200 "andl $0x7ff00000,%0\n\t" 201 "negl %0\n\t" 202 "shrl $31,%0\n\t" 203 : "+r" (ret) 204 : 205 : "cc"); 206 return (ret); 207 } 208 209 extern __GNU_INLINE double 210 floor(double d) 211 { 212 short rd = __swap87RD(fp_negative); 213 214 __asm__ __volatile__("frndint" : "+t" (d), "+r" (rd) : : "cc"); 215 __swap87RD(rd); 216 217 return (d); 218 } 219 220 /* 221 * branchless __isnan 222 * ((0x7ff00000-[((lx|-lx)>>31)&1]|ahx)>>31)&1 = 1 iff x is NaN 223 */ 224 extern __GNU_INLINE int 225 isnan(double d) 226 { 227 int ret; 228 229 __asm__ __volatile__( 230 "movl %1,%%ecx\n\t" 231 "negl %%ecx\n\t" /* ecx <-- -lo_32(x) */ 232 "orl %%ecx,%1\n\t" 233 "shrl $31,%1\n\t" /* 1 iff lx != 0 */ 234 "andl $0x7fffffff,%2\n\t" /* ecx <-- hi_32(abs(x)) */ 235 "orl %2,%1\n\t" 236 "subl $0x7ff00000,%1\n\t" 237 "negl %1\n\t" 238 "shrl $31,%1\n\t" 239 : "=r" (ret) 240 : "0" (_HI_WORD(d)), "r" (_LO_WORD(d)) 241 : "ecx"); 242 243 return (ret); 244 } 245 246 extern __GNU_INLINE int 247 isnanf(float f) 248 { 249 __asm__ __volatile__( 250 "andl $0x7fffffff,%0\n\t" 251 "negl %0\n\t" 252 "addl $0x7f800000,%0\n\t" 253 "shrl $31,%0\n\t" 254 : "+r" (f) 255 : 256 : "cc"); 257 258 return (f); 259 } 260 261 extern __GNU_INLINE double 262 rint(double a) 263 { 264 return (__inline_rint(a)); 265 } 266 267 extern __GNU_INLINE double 268 scalbn(double d, int n) 269 { 270 double dummy; 271 272 __asm__ __volatile__( 273 "fildl %2\n\t" /* Convert N to extended */ 274 "fxch\n\t" 275 "fscale\n\t" 276 : "+t" (d), "=u" (dummy) 277 : "m" (n) 278 : "cc"); 279 280 return (d); 281 } 282 283 extern __GNU_INLINE int 284 signbit(double d) 285 { 286 return (_HI_WORD(d) >> 31); 287 } 288 289 extern __GNU_INLINE int 290 signbitf(float f) 291 { 292 return ((*(uint32_t *)&f) >> 31); 293 } 294 295 extern __GNU_INLINE double 296 sqrt(double d) 297 { 298 return (__inline_sqrt(d)); 299 } 300 301 extern __GNU_INLINE float 302 sqrtf(float f) 303 { 304 return (__inline_sqrtf(f)); 305 } 306 307 extern __GNU_INLINE long double 308 sqrtl(long double ld) 309 { 310 __asm__ __volatile__("fsqrt" : "+t" (ld) : : "cc"); 311 return (ld); 312 } 313 314 extern __GNU_INLINE int 315 isnanl(long double ld) 316 { 317 int ret = _HIER_WORD(ld); 318 319 __asm__ __volatile__( 320 "andl $0x00007fff,%0\n\t" 321 "jz 1f\n\t" /* jump if exp is all 0 */ 322 "xorl $0x00007fff,%0\n\t" 323 "jz 2f\n\t" /* jump if exp is all 1 */ 324 "testl $0x80000000,%1\n\t" 325 "jz 3f\n\t" /* jump if leading bit is 0 */ 326 "xorl %0,%0\n\t" 327 "jmp 1f\n\t" 328 "2:\n\t" /* note that %0 = 0 from before */ 329 "cmpl $0x80000000,%1\n\t" /* what is first half of significand? */ 330 "jnz 3f\n\t" /* jump if not equal to 0x80000000 */ 331 "testl $0xffffffff,%2\n\t" /* is second half of significand 0? */ 332 "jnz 3f\n\t" /* jump if not equal to 0 */ 333 "jmp 1f\n\t" 334 "3:\n\t" 335 "movl $1,%0\n\t" 336 "1:\n\t" 337 : "+&r" (ret) 338 : "r" (_HI_WORD(ld)), "r" (_LO_WORD(ld)) 339 : "cc"); 340 341 return (ret); 342 } 343 344 #ifdef __cplusplus 345 } 346 #endif 347 348 #endif /* __GNUC__ */ 349 350 #endif /* _LIBM_INLINES_H */ 351