home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ldouble / powl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  13.0 KB  |  547 lines
  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      MAXNUM
  52.  * pow underflow   x**y < 1/MAXNUM       0.0
  53.  * pow domain      x<0 and y noninteger  0.0
  54.  *
  55.  */
  56.  
  57. /*
  58. Cephes Math Library Release 2.2:  January, 1991
  59. Copyright 1984, 1991 by Stephen L. Moshier
  60. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  61. */
  62.  
  63.  
  64. #include "mconf.h"
  65.  
  66. static char fname[] = {"powl"};
  67.  
  68. /* Table size */
  69. #define NXT 32
  70. /* log2(Table size) */
  71. #define LNXT 5
  72.  
  73. #ifdef UNK
  74. /* log(1+x) =  x - .5x^2 + x^3 *  P(z)/Q(z)
  75.  * on the domain  2^(-1/32) - 1  <=  x  <=  2^(1/32) - 1
  76.  */
  77. static long double P[] = {
  78.  8.3319510773868690346226E-4L,
  79.  4.9000050881978028599627E-1L,
  80.  1.7500123722550302671919E0L,
  81.  1.4000100839971580279335E0L,
  82. };
  83. static long double Q[] = {
  84. /* 1.0000000000000000000000E0L,*/
  85.  5.2500282295834889175431E0L,
  86.  8.4000598057587009834666E0L,
  87.  4.2000302519914740834728E0L,
  88. };
  89. /* A[i] = 2^(-i/32), rounded to IEEE long double precision.
  90.  * If i is even, A[i] + B[i/2] gives additional accuracy.
  91.  */
  92. static long double A[33] = {
  93.  1.0000000000000000000000E0L,
  94.  9.7857206208770013448287E-1L,
  95.  9.5760328069857364691013E-1L,
  96.  9.3708381705514995065011E-1L,
  97.  9.1700404320467123175367E-1L,
  98.  8.9735453750155359320742E-1L,
  99.  8.7812608018664974155474E-1L,
  100.  8.5930964906123895780165E-1L,
  101.  8.4089641525371454301892E-1L,
  102.  8.2287773907698242225554E-1L,
  103.  8.0524516597462715409607E-1L,
  104.  7.8799042255394324325455E-1L,
  105.  7.7110541270397041179298E-1L,
  106.  7.5458221379671136985669E-1L,
  107.  7.3841307296974965571198E-1L,
  108.  7.2259040348852331001267E-1L,
  109.  7.0710678118654752438189E-1L,
  110.  6.9195494098191597746178E-1L,
  111.  6.7712777346844636413344E-1L,
  112.  6.6261832157987064729696E-1L,
  113.  6.4841977732550483296079E-1L,
  114.  6.3452547859586661129850E-1L,
  115.  6.2092890603674202431705E-1L,
  116.  6.0762367999023443907803E-1L,
  117.  5.9460355750136053334378E-1L,
  118.  5.8186242938878875689693E-1L,
  119.  5.6939431737834582684856E-1L,
  120.  5.5719337129794626814472E-1L,
  121.  5.4525386633262882960438E-1L,
  122.  5.3357020033841180906486E-1L,
  123.  5.2213689121370692017331E-1L,
  124.  5.1094857432705833910408E-1L,
  125.  5.0000000000000000000000E-1L,
  126. };
  127. static long double B[17] = {
  128.  0.0000000000000000000000E0L,
  129.  2.6176170809902549338711E-20L,
  130. -1.0126791927256478897086E-20L,
  131.  1.3438228172316276937655E-21L,
  132.  1.2207982955417546912101E-20L,
  133. -6.3084814358060867200133E-21L,
  134.  1.3164426894366316434230E-20L,
  135. -1.8527916071632873716786E-20L,
  136.  1.8950325588932570796551E-20L,
  137.  1.5564775779538780478155E-20L,
  138.  6.0859793637556860974380E-21L,
  139. -2.0208749253662532228949E-20L,
  140.  1.4966292219224761844552E-20L,
  141.  3.3540909728056476875639E-21L,
  142. -8.6987564101742849540743E-22L,
  143. -1.2327176863327626135542E-20L,
  144.  0.0000000000000000000000E0L,
  145. };
  146.  
  147. /* 2^x = 1 + x P(x),
  148.  * on the interval -1/32 <= x <= 0
  149.  */
  150. static long double R[] = {
  151.  1.5089970579127659901157E-5L,
  152.  1.5402715328927013076125E-4L,
  153.  1.3333556028915671091390E-3L,
  154.  9.6181291046036762031786E-3L,
  155.  5.5504108664798463044015E-2L,
  156.  2.4022650695910062854352E-1L,
  157.  6.9314718055994530931447E-1L,
  158. };
  159.  
  160. #define douba(k) A[k]
  161. #define doubb(k) B[k]
  162. #define MEXP (NXT*16384.0L)
  163. /* The following if denormal numbers are supported, else -MEXP: */
  164. #define MNEXP (-NXT*(16384.0L-64.0L))
  165. /* log2(e) - 1 */
  166. #define LOG2EA 0.44269504088896340735992L
  167. #endif
  168.  
  169.  
  170. #ifdef IBMPC
  171. static short P[] = {
  172. 0xb804,0xa8b7,0xc6f4,0xda6a,0x3ff4,
  173. 0x7de9,0xcf02,0x58c0,0xfae1,0x3ffd,
  174. 0x405a,0x3722,0x67c9,0xe000,0x3fff,
  175. 0xcd99,0x6b43,0x87ca,0xb333,0x3fff,
  176. };
  177. static short Q[] = {
  178. /* 0x0000,0x0000,0x0000,0x8000,0x3fff, */
  179. 0x6307,0xa469,0x3b33,0xa800,0x4001,
  180. 0xfec2,0x62d7,0xa51c,0x8666,0x4002,
  181. 0xda32,0xd072,0xa5d7,0x8666,0x4001,
  182. };
  183. static short A[] = {
  184. 0x0000,0x0000,0x0000,0x8000,0x3fff,
  185. 0x033a,0x722a,0xb2db,0xfa83,0x3ffe,
  186. 0xcc2c,0x2486,0x7d15,0xf525,0x3ffe,
  187. 0xf5cb,0xdcda,0xb99b,0xefe4,0x3ffe,
  188. 0x392f,0xdd24,0xc6e7,0xeac0,0x3ffe,
  189. 0x48a8,0x7c83,0x06e7,0xe5b9,0x3ffe,
  190. 0xe111,0x2a94,0xdeec,0xe0cc,0x3ffe,
  191. 0x3755,0xdaf2,0xb797,0xdbfb,0x3ffe,
  192. 0x6af4,0xd69d,0xfcca,0xd744,0x3ffe,
  193. 0xe45a,0xf12a,0x1d91,0xd2a8,0x3ffe,
  194. 0x80e4,0x1f84,0x8c15,0xce24,0x3ffe,
  195. 0x27a3,0x6e2f,0xbd86,0xc9b9,0x3ffe,
  196. 0xdadd,0x5506,0x2a11,0xc567,0x3ffe,
  197. 0x9456,0x6670,0x4cca,0xc12c,0x3ffe,
  198. 0x36bf,0x580c,0xa39f,0xbd08,0x3ffe,
  199. 0x9ee9,0x62fb,0xaf47,0xb8fb,0x3ffe,
  200. 0x6484,0xf9de,0xf333,0xb504,0x3ffe,
  201. 0x2590,0xd2ac,0xf581,0xb123,0x3ffe,
  202. 0x4ac6,0x42a1,0x3eea,0xad58,0x3ffe,
  203. 0x0ef8,0xea7c,0x5ab4,0xa9a1,0x3ffe,
  204. 0x38ea,0xb151,0xd6a9,0xa5fe,0x3ffe,
  205. 0x6819,0x0c49,0x4303,0xa270,0x3ffe,
  206. 0x11ae,0x91a1,0x3260,0x9ef5,0x3ffe,
  207. 0x5539,0xd54e,0x39b9,0x9b8d,0x3ffe,
  208. 0xa96f,0x8db8,0xf051,0x9837,0x3ffe,
  209. 0x0961,0xfef7,0xefa8,0x94f4,0x3ffe,
  210. 0xc336,0xab11,0xd373,0x91c3,0x3ffe,
  211. 0x53c0,0x45cd,0x398b,0x8ea4,0x3ffe,
  212. 0xd6e7,0xea8b,0xc1e3,0x8b95,0x3ffe,
  213. 0x8527,0x92da,0x0e80,0x8898,0x3ffe,
  214. 0x7b15,0xcc48,0xc367,0x85aa,0x3ffe,
  215. 0xa1d7,0xac2b,0x8698,0x82cd,0x3ffe,
  216. 0x0000,0x0000,0x0000,0x8000,0x3ffe,
  217. };
  218. static short B[] = {
  219. 0x0000,0x0000,0x0000,0x0000,0x0000,
  220. 0x1f87,0xdb30,0x18f5,0xf73a,0x3fbd,
  221. 0xac15,0x3e46,0x2932,0xbf4a,0xbfbc,
  222. 0x7944,0xba66,0xa091,0xcb12,0x3fb9,
  223. 0xff78,0x40b4,0x2ee6,0xe69a,0x3fbc,
  224. 0xc895,0x5069,0xe383,0xee53,0xbfbb,
  225. 0x7cde,0x9376,0x4325,0xf8ab,0x3fbc,
  226. 0xa10c,0x25e0,0xc093,0xaefd,0xbfbd,
  227. 0x7d3e,0xea95,0x1366,0xb2fb,0x3fbd,
  228. 0x5d89,0xeb34,0x5191,0x9301,0x3fbd,
  229. 0x80d9,0xb883,0xfb10,0xe5eb,0x3fbb,
  230. 0x045d,0x288c,0xc1ec,0xbedd,0xbfbd,
  231. 0xeded,0x5c85,0x4630,0x8d5a,0x3fbd,
  232. 0x9d82,0xe5ac,0x8e0a,0xfd6d,0x3fba,
  233. 0x6dfd,0xeb58,0xaf14,0x8373,0xbfb9,
  234. 0xf938,0x7aac,0x91cf,0xe8da,0xbfbc,
  235. 0x0000,0x0000,0x0000,0x0000,0x0000,
  236. };
  237. static short R[] = {
  238. 0xa69b,0x530e,0xee1d,0xfd2a,0x3fee,
  239. 0xc746,0x8e7e,0x5960,0xa182,0x3ff2,
  240. 0x63b6,0xadda,0xfd6a,0xaec3,0x3ff5,
  241. 0xc104,0xfd99,0x5b7c,0x9d95,0x3ff8,
  242. 0xe05e,0x249d,0x46b8,0xe358,0x3ffa,
  243. 0x5d1d,0x162c,0xeffc,0xf5fd,0x3ffc,
  244. 0x79aa,0xd1cf,0x17f7,0xb172,0x3ffe,
  245. };
  246.  
  247. #define douba(k) (*(long double *)(&A[5*(k)]))
  248. #define doubb(k) (*(long double *)(&B[5*(k)]))
  249. #define MEXP (NXT*16384.0L)
  250. #define MNEXP (-NXT*16384.0L)
  251. static short L[5] = {0xc2ef,0x705f,0xeca5,0xe2a8,0x3ffd};
  252. #define LOG2EA (*(long double *)(&L[0]))
  253. #endif
  254.  
  255. #ifdef MIEEE
  256. static long P[] = {
  257. 0x3ff40000,0xda6ac6f4,0xa8b7b804,
  258. 0x3ffd0000,0xfae158c0,0xcf027de9,
  259. 0x3fff0000,0xe00067c9,0x3722405a,
  260. 0x3fff0000,0xb33387ca,0x6b43cd99,
  261. };
  262. static long Q[] = {
  263. /* 0x3fff0000,0x80000000,0x00000000, */
  264. 0x40010000,0xa8003b33,0xa4696307,
  265. 0x40020000,0x8666a51c,0x62d7fec2,
  266. 0x40010000,0x8666a5d7,0xd072da32,
  267. };
  268. static long A[] = {
  269. 0x3fff0000,0x80000000,0x00000000,
  270. 0x3ffe0000,0xfa83b2db,0x722a033a,
  271. 0x3ffe0000,0xf5257d15,0x2486cc2c,
  272. 0x3ffe0000,0xefe4b99b,0xdcdaf5cb,
  273. 0x3ffe0000,0xeac0c6e7,0xdd24392f,
  274. 0x3ffe0000,0xe5b906e7,0x7c8348a8,
  275. 0x3ffe0000,0xe0ccdeec,0x2a94e111,
  276. 0x3ffe0000,0xdbfbb797,0xdaf23755,
  277. 0x3ffe0000,0xd744fcca,0xd69d6af4,
  278. 0x3ffe0000,0xd2a81d91,0xf12ae45a,
  279. 0x3ffe0000,0xce248c15,0x1f8480e4,
  280. 0x3ffe0000,0xc9b9bd86,0x6e2f27a3,
  281. 0x3ffe0000,0xc5672a11,0x5506dadd,
  282. 0x3ffe0000,0xc12c4cca,0x66709456,
  283. 0x3ffe0000,0xbd08a39f,0x580c36bf,
  284. 0x3ffe0000,0xb8fbaf47,0x62fb9ee9,
  285. 0x3ffe0000,0xb504f333,0xf9de6484,
  286. 0x3ffe0000,0xb123f581,0xd2ac2590,
  287. 0x3ffe0000,0xad583eea,0x42a14ac6,
  288. 0x3ffe0000,0xa9a15ab4,0xea7c0ef8,
  289. 0x3ffe0000,0xa5fed6a9,0xb15138ea,
  290. 0x3ffe0000,0xa2704303,0x0c496819,
  291. 0x3ffe0000,0x9ef53260,0x91a111ae,
  292. 0x3ffe0000,0x9b8d39b9,0xd54e5539,
  293. 0x3ffe0000,0x9837f051,0x8db8a96f,
  294. 0x3ffe0000,0x94f4efa8,0xfef70961,
  295. 0x3ffe0000,0x91c3d373,0xab11c336,
  296. 0x3ffe0000,0x8ea4398b,0x45cd53c0,
  297. 0x3ffe0000,0x8b95c1e3,0xea8bd6e7,
  298. 0x3ffe0000,0x88980e80,0x92da8527,
  299. 0x3ffe0000,0x85aac367,0xcc487b15,
  300. 0x3ffe0000,0x82cd8698,0xac2ba1d7,
  301. 0x3ffe0000,0x80000000,0x00000000,
  302. };
  303. static long B[27] = {
  304. 0x00000000,0x00000000,0x00000000,
  305. 0x3fbd0000,0xf73a18f5,0xdb301f87,
  306. 0xbfbc0000,0xbf4a2932,0x3e46ac15,
  307. 0x3fb90000,0xcb12a091,0xba667944,
  308. 0x3fbc0000,0xe69a2ee6,0x40b4ff78,
  309. 0xbfbb0000,0xee53e383,0x5069c895,
  310. 0x3fbc0000,0xf8ab4325,0x93767cde,
  311. 0xbfbd0000,0xaefdc093,0x25e0a10c,
  312. 0x3fbd0000,0xb2fb1366,0xea957d3e,
  313. 0x3fbd0000,0x93015191,0xeb345d89,
  314. 0x3fbb0000,0xe5ebfb10,0xb88380d9,
  315. 0xbfbd0000,0xbeddc1ec,0x288c045d,
  316. 0x3fbd0000,0x8d5a4630,0x5c85eded,
  317. 0x3fba0000,0xfd6d8e0a,0xe5ac9d82,
  318. 0xbfb90000,0x8373af14,0xeb586dfd,
  319. 0xbfbc0000,0xe8da91cf,0x7aacf938,
  320. 0x00000000,0x00000000,0x00000000,
  321. };
  322. static long R[] = {
  323. 0x3fee0000,0xfd2aee1d,0x530ea69b,
  324. 0x3ff20000,0xa1825960,0x8e7ec746,
  325. 0x3ff50000,0xaec3fd6a,0xadda63b6,
  326. 0x3ff80000,0x9d955b7c,0xfd99c104,
  327. 0x3ffa0000,0xe35846b8,0x249de05e,
  328. 0x3ffc0000,0xf5fdeffc,0x162c5d1d,
  329. 0x3ffe0000,0xb17217f7,0xd1cf79aa,
  330. };
  331.  
  332. #define douba(k) (*(long double *)&A[3*(k)])
  333. #define doubb(k) (*(long double *)&B[3*(k)])
  334. #define MEXP (NXT*16384.0L)
  335. #define MNEXP (-NXT*16382.0L)
  336. static long L[3] = {0x3ffd0000,0xe2a8eca5,0x705fc2ef};
  337. #define LOG2EA (*(long double *)(&L[0]))
  338. #endif
  339.  
  340.  
  341. #define F W
  342. #define Fa Wa
  343. #define Fb Wb
  344. #define G W
  345. #define Ga Wa
  346. #define Gb u
  347. #define H W
  348. #define Ha Wb
  349. #define Hb Wb
  350.  
  351. extern long double MAXNUML;
  352. static VOLATILE long double z;
  353. static long double w, W, Wa, Wb, ya, yb, u, v;
  354. long double floorl(), fabsl(), frexpl(), ldexpl();
  355. long double reducl(), polevll(), p1evll(), powil();
  356.  
  357.  
  358.  
  359. long double powl( x, y )
  360. long double x, y;
  361. {
  362. /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
  363. int i, nflg;
  364. long e;
  365.  
  366. nflg = 0;    /* flag = 1 if x<0 raised to integer power */
  367. w = floorl(y);
  368. if( (w == y) && (fabsl(w) < 32768.0) )
  369.     {
  370.     i = w;
  371.     w = powil( x, i );
  372.     return( w );
  373.     }
  374.  
  375.  
  376. if( x <= 0.0L )
  377.     {
  378.     if( x == 0.0L )
  379.         {
  380.         if( y == 0.0L )
  381.             return( 1.0L );  /*   0**0   */
  382.         else  
  383.             return( 0.0L );  /*   0**y   */
  384.         }
  385.     else
  386.         {
  387.         if( w != y )
  388.             { /* noninteger power of negative number */
  389.             mtherr( fname, DOMAIN );
  390.             return(0.0L);
  391.             }
  392.         nflg = 1;
  393.         x = fabsl(x);
  394.         }
  395.     }
  396.  
  397. /* separate significand from exponent */
  398. x = frexpl( x, &i );
  399. e = i;
  400.  
  401. /* find significand in antilog table A[] */
  402. i = 1;
  403. if( x <= douba(17) )
  404.     i = 17;
  405. if( x <= douba(i+8) )
  406.     i += 8;
  407. if( x <= douba(i+4) )
  408.     i += 4;
  409. if( x <= douba(i+2) )
  410.     i += 2;
  411. if( x >= douba(1) )
  412.     i = -1;
  413. i += 1;
  414.  
  415.  
  416. /* Find (x - A[i])/A[i]
  417.  * in order to compute log(x/A[i]):
  418.  *
  419.  * log(x) = log( a x/a ) = log(a) + log(x/a)
  420.  *
  421.  * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
  422.  */
  423. x -= douba(i);
  424. x -= doubb(i/2);
  425. x /= douba(i);
  426.  
  427.  
  428. /* rational approximation for log(1+v):
  429.  *
  430.  * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
  431.  */
  432. z = x*x;
  433. w = x * ( z * polevll( x, P, 3 ) / p1evll( x, Q, 3 ) );
  434. /*w = (x *  z * polevll( x, P, 3 )) / p1evll( x, Q, 4 );*/
  435. w = w - ldexpl( z, -1 );   /*  w - 0.5 * z  */
  436.  
  437. /* Convert to base 2 logarithm:
  438.  * multiply by log2(e) = 1 + LOG2EA
  439.  */
  440. z = LOG2EA * w;
  441. z += w;
  442. z += LOG2EA * x;
  443. z += x;
  444.  
  445. /* Compute exponent term of the base 2 logarithm. */
  446. w = -i;
  447. w = ldexpl( w, -LNXT );    /* divide by NXT */
  448. w += e;
  449. /* Now base 2 log of x is w + z. */
  450.  
  451. /* Multiply base 2 log by y, in extended precision.
  452.  
  453. /* separate y into large part ya
  454.  * and small part yb less than 1/NXT
  455.  */
  456. ya = reducl(y);
  457. yb = y - ya;
  458.  
  459. /* (w+z)(ya+yb)
  460.  * = w*ya + w*yb + z*y
  461.  */
  462. F = z * y  +  w * yb;
  463. Fa = reducl(F);
  464. Fb = F - Fa;
  465.  
  466. G = Fa + w * ya;
  467. Ga = reducl(G);
  468. Gb = G - Ga;
  469.  
  470. H = Fb + Gb;
  471. Ha = reducl(H);
  472. w = ldexpl( Ga+Ha, LNXT );
  473.  
  474. /* Test the power of 2 for overflow */
  475. if( w > MEXP )
  476.     {
  477.     printf( "w = %.4Le ", w );
  478.     mtherr( fname, OVERFLOW );
  479.     return( MAXNUML );
  480.     }
  481.  
  482. if( w < MNEXP )
  483.     {
  484.     printf( "w = %.4Le ", w );
  485.     mtherr( fname, UNDERFLOW );
  486.     return( 0.0 );
  487.     }
  488.  
  489. e = w;
  490. Hb = H - Ha;
  491.  
  492. if( Hb > 0.0L )
  493.     {
  494.     e += 1;
  495.     Hb -= (1.0L/NXT);  /*0.0625L;*/
  496.     }
  497.  
  498. /* Now the product y * log2(x)  =  Hb + e/NXT.
  499.  *
  500.  * Compute base 2 exponential of Hb,
  501.  * where -0.0625 <= Hb <= 0.
  502.  */
  503. z = Hb * polevll( Hb, R, 6 );  /*    z  =  2**Hb - 1    */
  504.  
  505. /* Express e/NXT as an integer plus a negative number of (1/NXT)ths.
  506.  * Find lookup table entry for the fractional power of 2.
  507.  */
  508. if( e < 0 )
  509.     i = 0;
  510. else
  511.     i = 1;
  512. i = e/NXT + i;
  513. e = NXT*i - e;
  514. w = douba( e );
  515. z = w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
  516. z = z + w;
  517. z = ldexpl( z, i );  /* multiply by integer power of 2 */
  518.  
  519. if( nflg )
  520.     {
  521. /* For negative x,
  522.  * find out if the integer exponent
  523.  * is odd or even.
  524.  */
  525.     w = ldexpl( y, -1 );
  526.     w = floorl(w);
  527.     w = ldexpl( w, 1 );
  528.     if( w != y )
  529.         z = -z; /* odd exponent */
  530.     }
  531.  
  532. return( z );
  533. }
  534.  
  535.  
  536. /* Find a multiple of 1/NXT that is within 1/NXT of x. */
  537. static long double reducl(x)
  538. long double x;
  539. {
  540. long double t;
  541.  
  542. t = ldexpl( x, LNXT );
  543. t = floorl( t );
  544. t = ldexpl( t, -LNXT );
  545. return(t);
  546. }
  547.