s_cbrtf.c revision 8870
12116Sjkh/* s_cbrtf.c -- float version of s_cbrt.c. 22116Sjkh * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 32116Sjkh */ 42116Sjkh 52116Sjkh/* 62116Sjkh * ==================================================== 72116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 82116Sjkh * 92116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 102116Sjkh * Permission to use, copy, modify, and distribute this 118870Srgrimes * software is freely granted, provided that this notice 122116Sjkh * is preserved. 132116Sjkh * ==================================================== 142116Sjkh */ 152116Sjkh 162116Sjkh#ifndef lint 178870Srgrimesstatic char rcsid[] = "$Id: s_cbrtf.c,v 1.1.1.1 1994/08/19 09:39:57 jkh Exp $"; 182116Sjkh#endif 192116Sjkh 202116Sjkh#include "math.h" 212116Sjkh#include "math_private.h" 222116Sjkh 232116Sjkh/* cbrtf(x) 242116Sjkh * Return cube root of x 252116Sjkh */ 262116Sjkh#ifdef __STDC__ 278870Srgrimesstatic const unsigned 282116Sjkh#else 298870Srgrimesstatic unsigned 302116Sjkh#endif 312116Sjkh B1 = 709958130, /* B1 = (84+2/3-0.03306235651)*2**23 */ 322116Sjkh B2 = 642849266; /* B2 = (76+2/3-0.03306235651)*2**23 */ 332116Sjkh 342116Sjkh#ifdef __STDC__ 352116Sjkhstatic const float 362116Sjkh#else 372116Sjkhstatic float 382116Sjkh#endif 392116SjkhC = 5.4285717010e-01, /* 19/35 = 0x3f0af8b0 */ 402116SjkhD = -7.0530611277e-01, /* -864/1225 = 0xbf348ef1 */ 412116SjkhE = 1.4142856598e+00, /* 99/70 = 0x3fb50750 */ 422116SjkhF = 1.6071428061e+00, /* 45/28 = 0x3fcdb6db */ 432116SjkhG = 3.5714286566e-01; /* 5/14 = 0x3eb6db6e */ 442116Sjkh 452116Sjkh#ifdef __STDC__ 468870Srgrimes float cbrtf(float x) 472116Sjkh#else 488870Srgrimes float cbrtf(x) 492116Sjkh float x; 502116Sjkh#endif 512116Sjkh{ 522116Sjkh float r,s,t; 532116Sjkh int32_t hx; 542116Sjkh u_int32_t sign; 552116Sjkh u_int32_t high; 562116Sjkh 572116Sjkh GET_FLOAT_WORD(hx,x); 582116Sjkh sign=hx&0x80000000; /* sign= sign(x) */ 592116Sjkh hx ^=sign; 602116Sjkh if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */ 618870Srgrimes if(hx==0) 622116Sjkh return(x); /* cbrt(0) is itself */ 632116Sjkh 642116Sjkh SET_FLOAT_WORD(x,hx); /* x <- |x| */ 652116Sjkh /* rough cbrt to 5 bits */ 662116Sjkh if(hx<0x00800000) /* subnormal number */ 672116Sjkh {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */ 682116Sjkh t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2); 692116Sjkh } 702116Sjkh else 712116Sjkh SET_FLOAT_WORD(t,hx/3+B1); 722116Sjkh 732116Sjkh 742116Sjkh /* new cbrt to 23 bits */ 752116Sjkh r=t*t/x; 762116Sjkh s=C+r*t; 778870Srgrimes t*=G+F/(s+E+D/s); 782116Sjkh 792116Sjkh /* retore the sign bit */ 802116Sjkh GET_FLOAT_WORD(high,t); 812116Sjkh SET_FLOAT_WORD(t,high|sign); 822116Sjkh return(t); 832116Sjkh} 84