Open64: osprey/libm/mips/j1.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: j1.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.j1.c $
00049  *
00050  * Revision history:
00051  *  28-May-93 - Original Version
00052  *
00053  * Description: source code for bessel functions of first and second kind
00054  *    of order 1
00055  *
00056  * ====================================================================
00057  * ====================================================================
00058  */
00059 
00060 /* $Header: /home/bos/bk/kpro64-pending/libm/mips/j1.c 1.5 04/12/21 14:58:22-08:00 bos@eng-25.internal.keyresearch.com $ */
00061 
00062 /****************************  IMPORTANT NOTE  ****************************
00063  *
00064  *  - This file contains library routines which are "safe" for parallel
00065  *    processing.  This fact is indicated here and in 'math.h' by inclusion
00066  *    of the following directive(s):
00067  *
00068  *        #pragma no side effects (j1)
00069  *        #pragma no side effects (y1)
00070  *
00071  *  - To be considered safe a routine cannot do any I/O (including GL calls),
00072  *    modify any global location ('errno' for example) or have global state
00073  *    (as 'random' does).  Nor can it call any other routine which does any
00074  *    of these things.
00075  *
00076  *  - If in the future any of these routines are modified so that they are
00077  *    no longer "safe", then the appropriate 'pragma no side effects' must
00078  *    be manually removed from 'math.h' and this file.
00079  *
00080  *  - Note that this pragma only has meaning to Power C ('pca') and is
00081  *    properly passed on by 'cpp' & 'c++' and ignored by 'ccom'.
00082  *
00083  ******************************  END OF NOTE  ******************************/
00084 /*
00085   floating point Bessel's function
00086   of the first and second kinds
00087   of order one
00088 
00089   j1:
00090 
00091   Calls sin, cos, sqrt.
00092 
00093   Coefficients are from Hart & Cheney.
00094   #6050 (20.98D)
00095   #6750 (19.19D)
00096   #7150 (19.35D)
00097 
00098   y1:
00099 
00100   Calls sin, cos, sqrt, log, j1.
00101 
00102   Coefficients are from Hart & Cheney.
00103   #6447 (22.18D)
00104   #6750 (19.19D)
00105   #7150 (19.35D)
00106 */
00107 
00108 #ifdef _CALL_MATHERR
00109 #include <stdio.h>
00110 #include <math.h>
00111 #endif
00112 
00113 #include "libm.h"
00114 
00115 #if defined(mips) && !defined(__GNUC__)
00116 extern  double  j1(double);
00117 extern  double  y1(double);
00118 
00119 #pragma weak j1 = __j1
00120 #pragma weak y1 = __y1
00121 #endif
00122 
00123 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00124 extern double __j1(double);
00125 extern double __y1(double);
00126 #pragma weak j1
00127 #pragma weak y1
00128 double j1( double arg ) {
00129   return __j1( arg );
00130 }
00131 double y1( double arg ) {
00132   return __y1( arg );
00133 }
00134 #elif defined(__GNUC__)
00135 extern  double  __j1(double);
00136 
00137 double    j1() __attribute__ ((weak, alias ("__j1")));
00138 
00139 extern  double  __y1(double);
00140 
00141 double    y1() __attribute__ ((weak, alias ("__y1")));
00142 
00143 #endif
00144 
00145 extern  double  __sin(double);
00146 extern  double  __cos(double);
00147 extern  double  __log(double);
00148 
00149 static  double p1(double), q1(double);
00150 
00151 static const  du  Qnan =
00152 {D(QNANHI, QNANLO)};
00153 
00154 static const  du  Neginf =
00155 {D(0xfff00000, 0x00000000)};
00156 
00157 static  const double twobypi =
00158   0.63661977236758134;
00159 
00160 static  const double rsqrtpi =
00161   0.56418958354775629;
00162 
00163 static double s1[] = {
00164   0.581199354001606143928050809e21,
00165   -.6672106568924916298020941484e20,
00166   0.2316433580634002297931815435e19,
00167   -.3588817569910106050743641413e17,
00168   0.2908795263834775409737601689e15,
00169   -.1322983480332126453125473247e13,
00170   0.3413234182301700539091292655e10,
00171   -.4695753530642995859767162166e7,
00172   0.2701122710892323414856790990e4,
00173 };
00174 
00175 static double t1[] = {
00176   0.1162398708003212287858529400e22,
00177   0.1185770712190320999837113348e20,
00178   0.6092061398917521746105196863e17,
00179   0.2081661221307607351240184229e15,
00180   0.5243710262167649715406728642e12,
00181   0.1013863514358673989967045588e10,
00182   0.1501793594998585505921097578e7,
00183   0.1606931573481487801970916749e4,
00184   1.0,
00185 };
00186 
00187 static double s2[] = {
00188          1.0,
00189          2.2467834456732897,
00190          1.4955846127260149,
00191          3.4576449324751314e-1,
00192          2.4969642938750994e-2,
00193          3.6696523370953188e-4,
00194 };
00195 
00196 static double t2[] = {
00197          1.0,
00198          2.2449523909857897,
00199          1.4915091861206051,
00200          3.4310990865393919e-1,
00201          2.4388513111364209e-2,
00202          3.3134568440823793e-4,
00203         -2.2792852359037530e-7,
00204 };
00205 
00206 static double s3[] = {
00207          4.6875e-2,
00208          1.2043362532202404e-1,
00209          9.3819655587348968e-2,
00210          2.6273811593862773e-2,
00211          2.4287960814178163e-3,
00212          5.0281167606313395e-5,
00213 };
00214 
00215 static double t3[] = {
00216          1.0,
00217          2.5735231344740127,
00218          2.0123005490948824,
00219          5.6866066418403355e-1,
00220          5.3928104387696746e-2,
00221          1.2314671563610277e-3,
00222          1.4272352241007432e-6,
00223 };
00224 
00225 static double s4[] = {
00226   -.9963753424306922225996744354e23,
00227   0.2655473831434854326894248968e23,
00228   -.1212297555414509577913561535e22,
00229   0.2193107339917797592111427556e20,
00230   -.1965887462722140658820322248e18,
00231   0.9569930239921683481121552788e15,
00232   -.2580681702194450950541426399e13,
00233   0.3639488548124002058278999428e10,
00234   -.2108847540133123652824139923e7,
00235 };
00236 
00237 static double t4[] = {
00238   0.5082067366941243245314424152e24,
00239   0.5435310377188854170800653097e22,
00240   0.2954987935897148674290758119e20,
00241   0.1082258259408819552553850180e18,
00242   0.2976632125647276729292742282e15,
00243   0.6465340881265275571961681500e12,
00244   0.1128686837169442121732366891e10,
00245   0.1563282754899580604737366452e7,
00246   0.1612361029677000859332072312e4,
00247   1.0,
00248 };
00249 
00250 
00251 /* ====================================================================
00252  *
00253  * FunctionName   j1
00254  *
00255  * Description    computes bessel function of first kind of order 1
00256  *
00257  * ====================================================================
00258  */
00259 
00260 double
00261 __j1( arg )
00262 double arg;
00263 {
00264 #ifdef _32BIT_MACHINE
00265 int ix, xpt;
00266 #else
00267 long long ix, xpt;
00268 #endif
00269 
00270 double x, xsq;
00271 double num, denom;
00272 double c, s, t;
00273 double result;
00274 #ifdef _CALL_MATHERR
00275 struct exception  exstruct;
00276 #endif
00277 
00278   if( arg != arg )
00279   {
00280     /* arg is NaN */
00281 
00282 #ifdef _CALL_MATHERR
00283 
00284     exstruct.type = DOMAIN;
00285     exstruct.name = "j1";
00286     exstruct.arg1 = arg;
00287     exstruct.retval = Qnan.d;
00288 
00289     if ( matherr( &exstruct ) == 0 )
00290     {
00291       fprintf(stderr, "domain error in j1\n");
00292       SETERRNO(EDOM);
00293     }
00294 
00295     return ( exstruct.retval );
00296 #else
00297     NAN_SETERRNO(EDOM);
00298 
00299     return ( Qnan.d );
00300 #endif
00301   }
00302 
00303   /* extract exponent of x for some quick screening */
00304 
00305 #ifdef _32BIT_MACHINE
00306 
00307   DBLHI2INT(arg, ix); /* copy MSW of arg to ix  */
00308 #else
00309   DBL2LL(arg, ix);  /* copy arg to ix */
00310 #endif
00311 
00312   xpt = (ix >> DMANTWIDTH);
00313   xpt &= 0x7ff;
00314 
00315   x = fabs(arg);
00316 
00317   if ( xpt < 0x402 )
00318   {
00319     /*      |x| < 8.0 */
00320 
00321     if ( xpt >= 0x3e4 )
00322     {
00323       /* |x| >= 2^(-27) */
00324 
00325       xsq = x*x;
00326 
00327       num = (((((((s1[8]*xsq  + s1[7])*xsq + s1[6])*xsq +
00328         s1[5])*xsq + s1[4])*xsq + s1[3])*xsq +
00329         s1[2])*xsq + s1[1])*xsq + s1[0];
00330 
00331       denom = (((((((t1[8]*xsq  + t1[7])*xsq + t1[6])*xsq +
00332         t1[5])*xsq + t1[4])*xsq + t1[3])*xsq +
00333         t1[2])*xsq + t1[1])*xsq + t1[0];
00334 
00335       result = x*num/denom;
00336     }
00337     else
00338     {
00339       result = 0.5*x;
00340     }
00341   }
00342   else if ( xpt < 0x430 )
00343   {
00344 
00345     c = __cos(x);
00346     s = __sin(x);
00347     t = 8.0/x;
00348   
00349     /* j1(x) = 1/sqrt(pi)*(p1(t)*(s-c) + q1(t)*(s+c))/sqrt(x) */
00350   
00351     /* Use identity sin(x) +/- cos(x) =
00352        -cos(2*x)/(sin(x) -/+ cos(x))
00353        to rewrite s +/- c to avoid loss of significance due to
00354        subtraction.  Note that sin and cos are only defined
00355        out to 2**50 on this system.
00356     */
00357   
00358     if ( (c*s) >= 0.0 )
00359     {
00360       result = -p1(t)*(__cos(x + x)/(s + c)) + q1(t)*(s + c);
00361     }
00362     else
00363     {
00364       result = p1(t)*(s - c) - q1(t)*(__cos(x + x)/(s - c));
00365     }
00366   
00367     result = rsqrtpi*result/sqrt(x);
00368   }
00369   else
00370   {
00371 #ifdef _CALL_MATHERR
00372 
00373     exstruct.type = TLOSS;
00374     exstruct.name = "j1";
00375     exstruct.arg1 = arg;
00376     exstruct.retval = 0.0;
00377 
00378     if ( matherr( &exstruct ) == 0 )
00379     {
00380       fprintf(stderr, "total loss of significance \
00381 error in j1\n");
00382       SETERRNO(ERANGE);
00383     }
00384 
00385     return ( exstruct.retval );
00386 #else
00387     SETERRNO(ERANGE);
00388 
00389     return ( 0.0 );
00390 #endif
00391   }
00392 
00393   if ( ix >= 0 )
00394     return ( result );
00395   else
00396     return ( -result );
00397 
00398 }
00399 
00400 
00401 /* ====================================================================
00402  *
00403  * FunctionName   y1
00404  *
00405  * Description    computes bessel function of second kind of order 1
00406  *
00407  * ====================================================================
00408  */
00409 
00410 double
00411 __y1( arg )
00412 double arg;
00413 {
00414 #ifdef _32BIT_MACHINE
00415 int ix, xpt;
00416 #else
00417 long long ix, xpt;
00418 #endif
00419 
00420 double x, xsq;
00421 double num, denom, poly;
00422 double s, c, t;
00423 double result;
00424 #ifdef _CALL_MATHERR
00425 struct exception  exstruct;
00426 #endif
00427 
00428   if( arg != arg )
00429   {
00430     /* arg is NaN */
00431 
00432 #ifdef _CALL_MATHERR
00433 
00434     exstruct.type = DOMAIN;
00435     exstruct.name = "y1";
00436     exstruct.arg1 = arg;
00437     exstruct.retval = Qnan.d;
00438 
00439     if ( matherr( &exstruct ) == 0 )
00440     {
00441       fprintf(stderr, "domain error in y1\n");
00442       SETERRNO(EDOM);
00443     }
00444 
00445     return ( exstruct.retval );
00446 #else
00447     NAN_SETERRNO(EDOM);
00448 
00449     return ( Qnan.d );
00450 #endif
00451   }
00452 
00453   if ( arg == 0.0 )
00454   {
00455 #ifdef _CALL_MATHERR
00456 
00457     exstruct.type = OVERFLOW;
00458     exstruct.name = "y1";
00459     exstruct.arg1 = arg;
00460     exstruct.retval = Neginf.d;
00461 
00462     if ( matherr( &exstruct ) == 0 )
00463     {
00464       fprintf(stderr, "overflow range error in y1\n");
00465       SETERRNO(ERANGE);
00466     }
00467 
00468     return ( exstruct.retval );
00469 #else
00470     SETERRNO(ERANGE);
00471 
00472     return ( Neginf.d );
00473 #endif
00474   }
00475 
00476   if ( arg < 0.0 )
00477   {
00478 #ifdef _CALL_MATHERR
00479 
00480     exstruct.type = DOMAIN;
00481     exstruct.name = "y1";
00482     exstruct.arg1 = arg;
00483     exstruct.retval = Qnan.d;
00484 
00485     if ( matherr( &exstruct ) == 0 )
00486     {
00487       fprintf(stderr, "domain error in y1\n");
00488       SETERRNO(EDOM);
00489     }
00490 
00491     return ( exstruct.retval );
00492 #else
00493     SETERRNO(EDOM);
00494 
00495     return ( Qnan.d );
00496 #endif
00497   }
00498 
00499   x = arg;
00500 
00501   /* extract exponent of x for some quick screening */
00502 
00503 #ifdef _32BIT_MACHINE
00504 
00505   DBLHI2INT(x, ix); /* copy MSW of x to ix  */
00506 #else
00507   DBL2LL(x, ix);    /* copy x to ix */
00508 #endif
00509 
00510   xpt = (ix >> DMANTWIDTH);
00511   xpt &= 0x7ff;
00512 
00513   if ( xpt < 0x402 )
00514   {
00515     /*      |x| < 8.0 */
00516 
00517     if ( xpt >= 0x3c9 )
00518     {
00519       /* |x| >= 2^(-54) */
00520 
00521       xsq = x*x;
00522 
00523       num = (((((((s4[8]*xsq + s4[7])*xsq + s4[6])*xsq +
00524         s4[5])*xsq + s4[4])*xsq + s4[3])*xsq +
00525         s4[2])*xsq + s4[1])*xsq + s4[0];
00526 
00527       denom = ((((((((t4[9]*xsq + t4[8])*xsq + t4[7])*xsq +
00528         t4[6])*xsq + t4[5])*xsq + t4[4])*xsq + t4[3])*xsq +
00529         t4[2])*xsq + t4[1])*xsq + t4[0];
00530 
00531       poly = num/denom;
00532 
00533       result = x*poly + twobypi*(__j1(x)*__log(x) - 1.0/x);
00534 
00535       return ( result );
00536     }
00537     else
00538     {
00539       return ( -twobypi/x );
00540     }
00541 
00542   }
00543   else if ( xpt < 0x430 )
00544   {
00545     c = __cos(x);
00546     s = __sin(x);
00547     t = 8.0/x;
00548   
00549     /* y1(x) = 1/sqrt(pi)*(-p1(t)*(s+c) + q1(t)*(s-c))/sqrt(x) */
00550   
00551     /* Use identity sin(x) +/- cos(x) =
00552        -cos(2*x)/(sin(x) -/+ cos(x))
00553        to rewrite s +/- c to avoid loss of significance due to
00554        subtraction.  Note that sin and cos are only defined
00555        out to 2**50 on this system.
00556     */
00557   
00558     if ( (c*s) >= 0.0 )
00559     {
00560       result = -p1(t)*(s + c) - q1(t)*(__cos(x + x)/(s + c));
00561     }
00562     else
00563     {
00564       result = p1(t)*(__cos(x + x)/(s - c)) + q1(t)*(s - c);
00565     }
00566   
00567     return ( rsqrtpi*result/sqrt(arg) );
00568   }
00569   else
00570   {
00571     /* |x| exceeds 2**49; return zero  */
00572 
00573 #ifdef _CALL_MATHERR
00574 
00575     exstruct.type = TLOSS;
00576     exstruct.name = "y1";
00577     exstruct.arg1 = x;
00578     exstruct.retval = 0.0;
00579 
00580     if ( matherr( &exstruct ) == 0 )
00581     {
00582       fprintf(stderr, "total loss of significance \
00583 error in y1\n");
00584       SETERRNO(ERANGE);
00585     }
00586 
00587     return ( exstruct.retval );
00588 #else
00589     SETERRNO(ERANGE);
00590 
00591     return ( 0.0 );
00592 #endif
00593   }
00594 }
00595 
00596 static
00597 double
00598 p1( double t )
00599 {
00600 double  tsq;
00601 double  num, denom, result;
00602 
00603   tsq = t*t;
00604 
00605   num = ((((s2[5]*tsq + s2[4])*tsq + s2[3])*tsq +
00606     s2[2])*tsq + s2[1])*tsq + s2[0];
00607 
00608   denom = (((((t2[6]*tsq + t2[5])*tsq + t2[4])*tsq +
00609     t2[3])*tsq + t2[2])*tsq + t2[1])*tsq + t2[0];
00610 
00611   result = num/denom;
00612 
00613   return ( result );
00614 }
00615 
00616 static
00617 double
00618 q1( double t )
00619 {
00620 double  tsq;
00621 double  num, denom, result;
00622 
00623   tsq = t*t;
00624 
00625   num = ((((s3[5]*tsq + s3[4])*tsq + s3[3])*tsq +
00626     s3[2])*tsq + s3[1])*tsq + s3[0];
00627 
00628   denom = (((((t3[6]*tsq + t3[5])*tsq + t3[4])*tsq +
00629     t3[3])*tsq + t3[2])*tsq + t3[1])*tsq + t3[0];
00630 
00631   result = t*num/denom;
00632 
00633   return ( result );
00634 }
00635 
00636 #ifdef NO_LONG_DOUBLE
00637 
00638 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00639 extern long double __j1l(long double);
00640 extern long double __y1l(long double);
00641 #pragma weak j1l
00642 #pragma weak y1l
00643 long double j1l( long double x ) {
00644   return ( (long double)__j1((double)x) );
00645 }
00646 long double y1l( long double x ) {
00647   return ( (long double)__y1((double)x) );
00648 }
00649 #elif defined(__GNUC__)
00650 extern  long double  __j1l(long double);
00651 
00652 long double    j1l() __attribute__ ((weak, alias ("__j1l")));
00653 
00654 extern  long double  __y1l(long double);
00655 
00656 long double    y1l() __attribute__ ((weak, alias ("__y1l")));
00657 
00658 #endif
00659 
00660 long double
00661 __j1l( long double x )
00662 { 
00663   return ( (long double)__j1((double)x) );
00664 }
00665 
00666 long double
00667 __y1l( long double x )
00668 { 
00669   return ( (long double)__y1((double)x) );
00670 }
00671 
00672 #endif
00673