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