Open64: osprey/libm/mips/j0.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: j0.c
00045  * $Revision: 1.5 $
00046  * $Date: 04/12/21 14:58:22-08:00 $
00047  * $Author: bos@eng-25.internal.keyresearch.com $
00048  * $Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.j0.c $
00049  *
00050  * Revision history:
00051  *  20-Mar-98 - Original Version
00052  *
00053  * Description: source code for bessel functions of the first and second
00054  *  kinds of order zero
00055  *
00056  *  Based on the ideas in Hart & Cheney's book, Computer Approximations,
00057  *  John Wiley and Sons, 1968, Section 6.8.
00058  *
00059  *  j0:
00060  *
00061  *  Calls sin, cos, sqrt.
00062  *
00063  *  y0:
00064  *
00065  *  Calls sin, cos, sqrt, log, j0.
00066  *
00067  *
00068  * ====================================================================
00069  * ====================================================================
00070  */
00071 
00072 static char *rcs_id = "$Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.j0.c $ $Revision: 1.5 $";
00073 
00074 /****************************  IMPORTANT NOTE  ****************************
00075  *
00076  *  - This file contains library routines which are "safe" for parallel
00077  *    processing.  This fact is indicated here and in 'math.h' by inclusion
00078  *    of the following directive(s):
00079  *
00080  *        #pragma no side effects (j0)
00081  *        #pragma no side effects (y0)
00082  *
00083  *  - To be considered safe a routine cannot do any I/O (including GL calls),
00084  *    modify any global location ('errno' for example) or have global state
00085  *    (as 'random' does).  Nor can it call any other routine which does any
00086  *    of these things.
00087  *
00088  *  - If in the future any of these routines are modified so that they are
00089  *    no longer "safe", then the appropriate 'pragma no side effects' must
00090  *    be manually removed from 'math.h' and this file.
00091  *
00092  *  - Note that this pragma only has meaning to Power C ('pca') and is
00093  *    properly passed on by 'cpp' & 'c++' and ignored by 'ccom'.
00094  *
00095  ******************************  END OF NOTE  ******************************/
00096 
00097 #ifdef _CALL_MATHERR
00098 #include <stdio.h>
00099 #include <math.h>
00100 #include <errno.h>
00101 #endif
00102 
00103 #include "libm.h"
00104 
00105 #if defined(mips) && !defined(__GNUC__)
00106 extern double j0(double);
00107 extern double y0(double);
00108 
00109 #pragma weak j0 = __j0
00110 #pragma weak y0 = __y0
00111 #endif
00112 
00113 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00114 extern double __j0(double);
00115 extern double __y0(double);
00116 #pragma weak j0
00117 #pragma weak y0
00118 double j0( double arg ) {
00119   return __j0( arg );
00120 }
00121 double y0( double x ) {
00122   return __y0( x );
00123 }
00124 #elif defined(__GNUC__)
00125 extern  double  __j0(double);
00126 
00127 double    j0() __attribute__ ((weak, alias ("__j0")));
00128 
00129 extern  double  __y0(double);
00130 
00131 double    y0() __attribute__ ((weak, alias ("__y0")));
00132 
00133 #endif
00134 
00135 extern double __sin(double);
00136 extern double __cos(double);
00137 extern double __log(double);
00138 
00139 static double p0(double), q0(double);
00140 
00141 static const du twobypi = /* 2.0/pi */
00142 {
00143  D(0x3fe45f30, 0x6dc9c883)
00144 };
00145 
00146 static const du rsqrtpi = /* 1.0/sqrt(pi) */
00147 {
00148  D(0x3fe20dd7, 0x50429b6d)
00149 };
00150 
00151 static const du k0 =  /* 2.0/pi*(euler - log(2.0)) */
00152 {
00153  D(0xbfb2e4d6, 0x99cbd01f)
00154 };
00155 
00156 /* coefficients for rational approximation of j0
00157    on +/- 8.0
00158 */
00159 
00160 static const du s1[] =
00161 {
00162  D(0xbfd00000, 0x00000000),
00163  D(0x3f8a5dd2, 0x07400a00),
00164  D(0xbf322e84, 0x5625ccf7),
00165  D(0x3ec86226, 0xdcfbebcf),
00166  D(0xbe519953, 0xd3612c5b),
00167  D(0x3dcbf06d, 0xa8ff6623),
00168  D(0xbd371d0d, 0x96e69a82),
00169  D(0x3c8f6396, 0x070cb706),
00170 };
00171 
00172 static const du t1[] =
00173 {
00174  D(0x3ff00000, 0x00000000),
00175  D(0x3f8688b7, 0xe2ffd7ff),
00176  D(0x3f101192, 0xb33f1e37),
00177  D(0x3e8eac3a, 0x456bf6e2),
00178  D(0x3e05be02, 0xd273ac07),
00179  D(0x3d77e0cb, 0x01046ec0),
00180  D(0x3ce45590, 0x9323449e),
00181  D(0x3c497437, 0x5fe08079),
00182  D(0x3ba31c0e, 0x7b50a76b),
00183 };
00184 
00185 /* coefficients for rational approximation of p0 on (0.0, 1.0) */
00186 
00187 static const du s2[] =
00188 {
00189  D(0xbf51ffff, 0xffffffb8),
00190  D(0xbf6469e4, 0x36e62881),
00191  D(0xbf5b424a, 0x11c914b4),
00192  D(0xbf38b387, 0x1e57cade),
00193  D(0xbef99508, 0xd4953bc5),
00194  D(0xbe89282f, 0x2055d015),
00195 };
00196 
00197 static const du t2[] =
00198 {
00199  D(0x3ff00000, 0x00000000),
00200  D(0x40025847, 0x4d3e52b1),
00201  D(0x3ff91cde, 0x901a39c4),
00202  D(0x3fd83087, 0x0c521255),
00203  D(0x3f9dde76, 0xb8de94f3),
00204  D(0x3f400cef, 0x4c143276),
00205  D(0x3e89282f, 0x2055d015),
00206 };
00207 
00208 /* coefficients for rational approximation of q0 on (0.0, 1.0) */
00209 
00210 static const du s3[] = 
00211 {
00212  D(0xbf900000, 0x00000000),
00213  D(0x3f22bfff, 0xfffffe5d),
00214  D(0x3f381dc6, 0x0bdd69d3),
00215  D(0x3f328f71, 0xc888d5b2),
00216  D(0x3f13cc8a, 0x06fa2970),
00217  D(0x3ed8e553, 0xf75ae5b9),
00218  D(0x3e70a26b, 0xd3f48e4f),
00219 };
00220 
00221 static const du t3[] = 
00222 {
00223  D(0x3ff00000, 0x00000000),
00224  D(0x4004f783, 0x06b5e18d),
00225  D(0x4000ced3, 0x6118ee31),
00226  D(0x3fe3b5ca, 0x55784db2),
00227  D(0x3fafcd81, 0xa67601cb),
00228  D(0x3f5a87e4, 0x0ce687ca),
00229  D(0x3ed0a26b, 0xd3f48e4f),
00230 };
00231 
00232 /* coefficients for rational approximation of y0 on (0.0, 8.0) */
00233 
00234 static const du s4[] =
00235 {
00236  D(0xc0034066, 0x8510b96f),
00237  D(0x3fc86981, 0x11fdbb9b),
00238  D(0xbf747d16, 0xe102a52f),
00239  D(0x3f0f6db4, 0x4d204dbf),
00240  D(0xbe99068d, 0xa7ea12f5),
00241  D(0x3e1577e8, 0xd5bafc4a),
00242  D(0xbd82ec16, 0x2f225a36),
00243  D(0x3cdb1b0b, 0x6c893369),
00244 };
00245 
00246 static const du t4[] =
00247 {
00248  D(0x3ff00000, 0x00000000),
00249  D(0x3f870a22, 0xf33840f5),
00250  D(0x3f10d4c2, 0x2520436c),
00251  D(0x3e907f07, 0xa7ebc30f),
00252  D(0x3e0815f6, 0x8930073b),
00253  D(0x3d7b583a, 0x90523001),
00254  D(0x3ce831c7, 0xe53c886c),
00255  D(0x3c4fb50f, 0x5f529cf6),
00256  D(0x3ba94ca1, 0x99c520de),
00257 };
00258 
00259 static const du Qnan =
00260 {D(QNANHI, QNANLO)};
00261 
00262 static const du Neginf =
00263 {D(0xfff00000, 0x00000000)};
00264 
00265 
00266 /* ====================================================================
00267  *
00268  * FunctionName   j0
00269  *
00270  * Description    computes bessel function of first kind of order 0
00271  *
00272  * ====================================================================
00273  */
00274 
00275 double
00276 __j0( double arg )
00277 {
00278 #ifdef _32BIT_MACHINE
00279 int ix, xpt;
00280 #else
00281 long long ix, xpt;
00282 #endif
00283 double x;
00284 double xsq, rx;
00285 double num, denom;
00286 double  s, c, t;
00287 double result;
00288 #ifdef _CALL_MATHERR
00289 struct exception  exstruct;
00290 #endif
00291 
00292   if( arg != arg )
00293   {
00294     /* arg is NaN */
00295 
00296 #ifdef _CALL_MATHERR
00297 
00298     exstruct.type = DOMAIN;
00299     exstruct.name = "j0";
00300     exstruct.arg1 = arg;
00301     exstruct.retval = Qnan.d;
00302 
00303     if ( matherr( &exstruct ) == 0 )
00304     {
00305       fprintf(stderr, "domain error in j0\n");
00306       SETERRNO(EDOM);
00307     }
00308 
00309     return ( exstruct.retval );
00310 #else
00311     NAN_SETERRNO(EDOM);
00312 
00313     return ( Qnan.d );
00314 #endif
00315   }
00316 
00317   x = fabs(arg);
00318 
00319   /* extract exponent of x for some quick screening */
00320 
00321 #ifdef _32BIT_MACHINE
00322   DBLHI2INT(x, ix); /* copy hi part of x to ix */
00323 #else
00324   DBL2LL(x, ix);    /* copy x to ix */
00325 #endif
00326 
00327   xpt = (ix >> DMANTWIDTH);
00328   xpt &= 0x7ff;
00329 
00330   if ( xpt < 0x402 )
00331   {
00332     /*      |x| < 8.0 */
00333 
00334     if ( xpt >= 0x3e4 )
00335     {
00336       /* |x| >= 2^(-27) */
00337 
00338       rx = x;
00339       xsq = rx*rx;
00340 
00341       num = ((((((s1[7].d*xsq + s1[6].d)*xsq +
00342         s1[5].d)*xsq + s1[4].d)*xsq + s1[3].d)*xsq +
00343         s1[2].d)*xsq + s1[1].d)*xsq + s1[0].d;
00344 
00345       denom = (((((((t1[8].d*xsq + t1[7].d)*xsq + t1[6].d)*xsq +
00346         t1[5].d)*xsq + t1[4].d)*xsq + t1[3].d)*xsq +
00347         t1[2].d)*xsq + t1[1].d)*xsq + t1[0].d;
00348 
00349       result = 1.0 + xsq*(num/denom);
00350 
00351       return ( (double)result );
00352     }
00353     else
00354     {
00355       return ( 1.0 );
00356     }
00357   }
00358   else if ( xpt < 0x430 )
00359   {
00360     rx = x;
00361     c = __cos(x);
00362     s = __sin(x);
00363     t = 8.0/rx;
00364   
00365     /* j0(x) = 1/sqrt(pi)*(p0(t)*(s+c) - q0(t)*(s-c))/sqrt(x) */
00366   
00367     /* Use identity sin(x) +/- cos(x) =
00368        -cos(2*x)/(sin(x) -/+ cos(x))
00369        to rewrite s +/- c to avoid loss of significance due to
00370        subtraction.  Note that sin and cos are only defined
00371        out to 2**50 on this system.
00372     */
00373   
00374     if ( (c*s) >= 0.0 )
00375     {
00376       result = 
00377       p0(t)*(s + c) + q0(t)*(__cos(x + x)/(s + c));
00378     }
00379     else
00380     {
00381       result = 
00382       -p0(t)*(__cos(x + x)/(s - c)) - q0(t)*(s - c);
00383     }
00384   
00385     return ( rsqrtpi.d*result/sqrt(rx) );
00386   }
00387   else
00388   {
00389     /* |x| exceeds 2**49; return zero  */
00390 
00391 #ifdef _CALL_MATHERR
00392 
00393     exstruct.type = TLOSS;
00394     exstruct.name = "j0";
00395     exstruct.arg1 = arg;
00396     exstruct.retval = 0.0;
00397 
00398     if ( matherr( &exstruct ) == 0 )
00399     {
00400       fprintf(stderr, "total loss of significance \
00401 error in j0\n");
00402       SETERRNO(ERANGE);
00403     }
00404 
00405     return ( exstruct.retval );
00406 #else
00407     SETERRNO(ERANGE);
00408 
00409     return ( 0.0 );
00410 #endif
00411   }
00412 
00413 }
00414 
00415 
00416 /* ====================================================================
00417  *
00418  * FunctionName   y0
00419  *
00420  * Description    computes bessel function of second kind of order 0
00421  *
00422  * ====================================================================
00423  */
00424 
00425 double
00426 __y0( double x )
00427 {
00428 #ifdef _32BIT_MACHINE
00429 int ix, xpt;
00430 #else
00431 long long ix, xpt;
00432 #endif
00433 double  rx;
00434 double  xsq, t;
00435 double  num, denom;
00436 double  s, c;
00437 double  result;
00438 #ifdef _CALL_MATHERR
00439 struct exception  exstruct;
00440 #endif
00441 
00442   if( x != x )
00443   {
00444     /* arg is NaN */
00445 
00446 #ifdef _CALL_MATHERR
00447 
00448     exstruct.type = DOMAIN;
00449     exstruct.name = "y0";
00450     exstruct.arg1 = x;
00451     exstruct.retval = Qnan.d;
00452 
00453     if ( matherr( &exstruct ) == 0 )
00454     {
00455       fprintf(stderr, "domain error in y0\n");
00456       SETERRNO(EDOM);
00457     }
00458 
00459     return ( exstruct.retval );
00460 #else
00461     NAN_SETERRNO(EDOM);
00462 
00463     return ( Qnan.d );
00464 #endif
00465   }
00466 
00467   if ( x == 0.0 )
00468   {
00469 #ifdef _CALL_MATHERR
00470 
00471     exstruct.type = OVERFLOW;
00472     exstruct.name = "y0";
00473     exstruct.arg1 = x;
00474     exstruct.retval = Neginf.d;
00475 
00476     if ( matherr( &exstruct ) == 0 )
00477     {
00478       fprintf(stderr, "overflow range error in y0\n");
00479       SETERRNO(ERANGE);
00480     }
00481 
00482     return ( exstruct.retval );
00483 #else
00484     SETERRNO(ERANGE);
00485 
00486     return ( Neginf.d );
00487 #endif
00488   }
00489 
00490   if ( x < 0.0 )
00491   {
00492 #ifdef _CALL_MATHERR
00493 
00494     exstruct.type = DOMAIN;
00495     exstruct.name = "y0";
00496     exstruct.arg1 = x;
00497     exstruct.retval = Qnan.d;
00498 
00499     if ( matherr( &exstruct ) == 0 )
00500     {
00501       fprintf(stderr, "domain error in y0\n");
00502       SETERRNO(EDOM);
00503     }
00504 
00505     return ( exstruct.retval );
00506 #else
00507     SETERRNO(EDOM);
00508 
00509     return ( Qnan.d );
00510 #endif
00511   }
00512 
00513   /* extract exponent of x for some quick screening */
00514 
00515 #ifdef _32BIT_MACHINE
00516   DBLHI2INT(x, ix); /* copy hi part of x to ix */
00517 #else
00518   DBL2LL(x, ix);    /* copy x to ix */
00519 #endif
00520 
00521   xpt = (ix >> DMANTWIDTH);
00522   xpt &= 0x7ff;
00523 
00524   if ( xpt < 0x402 )
00525   {
00526     /*      |x| < 8.0 */
00527 
00528     if ( xpt >= 0x3e4 )
00529     {
00530       /* |x| >= 2^(-27) */
00531 
00532       rx = x;
00533       xsq = rx*rx;
00534 
00535       num = ((((((s4[7].d*xsq + s4[6].d)*xsq + s4[5].d)*xsq +
00536         s4[4].d)*xsq + s4[3].d)*xsq + s4[2].d)*xsq +
00537         s4[1].d)*xsq + s4[0].d;
00538 
00539       denom = (((((((t4[8].d*xsq + t4[7].d)*xsq + t4[6].d)*xsq +
00540         t4[5].d)*xsq + t4[4].d)*xsq + t4[3].d)*xsq +
00541         t4[2].d)*xsq + t4[1].d)*xsq + t4[0].d;
00542 
00543       result = k0.d + k0.d*(xsq*num/denom) +
00544         twobypi.d*(__j0(x)*__log(x));
00545 
00546       return ( (double)result );
00547     }
00548     else
00549     {
00550       result = k0.d + twobypi.d*__log(x);
00551 
00552       return ( (double)result );
00553     }
00554   }
00555   else if ( xpt < 0x430 )
00556   {
00557     rx = x;
00558     c = __cos(x);
00559     s = __sin(x);
00560     t = 8.0/rx;
00561   
00562     /* y0(x) = 1/sqrt(pi)*(p0(t)*(s-c) + q0(t)*(s+c))/sqrt(x) */
00563   
00564     /* Use identity sin(x) +/- cos(x) =
00565        -cos(2*x)/(sin(x) -/+ cos(x))
00566        to rewrite s +/- c to avoid loss of significance due to
00567        subtraction.  Note that sin and cos are only defined
00568        out to 2**50 on this system.
00569     */
00570   
00571     if ( (c*s) >= 0.0 )
00572     {
00573       result = 
00574       -p0(t)*(__cos(x + x)/(s + c)) + q0(t)*(s + c);
00575     }
00576     else
00577     {
00578       result = 
00579       p0(t)*(s - c) - q0(t)*(__cos(x + x)/(s - c));
00580     }
00581   
00582     return ( rsqrtpi.d*result/sqrt(rx) );
00583   }
00584   else
00585   {
00586     /* |x| exceeds 2**49; return zero  */
00587 
00588 #ifdef _CALL_MATHERR
00589 
00590     exstruct.type = TLOSS;
00591     exstruct.name = "y0";
00592     exstruct.arg1 = x;
00593     exstruct.retval = 0.0;
00594 
00595     if ( matherr( &exstruct ) == 0 )
00596     {
00597       fprintf(stderr, "total loss of significance \
00598 error in y0\n");
00599       SETERRNO(ERANGE);
00600     }
00601 
00602     return ( exstruct.retval );
00603 #else
00604     SETERRNO(ERANGE);
00605 
00606     return ( 0.0 );
00607 #endif
00608   }
00609 }
00610 
00611 static
00612 double
00613 p0( double t )
00614 {
00615 double  tsq;
00616 double  num, denom, result;
00617 
00618   tsq = t*t;
00619 
00620   num = ((((s2[5].d*tsq + s2[4].d)*tsq + s2[3].d)*tsq +
00621       s2[2].d)*tsq + s2[1].d)*tsq + s2[0].d;
00622 
00623   denom = (((((t2[6].d*tsq + t2[5].d)*tsq + t2[4].d)*tsq +
00624       t2[3].d)*tsq  + t2[2].d)*tsq + t2[1].d)*tsq + 
00625       t2[0].d;
00626 
00627 
00628   result = 1.0 + tsq*(num/denom);
00629 
00630   return ( result );
00631 }
00632 
00633 static
00634 double
00635 q0( double t )
00636 {
00637 double  tsq;
00638 double  num, denom, result;
00639 
00640   tsq = t*t;
00641 
00642   num = ((((s3[6].d*tsq + s3[5].d)*tsq + s3[4].d)*tsq +
00643       s3[3].d)*tsq + s3[2].d)*tsq + s3[1].d;
00644 
00645   denom = (((((t3[6].d*tsq + t3[5].d)*tsq + t3[4].d)*tsq +
00646       t3[3].d)*tsq  + t3[2].d)*tsq + t3[1].d)*tsq + 
00647       t3[0].d;
00648 
00649 
00650   result = s3[0].d*t + tsq*t*num/denom;
00651 
00652   return ( result );
00653 }
00654 
00655 #ifdef NO_LONG_DOUBLE
00656 
00657 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00658 extern  long double  __j0l(long double);
00659 extern  long double  __y0l(long double);
00660 long double j0l( long double x ) {  
00661   return ( (long double)__j0((double)x) );
00662 }
00663 
00664 long double y0l( long double x ) {  
00665   return ( (long double)__y0((double)x) );
00666 }
00667 #elif defined(__GNUC__)
00668 extern  long double  __j0l(long double);
00669 
00670 long double    j0l() __attribute__ ((weak, alias ("__j0l")));
00671 
00672 extern  long double  __y0l(long double);
00673 
00674 long double    y0l() __attribute__ ((weak, alias ("__y0l")));
00675 
00676 #endif
00677 
00678 long double
00679 __j0l( long double x )
00680 { 
00681   return ( (long double)__j0((double)x) );
00682 }
00683 
00684 long double
00685 __y0l( long double x )
00686 { 
00687   return ( (long double)__y0((double)x) );
00688 }
00689 
00690 #endif
00691