Subversion Repositories Kolibri OS

Rev

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

  1. /*                                                      powl.c
  2.  *
  3.  *      Power function, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, z, powl();
  10.  *
  11.  * z = powl( 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/32 and pseudo extended precision arithmetic to
  23.  * obtain several extra bits of accuracy in both the logarithm
  24.  * and the exponential.
  25.  *
  26.  *
  27.  *
  28.  * ACCURACY:
  29.  *
  30.  * The relative error of pow(x,y) can be estimated
  31.  * by   y dl ln(2),   where dl is the absolute error of
  32.  * the internally computed base 2 logarithm.  At the ends
  33.  * of the approximation interval the logarithm equal 1/32
  34.  * and its relative error is about 1 lsb = 1.1e-19.  Hence
  35.  * the predicted relative error in the result is 2.3e-21 y .
  36.  *
  37.  *                      Relative error:
  38.  * arithmetic   domain     # trials      peak         rms
  39.  *
  40.  *    IEEE     +-1000       40000      2.8e-18      3.7e-19
  41.  * .001 < x < 1000, with log(x) uniformly distributed.
  42.  * -1000 < y < 1000, y uniformly distributed.
  43.  *
  44.  *    IEEE     0,8700       60000      6.5e-18      1.0e-18
  45.  * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
  46.  *
  47.  *
  48.  * ERROR MESSAGES:
  49.  *
  50.  *   message         condition      value returned
  51.  * pow overflow     x**y > MAXNUM      INFINITY
  52.  * pow underflow   x**y < 1/MAXNUM       0.0
  53.  * pow domain      x<0 and y noninteger  0.0
  54.  *
  55.  */
  57. /*
  58. Cephes Math Library Release 2.7:  May, 1998
  59. Copyright 1984, 1991, 1998 by Stephen L. Moshier
  60. */
  61.  
  62. /*
  63. Modified for mingw
  64. 2002-07-22 Danny Smith <dannysmith@users.sourceforge.net>
  65. */
  66.  
  67. #ifdef __MINGW32__
  68. #include "cephes_mconf.h"
  69. #else
  70. #include "mconf.h"
  71.  
  72. static char fname[] = {"powl"};
  73. #endif
  74.  
  75. #ifndef _SET_ERRNO
  76. #define _SET_ERRNO(x)
  77. #endif
  78.  
  79.  
  80. /* Table size */
  81. #define NXT 32
  82. /* log2(Table size) */
  83. #define LNXT 5
  84.  
  85. #ifdef UNK
  86. /* log(1+x) =  x - .5x^2 + x^3 *  P(z)/Q(z)
  87.  * on the domain  2^(-1/32) - 1  <=  x  <=  2^(1/32) - 1
  88.  */
  89. static long double P[] = {
  90.  8.3319510773868690346226E-4L,
  91.  4.9000050881978028599627E-1L,
  92.  1.7500123722550302671919E0L,
  93.  1.4000100839971580279335E0L,
  94. };
  95. static long double Q[] = {
  96. /* 1.0000000000000000000000E0L,*/
  97.  5.2500282295834889175431E0L,
  98.  8.4000598057587009834666E0L,
  99.  4.2000302519914740834728E0L,
  100. };
  101. /* A[i] = 2^(-i/32), rounded to IEEE long double precision.
  102.  * If i is even, A[i] + B[i/2] gives additional accuracy.
  103.  */
  104. static long double A[33] = {
  105.  1.0000000000000000000000E0L,
  106.  9.7857206208770013448287E-1L,
  107.  9.5760328069857364691013E-1L,
  108.  9.3708381705514995065011E-1L,
  109.  9.1700404320467123175367E-1L,
  110.  8.9735453750155359320742E-1L,
  111.  8.7812608018664974155474E-1L,
  112.  8.5930964906123895780165E-1L,
  113.  8.4089641525371454301892E-1L,
  114.  8.2287773907698242225554E-1L,
  115.  8.0524516597462715409607E-1L,
  116.  7.8799042255394324325455E-1L,
  117.  7.7110541270397041179298E-1L,
  118.  7.5458221379671136985669E-1L,
  119.  7.3841307296974965571198E-1L,
  120.  7.2259040348852331001267E-1L,
  121.  7.0710678118654752438189E-1L,
  122.  6.9195494098191597746178E-1L,
  123.  6.7712777346844636413344E-1L,
  124.  6.6261832157987064729696E-1L,
  125.  6.4841977732550483296079E-1L,
  126.  6.3452547859586661129850E-1L,
  127.  6.2092890603674202431705E-1L,
  128.  6.0762367999023443907803E-1L,
  129.  5.9460355750136053334378E-1L,
  130.  5.8186242938878875689693E-1L,
  131.  5.6939431737834582684856E-1L,
  132.  5.5719337129794626814472E-1L,
  133.  5.4525386633262882960438E-1L,
  134.  5.3357020033841180906486E-1L,
  135.  5.2213689121370692017331E-1L,
  136.  5.1094857432705833910408E-1L,
  137.  5.0000000000000000000000E-1L,
  138. };
  139. static long double B[17] = {
  140.  0.0000000000000000000000E0L,
  141.  2.6176170809902549338711E-20L,
  142. -1.0126791927256478897086E-20L,
  143.  1.3438228172316276937655E-21L,
  144.  1.2207982955417546912101E-20L,
  145. -6.3084814358060867200133E-21L,
  146.  1.3164426894366316434230E-20L,
  147. -1.8527916071632873716786E-20L,
  148.  1.8950325588932570796551E-20L,
  149.  1.5564775779538780478155E-20L,
  150.  6.0859793637556860974380E-21L,
  151. -2.0208749253662532228949E-20L,
  152.  1.4966292219224761844552E-20L,
  153.  3.3540909728056476875639E-21L,
  154. -8.6987564101742849540743E-22L,
  155. -1.2327176863327626135542E-20L,
  156.  0.0000000000000000000000E0L,
  157. };
  158.  
  159. /* 2^x = 1 + x P(x),
  160.  * on the interval -1/32 <= x <= 0
  161.  */
  162. static long double R[] = {
  163.  1.5089970579127659901157E-5L,
  164.  1.5402715328927013076125E-4L,
  165.  1.3333556028915671091390E-3L,
  166.  9.6181291046036762031786E-3L,
  167.  5.5504108664798463044015E-2L,
  168.  2.4022650695910062854352E-1L,
  169.  6.9314718055994530931447E-1L,
  170. };
  171.  
  172. #define douba(k) A[k]
  173. #define doubb(k) B[k]
  174. #define MEXP (NXT*16384.0L)
  175. /* The following if denormal numbers are supported, else -MEXP: */
  176. #ifdef DENORMAL
  177. #define MNEXP (-NXT*(16384.0L+64.0L))
  178. #else
  179. #define MNEXP (-NXT*16384.0L)
  180. #endif
  181. /* log2(e) - 1 */
  182. #define LOG2EA 0.44269504088896340735992L
  183. #endif
  184.  
  185.  
  186. #ifdef IBMPC
  187. static const unsigned short P[] = {
  188. 0xb804,0xa8b7,0xc6f4,0xda6a,0x3ff4, XPD
  189. 0x7de9,0xcf02,0x58c0,0xfae1,0x3ffd, XPD
  190. 0x405a,0x3722,0x67c9,0xe000,0x3fff, XPD
  191. 0xcd99,0x6b43,0x87ca,0xb333,0x3fff, XPD
  192. };
  193. static const unsigned short Q[] = {
  194. /* 0x0000,0x0000,0x0000,0x8000,0x3fff, */
  195. 0x6307,0xa469,0x3b33,0xa800,0x4001, XPD
  196. 0xfec2,0x62d7,0xa51c,0x8666,0x4002, XPD
  197. 0xda32,0xd072,0xa5d7,0x8666,0x4001, XPD
  198. };
  199. static const unsigned short A[] = {
  200. 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
  201. 0x033a,0x722a,0xb2db,0xfa83,0x3ffe, XPD
  202. 0xcc2c,0x2486,0x7d15,0xf525,0x3ffe, XPD
  203. 0xf5cb,0xdcda,0xb99b,0xefe4,0x3ffe, XPD
  204. 0x392f,0xdd24,0xc6e7,0xeac0,0x3ffe, XPD
  205. 0x48a8,0x7c83,0x06e7,0xe5b9,0x3ffe, XPD
  206. 0xe111,0x2a94,0xdeec,0xe0cc,0x3ffe, XPD
  207. 0x3755,0xdaf2,0xb797,0xdbfb,0x3ffe, XPD
  208. 0x6af4,0xd69d,0xfcca,0xd744,0x3ffe, XPD
  209. 0xe45a,0xf12a,0x1d91,0xd2a8,0x3ffe, XPD
  210. 0x80e4,0x1f84,0x8c15,0xce24,0x3ffe, XPD
  211. 0x27a3,0x6e2f,0xbd86,0xc9b9,0x3ffe, XPD
  212. 0xdadd,0x5506,0x2a11,0xc567,0x3ffe, XPD
  213. 0x9456,0x6670,0x4cca,0xc12c,0x3ffe, XPD
  214. 0x36bf,0x580c,0xa39f,0xbd08,0x3ffe, XPD
  215. 0x9ee9,0x62fb,0xaf47,0xb8fb,0x3ffe, XPD
  216. 0x6484,0xf9de,0xf333,0xb504,0x3ffe, XPD
  217. 0x2590,0xd2ac,0xf581,0xb123,0x3ffe, XPD
  218. 0x4ac6,0x42a1,0x3eea,0xad58,0x3ffe, XPD
  219. 0x0ef8,0xea7c,0x5ab4,0xa9a1,0x3ffe, XPD
  220. 0x38ea,0xb151,0xd6a9,0xa5fe,0x3ffe, XPD
  221. 0x6819,0x0c49,0x4303,0xa270,0x3ffe, XPD
  222. 0x11ae,0x91a1,0x3260,0x9ef5,0x3ffe, XPD
  223. 0x5539,0xd54e,0x39b9,0x9b8d,0x3ffe, XPD
  224. 0xa96f,0x8db8,0xf051,0x9837,0x3ffe, XPD
  225. 0x0961,0xfef7,0xefa8,0x94f4,0x3ffe, XPD
  226. 0xc336,0xab11,0xd373,0x91c3,0x3ffe, XPD
  227. 0x53c0,0x45cd,0x398b,0x8ea4,0x3ffe, XPD
  228. 0xd6e7,0xea8b,0xc1e3,0x8b95,0x3ffe, XPD
  229. 0x8527,0x92da,0x0e80,0x8898,0x3ffe, XPD
  230. 0x7b15,0xcc48,0xc367,0x85aa,0x3ffe, XPD
  231. 0xa1d7,0xac2b,0x8698,0x82cd,0x3ffe, XPD
  232. 0x0000,0x0000,0x0000,0x8000,0x3ffe, XPD
  233. };
  234. static const unsigned short B[] = {
  235. 0x0000,0x0000,0x0000,0x0000,0x0000, XPD
  236. 0x1f87,0xdb30,0x18f5,0xf73a,0x3fbd, XPD
  237. 0xac15,0x3e46,0x2932,0xbf4a,0xbfbc, XPD
  238. 0x7944,0xba66,0xa091,0xcb12,0x3fb9, XPD
  239. 0xff78,0x40b4,0x2ee6,0xe69a,0x3fbc, XPD
  240. 0xc895,0x5069,0xe383,0xee53,0xbfbb, XPD
  241. 0x7cde,0x9376,0x4325,0xf8ab,0x3fbc, XPD
  242. 0xa10c,0x25e0,0xc093,0xaefd,0xbfbd, XPD
  243. 0x7d3e,0xea95,0x1366,0xb2fb,0x3fbd, XPD
  244. 0x5d89,0xeb34,0x5191,0x9301,0x3fbd, XPD
  245. 0x80d9,0xb883,0xfb10,0xe5eb,0x3fbb, XPD
  246. 0x045d,0x288c,0xc1ec,0xbedd,0xbfbd, XPD
  247. 0xeded,0x5c85,0x4630,0x8d5a,0x3fbd, XPD
  248. 0x9d82,0xe5ac,0x8e0a,0xfd6d,0x3fba, XPD
  249. 0x6dfd,0xeb58,0xaf14,0x8373,0xbfb9, XPD
  250. 0xf938,0x7aac,0x91cf,0xe8da,0xbfbc, XPD
  251. 0x0000,0x0000,0x0000,0x0000,0x0000, XPD
  252. };
  253. static const unsigned short R[] = {
  254. 0xa69b,0x530e,0xee1d,0xfd2a,0x3fee, XPD
  255. 0xc746,0x8e7e,0x5960,0xa182,0x3ff2, XPD
  256. 0x63b6,0xadda,0xfd6a,0xaec3,0x3ff5, XPD
  257. 0xc104,0xfd99,0x5b7c,0x9d95,0x3ff8, XPD
  258. 0xe05e,0x249d,0x46b8,0xe358,0x3ffa, XPD
  259. 0x5d1d,0x162c,0xeffc,0xf5fd,0x3ffc, XPD
  260. 0x79aa,0xd1cf,0x17f7,0xb172,0x3ffe, XPD
  261. };
  262.  
  263. /* 10 byte sizes versus 12 byte */
  264. #define douba(k) (*(long double *)(&A[(sizeof( long double )/2)*(k)]))
  265. #define doubb(k) (*(long double *)(&B[(sizeof( long double )/2)*(k)]))
  266. #define MEXP (NXT*16384.0L)
  267. #ifdef DENORMAL
  268. #define MNEXP (-NXT*(16384.0L+64.0L))
  269. #else
  270. #define MNEXP (-NXT*16384.0L)
  271. #endif
  272. static const
  273. union
  274. {
  275.   unsigned short L[6];
  276.   long double ld;
  277. } log2ea = {{0xc2ef,0x705f,0xeca5,0xe2a8,0x3ffd, XPD}};
  278.  
  279. #define LOG2EA (log2ea.ld)
  280. /*
  281. #define LOG2EA 0.44269504088896340735992L
  282. */
  283. #endif
  284.  
  285. #ifdef MIEEE
  286. static long P[] = {
  287. 0x3ff40000,0xda6ac6f4,0xa8b7b804,
  288. 0x3ffd0000,0xfae158c0,0xcf027de9,
  289. 0x3fff0000,0xe00067c9,0x3722405a,
  290. 0x3fff0000,0xb33387ca,0x6b43cd99,
  291. };
  292. static long Q[] = {
  293. /* 0x3fff0000,0x80000000,0x00000000, */
  294. 0x40010000,0xa8003b33,0xa4696307,
  295. 0x40020000,0x8666a51c,0x62d7fec2,
  296. 0x40010000,0x8666a5d7,0xd072da32,
  297. };
  298. static long A[] = {
  299. 0x3fff0000,0x80000000,0x00000000,
  300. 0x3ffe0000,0xfa83b2db,0x722a033a,
  301. 0x3ffe0000,0xf5257d15,0x2486cc2c,
  302. 0x3ffe0000,0xefe4b99b,0xdcdaf5cb,
  303. 0x3ffe0000,0xeac0c6e7,0xdd24392f,
  304. 0x3ffe0000,0xe5b906e7,0x7c8348a8,
  305. 0x3ffe0000,0xe0ccdeec,0x2a94e111,
  306. 0x3ffe0000,0xdbfbb797,0xdaf23755,
  307. 0x3ffe0000,0xd744fcca,0xd69d6af4,
  308. 0x3ffe0000,0xd2a81d91,0xf12ae45a,
  309. 0x3ffe0000,0xce248c15,0x1f8480e4,
  310. 0x3ffe0000,0xc9b9bd86,0x6e2f27a3,
  311. 0x3ffe0000,0xc5672a11,0x5506dadd,
  312. 0x3ffe0000,0xc12c4cca,0x66709456,
  313. 0x3ffe0000,0xbd08a39f,0x580c36bf,
  314. 0x3ffe0000,0xb8fbaf47,0x62fb9ee9,
  315. 0x3ffe0000,0xb504f333,0xf9de6484,
  316. 0x3ffe0000,0xb123f581,0xd2ac2590,
  317. 0x3ffe0000,0xad583eea,0x42a14ac6,
  318. 0x3ffe0000,0xa9a15ab4,0xea7c0ef8,
  319. 0x3ffe0000,0xa5fed6a9,0xb15138ea,
  320. 0x3ffe0000,0xa2704303,0x0c496819,
  321. 0x3ffe0000,0x9ef53260,0x91a111ae,
  322. 0x3ffe0000,0x9b8d39b9,0xd54e5539,
  323. 0x3ffe0000,0x9837f051,0x8db8a96f,
  324. 0x3ffe0000,0x94f4efa8,0xfef70961,
  325. 0x3ffe0000,0x91c3d373,0xab11c336,
  326. 0x3ffe0000,0x8ea4398b,0x45cd53c0,
  327. 0x3ffe0000,0x8b95c1e3,0xea8bd6e7,
  328. 0x3ffe0000,0x88980e80,0x92da8527,
  329. 0x3ffe0000,0x85aac367,0xcc487b15,
  330. 0x3ffe0000,0x82cd8698,0xac2ba1d7,
  331. 0x3ffe0000,0x80000000,0x00000000,
  332. };
  333. static long B[51] = {
  334. 0x00000000,0x00000000,0x00000000,
  335. 0x3fbd0000,0xf73a18f5,0xdb301f87,
  336. 0xbfbc0000,0xbf4a2932,0x3e46ac15,
  337. 0x3fb90000,0xcb12a091,0xba667944,
  338. 0x3fbc0000,0xe69a2ee6,0x40b4ff78,
  339. 0xbfbb0000,0xee53e383,0x5069c895,
  340. 0x3fbc0000,0xf8ab4325,0x93767cde,
  341. 0xbfbd0000,0xaefdc093,0x25e0a10c,
  342. 0x3fbd0000,0xb2fb1366,0xea957d3e,
  343. 0x3fbd0000,0x93015191,0xeb345d89,
  344. 0x3fbb0000,0xe5ebfb10,0xb88380d9,
  345. 0xbfbd0000,0xbeddc1ec,0x288c045d,
  346. 0x3fbd0000,0x8d5a4630,0x5c85eded,
  347. 0x3fba0000,0xfd6d8e0a,0xe5ac9d82,
  348. 0xbfb90000,0x8373af14,0xeb586dfd,
  349. 0xbfbc0000,0xe8da91cf,0x7aacf938,
  350. 0x00000000,0x00000000,0x00000000,
  351. };
  352. static long R[] = {
  353. 0x3fee0000,0xfd2aee1d,0x530ea69b,
  354. 0x3ff20000,0xa1825960,0x8e7ec746,
  355. 0x3ff50000,0xaec3fd6a,0xadda63b6,
  356. 0x3ff80000,0x9d955b7c,0xfd99c104,
  357. 0x3ffa0000,0xe35846b8,0x249de05e,
  358. 0x3ffc0000,0xf5fdeffc,0x162c5d1d,
  359. 0x3ffe0000,0xb17217f7,0xd1cf79aa,
  360. };
  361.  
  362. #define douba(k) (*(long double *)&A[3*(k)])
  363. #define doubb(k) (*(long double *)&B[3*(k)])
  364. #define MEXP (NXT*16384.0L)
  365. #ifdef DENORMAL
  366. #define MNEXP (-NXT*(16384.0L+64.0L))
  367. #else
  368. #define MNEXP (-NXT*16382.0L)
  369. #endif
  370. static long L[3] = {0x3ffd0000,0xe2a8eca5,0x705fc2ef};
  371. #define LOG2EA (*(long double *)(&L[0]))
  372. #endif
  373.  
  374.  
  375. #define F W
  376. #define Fa Wa
  377. #define Fb Wb
  378. #define G W
  379. #define Ga Wa
  380. #define Gb u
  381. #define H W
  382. #define Ha Wb
  383. #define Hb Wb
  384.  
  385. #ifndef __MINGW32__
  386. extern long double MAXNUML;
  387. #endif
  388.  
  389. static VOLATILE long double z;
  390. static long double w, W, Wa, Wb, ya, yb, u;
  391.  
  392. #ifdef __MINGW32__
  393. static __inline__ long double reducl( long double );
  394. extern long double __powil ( long double, int );
  395. extern long double powl ( long double x, long double y);
  396. #else
  397. #ifdef ANSIPROT
  398. extern long double floorl ( long double );
  399. extern long double fabsl ( long double );
  400. extern long double frexpl ( long double, int * );
  401. extern long double ldexpl ( long double, int );
  402. extern long double polevll ( long double, void *, int );
  403. extern long double p1evll ( long double, void *, int );
  404. extern long double __powil ( long double, int );
  405. extern int isnanl ( long double );
  406. extern int isfinitel ( long double );
  407. static long double reducl( long double );
  408. extern int signbitl ( long double );
  409. #else
  410. long double floorl(), fabsl(), frexpl(), ldexpl();
  411. long double polevll(), p1evll(), __powil();
  412. static long double reducl();
  413. int isnanl(), isfinitel(), signbitl();
  414. #endif  /* __MINGW32__ */
  415.  
  416. #ifdef INFINITIES
  417. extern long double INFINITYL;
  418. #else
  419. #define INFINITYL MAXNUML
  420. #endif
  421.  
  422. #ifdef NANS
  423. extern long double NANL;
  424. #endif
  425. #ifdef MINUSZERO
  426. extern long double NEGZEROL;
  427. #endif
  428.  
  429. #endif /* __MINGW32__ */
  430.  
  431. #ifdef __MINGW32__
  432.  
  433. /* No error checking. We handle Infs and zeros ourselves.  */
  434. static __inline__ long double
  435. __fast_ldexpl (long double x, int expn)
  436. {
  437.   long double res;
  438.   __asm__ ("fscale"
  439.             : "=t" (res)
  440.             : "0" (x), "u" ((long double) expn));
  441.   return res;
  442. }
  443.  
  444. #define ldexpl  __fast_ldexpl
  445.  
  446. #endif
  447.  
  448.  
  449. long double powl( x, y )
  450. long double x, y;
  451. {
  452. /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
  453. int i, nflg, iyflg, yoddint;
  454. long e;
  455.  
  456. if( y == 0.0L )
  457.         return( 1.0L );
  458.  
  459. #ifdef NANS
  460. if( isnanl(x) )
  461.         {
  462.         _SET_ERRNO (EDOM);
  463.         return( x );
  464.         }
  465. if( isnanl(y) )
  466.         {
  467.         _SET_ERRNO (EDOM);
  468.         return( y );
  469.         }
  470. #endif
  471.  
  472. if( y == 1.0L )
  473.         return( x );
  474.  
  475. if( isinfl(y) && (x == -1.0L || x == 1.0L) )
  476.         return( y );
  477.  
  478. if( x == 1.0L )
  479.         return( 1.0L );
  480.  
  481. if( y >= MAXNUML )
  482.         {
  483.         _SET_ERRNO (ERANGE);
  484. #ifdef INFINITIES
  485.         if( x > 1.0L )
  486.                 return( INFINITYL );
  487. #else
  488.         if( x > 1.0L )
  489.                 return( MAXNUML );
  490. #endif
  491.         if( x > 0.0L && x < 1.0L )
  492.                 return( 0.0L );
  493. #ifdef INFINITIES
  494.         if( x < -1.0L )
  495.                 return( INFINITYL );
  496. #else
  497.         if( x < -1.0L )
  498.                 return( MAXNUML );
  499. #endif
  500.         if( x > -1.0L && x < 0.0L )
  501.                 return( 0.0L );
  502.         }
  503. if( y <= -MAXNUML )
  504.         {
  505.         _SET_ERRNO (ERANGE);
  506.         if( x > 1.0L )
  507.                 return( 0.0L );
  508. #ifdef INFINITIES
  509.         if( x > 0.0L && x < 1.0L )
  510.                 return( INFINITYL );
  511. #else
  512.         if( x > 0.0L && x < 1.0L )
  513.                 return( MAXNUML );
  514. #endif
  515.         if( x < -1.0L )
  516.                 return( 0.0L );
  517. #ifdef INFINITIES
  518.         if( x > -1.0L && x < 0.0L )
  519.                 return( INFINITYL );
  520. #else
  521.         if( x > -1.0L && x < 0.0L )
  522.                 return( MAXNUML );
  523. #endif
  524.         }
  525. if( x >= MAXNUML )
  526.         {
  527. #if INFINITIES
  528.         if( y > 0.0L )
  529.                 return( INFINITYL );
  530. #else
  531.         if( y > 0.0L )
  532.                 return( MAXNUML );
  533. #endif
  534.         return( 0.0L );
  535.         }
  536.  
  537. w = floorl(y);
  538. /* Set iyflg to 1 if y is an integer.  */
  539. iyflg = 0;
  540. if( w == y )
  541.         iyflg = 1;
  542.  
  543. /* Test for odd integer y.  */
  544. yoddint = 0;
  545. if( iyflg )
  546.         {
  547.         ya = fabsl(y);
  548.         ya = floorl(0.5L * ya);
  549.         yb = 0.5L * fabsl(w);
  550.         if( ya != yb )
  551.                 yoddint = 1;
  552.         }
  553.  
  554. if( x <= -MAXNUML )
  555.         {
  556.         if( y > 0.0L )
  557.                 {
  558. #ifdef INFINITIES
  559.                 if( yoddint )
  560.                         return( -INFINITYL );
  561.                 return( INFINITYL );
  562. #else
  563.                 if( yoddint )
  564.                         return( -MAXNUML );
  565.                 return( MAXNUML );
  566. #endif
  567.                 }
  568.         if( y < 0.0L )
  569.                 {
  570. #ifdef MINUSZERO
  571.                 if( yoddint )
  572.                         return( NEGZEROL );
  573. #endif
  574.                 return( 0.0 );
  575.                 }
  576.         }
  577.  
  578.  
  579. nflg = 0;       /* flag = 1 if x<0 raised to integer power */
  580. if( x <= 0.0L )
  581.         {
  582.         if( x == 0.0L )
  583.                 {
  584.                 if( y < 0.0 )
  585.                         {
  586. #ifdef MINUSZERO
  587.                         if( signbitl(x) && yoddint )
  588.                                 return( -INFINITYL );
  589. #endif
  590. #ifdef INFINITIES
  591.                         return( INFINITYL );
  592. #else
  593.                         return( MAXNUML );
  594. #endif
  595.                         }
  596.                 if( y > 0.0 )
  597.                         {
  598. #ifdef MINUSZERO
  599.                         if( signbitl(x) && yoddint )
  600.                                 return( NEGZEROL );
  601. #endif
  602.                         return( 0.0 );
  603.                         }
  604.                 if( y == 0.0L )
  605.                         return( 1.0L );  /*   0**0   */
  606.                 else  
  607.                         return( 0.0L );  /*   0**y   */
  608.                 }
  609.         else
  610.                 {
  611.                 if( iyflg == 0 )
  612.                         { /* noninteger power of negative number */
  613.                         mtherr( fname, DOMAIN );
  614.                         _SET_ERRNO (EDOM);
  615. #ifdef NANS
  616.                         return(NANL);
  617. #else
  618.                         return(0.0L);
  619. #endif
  620.                         }
  621.                 nflg = 1;
  622.                 }
  623.         }
  624.  
  625. /* Integer power of an integer.  */
  626.  
  627. if( iyflg )
  628.         {
  629.         i = w;
  630.         w = floorl(x);
  631.         if( (w == x) && (fabsl(y) < 32768.0) )
  632.                 {
  633.                 w = __powil( x, (int) y );
  634.                 return( w );
  635.                 }
  636.         }
  637.  
  638.  
  639. if( nflg )
  640.         x = fabsl(x);
  641.  
  642. /* separate significand from exponent */
  643. x = frexpl( x, &i );
  644. e = i;
  645.  
  646. /* find significand in antilog table A[] */
  647. i = 1;
  648. if( x <= douba(17) )
  649.         i = 17;
  650. if( x <= douba(i+8) )
  651.         i += 8;
  652. if( x <= douba(i+4) )
  653.         i += 4;
  654. if( x <= douba(i+2) )
  655.         i += 2;
  656. if( x >= douba(1) )
  657.         i = -1;
  658. i += 1;
  659.  
  660.  
  661. /* Find (x - A[i])/A[i]
  662.  * in order to compute log(x/A[i]):
  663.  *
  664.  * log(x) = log( a x/a ) = log(a) + log(x/a)
  665.  *
  666.  * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
  667.  */
  668. x -= douba(i);
  669. x -= doubb(i/2);
  670. x /= douba(i);
  671.  
  672.  
  673. /* rational approximation for log(1+v):
  674.  *
  675.  * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
  676.  */
  677. z = x*x;
  678. w = x * ( z * polevll( x, P, 3 ) / p1evll( x, Q, 3 ) );
  679. w = w - ldexpl( z, -1 );   /*  w - 0.5 * z  */
  680.  
  681. /* Convert to base 2 logarithm:
  682.  * multiply by log2(e) = 1 + LOG2EA
  683.  */
  684. z = LOG2EA * w;
  685. z += w;
  686. z += LOG2EA * x;
  687. z += x;
  688.  
  689. /* Compute exponent term of the base 2 logarithm. */
  690. w = -i;
  691. w = ldexpl( w, -LNXT ); /* divide by NXT */
  692. w += e;
  693. /* Now base 2 log of x is w + z. */
  694.  
  695. /* Multiply base 2 log by y, in extended precision. */
  696.  
  697. /* separate y into large part ya
  698.  * and small part yb less than 1/NXT
  699.  */
  700. ya = reducl(y);
  701. yb = y - ya;
  702.  
  703. /* (w+z)(ya+yb)
  704.  * = w*ya + w*yb + z*y
  705.  */
  706. F = z * y  +  w * yb;
  707. Fa = reducl(F);
  708. Fb = F - Fa;
  709.  
  710. G = Fa + w * ya;
  711. Ga = reducl(G);
  712. Gb = G - Ga;
  713.  
  714. H = Fb + Gb;
  715. Ha = reducl(H);
  716. w = ldexpl( Ga + Ha, LNXT );
  717.  
  718. /* Test the power of 2 for overflow */
  719. if( w > MEXP )
  720.         {
  721.         _SET_ERRNO (ERANGE);
  722.         mtherr( fname, OVERFLOW );
  723.         return( MAXNUML );
  724.         }
  725.  
  726. if( w < MNEXP )
  727.         {
  728.         _SET_ERRNO (ERANGE);
  729.         mtherr( fname, UNDERFLOW );
  730.         return( 0.0L );
  731.         }
  732.  
  733. e = w;
  734. Hb = H - Ha;
  735.  
  736. if( Hb > 0.0L )
  737.         {
  738.         e += 1;
  739.         Hb -= (1.0L/NXT);  /*0.0625L;*/
  740.         }
  741.  
  742. /* Now the product y * log2(x)  =  Hb + e/NXT.
  743.  *
  744.  * Compute base 2 exponential of Hb,
  745.  * where -0.0625 <= Hb <= 0.
  746.  */
  747. z = Hb * polevll( Hb, R, 6 );  /*    z  =  2**Hb - 1    */
  748.  
  749. /* Express e/NXT as an integer plus a negative number of (1/NXT)ths.
  750.  * Find lookup table entry for the fractional power of 2.
  751.  */
  752. if( e < 0 )
  753.         i = 0;
  754. else
  755.         i = 1;
  756. i = e/NXT + i;
  757. e = NXT*i - e;
  758. w = douba( e );
  759. z = w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
  760. z = z + w;
  761. z = ldexpl( z, i );  /* multiply by integer power of 2 */
  762.  
  763. if( nflg )
  764.         {
  765. /* For negative x,
  766.  * find out if the integer exponent
  767.  * is odd or even.
  768.  */
  769.         w = ldexpl( y, -1 );
  770.         w = floorl(w);
  771.         w = ldexpl( w, 1 );
  772.         if( w != y )
  773.                 z = -z; /* odd exponent */
  774.         }
  775.  
  776. return( z );
  777. }
  778.  
  779. static __inline__ long double
  780. __convert_inf_to_maxnum(long double x)
  781. {
  782.   if (isinf(x))
  783.     return (x > 0.0L ? MAXNUML : -MAXNUML);
  784.   else
  785.     return x;
  786. }
  787.  
  788.  
  789. /* Find a multiple of 1/NXT that is within 1/NXT of x. */
  790. static __inline__ long double reducl(x)
  791. long double x;
  792. {
  793. long double t;
  794.  
  795. /* If the call to ldexpl overflows, set it to MAXNUML.
  796.    This avoids Inf - Inf = Nan result when calculating the 'small'
  797.    part of a reduction.  Instead, the small part becomes Inf,
  798.    causing under/overflow when adding it to the 'large' part.
  799.    There must be a cleaner way of doing this. */
  800. t =  __convert_inf_to_maxnum (ldexpl( x, LNXT ));
  801. t = floorl( t );
  802. t = ldexpl( t, -LNXT );
  803. return(t);
  804. }
  805.