Open64: osprey/libm/mips/acos.c Source File

00001 /*
00002  * Copyright 2003, 2004 PathScale, Inc.  All Rights Reserved.
00003  */
00004 
00005 
00006 /*
00007 
00008   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00009 
00010   This program is free software; you can redistribute it and/or modify it
00011   under the terms of version 2.1 of the GNU Lesser General Public License 
00012   as published by the Free Software Foundation.
00013 
00014   This program is distributed in the hope that it would be useful, but
00015   WITHOUT ANY WARRANTY; without even the implied warranty of
00016   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00017 
00018   Further, this software is distributed without any warranty that it is
00019   free of the rightful claim of any third person regarding infringement 
00020   or the like.  Any license provided herein, whether implied or 
00021   otherwise, applies only to this software file.  Patent licenses, if
00022   any, provided herein do not apply to combinations of this program with 
00023   other software, or any other product whatsoever.  
00024 
00025   You should have received a copy of the GNU Lesser General Public 
00026   License along with this program; if not, write the Free Software 
00027   Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307, 
00028   USA.
00029 
00030   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00031   Mountain View, CA 94043, or:
00032 
00033   http://www.sgi.com
00034 
00035   For further information regarding this notice, see:
00036 
00037   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00038 
00039 */
00040 
00041 /* ====================================================================
00042  * ====================================================================
00043  *
00044  * Module: acos.c
00045  * $Revision: 1.1.1.1 $
00046  * $Date: 2005/10/21 19:00:00 $
00047  * $Author: marcel $
00048  * $Source: /proj/osprey/CVS/open64/osprey1.0/libm/mips/acos.c,v $
00049  *
00050  * Revision history:
00051  *  09-Jun-93 - Original Version
00052  *
00053  * Description: source code for acos function
00054  *
00055  * ====================================================================
00056  * ====================================================================
00057  */
00058 
00059 static char *rcs_id = "$Source: /proj/osprey/CVS/open64/osprey1.0/libm/mips/acos.c,v $ $Revision: 1.1.1.1 $";
00060 
00061 #ifdef _CALL_MATHERR
00062 #include <stdio.h>
00063 #include <math.h>
00064 #include <errno.h>
00065 #endif
00066 
00067 #include "libm.h"
00068 
00069 #if defined(mips) && !defined(__GNUC__)
00070 extern  double  acos(double);
00071 
00072 #pragma weak acos = __acos
00073 #endif
00074 
00075 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00076 extern double __acos(double);
00077 #pragma weak acos
00078 double acos(double x) {
00079   return __acos(x);
00080 }
00081 #elif defined(__GNUC__)
00082 extern  double  __acos(double);
00083 
00084 double    acos() __attribute__ ((weak, alias ("__acos")));
00085 
00086 #endif
00087 
00088 /* coefficients for rational approximation of asin on +/- 0.5    */
00089 
00090 static const du P[] =
00091 {
00092 {D(0x00000000, 0x00000000)},
00093 {D(0x400cf4d7, 0x1166375d)},
00094 {D(0xc019dde2, 0xfd680d32)},
00095 {D(0x400c15ee, 0x4cc68a6a)},
00096 {D(0xbfe20f7e, 0xdc1c40fe)},
00097 {D(0x3f780cd5, 0x52bc78fd)},
00098 };
00099 
00100 static const du Q[] =
00101 {
00102 {D(0x4035b7a1, 0x4d0ca925)},
00103 {D(0xc0484954, 0xef63fb64)},
00104 {D(0x40428d6e, 0x183c02d2)},
00105 {D(0xc026104e, 0x748fbc12)},
00106 {D(0x3ff00000, 0x00000000)},
00107 };
00108 
00109 /* coefficients for rational approximation of asin on +/- sqrt(2 - sqrt(3))/2 */
00110 
00111 static const du P2[] =
00112 {
00113 {D(0x00000000, 0x00000000)},
00114 {D(0xc0097a02, 0x2e7a0e13)},
00115 {D(0x4009aeb0, 0x0736a7c7)},
00116 {D(0xbfe37279, 0x9e195a2e)},
00117 };
00118 
00119 static const du Q2[] =
00120 {
00121 {D(0xc0331b81, 0xa2db8f8a)},
00122 {D(0x403bdc31, 0x8eb7aeb5)},
00123 {D(0xc0262173, 0xdc9ece8a)},
00124 {D(0x3ff00000, 0x00000000)},
00125 };
00126 
00127 static const  du  Qnan =
00128 {D(QNANHI, QNANLO)};
00129 
00130 static const du one =
00131 {D(0x3ff00000, 0x00000000)};
00132 
00133 static const du m_one =
00134 {D(0xbff00000, 0x00000000)};
00135 
00136 static const du piby4 =
00137 {D(0x3fe921fb, 0x54442d18)};
00138 
00139 static const du piby2 =
00140 {D(0x3ff921fb, 0x54442d18)};
00141 
00142 static const du pi =
00143 {D(0x400921fb, 0x54442d18)};
00144 
00145 #ifdef _32BIT_MACHINE
00146 
00147 static const int  root3by2 = 0x3febb67a;
00148 
00149 #else
00150 
00151 static const long long  root3by2 = 0x3febb67ae8584caall;
00152 
00153 #endif
00154 
00155 
00156 /* ====================================================================
00157  *
00158  * FunctionName   acos  
00159  *
00160  * Description    computes arccosine of arg
00161  *
00162  * ====================================================================
00163  */
00164 
00165 double
00166 __acos( double x )
00167 {
00168 #ifdef _32BIT_MACHINE
00169 
00170 int ix, xpt;
00171 unsigned int  iabsx;
00172 
00173 #else
00174 
00175 long long ix, xpt;
00176 unsigned long long iabsx;
00177 
00178 #endif
00179 
00180 double  xsq, num, denom, result;
00181 double  y, ysq;
00182 double  absx;
00183 double  q;
00184 #ifdef _CALL_MATHERR
00185 struct exception  exstruct;
00186 #endif
00187 
00188   /* Note that asin(x) + acos(x) = pi/2, so we can use
00189    * approximations for asin(x) to compute acos(x).
00190    */
00191 
00192   /* extract exponent of x for some quick screening */
00193 
00194 #ifdef _32BIT_MACHINE
00195 
00196   DBLHI2INT(x, ix);
00197 #else
00198   DBL2LL(x, ix);
00199 #endif
00200   iabsx = (ix << 1);
00201   iabsx >>= 1;
00202   xpt = (iabsx >> DMANTWIDTH);
00203 
00204   if ( xpt < 0x3fe )
00205   {
00206     /*   |x| < 0.5    */
00207 
00208     if ( xpt >= 0x3c8 )
00209     {
00210       /* |x| >= 2^(-55)   */
00211 
00212       xsq = x*x;
00213       num = ((((P[5].d*xsq + P[4].d)*xsq + P[3].d)*xsq + P[2].d)*xsq +
00214         P[1].d);
00215 
00216       denom = ((((xsq + Q[3].d)*xsq + Q[2].d)*xsq + Q[1].d)*xsq +
00217         Q[0].d);
00218 
00219       result = x + x*(xsq*num)/denom;
00220 
00221       return ( piby4.d - result + piby4.d );
00222     }
00223 
00224     return ( piby2.d );
00225   }
00226 
00227   if ( xpt < 0x3ff )
00228   {
00229     /*  |x| < 1.0   */
00230 
00231     if ( iabsx < root3by2 )
00232     {
00233       /* uses the identity asin(x) = pi/4 + 0.5*asin(2*x**2 - 1) */
00234 
00235       xsq = x*x;
00236       y = xsq + xsq - one.d;
00237       xsq = y*y;
00238 
00239       num = ((((P[5].d*xsq + P[4].d)*xsq + P[3].d)*xsq + P[2].d)*xsq +
00240         P[1].d);
00241 
00242       denom = ((((xsq + Q[3].d)*xsq + Q[2].d)*xsq + Q[1].d)*xsq +
00243         Q[0].d);
00244 
00245       result = y + y*(xsq*num)/denom;
00246 
00247       if ( x > 0.0 )
00248         result = piby4.d - 0.5*result;
00249       else
00250         result = piby2.d + (piby4.d + 0.5*result);
00251 
00252       return ( result );
00253     }
00254 
00255     /* uses the identity asin(x) = pi/2 - 2*asin(sqrt(.5*(1 - x)) */
00256 
00257     absx = fabs(x);
00258 
00259     ysq = 0.5*(one.d - absx);
00260     y = sqrt(ysq);
00261 
00262     num = ((P2[3].d*ysq + P2[2].d)*ysq + P2[1].d);
00263     denom = (((ysq + Q2[2].d)*ysq + Q2[1].d)*ysq + Q2[0].d);
00264 
00265     q = y + y*(ysq*num)/denom;
00266 
00267     if ( x > 0.0 )
00268       result = q + q;
00269     else
00270       result = pi.d - (q + q);
00271 
00272     return ( result );
00273   }
00274 
00275   if ( x != x )
00276   {
00277     /* x is a NaN; return a quiet NaN */
00278 
00279 #ifdef _CALL_MATHERR
00280 
00281     exstruct.type = DOMAIN;
00282     exstruct.name = "acos";
00283     exstruct.arg1 = x;
00284     exstruct.retval = Qnan.d;
00285 
00286     if ( matherr( &exstruct ) == 0 )
00287     {
00288       fprintf(stderr, "domain error in acos\n");
00289       SETERRNO(EDOM);
00290     }
00291 
00292     return ( exstruct.retval );
00293 #else
00294     NAN_SETERRNO(EDOM);
00295 
00296     return ( Qnan.d );
00297 #endif
00298   }
00299 
00300   if ( x == one.d )
00301     return ( 0.0 );
00302 
00303   if ( x == m_one.d )
00304     return ( pi.d );
00305 
00306   /*  |x| > 1.0; return a NaN */
00307 
00308 #ifdef _CALL_MATHERR
00309 
00310   exstruct.type = DOMAIN;
00311   exstruct.name = "acos";
00312   exstruct.arg1 = x;
00313   exstruct.retval = Qnan.d;
00314 
00315   if ( matherr( &exstruct ) == 0 )
00316   {
00317     fprintf(stderr, "domain error in acos\n");
00318     SETERRNO(EDOM);
00319   }
00320 
00321   return ( exstruct.retval );
00322 #else
00323   SETERRNO(EDOM);
00324 
00325   return ( Qnan.d );
00326 #endif
00327 }
00328 
00329 #ifdef NO_LONG_DOUBLE
00330 
00331 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00332 #pragma weak acosl
00333 long double acosl( long double x ) {  
00334   return ( (long double)__acos((double)x) );
00335 }
00336 #elif defined(__GNUC__)
00337 extern  long double  __acosl(long double);
00338 
00339 long double    acosl() __attribute__ ((weak, alias ("__acosl")));
00340 
00341 #endif
00342 
00343 long double
00344 __acosl( long double x )
00345 { 
00346   return ( (long double)__acos((double)x) );
00347 }
00348 
00349 #endif
00350