lib/libm/sqrtf.c

   1 /*                                                      sqrtf.c
   2  *
   3  *      Square root
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, sqrtf();
  10  *
  11  * y = sqrtf( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the square root of x.
  18  *
  19  * Range reduction involves isolating the power of two of the
  20  * argument and using a polynomial approximation to obtain
  21  * a rough value for the square root.  Then Heron's iteration
  22  * is used three times to converge to an accurate value.
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *
  29  *                      Relative error:
  30  * arithmetic   domain     # trials      peak         rms
  31  *    IEEE      0,1.e38     100000       8.7e-8     2.9e-8
  32  *
  33  *
  34  * ERROR MESSAGES:
  35  *
  36  *   message         condition      value returned
  37  * sqrtf domain        x < 0            0.0
  38  *
  39  */
  40 \f
  41 /*
  42 Cephes Math Library Release 2.2:  June, 1992
  43 Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
  44 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  45 */
  46
  47 /* Single precision square root
  48  * test interval: [sqrt(2)/2, sqrt(2)]
  49  * trials: 30000
  50  * peak relative error: 8.8e-8
  51  * rms relative error: 3.3e-8
  52  *
  53  * test interval: [0.01, 100.0]
  54  * trials: 50000
  55  * peak relative error: 8.7e-8
  56  * rms relative error: 3.3e-8
  57  *
  58  * Copyright (C) 1989 by Stephen L. Moshier.  All rights reserved.
  59  */
  60 #include "mconf.h"
  61
  62 #ifdef ANSIC
  63 float frexpf( float, int * );
  64 float ldexpf( float, int );
  65
  66 float sqrtf( float xx )
  67 #else
  68 float frexpf(), ldexpf();
  69
  70 float sqrtf(xx)
  71 float xx;
  72 #endif
  73 {
  74 float f, x, y;
  75 int e;
  76
  77 f = xx;
  78 if( f <= 0.0 )
  79         {
  80         if( f < 0.0 )
  81                 mtherr( "sqrtf", DOMAIN );
  82         return( 0.0 );
  83         }
  84
  85 x = frexpf( f, &e );    /* f = x * 2**e,   0.5 <= x < 1.0 */
  86 /* If power of 2 is odd, double x and decrement the power of 2. */
  87 if( e & 1 )
  88         {
  89         x = x + x;
  90         e -= 1;
  91         }
  92
  93 e >>= 1;        /* The power of 2 of the square root. */
  94
  95 if( x > 1.41421356237 )
  96         {
  97 /* x is between sqrt(2) and 2. */
  98         x = x - 2.0;
  99         y =
 100         ((((( -9.8843065718E-4 * x
 101           + 7.9479950957E-4) * x
 102           - 3.5890535377E-3) * x
 103           + 1.1028809744E-2) * x
 104           - 4.4195203560E-2) * x
 105           + 3.5355338194E-1) * x
 106           + 1.41421356237E0;
 107         goto sqdon;
 108         }
 109
 110 if( x > 0.707106781187 )
 111         {
 112 /* x is between sqrt(2)/2 and sqrt(2). */
 113         x = x - 1.0;
 114         y =
 115         ((((( 1.35199291026E-2 * x
 116           - 2.26657767832E-2) * x
 117           + 2.78720776889E-2) * x
 118           - 3.89582788321E-2) * x
 119           + 6.24811144548E-2) * x
 120           - 1.25001503933E-1) * x * x
 121           + 0.5 * x
 122           + 1.0;
 123         goto sqdon;
 124         }
 125
 126 /* x is between 0.5 and sqrt(2)/2. */
 127 x = x - 0.5;
 128 y =
 129 ((((( -3.9495006054E-1 * x
 130   + 5.1743034569E-1) * x
 131   - 4.3214437330E-1) * x
 132   + 3.5310730460E-1) * x
 133   - 3.5354581892E-1) * x
 134   + 7.0710676017E-1) * x
 135   + 7.07106781187E-1;
 136
 137 sqdon:
 138 y = ldexpf( y, e );  /* y = y * 2**e */
 139 return( y);
 140 }