Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /*                                                      pow.c
  2.  *
  3.  *      Power function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, z, pow();
  10.  *
  11.  * z = pow( x, y );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Computes x raised to the yth power.  Analytically,
  18.  *
  19.  *      x**y  =  exp( y log(x) ).
  20.  *
  21.  * Following Cody and Waite, this program uses a lookup table
  22.  * of 2**-i/16 and pseudo extended precision arithmetic to
  23.  * obtain an extra three bits of accuracy in both the logarithm
  24.  * and the exponential.
  25.  *
  26.  *
  27.  *
  28.  * ACCURACY:
  29.  *
  30.  *                      Relative error:
  31.  * arithmetic   domain     # trials      peak         rms
  32.  *    IEEE     -26,26       30000      4.2e-16      7.7e-17
  33.  *    DEC      -26,26       60000      4.8e-17      9.1e-18
  34.  * 1/26 < x < 26, with log(x) uniformly distributed.
  35.  * -26 < y < 26, y uniformly distributed.
  36.  *    IEEE     0,8700       30000      1.5e-14      2.1e-15
  37.  * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
  38.  *
  39.  *
  40.  * ERROR MESSAGES:
  41.  *
  42.  *   message         condition      value returned
  43.  * pow overflow     x**y > MAXNUM      INFINITY
  44.  * pow underflow   x**y < 1/MAXNUM       0.0
  45.  * pow domain      x<0 and y noninteger  0.0
  46.  *
  47.  */
  49. /*
  50. Cephes Math Library Release 2.8:  June, 2000
  51. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  52. */
  53.  
  54. /*
  55. Modified for mingw
  56. 2002-09-27 Danny Smith <dannysmith@users.sourceforge.net>
  57. */
  58.  
  59. #ifdef __MINGW32__
  60. #include "cephes_mconf.h"
  61. #else
  62. #include "mconf.h"
  63. static char fname[] = {"pow"};
  64. #endif
  65.  
  66. #ifndef _SET_ERRNO
  67. #define _SET_ERRNO(x)
  68. #endif
  69.  
  70. #define SQRTH 0.70710678118654752440
  71.  
  72. #ifdef UNK
  73. static double P[] = {
  74.   4.97778295871696322025E-1,
  75.   3.73336776063286838734E0,
  76.   7.69994162726912503298E0,
  77.   4.66651806774358464979E0
  78. };
  79. static double Q[] = {
  80. /* 1.00000000000000000000E0, */
  81.   9.33340916416696166113E0,
  82.   2.79999886606328401649E1,
  83.   3.35994905342304405431E1,
  84.   1.39995542032307539578E1
  85. };
  86. /* 2^(-i/16), IEEE precision */
  87. static double A[] = {
  88.   1.00000000000000000000E0,
  89.   9.57603280698573700036E-1,
  90.   9.17004043204671215328E-1,
  91.   8.78126080186649726755E-1,
  92.   8.40896415253714502036E-1,
  93.   8.05245165974627141736E-1,
  94.   7.71105412703970372057E-1,
  95.   7.38413072969749673113E-1,
  96.   7.07106781186547572737E-1,
  97.   6.77127773468446325644E-1,
  98.   6.48419777325504820276E-1,
  99.   6.20928906036742001007E-1,
  100.   5.94603557501360513449E-1,
  101.   5.69394317378345782288E-1,
  102.   5.45253866332628844837E-1,
  103.   5.22136891213706877402E-1,
  104.   5.00000000000000000000E-1
  105. };
  106. static double B[] = {
  107.  0.00000000000000000000E0,
  108.  1.64155361212281360176E-17,
  109.  4.09950501029074826006E-17,
  110.  3.97491740484881042808E-17,
  111. -4.83364665672645672553E-17,
  112.  1.26912513974441574796E-17,
  113.  1.99100761573282305549E-17,
  114. -1.52339103990623557348E-17,
  115.  0.00000000000000000000E0
  116. };
  117. static double R[] = {
  118.  1.49664108433729301083E-5,
  119.  1.54010762792771901396E-4,
  120.  1.33335476964097721140E-3,
  121.  9.61812908476554225149E-3,
  122.  5.55041086645832347466E-2,
  123.  2.40226506959099779976E-1,
  124.  6.93147180559945308821E-1
  125. };
  126.  
  127. #define douba(k) A[k]
  128. #define doubb(k) B[k]
  129. #define MEXP 16383.0
  130. #ifdef DENORMAL
  131. #define MNEXP -17183.0
  132. #else
  133. #define MNEXP -16383.0
  134. #endif
  135. #endif
  136.  
  137. #ifdef DEC
  138. static unsigned short P[] = {
  139. 0037776,0156313,0175332,0163602,
  140. 0040556,0167577,0052366,0174245,
  141. 0040766,0062753,0175707,0055564,
  142. 0040625,0052035,0131344,0155636,
  143. };
  144. static unsigned short Q[] = {
  145. /*0040200,0000000,0000000,0000000,*/
  146. 0041025,0052644,0154404,0105155,
  147. 0041337,0177772,0007016,0047646,
  148. 0041406,0062740,0154273,0020020,
  149. 0041137,0177054,0106127,0044555,
  150. };
  151. static unsigned short A[] = {
  152. 0040200,0000000,0000000,0000000,
  153. 0040165,0022575,0012444,0103314,
  154. 0040152,0140306,0163735,0022071,
  155. 0040140,0146336,0166052,0112341,
  156. 0040127,0042374,0145326,0116553,
  157. 0040116,0022214,0012437,0102201,
  158. 0040105,0063452,0010525,0003333,
  159. 0040075,0004243,0117530,0006067,
  160. 0040065,0002363,0031771,0157145,
  161. 0040055,0054076,0165102,0120513,
  162. 0040045,0177326,0124661,0050471,
  163. 0040036,0172462,0060221,0120422,
  164. 0040030,0033760,0050615,0134251,
  165. 0040021,0141723,0071653,0010703,
  166. 0040013,0112701,0161752,0105727,
  167. 0040005,0125303,0063714,0044173,
  168. 0040000,0000000,0000000,0000000
  169. };
  170. static unsigned short B[] = {
  171. 0000000,0000000,0000000,0000000,
  172. 0021473,0040265,0153315,0140671,
  173. 0121074,0062627,0042146,0176454,
  174. 0121413,0003524,0136332,0066212,
  175. 0121767,0046404,0166231,0012553,
  176. 0121257,0015024,0002357,0043574,
  177. 0021736,0106532,0043060,0056206,
  178. 0121310,0020334,0165705,0035326,
  179. 0000000,0000000,0000000,0000000
  180. };
  181.  
  182. static unsigned short R[] = {
  183. 0034173,0014076,0137624,0115771,
  184. 0035041,0076763,0003744,0111311,
  185. 0035656,0141766,0041127,0074351,
  186. 0036435,0112533,0073611,0116664,
  187. 0037143,0054106,0134040,0152223,
  188. 0037565,0176757,0176026,0025551,
  189. 0040061,0071027,0173721,0147572
  190. };
  191.  
  192. /*
  193. static double R[] = {
  194. 0.14928852680595608186e-4,
  195. 0.15400290440989764601e-3,
  196. 0.13333541313585784703e-2,
  197. 0.96181290595172416964e-2,
  198. 0.55504108664085595326e-1,
  199. 0.24022650695909537056e0,
  200. 0.69314718055994529629e0
  201. };
  202. */
  203. #define douba(k) (*(double *)&A[(k)<<2])
  204. #define doubb(k) (*(double *)&B[(k)<<2])
  205. #define MEXP 2031.0
  206. #define MNEXP -2031.0
  207. #endif
  208.  
  209. #ifdef IBMPC
  210. static const unsigned short P[] = {
  211. 0x5cf0,0x7f5b,0xdb99,0x3fdf,
  212. 0xdf15,0xea9e,0xddef,0x400d,
  213. 0xeb6f,0x7f78,0xccbd,0x401e,
  214. 0x9b74,0xb65c,0xaa83,0x4012,
  215. };
  216. static const unsigned short Q[] = {
  217. /*0x0000,0x0000,0x0000,0x3ff0,*/
  218. 0x914e,0x9b20,0xaab4,0x4022,
  219. 0xc9f5,0x41c1,0xffff,0x403b,
  220. 0x6402,0x1b17,0xccbc,0x4040,
  221. 0xe92e,0x918a,0xffc5,0x402b,
  222. };
  223. static const unsigned short A[] = {
  224. 0x0000,0x0000,0x0000,0x3ff0,
  225. 0x90da,0xa2a4,0xa4af,0x3fee,
  226. 0xa487,0xdcfb,0x5818,0x3fed,
  227. 0x529c,0xdd85,0x199b,0x3fec,
  228. 0xd3ad,0x995a,0xe89f,0x3fea,
  229. 0xf090,0x82a3,0xc491,0x3fe9,
  230. 0xa0db,0x422a,0xace5,0x3fe8,
  231. 0x0187,0x73eb,0xa114,0x3fe7,
  232. 0x3bcd,0x667f,0xa09e,0x3fe6,
  233. 0x5429,0xdd48,0xab07,0x3fe5,
  234. 0x2a27,0xd536,0xbfda,0x3fe4,
  235. 0x3422,0x4c12,0xdea6,0x3fe3,
  236. 0xb715,0x0a31,0x06fe,0x3fe3,
  237. 0x6238,0x6e75,0x387a,0x3fe2,
  238. 0x517b,0x3c7d,0x72b8,0x3fe1,
  239. 0x890f,0x6cf9,0xb558,0x3fe0,
  240. 0x0000,0x0000,0x0000,0x3fe0
  241. };
  242. static const unsigned short B[] = {
  243. 0x0000,0x0000,0x0000,0x0000,
  244. 0x3707,0xd75b,0xed02,0x3c72,
  245. 0xcc81,0x345d,0xa1cd,0x3c87,
  246. 0x4b27,0x5686,0xe9f1,0x3c86,
  247. 0x6456,0x13b2,0xdd34,0xbc8b,
  248. 0x42e2,0xafec,0x4397,0x3c6d,
  249. 0x82e4,0xd231,0xf46a,0x3c76,
  250. 0x8a76,0xb9d7,0x9041,0xbc71,
  251. 0x0000,0x0000,0x0000,0x0000
  252. };
  253. static const unsigned short R[] = {
  254. 0x937f,0xd7f2,0x6307,0x3eef,
  255. 0x9259,0x60fc,0x2fbe,0x3f24,
  256. 0xef1d,0xc84a,0xd87e,0x3f55,
  257. 0x33b7,0x6ef1,0xb2ab,0x3f83,
  258. 0x1a92,0xd704,0x6b08,0x3fac,
  259. 0xc56d,0xff82,0xbfbd,0x3fce,
  260. 0x39ef,0xfefa,0x2e42,0x3fe6
  261. };
  262.  
  263. #define douba(k) (*(double *)&A[(k)<<2])
  264. #define doubb(k) (*(double *)&B[(k)<<2])
  265. #define MEXP 16383.0
  266. #ifdef DENORMAL
  267. #define MNEXP -17183.0
  268. #else
  269. #define MNEXP -16383.0
  270. #endif
  271. #endif
  272.  
  273. #ifdef MIEEE
  274. static unsigned short P[] = {
  275. 0x3fdf,0xdb99,0x7f5b,0x5cf0,
  276. 0x400d,0xddef,0xea9e,0xdf15,
  277. 0x401e,0xccbd,0x7f78,0xeb6f,
  278. 0x4012,0xaa83,0xb65c,0x9b74
  279. };
  280. static unsigned short Q[] = {
  281. 0x4022,0xaab4,0x9b20,0x914e,
  282. 0x403b,0xffff,0x41c1,0xc9f5,
  283. 0x4040,0xccbc,0x1b17,0x6402,
  284. 0x402b,0xffc5,0x918a,0xe92e
  285. };
  286. static unsigned short A[] = {
  287. 0x3ff0,0x0000,0x0000,0x0000,
  288. 0x3fee,0xa4af,0xa2a4,0x90da,
  289. 0x3fed,0x5818,0xdcfb,0xa487,
  290. 0x3fec,0x199b,0xdd85,0x529c,
  291. 0x3fea,0xe89f,0x995a,0xd3ad,
  292. 0x3fe9,0xc491,0x82a3,0xf090,
  293. 0x3fe8,0xace5,0x422a,0xa0db,
  294. 0x3fe7,0xa114,0x73eb,0x0187,
  295. 0x3fe6,0xa09e,0x667f,0x3bcd,
  296. 0x3fe5,0xab07,0xdd48,0x5429,
  297. 0x3fe4,0xbfda,0xd536,0x2a27,
  298. 0x3fe3,0xdea6,0x4c12,0x3422,
  299. 0x3fe3,0x06fe,0x0a31,0xb715,
  300. 0x3fe2,0x387a,0x6e75,0x6238,
  301. 0x3fe1,0x72b8,0x3c7d,0x517b,
  302. 0x3fe0,0xb558,0x6cf9,0x890f,
  303. 0x3fe0,0x0000,0x0000,0x0000
  304. };
  305. static unsigned short B[] = {
  306. 0x0000,0x0000,0x0000,0x0000,
  307. 0x3c72,0xed02,0xd75b,0x3707,
  308. 0x3c87,0xa1cd,0x345d,0xcc81,
  309. 0x3c86,0xe9f1,0x5686,0x4b27,
  310. 0xbc8b,0xdd34,0x13b2,0x6456,
  311. 0x3c6d,0x4397,0xafec,0x42e2,
  312. 0x3c76,0xf46a,0xd231,0x82e4,
  313. 0xbc71,0x9041,0xb9d7,0x8a76,
  314. 0x0000,0x0000,0x0000,0x0000
  315. };
  316. static unsigned short R[] = {
  317. 0x3eef,0x6307,0xd7f2,0x937f,
  318. 0x3f24,0x2fbe,0x60fc,0x9259,
  319. 0x3f55,0xd87e,0xc84a,0xef1d,
  320. 0x3f83,0xb2ab,0x6ef1,0x33b7,
  321. 0x3fac,0x6b08,0xd704,0x1a92,
  322. 0x3fce,0xbfbd,0xff82,0xc56d,
  323. 0x3fe6,0x2e42,0xfefa,0x39ef
  324. };
  325.  
  326. #define douba(k) (*(double *)&A[(k)<<2])
  327. #define doubb(k) (*(double *)&B[(k)<<2])
  328. #define MEXP 16383.0
  329. #ifdef DENORMAL
  330. #define MNEXP -17183.0
  331. #else
  332. #define MNEXP -16383.0
  333. #endif
  334. #endif
  335.  
  336. /* log2(e) - 1 */
  337. #define LOG2EA 0.44269504088896340736
  338.  
  339. #define F W
  340. #define Fa Wa
  341. #define Fb Wb
  342. #define G W
  343. #define Ga Wa
  344. #define Gb u
  345. #define H W
  346. #define Ha Wb
  347. #define Hb Wb
  348.  
  349. #ifdef __MINGW32__
  350. static __inline__ double reduc( double );
  351. extern double __powi ( double, int );
  352. extern double pow ( double x,  double y);
  353.  
  354. #else /* __MINGW32__ */
  355.  
  356. #ifdef ANSIPROT
  357. extern double floor ( double );
  358. extern double fabs ( double );
  359. extern double frexp ( double, int * );
  360. extern double ldexp ( double, int );
  361. extern double polevl ( double, void *, int );
  362. extern double p1evl ( double, void *, int );
  363. extern double __powi ( double, int );
  364. extern int signbit ( double );
  365. extern int isnan ( double );
  366. extern int isfinite ( double );
  367. static double reduc ( double );
  368. #else
  369. double floor(), fabs(), frexp(), ldexp();
  370. double polevl(), p1evl(), __powi();
  371. int signbit(), isnan(), isfinite();
  372. static double reduc();
  373. #endif
  374. extern double MAXNUM;
  375. #ifdef INFINITIES
  376. extern double INFINITY;
  377. #endif
  378. #ifdef NANS
  379. extern double NAN;
  380. #endif
  381. #ifdef MINUSZERO
  382. extern double NEGZERO;
  383. #endif
  384.  
  385. #endif /* __MINGW32__ */
  386.  
  387. double pow( x, y )
  388. double x, y;
  389. {
  390. double w, z, W, Wa, Wb, ya, yb, u;
  391. /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
  392. double aw, ay, wy;
  393. int e, i, nflg, iyflg, yoddint;
  394.  
  395. if( y == 0.0 )
  396.         return( 1.0 );
  397. #ifdef NANS
  398. if( isnan(x) || isnan(y) )
  399.         {
  400.         _SET_ERRNO (EDOM);
  401.         return( x + y );
  402.         }
  403. #endif
  404. if( y == 1.0 )
  405.         return( x );
  406.  
  407.  
  408. #ifdef INFINITIES
  409. if( !isfinite(y) && (x == 1.0 || x == -1.0) )
  410.         {
  411.         mtherr( "pow", DOMAIN );
  412. #ifdef NANS
  413.         return( NAN );
  414. #else
  415.         return( INFINITY );
  416. #endif
  417.         }
  418. #endif
  419.  
  420. if( x == 1.0 )
  421.         return( 1.0 );
  422.  
  423. if( y >= MAXNUM )
  424.         {
  425.         _SET_ERRNO (ERANGE);
  426. #ifdef INFINITIES
  427.         if( x > 1.0 )
  428.                 return( INFINITY );
  429. #else
  430.         if( x > 1.0 )
  431.                 return( MAXNUM );
  432. #endif
  433.         if( x > 0.0 && x < 1.0 )
  434.                 return( 0.0);
  435.         if( x < -1.0 )
  436.                 {
  437. #ifdef INFINITIES
  438.                 return( INFINITY );
  439. #else
  440.                 return( MAXNUM );
  441. #endif
  442.                 }
  443.         if( x > -1.0 && x < 0.0 )
  444.                 return( 0.0 );
  445.         }
  446. if( y <= -MAXNUM )
  447.         {
  448.         _SET_ERRNO (ERANGE);
  449.         if( x > 1.0 )
  450.                 return( 0.0 );
  451. #ifdef INFINITIES
  452.         if( x > 0.0 && x < 1.0 )
  453.                 return( INFINITY );
  454. #else
  455.         if( x > 0.0 && x < 1.0 )
  456.                 return( MAXNUM );
  457. #endif
  458.         if( x < -1.0 )
  459.                 return( 0.0 );
  460. #ifdef INFINITIES
  461.         if( x > -1.0 && x < 0.0 )
  462.                 return( INFINITY );
  463. #else
  464.         if( x > -1.0 && x < 0.0 )
  465.                 return( MAXNUM );
  466. #endif
  467.         }
  468. if( x >= MAXNUM )
  469.         {
  470. #if INFINITIES
  471.         if( y > 0.0 )
  472.                 return( INFINITY );
  473. #else
  474.         if( y > 0.0 )
  475.                 return( MAXNUM );
  476. #endif
  477.         return(0.0);
  478.         }
  479. /* Set iyflg to 1 if y is an integer.  */
  480. iyflg = 0;
  481. w = floor(y);
  482. if( w == y )
  483.         iyflg = 1;
  484.  
  485. /* Test for odd integer y.  */
  486. yoddint = 0;
  487. if( iyflg )
  488.         {
  489.         ya = fabs(y);
  490.         ya = floor(0.5 * ya);
  491.         yb = 0.5 * fabs(w);
  492.         if( ya != yb )
  493.                 yoddint = 1;
  494.         }
  495.  
  496. if( x <= -MAXNUM )
  497.         {
  498.         if( y > 0.0 )
  499.                 {
  500. #ifdef INFINITIES
  501.                 if( yoddint )
  502.                         return( -INFINITY );
  503.                 return( INFINITY );
  504. #else
  505.                 if( yoddint )
  506.                         return( -MAXNUM );
  507.                 return( MAXNUM );
  508. #endif
  509.                 }
  510.         if( y < 0.0 )
  511.                 {
  512. #ifdef MINUSZERO
  513.                 if( yoddint )
  514.                         return( NEGZERO );
  515. #endif
  516.                 return( 0.0 );
  517.                 }
  518.         }
  519.  
  520. nflg = 0;       /* flag = 1 if x<0 raised to integer power */
  521. if( x <= 0.0 )
  522.         {
  523.         if( x == 0.0 )
  524.                 {
  525.                 if( y < 0.0 )
  526.                         {
  527. #ifdef MINUSZERO
  528.                         if( signbit(x) && yoddint )
  529.                                 return( -INFINITY );
  530. #endif
  531. #ifdef INFINITIES
  532.                         return( INFINITY );
  533. #else
  534.                         return( MAXNUM );
  535. #endif
  536.                         }
  537.                 if( y > 0.0 )
  538.                         {
  539. #ifdef MINUSZERO
  540.                         if( signbit(x) && yoddint )
  541.                                 return( NEGZERO );
  542. #endif
  543.                         return( 0.0 );
  544.                         }
  545.                 return( 1.0 );
  546.                 }
  547.         else
  548.                 {
  549.                 if( iyflg == 0 )
  550.                         { /* noninteger power of negative number */
  551.                         mtherr( fname, DOMAIN );
  552.                         _SET_ERRNO (EDOM);
  553. #ifdef NANS
  554.                         return(NAN);
  555. #else
  556.                         return(0.0L);
  557. #endif
  558.                         }
  559.                 nflg = 1;
  560.                 }
  561.         }
  562.  
  563. /* Integer power of an integer.  */
  564.  
  565. if( iyflg )
  566.         {
  567.         i = w;
  568.         w = floor(x);
  569.         if( (w == x) && (fabs(y) < 32768.0) )
  570.                 {
  571.                 w = __powi( x, (int) y );
  572.                 return( w );
  573.                 }
  574.         }
  575.  
  576. if( nflg )
  577.         x = fabs(x);
  578.  
  579. /* For results close to 1, use a series expansion.  */
  580. w = x - 1.0;
  581. aw = fabs(w);
  582. ay = fabs(y);
  583. wy = w * y;
  584. ya = fabs(wy);
  585. if((aw <= 1.0e-3 && ay <= 1.0)
  586.    || (ya <= 1.0e-3 && ay >= 1.0))
  587.         {
  588.         z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
  589.                 + 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
  590.         goto done;
  591.         }
  592. /* These are probably too much trouble.  */
  593. #if 0
  594. w = y * log(x);
  595. if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
  596.   {
  597.     z = ((((((
  598.     w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
  599.     goto done;
  600.   }
  601.  
  602. if(ya <= 1.0e-3 && aw <= 1.0e-4)
  603.   {
  604.     z = (((((
  605.              wy*1./720.
  606.              + (-w*1./48. + 1./120.) )*wy
  607.             + ((w*17./144. - 1./12.)*w + 1./24.) )*wy
  608.            + (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
  609.           + ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
  610.          + (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
  611.            + wy + 1.0;
  612.     goto done;
  613.   }
  614. #endif
  615.  
  616. /* separate significand from exponent */
  617. x = frexp( x, &e );
  618.  
  619. #if 0
  620. /* For debugging, check for gross overflow. */
  621. if( (e * y)  > (MEXP + 1024) )
  622.         goto overflow;
  623. #endif
  624.  
  625. /* Find significand of x in antilog table A[]. */
  626. i = 1;
  627. if( x <= douba(9) )
  628.         i = 9;
  629. if( x <= douba(i+4) )
  630.         i += 4;
  631. if( x <= douba(i+2) )
  632.         i += 2;
  633. if( x >= douba(1) )
  634.         i = -1;
  635. i += 1;
  636.  
  637.  
  638. /* Find (x - A[i])/A[i]
  639.  * in order to compute log(x/A[i]):
  640.  *
  641.  * log(x) = log( a x/a ) = log(a) + log(x/a)
  642.  *
  643.  * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
  644.  */
  645. x -= douba(i);
  646. x -= doubb(i/2);
  647. x /= douba(i);
  648.  
  649.  
  650. /* rational approximation for log(1+v):
  651.  *
  652.  * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
  653.  */
  654. z = x*x;
  655. w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
  656. w = w - ldexp( z, -1 );   /*  w - 0.5 * z  */
  657.  
  658. /* Convert to base 2 logarithm:
  659.  * multiply by log2(e)
  660.  */
  661. w = w + LOG2EA * w;
  662. /* Note x was not yet added in
  663.  * to above rational approximation,
  664.  * so do it now, while multiplying
  665.  * by log2(e).
  666.  */
  667. z = w + LOG2EA * x;
  668. z = z + x;
  669.  
  670. /* Compute exponent term of the base 2 logarithm. */
  671. w = -i;
  672. w = ldexp( w, -4 );     /* divide by 16 */
  673. w += e;
  674. /* Now base 2 log of x is w + z. */
  675.  
  676. /* Multiply base 2 log by y, in extended precision. */
  677.  
  678. /* separate y into large part ya
  679.  * and small part yb less than 1/16
  680.  */
  681. ya = reduc(y);
  682. yb = y - ya;
  683.  
  684.  
  685. F = z * y  +  w * yb;
  686. Fa = reduc(F);
  687. Fb = F - Fa;
  688.  
  689. G = Fa + w * ya;
  690. Ga = reduc(G);
  691. Gb = G - Ga;
  692.  
  693. H = Fb + Gb;
  694. Ha = reduc(H);
  695. w = ldexp( Ga+Ha, 4 );
  696.  
  697. /* Test the power of 2 for overflow */
  698. if( w > MEXP )
  699.         {
  700. #ifndef INFINITIES
  701.         mtherr( fname, OVERFLOW );
  702. #endif
  703. #ifdef INFINITIES
  704.         if( nflg && yoddint )
  705.           return( -INFINITY );
  706.         return( INFINITY );
  707. #else
  708.         if( nflg && yoddint )
  709.           return( -MAXNUM );
  710.         return( MAXNUM );
  711. #endif
  712.         }
  713.  
  714. if( w < (MNEXP - 1) )
  715.         {
  716. #ifndef DENORMAL
  717.         mtherr( fname, UNDERFLOW );
  718. #endif
  719. #ifdef MINUSZERO
  720.         if( nflg && yoddint )
  721.           return( NEGZERO );
  722. #endif
  723.         return( 0.0 );
  724.         }
  725.  
  726. e = w;
  727. Hb = H - Ha;
  728.  
  729. if( Hb > 0.0 )
  730.         {
  731.         e += 1;
  732.         Hb -= 0.0625;
  733.         }
  734.  
  735. /* Now the product y * log2(x)  =  Hb + e/16.0.
  736.  *
  737.  * Compute base 2 exponential of Hb,
  738.  * where -0.0625 <= Hb <= 0.
  739.  */
  740. z = Hb * polevl( Hb, R, 6 );  /*    z  =  2**Hb - 1    */
  741.  
  742. /* Express e/16 as an integer plus a negative number of 16ths.
  743.  * Find lookup table entry for the fractional power of 2.
  744.  */
  745. if( e < 0 )
  746.         i = 0;
  747. else
  748.         i = 1;
  749. i = e/16 + i;
  750. e = 16*i - e;
  751. w = douba( e );
  752. z = w + w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
  753. z = ldexp( z, i );  /* multiply by integer power of 2 */
  754.  
  755. done:
  756.  
  757. /* Negate if odd integer power of negative number */
  758. if( nflg && yoddint )
  759.         {
  760. #ifdef MINUSZERO
  761.         if( z == 0.0 )
  762.                 z = NEGZERO;
  763.         else
  764. #endif
  765.                 z = -z;
  766.         }
  767. return( z );
  768. }
  769.  
  770.  
  771. /* Find a multiple of 1/16 that is within 1/16 of x. */
  772. static __inline__ double reduc(x)
  773. double x;
  774. {
  775. double t;
  776.  
  777. t = ldexp( x, 4 );
  778. t = floor( t );
  779. t = ldexp( t, -4 );
  780. return(t);
  781. }
  782.