/* cbrt.c - software cube-root function *//* Copyright 1992 Wind River Systems, Inc. *//*modification history--------------------01a,30jul92,smb  documentation.*//*DESCRIPTION* Copyright (c) 1985 Regents of the University of California.* All rights reserved.** Redistribution and use in source and binary forms are permitted* provided that the above copyright notice and this paragraph are* duplicated in all such forms and that any documentation,* advertising materials, and other materials related to such* distribution and use acknowledge that the software was developed* by the University of California, Berkeley.  The name of the* University may not be used to endorse or promote products derived* from this software without specific prior written permission.* THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.** All recipients should regard themselves as participants in an ongoing* research project and hence should feel obligated to report their* experiences (good or bad) with these elementary function codes, using* the sendbug(8) program, to the authors.*#ifndef lintSEE ALSO: American National Standard X3.159-1989*/#include "vxWorks.h"#include "math.h"#if !defined(vax)&&!defined(tahoe)static unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */	             B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */static double	    C= 19./35.,	    D= -864./1225.,	    E= 99./70.,	    F= 45./28.,	    G= 5./14.;/****************************************************************************** cbrt	-** kahan's cube root (53 bits IEEE double precision)* for IEEE machines only* coded in C by K.C. Ng, 4/30/85** Accuracy:*	better than 0.667 ulps according to an error analysis. Maximum* error observed was 0.666 ulps in an 1,000,000 random arguments test.** Warning: this code is semi machine dependent; the ordering of words in* a floating point number must be known in advance. I assume that the* long interger at the address of a floating point number will be the* leading 32 bits of that floating point number (i.e., sign, exponent,* and the 20 most significant bits).* On a National machine, it has different ordering; therefore, this code* must be compiled with flag -DNATIONAL.*/double cbrt(x)double x;{	double r,s,t=0.0,w;	unsigned long *px = (unsigned long *) &x,	              *pt = (unsigned long *) &t,		      mexp,sign;#ifdef national /* ordering of words in a floating points number */	int n0=1,n1=0;#else	/* national */	int n0=0,n1=1;#endif	/* national */	mexp=px[n0]&0x7ff00000;	if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */	if(x==0.0) return(x);		/* cbrt(0) is itself */	sign=px[n0]&0x80000000; /* sign= sign(x) */	px[n0] ^= sign;		/* x=|x| */    /* rough cbrt to 5 bits */	if(mexp==0) 		/* subnormal number */	  {pt[n0]=0x43500000; 	/* set t= 2**54 */	   t*=x; pt[n0]=pt[n0]/3+B2;	  }	else	  pt[n0]=px[n0]/3+B1;    /* new cbrt to 23 bits, may be implemented in single precision */	r=t*t/x;	s=C+r*t;	t*=G+F/(s+E+D/s);    /* chopped to 20 bits and make it larger than cbrt(x) */	pt[n1]=0; pt[n0]+=0x00000001;    /* one step newton iteration to 53 bits with error less than 0.667 ulps */	s=t*t;		/* t*t is exact */	r=x/s;	w=t+t;	r=(r-t)/(w+r);	/* r-t is exact */	t=t+t*r;    /* retore the sign bit */	pt[n0] |= sign;	return(t);}#endif