1b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper 2903924f0541ae61edddee994437549f1db343047Mark Wielaard/* @(#)e_j1.c 1.3 95/01/18 */ 3de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard/* 4b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * ==================================================== 5de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * 7b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * Developed at SunSoft, a Sun Microsystems, Inc. business. 8de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * Permission to use, copy, modify, and distribute this 9de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * software is freely granted, provided that this notice 10de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * is preserved. 11de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * ==================================================== 12de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard */ 13de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard 14de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard/* __ieee754_j1(x), __ieee754_y1(x) 15de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * Bessel function of the first and second kinds of order zero. 16de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * Method -- ieee_j1(x): 17de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * 1. For tiny x, we use ieee_j1(x) = x/2 - x^3/16 + x^5/384 - ... 18de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * 2. Reduce x to |x| since ieee_j1(x)=-ieee_j1(-x), and 19de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * for x in (0,2) 20de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * j1(x) = x/2 + x*z*R0/S0, where z = x*x; 21361df7da6dfecd817b27e62b91752ac316d7cdd4Ulrich Drepper * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) 22361df7da6dfecd817b27e62b91752ac316d7cdd4Ulrich Drepper * for x in (2,inf) 23361df7da6dfecd817b27e62b91752ac316d7cdd4Ulrich Drepper * j1(x) = ieee_sqrt(2/(pi*x))*(p1(x)*ieee_cos(x1)-q1(x)*ieee_sin(x1)) 24361df7da6dfecd817b27e62b91752ac316d7cdd4Ulrich Drepper * y1(x) = ieee_sqrt(2/(pi*x))*(p1(x)*ieee_sin(x1)+q1(x)*ieee_cos(x1)) 25de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * where x1 = x-3*pi/4. It is better to compute ieee_sin(x1),cos(x1) 26de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * as follow: 27de2ed97f33139af5c7a0811e4ec66fc896a13cf2Mark Wielaard * cos(x1) = ieee_cos(x)cos(3pi/4)+ieee_sin(x)sin(3pi/4) 28b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * = 1/ieee_sqrt(2) * (ieee_sin(x) - ieee_cos(x)) 29b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * sin(x1) = ieee_sin(x)cos(3pi/4)-ieee_cos(x)sin(3pi/4) 30920f03dcd3a6dd46b3d88c45be8dbfa0adc780ddFlorian Weimer * = -1/ieee_sqrt(2) * (ieee_sin(x) + ieee_cos(x)) 31b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * (To avoid cancellation, use 32b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * sin(x) +- ieee_cos(x) = -ieee_cos(2x)/(ieee_sin(x) -+ ieee_cos(x)) 33d17fac7e89666b47811581b10b5ca0d253a3a82dRoland McGrath * to compute the worse one.) 34b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * 355083a70d3b64946fa47ea5766943a15a3ecc6891Mark Wielaard * 3 Special cases 36b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * j1(nan)= nan 37b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * j1(0) = 0 38df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard * j1(inf) = 0 39b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * 40df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard * Method -- ieee_y1(x): 41b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * 1. screen out x<=0 cases: ieee_y1(0)=-inf, ieee_y1(x<0)=NaN 42e07da4f30d2fabc6e2d4a4bfcc05e02091c99f0aRoland McGrath * 2. For x<2. 43e07da4f30d2fabc6e2d4a4bfcc05e02091c99f0aRoland McGrath * Since 44e07da4f30d2fabc6e2d4a4bfcc05e02091c99f0aRoland McGrath * y1(x) = 2/pi*(ieee_j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) 45e07da4f30d2fabc6e2d4a4bfcc05e02091c99f0aRoland McGrath * therefore ieee_y1(x)-2/pi*ieee_j1(x)*ln(x)-1/x is an odd function. 46059c83e5db89955913a39fe6705acca571c32c3fRoland McGrath * We use the following function to approximate y1, 47059c83e5db89955913a39fe6705acca571c32c3fRoland McGrath * y1(x) = x*U(z)/V(z) + (2/pi)*(ieee_j1(x)*ln(x)-1/x), z= x^2 48df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard * where for x in [0,2] (abs err less than 2**-65.89) 493425454a10d307fae891fb667cf7969e945cde79Josh Stone * U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4 50059c83e5db89955913a39fe6705acca571c32c3fRoland McGrath * V(z) = 1 + v0[0]*z + ... + v0[4]*z^5 51df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard * Note: For tiny x, 1/x dominate y1 and hence 52b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) 53b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * 3. For x>=2. 54df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard * y1(x) = ieee_sqrt(2/(pi*x))*(p1(x)*ieee_sin(x1)+q1(x)*ieee_cos(x1)) 55b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * where x1 = x-3*pi/4. It is better to compute ieee_sin(x1),cos(x1) 56df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard * by method mentioned above. 57122f388b300a489b54829d80f5e5fffb3f743ab7Roland McGrath */ 58df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard 59df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard#include "fdlibm.h" 60df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard 61df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard#ifdef __STDC__ 62122f388b300a489b54829d80f5e5fffb3f743ab7Roland McGrathstatic double pone(double), qone(double); 63122f388b300a489b54829d80f5e5fffb3f743ab7Roland McGrath#else 64122f388b300a489b54829d80f5e5fffb3f743ab7Roland McGrathstatic double pone(), qone(); 65df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard#endif 66df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard 67df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard#ifdef __STDC__ 68df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaardstatic const double 69df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard#else 70df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaardstatic double 71df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard#endif 72df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaardhuge = 1e300, 73df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaardone = 1.0, 74df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaardinvsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ 75df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaardtpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 76df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaard /* R0/S0 on [0,2] */ 77df85bf99021119fcbb2ced66dd69f1cceafb180cMark Wielaardr00 = -6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */ 78b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepperr01 = 1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */ 79b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepperr02 = -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */ 806258e7486eb3eed6e50005946795c5fbf73aa106Ulrich Drepperr03 = 4.96727999609584448412e-08, /* 0x3E6AAAFA, 0x46CA0BD9 */ 816258e7486eb3eed6e50005946795c5fbf73aa106Ulrich Dreppers01 = 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */ 82ed431ddb74331f24add8c6d932ebed129c4385d8Roland McGraths02 = 1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */ 83ed431ddb74331f24add8c6d932ebed129c4385d8Roland McGraths03 = 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */ 846258e7486eb3eed6e50005946795c5fbf73aa106Ulrich Dreppers04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */ 856258e7486eb3eed6e50005946795c5fbf73aa106Ulrich Dreppers05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */ 86122f388b300a489b54829d80f5e5fffb3f743ab7Roland McGrath 876258e7486eb3eed6e50005946795c5fbf73aa106Ulrich Drepperstatic double zero = 0.0; 88b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper 8965cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone#ifdef __STDC__ 90f95760aff004850544f83626404c134d6a07c630Roland McGrath double __ieee754_j1(double x) 912c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath#else 922c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath double __ieee754_j1(x) 932c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath double x; 942c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath#endif 952c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath{ 962c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath double z, s,c,ss,cc,r,u,v,y; 972c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath int hx,ix; 982c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath 992c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath hx = __HI(x); 1002c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath ix = hx&0x7fffffff; 1012c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath if(ix>=0x7ff00000) return one/x; 1022c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath y = ieee_fabs(x); 1032c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath if(ix >= 0x40000000) { /* |x| >= 2.0 */ 1042c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath s = ieee_sin(y); 1052c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath c = ieee_cos(y); 1062c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath ss = -s-c; 1072c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath cc = s-c; 1082c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath if(ix<0x7fe00000) { /* make sure y+y not overflow */ 1092c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath z = ieee_cos(y+y); 1102c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath if ((s*c)>zero) cc = z/ss; 1112c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath else ss = z/cc; 1122c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath } 1132c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath /* 1142c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath * ieee_j1(x) = 1/ieee_sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / ieee_sqrt(x) 1155ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * ieee_y1(x) = 1/ieee_sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / ieee_sqrt(x) 1165ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath */ 1175ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath if(ix>0x48000000) z = (invsqrtpi*cc)/ieee_sqrt(y); 1185ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath else { 1195ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath u = pone(y); v = qone(y); 1205ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath z = invsqrtpi*(u*cc-v*ss)/ieee_sqrt(y); 1215ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath } 1222c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrath if(hx<0) return -z; 123f95760aff004850544f83626404c134d6a07c630Roland McGrath else return z; 124f95760aff004850544f83626404c134d6a07c630Roland McGrath } 125f95760aff004850544f83626404c134d6a07c630Roland McGrath if(ix<0x3e400000) { /* |x|<2**-27 */ 126f95760aff004850544f83626404c134d6a07c630Roland McGrath if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */ 1271743d7f010bead5e869d097e23ce840583913381Roland McGrath } 128f95760aff004850544f83626404c134d6a07c630Roland McGrath z = x*x; 129f95760aff004850544f83626404c134d6a07c630Roland McGrath r = z*(r00+z*(r01+z*(r02+z*r03))); 130f95760aff004850544f83626404c134d6a07c630Roland McGrath s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); 131f95760aff004850544f83626404c134d6a07c630Roland McGrath r *= x; 1321743d7f010bead5e869d097e23ce840583913381Roland McGrath return(x*0.5+r/s); 133f95760aff004850544f83626404c134d6a07c630Roland McGrath} 1345ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 135f95760aff004850544f83626404c134d6a07c630Roland McGrath#ifdef __STDC__ 1365ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic const double U0[5] = { 1375ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 1382c7d0ddfb8a08b8c450d0d3a5e297ec5d624661dRoland McGrathstatic double U0[5] = { 139f95760aff004850544f83626404c134d6a07c630Roland McGrath#endif 140f95760aff004850544f83626404c134d6a07c630Roland McGrath -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */ 141f95760aff004850544f83626404c134d6a07c630Roland McGrath 5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */ 142b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper -1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */ 14365cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone 2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */ 14465cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone -9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */ 14565cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone}; 14665cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone#ifdef __STDC__ 14765cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stonestatic const double V0[5] = { 14865cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone#else 14965cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stonestatic double V0[5] = { 15065cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone#endif 15165cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone 1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */ 15265cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone 2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */ 15365cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone 1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */ 1546258e7486eb3eed6e50005946795c5fbf73aa106Ulrich Drepper 6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */ 15565cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone 1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */ 15665cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone}; 15765cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone 15865cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone#ifdef __STDC__ 15965cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone double __ieee754_y1(double x) 16065cefbd0793c0f9e90a326d7bebf0a47c93294adJosh Stone#else 1616258e7486eb3eed6e50005946795c5fbf73aa106Ulrich Drepper double __ieee754_y1(x) 162b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper double x; 163b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper#endif 164b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper{ 1653fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil double z, s,c,ss,cc,u,v; 1663fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil int hx,ix,lx; 1673fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil 1683fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil hx = __HI(x); 1693fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil ix = 0x7fffffff&hx; 1703fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil lx = __LO(x); 1713fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ 1723fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil if(ix>=0x7ff00000) return one/(x+x*x); 1733fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil if((ix|lx)==0) return -one/zero; 1743fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil if(hx<0) return zero/zero; 1753fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil if(ix >= 0x40000000) { /* |x| >= 2.0 */ 1763fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil s = ieee_sin(x); 1773fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil c = ieee_cos(x); 1783fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil ss = -s-c; 1793fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil cc = s-c; 1803fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil if(ix<0x7fe00000) { /* make sure x+x not overflow */ 1813fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil z = ieee_cos(x+x); 1823fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil if ((s*c)>zero) cc = z/ss; 1833fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil else ss = z/cc; 1843fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil } 1853fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil /* ieee_y1(x) = ieee_sqrt(2/(pi*x))*(p1(x)*ieee_sin(x0)+q1(x)*ieee_cos(x0)) 1863fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * where x0 = x-3pi/4 1873fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * Better formula: 1883fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * ieee_cos(x0) = ieee_cos(x)cos(3pi/4)+ieee_sin(x)sin(3pi/4) 1893fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * = 1/ieee_sqrt(2) * (ieee_sin(x) - ieee_cos(x)) 1903fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * ieee_sin(x0) = ieee_sin(x)cos(3pi/4)-ieee_cos(x)sin(3pi/4) 1913fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * = -1/ieee_sqrt(2) * (ieee_cos(x) + ieee_sin(x)) 1923fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * To avoid cancellation, use 1933fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * ieee_sin(x) +- ieee_cos(x) = -ieee_cos(2x)/(ieee_sin(x) -+ ieee_cos(x)) 1943fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil * to compute the worse one. 1953fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil */ 1963fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil if(ix>0x48000000) z = (invsqrtpi*ss)/ieee_sqrt(x); 1973fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil else { 1983fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil u = pone(x); v = qone(x); 1993fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil z = invsqrtpi*(u*ss+v*cc)/ieee_sqrt(x); 2003fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil } 2013fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil return z; 202b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper } 203b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper if(ix<=0x3c900000) { /* x < 2**-54 */ 204a605a3c8f354895f4cca56aca8f8aa468b506fdfRoland McGrath return(-tpi/x); 2051036157150672569a79d30c9b59cc7623e3e8e60Roland McGrath } 206a605a3c8f354895f4cca56aca8f8aa468b506fdfRoland McGrath z = x*x; 207b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); 208b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); 209b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x)); 210b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper} 211b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper 212b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper/* For x >= 8, the asymptotic expansions of pone is 213b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. 214b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * We approximate pone by 215ed431ddb74331f24add8c6d932ebed129c4385d8Roland McGrath * pone(x) = 1 + (R/S) 216b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 217ed431ddb74331f24add8c6d932ebed129c4385d8Roland McGrath * S = 1 + ps0*s^2 + ... + ps4*s^10 218ed431ddb74331f24add8c6d932ebed129c4385d8Roland McGrath * and 21959ea7f33f781e6e3f8c9d81d457e5d99eee8f1ceRoland McGrath * | pone(x)-1-R/S | <= 2 ** ( -60.06) 220ed431ddb74331f24add8c6d932ebed129c4385d8Roland McGrath */ 22159ea7f33f781e6e3f8c9d81d457e5d99eee8f1ceRoland McGrath 22259ea7f33f781e6e3f8c9d81d457e5d99eee8f1ceRoland McGrath#ifdef __STDC__ 22359ea7f33f781e6e3f8c9d81d457e5d99eee8f1ceRoland McGrathstatic const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 224bfc01703fe30815b450526ebcfbb37ec9f5fe9c7Roland McGrath#else 225bfc01703fe30815b450526ebcfbb37ec9f5fe9c7Roland McGrathstatic double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 22659ea7f33f781e6e3f8c9d81d457e5d99eee8f1ceRoland McGrath#endif 22759ea7f33f781e6e3f8c9d81d457e5d99eee8f1ceRoland McGrath 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 228ed431ddb74331f24add8c6d932ebed129c4385d8Roland McGrath 1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */ 2293fdcdcbcb424c45645aff1200adfb34450caf0b2Jan Kratochvil 1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */ 2301743d7f010bead5e869d097e23ce840583913381Roland McGrath 4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */ 2311743d7f010bead5e869d097e23ce840583913381Roland McGrath 3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */ 232b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper 7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */ 233b08d5a8fb42f4586d756068065186b5af7e48daUlrich Drepper}; 2344e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh#ifdef __STDC__ 2354e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsiehstatic const double ps8[5] = { 2364e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh#else 2374e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsiehstatic double ps8[5] = { 2384e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh#endif 2394e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh 1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */ 2404e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh 3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */ 2414e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh 3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */ 2424e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh 9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */ 2434e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh 3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */ 2444e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh}; 2454e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh 2464e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh#ifdef __STDC__ 2474e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsiehstatic const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 2484e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh#else 2494e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsiehstatic double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 2504e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh#endif 2514e249605bb605b9bcc7e86e48d91247abd82ff5eChih-Hung Hsieh 1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */ 2525ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */ 2535ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */ 2545ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */ 2555ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */ 2565ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */ 2575ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 2585ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 2595ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic const double ps5[5] = { 2605ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 2615ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double ps5[5] = { 2625ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 2635ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */ 2645ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */ 2655ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */ 2665ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */ 2675ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */ 2685ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 2695083a70d3b64946fa47ea5766943a15a3ecc6891Mark Wielaard 2705ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 2715ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic const double pr3[6] = { 2725ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 2735ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 2745ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 2755ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */ 2765ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */ 2775ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */ 2785ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */ 2795ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */ 2805ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */ 2815ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 2825ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 2835ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic const double ps3[5] = { 2845ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 2855ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double ps3[5] = { 2865ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 2875ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */ 2885ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */ 2895ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */ 2905ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */ 2915ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */ 2925ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 2935ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 2945ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 2955ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 2965ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 2975ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 2985ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 2995ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */ 3005ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */ 3015ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */ 3025ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */ 3035ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */ 3045ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */ 3055ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 3065ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 3075ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic const double ps2[5] = { 3085ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 3095ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double ps2[5] = { 3105ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 3115ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */ 3125ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */ 3135ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */ 3145ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */ 3155ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */ 3165ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 3175ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3185ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 3195ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath static double pone(double x) 3205ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 3215ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath static double pone(x) 3225ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath double x; 3235ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 3245ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath{ 3255ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 3265ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath const double *p,*q; 3275ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 3285ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath double *p,*q; 3295ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 3305ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath double z,r,s; 3315ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath int ix; 3325ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath ix = 0x7fffffff&__HI(x); 3335ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath if(ix>=0x40200000) {p = pr8; q= ps8;} 3345ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath else if(ix>=0x40122E8B){p = pr5; q= ps5;} 3355ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath else if(ix>=0x4006DB6D){p = pr3; q= ps3;} 3365ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath else if(ix>=0x40000000){p = pr2; q= ps2;} 3375ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath z = one/(x*x); 3385ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 3395ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 3405ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath return one+ r/s; 3415ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath} 3425ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3435ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3445ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath/* For x >= 8, the asymptotic expansions of qone is 3455ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. 3465ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * We approximate pone by 3475ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * qone(x) = s*(0.375 + (R/S)) 3485ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 3495ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * S = 1 + qs1*s^2 + ... + qs6*s^12 3505ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * and 3515ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) 3525ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath */ 3535ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3545ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 3555ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 3565ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 3575ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 3585ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 3595ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 3605ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */ 3615ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */ 3625ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */ 3635ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */ 3645ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */ 3650ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath}; 3660ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath#ifdef __STDC__ 3670ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrathstatic const double qs8[6] = { 3680ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath#else 3690ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrathstatic double qs8[6] = { 3700ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath#endif 3710ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath 1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */ 3720ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath 7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */ 3730ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath 1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */ 3740ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath 7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */ 3750ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath 6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */ 3760ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath -2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */ 3770ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath}; 3780ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath 3790ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath#ifdef __STDC__ 3800ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrathstatic const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 3810ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath#else 3820ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrathstatic double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 3835ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 3845ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */ 3855ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */ 3865ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */ 3870ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath -1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */ 38886ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaard -1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */ 3895ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */ 39086ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaard}; 39186ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaard#ifdef __STDC__ 3927eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsiehstatic const double qs5[6] = { 39386ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaard#else 3947eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsiehstatic double qs5[6] = { 395903924f0541ae61edddee994437549f1db343047Mark Wielaard#endif 396903924f0541ae61edddee994437549f1db343047Mark Wielaard 8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */ 397903924f0541ae61edddee994437549f1db343047Mark Wielaard 1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */ 3987eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsieh 1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */ 3995ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */ 4005ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */ 401903924f0541ae61edddee994437549f1db343047Mark Wielaard -4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */ 402903924f0541ae61edddee994437549f1db343047Mark Wielaard}; 403903924f0541ae61edddee994437549f1db343047Mark Wielaard 404903924f0541ae61edddee994437549f1db343047Mark Wielaard#ifdef __STDC__ 40586ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaardstatic const double qr3[6] = { 4065ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 40786ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaardstatic double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 4085ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 4097eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsieh -5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */ 4105ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */ 4117eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsieh -4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */ 4125ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */ 4135ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */ 4145ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */ 4155ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 4165ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 41786ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaardstatic const double qs3[6] = { 4185ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 4197eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsiehstatic double qs3[6] = { 4205ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 4217eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsieh 4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */ 4225ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */ 4235ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */ 4245ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */ 425903924f0541ae61edddee994437549f1db343047Mark Wielaard 1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */ 4265ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */ 4275ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 4280ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath 4290ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrath#ifdef __STDC__ 4300ef402035d71bf68c818e7b070d4c6d9744f951dRoland McGrathstatic const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 4315ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#else 4325ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 4335ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 4345ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */ 43586ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaard -1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */ 43686ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaard -2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */ 4377eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsieh -1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */ 43886ed7f7f53179d7a893329e6b9851dbb75aba405Mark Wielaard -4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */ 4397eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsieh -2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */ 440903924f0541ae61edddee994437549f1db343047Mark Wielaard}; 441903924f0541ae61edddee994437549f1db343047Mark Wielaard#ifdef __STDC__ 442903924f0541ae61edddee994437549f1db343047Mark Wielaardstatic const double qs2[6] = { 4437eff36d5daa6ebca5e6399638a7643af105ae5b0Chih-Hung Hsieh#else 4445ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrathstatic double qs2[6] = { 4455ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#endif 446903924f0541ae61edddee994437549f1db343047Mark Wielaard 2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */ 447903924f0541ae61edddee994437549f1db343047Mark Wielaard 2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */ 448903924f0541ae61edddee994437549f1db343047Mark Wielaard 7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */ 449903924f0541ae61edddee994437549f1db343047Mark Wielaard 7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */ 4505ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */ 4515ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */ 4525ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath}; 4535ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath 4545ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath#ifdef __STDC__ 4555ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath static double qone(double x) 4563a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath#else 4573a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath static double qone(x) 4583a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath double x; 4593a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath#endif 4603a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath{ 4613a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath#ifdef __STDC__ 4623a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath const double *p,*q; 4633a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath#else 4643a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath double *p,*q; 4653a44c9a73887acf982ea2752fa16ce66e42105b7Roland McGrath#endif 4665ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath double s,r,z; 4675ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath int ix; 4685ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath ix = 0x7fffffff&__HI(x); 4695ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath if(ix>=0x40200000) {p = qr8; q= qs8;} 4705ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath else if(ix>=0x40122E8B){p = qr5; q= qs5;} 4715ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath else if(ix>=0x4006DB6D){p = qr3; q= qs3;} 4725ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath else if(ix>=0x40000000){p = qr2; q= qs2;} 4735ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath z = one/(x*x); 4745ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 4755ad466d544bdca479626d07d65288b2cc2638b4aRoland McGrath s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 476903924f0541ae61edddee994437549f1db343047Mark Wielaard return (.375 + r/s)/x; 477903924f0541ae61edddee994437549f1db343047Mark Wielaard} 478903924f0541ae61edddee994437549f1db343047Mark Wielaard