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

00001 /*
00002  * Copyright 2003, 2004, 2005, 2006 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: sinhf.c
00045  * $Revision: 1.5 $
00046  * $Date: 04/12/21 14:58:23-08:00 $
00047  * $Author: bos@eng-25.internal.keyresearch.com $
00048  * $Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.sinhf.c $
00049  *
00050  * Revision history:
00051  *  09-Jun-93 - Original Version
00052  *
00053  * Description: source code for sinhf function
00054  *
00055  * ====================================================================
00056  * ====================================================================
00057  */
00058 
00059 static char *rcs_id = "$Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.sinhf.c $ $Revision: 1.5 $";
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  float fsinh(float);
00071 extern  float sinhf(float);
00072 extern  float fcosh(float);
00073 extern  float coshf(float);
00074 
00075 #pragma weak fsinh = __sinhf
00076 #pragma weak sinhf = __sinhf
00077 
00078 #pragma weak fcosh = __coshf
00079 #pragma weak coshf = __coshf
00080 #endif
00081 
00082 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00083 extern  float  __sinhf(float);
00084 extern  float  __coshf(float);
00085 #pragma weak sinhf
00086 #pragma weak coshf
00087 float sinhf( float x ) {
00088   return __sinhf( x );
00089 }
00090 float coshf( float x ) {
00091   return __coshf( x );
00092 }
00093 #elif defined(__GNUC__)
00094 extern  float  __sinhf(float);
00095 float    sinhf(float) __attribute__ ((weak, alias ("__sinhf")));
00096 
00097 extern  float  __coshf(float);
00098 float    coshf(float) __attribute__ ((weak, alias ("__coshf")));
00099 #endif
00100 
00101 static const du rln2 =
00102 {D(0x3ff71547, 0x652b82fe)};
00103 
00104 static const du ln2 =
00105 {D(0x3fe62e42, 0xfefa39ef)};
00106 
00107 
00108 /* coefficients for polynomial approximation of cosh on +/- log(2)/2     */
00109 
00110 static const du C[] =
00111 {
00112 {D(0x3ff00000, 0x00000000)},
00113 {D(0x3fe00000, 0x042dcab5)},
00114 {D(0x3fa55548, 0x0b2e33ce)},
00115 {D(0x3f56d962, 0x2f30cebb)},
00116 };
00117 
00118 /* coefficients for polynomial approximation of sinh on +/- log(2)/2     */
00119 
00120 static const du S[] =
00121 {
00122 {D(0x3ff00000, 0x00000000)},
00123 {D(0x3fc55555, 0x5735979e)},
00124 {D(0x3f81110b, 0x21321a53)},
00125 {D(0x3f2a16f7, 0x7ea170a1)},
00126 };
00127 
00128 static const fu f_one = {0x3f800000};
00129 
00130 static const fu Llimit = {0xc2b2d4fc};
00131 
00132 static const fu Ulimit = {0x42b2d4fc};
00133 
00134 static const  du  sinhtab[] =
00135 {
00136 D(0x00000000, 0x00000000),
00137 D(0x3fe80000, 0x00000000),
00138 D(0x3ffe0000, 0x00000000),
00139 D(0x400f8000, 0x00000000),
00140 D(0x401fe000, 0x00000000),
00141 D(0x402ff800, 0x00000000),
00142 D(0x403ffe00, 0x00000000),
00143 D(0x404fff80, 0x00000000),
00144 D(0x405fffe0, 0x00000000),
00145 D(0x406ffff8, 0x00000000),
00146 D(0x407ffffe, 0x00000000),
00147 D(0x408fffff, 0x80000000),
00148 D(0x409fffff, 0xe0000000),
00149 D(0x40afffff, 0xf8000000),
00150 };
00151 
00152 static const  du  coshtab[] =
00153 {
00154 D(0x3ff00000, 0x00000000),
00155 D(0x3ff40000, 0x00000000),
00156 D(0x40010000, 0x00000000),
00157 D(0x40104000, 0x00000000),
00158 D(0x40201000, 0x00000000),
00159 D(0x40300400, 0x00000000),
00160 D(0x40400100, 0x00000000),
00161 D(0x40500040, 0x00000000),
00162 D(0x40600010, 0x00000000),
00163 D(0x40700004, 0x00000000),
00164 D(0x40800001, 0x00000000),
00165 D(0x40900000, 0x40000000),
00166 D(0x40a00000, 0x10000000),
00167 D(0x40b00000, 0x04000000),
00168 };
00169 
00170 static const  fu  Qnan = {QNANF};
00171 
00172 static const fu Neginf = {0x7f800000};
00173 
00174 static const fu Inf = {0x7f800000};
00175 
00176 
00177 /* ====================================================================
00178  *
00179  * FunctionName   sinhf
00180  *
00181  * Description    computes hyperbolic sine of arg
00182  *
00183  * ====================================================================
00184  */
00185 
00186 float
00187 __sinhf( float x )
00188 {
00189 #ifdef _32BIT_MACHINE
00190 
00191 int m, n;
00192 
00193 #else
00194 
00195 long long m, n;
00196 
00197 #endif
00198 
00199 int ix, xpt;
00200 double  result;
00201 double  nd;
00202 double  dx, absdx;
00203 double  y, ysq;
00204 double  sy, cy;
00205 #ifdef _CALL_MATHERR
00206 struct exception  exstruct;
00207 #endif
00208 
00209   /* extract exponent of x for some quick screening */
00210 
00211   FLT2INT(x, ix); /* copy arg to an integer */
00212   xpt = (ix >> MANTWIDTH);  /* shift off mantissa bits  */
00213   xpt &= 0xff;
00214 
00215   if ( xpt < 0x85 )
00216   {
00217     /* |x| < 64.0 */
00218 
00219     if ( xpt >= 0x73 )
00220     {
00221       /* |x| >= 2^(-12) */
00222 
00223       dx = x;
00224       absdx = fabs(dx);
00225 
00226       nd = absdx*rln2.d;
00227       n = ROUND(nd);
00228 
00229       /* n = NINT(|x|/log(2)) */
00230 
00231       /* reduce argument to +/- log(2)/2 */
00232 
00233       nd = n;
00234 
00235       y = absdx - nd*ln2.d;
00236 
00237       ysq = y*y;
00238 
00239       /* compute sinh(y) and cosh(y) - 1.0 */
00240 
00241       sy = ((S[3].d*ysq + S[2].d)*ysq + S[1].d)*(ysq*y) + y;
00242 
00243       cy = ((C[3].d*ysq + C[2].d)*ysq + C[1].d)*ysq + C[0].d;
00244 
00245       if ( n < 14 )
00246       {
00247         /* sinh(|x|) = sinh(y)*coshtab[n] + cosh(y)*sinhtab[n] */
00248 
00249         result = sy*coshtab[n].d + cy*sinhtab[n].d;
00250 
00251         if ( x < 0.0f )
00252           result = -result;
00253 
00254         return ( (float)result );
00255       }
00256 
00257 L:
00258       /* sinh(|x|) = 0.5*2^n*e^y */
00259 
00260       result = cy + sy;
00261 
00262       /* multiply result by 2**(n-1)  */
00263 
00264 #ifdef _32BIT_MACHINE
00265 
00266       DBLHI2INT(result, m);
00267       m += ((n - 1) << DMANTWIDTH);
00268       INT2DBLHI(m, result);
00269 #else
00270       DBL2LL(result, m);
00271       m += ((n - 1) << DMANTWIDTH);
00272       LL2DBL(m, result);
00273 #endif
00274       if ( x < 0.0f )
00275         result = -result;
00276 
00277       return ( (float)result );
00278     }
00279 
00280     return ( x );
00281   }
00282 
00283   if ( x != x )
00284   {
00285     /* x is a NaN; return a quiet NaN */
00286 
00287 #ifdef _CALL_MATHERR
00288 
00289                 exstruct.type = DOMAIN;
00290                 exstruct.name = "sinhf";
00291                 exstruct.arg1 = x;
00292                 exstruct.retval = Qnan.f;
00293 
00294                 if ( matherr( &exstruct ) == 0 )
00295                 {
00296                         fprintf(stderr, "domain error in sinhf\n");
00297                         SETERRNO(EDOM);
00298                 }
00299 
00300                 return ( exstruct.retval );
00301 #else
00302     NAN_SETERRNO(EDOM);
00303     
00304     return ( Qnan.f );
00305 #endif
00306   }
00307 
00308   if ( x > Ulimit.f )
00309   {
00310 #ifdef _CALL_MATHERR
00311 
00312                 exstruct.type = OVERFLOW;
00313                 exstruct.name = "sinhf";
00314                 exstruct.arg1 = x;
00315                 exstruct.retval = Inf.f;
00316 
00317                 if ( matherr( &exstruct ) == 0 )
00318                 {
00319                         SETERRNO(ERANGE);
00320                 }
00321 
00322                 return ( exstruct.retval );
00323 #else
00324     SETERRNO(ERANGE);
00325 
00326     return ( Inf.f );
00327 #endif
00328   }
00329 
00330   if ( x < Llimit.f )
00331   {
00332 #ifdef _CALL_MATHERR
00333 
00334                 exstruct.type = OVERFLOW;
00335                 exstruct.name = "sinhf";
00336                 exstruct.arg1 = x;
00337                 exstruct.retval = Neginf.f;
00338 
00339                 if ( matherr( &exstruct ) == 0 )
00340                 {
00341                         SETERRNO(ERANGE);
00342                 }
00343 
00344                 return ( exstruct.retval );
00345 #else
00346     SETERRNO(ERANGE);
00347 
00348     return ( Neginf.f );
00349 #endif
00350   }
00351 
00352   /* essentially the same algorithm as above, however here n > 14,
00353      so we can skip a couple of the tests 
00354   */
00355 
00356   dx = x;
00357   absdx = fabs(dx);
00358 
00359   nd = absdx*rln2.d;
00360   n = ROUND(nd);
00361   nd = n;
00362 
00363   y = absdx - nd*ln2.d;
00364 
00365   ysq = y*y;
00366 
00367   sy = ((S[3].d*ysq + S[2].d)*ysq + S[1].d)*(ysq*y) + y;
00368 
00369   cy = ((C[3].d*ysq + C[2].d)*ysq + C[1].d)*ysq + C[0].d;
00370 
00371   goto L;
00372 }
00373 
00374 
00375 /* ====================================================================
00376  *
00377  * FunctionName   coshf
00378  *
00379  * Description    computes hyperbolic cosine of arg
00380  *
00381  * ====================================================================
00382  */
00383 
00384 float
00385 __coshf( float x )
00386 {
00387 #ifdef _32BIT_MACHINE
00388 
00389 int m, n;
00390 
00391 #else
00392 
00393 long long m, n;
00394 
00395 #endif
00396 
00397 int ix, xpt;
00398 double  result;
00399 double  md;
00400 double  dx, absdx;
00401 double  y, ysq;
00402 double  sy, cy;
00403 #ifdef _CALL_MATHERR
00404 struct exception  exstruct;
00405 #endif
00406 
00407   /* extract exponent of x for some quick screening */
00408 
00409   FLT2INT(x, ix); /* copy arg to an integer */
00410   xpt = (ix >> MANTWIDTH);  /* shift off mantissa */
00411   xpt &= 0xff;
00412 
00413   if ( xpt < 0x85 )
00414   {
00415     /* |x| < 64.0 */
00416 
00417     if ( xpt >= 0x66 )
00418     {
00419       /* |x| >= 2^(-25) */
00420 
00421       dx = x;
00422 
00423       absdx = fabs(dx);
00424 
00425       md = absdx*rln2.d;
00426 
00427       n = ROUND(md);
00428 
00429       /* n = NINT(x/log(2)) */
00430 
00431       /* reduce argument to +/- log(2)/2 */
00432 
00433       md = n;
00434 
00435       y = absdx - md*ln2.d;
00436 
00437       ysq = y*y;
00438 
00439       /* compute sinh(y) and cosh(y) - 1.0 */
00440 
00441       sy = ((S[3].d*ysq + S[2].d)*ysq + S[1].d)*(ysq*y) + y;
00442 
00443       cy = ((C[3].d*ysq + C[2].d)*ysq + C[1].d)*ysq + C[0].d;
00444 
00445       if ( n < 14 )
00446       {
00447         /* cosh(x) = cosh(y)*coshtab[n] + sinh(y)*sinhtab[n] */
00448 
00449         result = sy*sinhtab[n].d + cy*coshtab[n].d;
00450 
00451         return ( (float)result );
00452       }
00453 
00454 L:
00455       /* cosh(x) = 0.5*2^n*e^y */
00456 
00457       result = cy + sy;
00458 
00459       /* multiply result by 2**(n-1)  */
00460 
00461 #ifdef _32BIT_MACHINE
00462 
00463       DBLHI2INT(result, m);
00464       m += ((n - 1) << DMANTWIDTH);
00465       INT2DBLHI(m, result);
00466 #else
00467       DBL2LL(result, m);
00468       m += ((n - 1) << DMANTWIDTH);
00469       LL2DBL(m, result);
00470 #endif
00471       return ( (float)result );
00472     }
00473 
00474     return ( f_one.f );
00475   }
00476 
00477   if ( x != x )
00478   {
00479     /* x is a NaN; return a quiet NaN */
00480 
00481 #ifdef _CALL_MATHERR
00482 
00483                 exstruct.type = DOMAIN;
00484                 exstruct.name = "coshf";
00485                 exstruct.arg1 = x;
00486                 exstruct.retval = Qnan.f;
00487 
00488                 if ( matherr( &exstruct ) == 0 )
00489                 {
00490                         fprintf(stderr, "domain error in coshf\n");
00491                         SETERRNO(EDOM);
00492                 }
00493 
00494                 return ( exstruct.retval );
00495 #else
00496     NAN_SETERRNO(EDOM);
00497     
00498     return ( Qnan.f );
00499 #endif
00500   }
00501 
00502   if ( fabsf(x) > Ulimit.f )
00503   {
00504 #ifdef _CALL_MATHERR
00505 
00506                 exstruct.type = OVERFLOW;
00507                 exstruct.name = "coshf";
00508                 exstruct.arg1 = x;
00509                 exstruct.retval = Inf.f;
00510 
00511                 if ( matherr( &exstruct ) == 0 )
00512                 {
00513                         SETERRNO(ERANGE);
00514                 }
00515 
00516                 return ( exstruct.retval );
00517 #else
00518     SETERRNO(ERANGE);
00519 
00520     return ( Inf.f );
00521 #endif
00522   }
00523 
00524   /* essentially the same algorithm as above, however here n > 14,
00525      so we can skip a couple of the tests 
00526   */
00527 
00528   dx = x;
00529 
00530   absdx = fabs(dx);
00531 
00532   md = absdx*rln2.d;
00533   n = ROUND(md);
00534   md = n;
00535 
00536   y = absdx - md*ln2.d;
00537 
00538   ysq = y*y;
00539 
00540   sy = ((S[3].d*ysq + S[2].d)*ysq + S[1].d)*(ysq*y) + y;
00541 
00542   cy = ((C[3].d*ysq + C[2].d)*ysq + C[1].d)*ysq + C[0].d;
00543 
00544   goto L;
00545 }
00546