11.1Sjtc/* @(#)s_cbrt.c 5.1 93/09/24 */
21.1Sjtc/*
31.1Sjtc * ====================================================
41.1Sjtc * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
51.1Sjtc *
61.1Sjtc * Developed at SunPro, a Sun Microsystems, Inc. business.
71.1Sjtc * Permission to use, copy, modify, and distribute this
81.10Ssimonb * software is freely granted, provided that this notice
91.1Sjtc * is preserved.
101.1Sjtc * ====================================================
111.1Sjtc */
121.3Sjtc
131.9Slukem#include <sys/cdefs.h>
141.7Sjtc#if defined(LIBM_SCCS) && !defined(lint)
151.13Sjoerg__RCSID("$NetBSD: s_cbrt.c,v 1.13 2013/11/20 12:29:13 joerg Exp $");
161.3Sjtc#endif
171.1Sjtc
181.12Sjoerg#include "namespace.h"
191.5Sjtc#include "math.h"
201.5Sjtc#include "math_private.h"
211.1Sjtc
221.12Sjoerg#ifndef __HAVE_LONG_DOUBLE
231.13Sjoerg__strong_alias(_cbrtl, cbrt)
241.12Sjoerg__weak_alias(cbrtl, _cbrtl)
251.12Sjoerg#endif
261.12Sjoerg
271.1Sjtc/* cbrt(x)
281.1Sjtc * Return cube root of x
291.1Sjtc */
301.6Sjtcstatic const u_int32_t
311.1Sjtc	B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
321.1Sjtc	B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
331.1Sjtc
341.1Sjtcstatic const double
351.1SjtcC =  5.42857142857142815906e-01, /* 19/35     = 0x3FE15F15, 0xF15F15F1 */
361.1SjtcD = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */
371.1SjtcE =  1.41428571428571436819e+00, /* 99/70     = 0x3FF6A0EA, 0x0EA0EA0F */
381.1SjtcF =  1.60714285714285720630e+00, /* 45/28     = 0x3FF9B6DB, 0x6DB6DB6E */
391.1SjtcG =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
401.1Sjtc
411.11Swizdouble
421.11Swizcbrt(double x)
431.1Sjtc{
441.6Sjtc	int32_t	hx;
451.1Sjtc	double r,s,t=0.0,w;
461.6Sjtc	u_int32_t sign;
471.6Sjtc	u_int32_t high,low;
481.1Sjtc
491.5Sjtc	GET_HIGH_WORD(hx,x);
501.1Sjtc	sign=hx&0x80000000; 		/* sign= sign(x) */
511.1Sjtc	hx  ^=sign;
521.1Sjtc	if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
531.5Sjtc	GET_LOW_WORD(low,x);
541.10Ssimonb	if((hx|low)==0)
551.1Sjtc	    return(x);		/* cbrt(0) is itself */
561.1Sjtc
571.5Sjtc	SET_HIGH_WORD(x,hx);	/* x <- |x| */
581.1Sjtc    /* rough cbrt to 5 bits */
591.1Sjtc	if(hx<0x00100000) 		/* subnormal number */
601.5Sjtc	  {SET_HIGH_WORD(t,0x43500000);	/* set t= 2**54 */
611.5Sjtc	   t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2);
621.1Sjtc	  }
631.1Sjtc	else
641.5Sjtc	  SET_HIGH_WORD(t,hx/3+B1);
651.1Sjtc
661.1Sjtc
671.1Sjtc    /* new cbrt to 23 bits, may be implemented in single precision */
681.1Sjtc	r=t*t/x;
691.1Sjtc	s=C+r*t;
701.10Ssimonb	t*=G+F/(s+E+D/s);
711.1Sjtc
721.10Ssimonb    /* chopped to 20 bits and make it larger than cbrt(x) */
731.5Sjtc	GET_HIGH_WORD(high,t);
741.5Sjtc	INSERT_WORDS(t,high+0x00000001,0);
751.1Sjtc
761.1Sjtc
771.1Sjtc    /* one step newton iteration to 53 bits with error less than 0.667 ulps */
781.1Sjtc	s=t*t;		/* t*t is exact */
791.1Sjtc	r=x/s;
801.1Sjtc	w=t+t;
811.1Sjtc	r=(r-t)/(w+r);	/* r-s is exact */
821.1Sjtc	t=t+t*r;
831.1Sjtc
841.1Sjtc    /* retore the sign bit */
851.5Sjtc	GET_HIGH_WORD(high,t);
861.5Sjtc	SET_HIGH_WORD(t,high|sign);
871.1Sjtc	return(t);
881.1Sjtc}
89