gcc/testsuite/gcc.c-torture/execute/cbrt.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10 */
  11
  12 void abort (void);
  13 void exit (int);
  14
  15 #ifndef __vax__
  16 static const unsigned long
  17         B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
  18         B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
  19
  20 static const double
  21         C =  5.42857142857142815906e-01, /* 19/35     = 0x3FE15F15, 0xF15F15F1 */
  22         D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */
  23         E =  1.41428571428571436819e+00, /* 99/70     = 0x3FF6A0EA, 0x0EA0EA0F */
  24         F =  1.60714285714285720630e+00, /* 45/28     = 0x3FF9B6DB, 0x6DB6DB6E */
  25         G =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
  26
  27 double
  28 cbrtl (double x)
  29 {
  30   long hx;
  31   double r,s,w;
  32   double lt;
  33   unsigned sign;
  34   typedef unsigned unsigned32 __attribute__((mode(SI)));
  35   union {
  36     double t;
  37     unsigned32 pt[2];
  38   } ut, ux;
  39   int n0;
  40
  41   ut.t = 1.0;
  42   n0 = (ut.pt[0] == 0);
  43
  44   ut.t = 0.0;
  45   ux.t = x;
  46
  47   hx = ux.pt[n0];                       /* high word of x */
  48   sign=hx&0x80000000;                   /* sign= sign(x) */
  49   hx  ^=sign;
  50   if(hx>=0x7ff00000) return(x+x);       /* cbrt(NaN,INF) is itself */
  51   if((hx| ux.pt[1-n0])==0)
  52     return(ux.t);                       /* cbrt(0) is itself */
  53
  54   ux.pt[n0] = hx;
  55   /* rough cbrt to 5 bits */
  56   if(hx<0x00100000)                     /* subnormal number */
  57     {ut.pt[n0]=0x43500000;              /* set t= 2**54 */
  58      ut.t*=x; ut.pt[n0]=ut.pt[n0]/3+B2;
  59    }
  60   else
  61     ut.pt[n0]=hx/3+B1;
  62
  63   /* new cbrt to 23 bits, may be implemented in single precision */
  64   r=ut.t*ut.t/ux.t;
  65   s=C+r*ut.t;
  66   ut.t*=G+F/(s+E+D/s);
  67
  68   /* chopped to 20 bits and make it larger than cbrt(x) */
  69   ut.pt[1-n0]=0; ut.pt[n0]+=0x00000001;
  70
  71   /* one step newton iteration to 53 bits with error less than 0.667 ulps */
  72   s=ut.t*ut.t;                          /* t*t is exact */
  73   r=ux.t/s;
  74   w=ut.t+ut.t;
  75   r=(r-ut.t)/(w+r);                     /* r-s is exact */
  76   ut.t=ut.t+ut.t*r;
  77
  78   /* restore the sign bit */
  79   ut.pt[n0] |= sign;
  80
  81   lt = ut.t;
  82   lt -= (lt - (x/(lt*lt))) * 0.333333333333333333333;
  83   return lt;
  84 }
  85
  86 int
  87 main (void)
  88 {
  89   if ((int) (cbrtl (27.0) + 0.5) != 3)
  90     abort ();
  91
  92   exit (0);
  93 }
  94 #else
  95 int main (void) { exit (0); }
  96 #endif