home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c082_122 / 2.ddi / MATHSRC.ZIP / POWL.CAS < prev    next >
Encoding:
Text File  |  1992-06-10  |  12.8 KB  |  393 lines
  1. /*------------------------------------------------------------------------
  2.  * filename - powl.cas
  3.  *
  4.  * function(s)
  5.  *        powl - long double power function, x^y
  6.  *-----------------------------------------------------------------------*/
  7.  
  8. /*
  9.  *      C/C++ Run Time Library - Version 5.0
  10.  *
  11.  *      Copyright (c) 1987, 1992 by Borland International
  12.  *      All Rights Reserved.
  13.  *
  14.  */
  15.  
  16.  
  17. #pragma inline
  18. #include <asmrules.h>
  19.  
  20. #include <_math.h>
  21. #include <math.h>
  22. #include <errno.h>
  23. #include <stddef.h>
  24.  
  25. #define  EXTPROC1(x)  asm call near ptr (x)
  26.  
  27. static  unsigned short   NANLOGL [5] = {0,0,0,0xC024, 0xFFFF};
  28.  
  29. /*--------------------------------------------------------------------------*
  30.  
  31. Name            powl - power function, x^y
  32.  
  33. Usage           long double powl(long double x, long double y);
  34.  
  35. Prototype in    math.h
  36.  
  37. Description     Return the value of x to the  power of y. If x is zero then
  38.                 y must be  positive, and if  x is negative  then y must  be
  39.                 integral.
  40.  
  41.                 First the special cases  Y == 0 or X == 0  or X == infinity
  42.                 are detected and given standard answers.
  43.  
  44.                 Two methods of calculation are used, depending upon whether
  45.                 Y is an integer of less than 64 bits. If not, then
  46.  
  47.                         X^Y = 2^(Log2(X) * Y)
  48.                             = DoExps ( DoLogs (X, Y), not scaled)
  49.  
  50.                 If Y is  an integer then it can be  represented as a binary
  51.                 number
  52.  
  53.                         Y = B0 + B1.2 + B2.+ .. Bn.2^n
  54.  
  55.                 where the  B coefficients are 0  or 1, and Bn  is always 1,
  56.                 for some n. The power of X is then calculated as:
  57.  
  58.                         Z = X;
  59.                         while (n-- > 0)
  60.                                 Z *= Z;
  61.                                 if (Bn) Z *= X;
  62.  
  63.                 which is  the standard trick  for fast integral  powers. It
  64.                 works for any X, positive or negative, if Y is not zero. In
  65.                 practice  it will  run faster than the  DoExps (DoLogs())
  66.                 method for  |Y|  <  100,  roughly,  and  slower for larger
  67.                 powers. Such large powers are very rare in actual usage.
  68.  
  69.                 Powers greater  than  2^64  in  theory  may  be  integers.
  70.                 However,  it is  also likely  that such  large numbers have
  71.                 lost  precision (especially  when you  consider that  the C
  72.                 data  type "long double"  is 53   bits precise).  These will  be
  73.                 treated as if  fractional. If X is positive  and very close
  74.                 to  1.0,  then an  answer  may  be  possible, but if X is
  75.                 negative an  exception is generated. The  rationale for the
  76.                 exception is  that if the  least bits of  Y have been  lost
  77.                 then  it is  not possible   to be  sure whether  the result
  78.                 should be positive or negative, so there is a total loss of
  79.                 precision.
  80.  
  81. Return value    Return the value of x to the power of y.
  82.                 When  the correct  value  would  overflow, powl returns the
  83.                 value _LHUGE_VAL.
  84.  
  85.                 If the argument x passed to powl  is less than or equal to 0
  86.                 and y  is not a whole  number, then errno is set to
  87.                         EDOM    Domain error
  88.  
  89.                 If x and y are both zero, then the return value is 1.0 and
  90.                 there is no error.  Many C compilers consider this a DOMAIN
  91.                 error as technically 0^0 is undefined.  There are continuous
  92.                 function f(x) and g(x) such that f(0) = 0 and g(0) = 0, but
  93.                 with the limit of f(x)/g(x) as x tends to 0 being any real
  94.                 number.  However, there is an elementary theorem that states
  95.                 that f(x) and g(x) are analytic and nonzero, then the limit
  96.                 is always 1.  Thus in a finite precision computing
  97.                 environment, it is hard to imagine a situation where a
  98.                 number other than 1 is desirable.
  99.  
  100. *---------------------------------------------------------------------------*/
  101.  
  102. #pragma warn -rvl       /* Function should return a value */
  103.  
  104. /*
  105. Identical to expl(), but calls __matherrl with different arguments.
  106. */
  107. static long double near pascal __expl (long double x)
  108. {
  109. asm     FLD     LONGDOUBLE (x)
  110. asm     mov     ax, x [8]       /* select exponent */
  111. asm     and     ah, 7Fh         /* remove sign bit */
  112. asm     cmp     ax, 3fffh+13
  113. asm     jb      exp_OK          /* expl (+-2^13) is the limit for long double */
  114.  
  115. exp_tooBig:
  116. asm     mov     ax, 0FFFFh      /* force extreme */
  117. asm     ja      exp_excess
  118. asm     mov     ax, x [6]
  119.  
  120. exp_excess:
  121. asm     test    BY0 (x [9]), 80h
  122. asm     jnz     exp_tooTiny
  123. asm     cmp     ax, 0B172h
  124. asm     jb      exp_OK
  125. asm     mov     si, OVERFLOW
  126. asm     jmp     short   exp_err
  127.  
  128. exp_tooTiny:
  129. asm     cmp     ax, 0B16Ch
  130. asm     jb      exp_OK
  131. asm     mov     si, UNDERFLOW
  132.  
  133. exp_err:
  134. asm     FSTP    ST(0)           /* discard ST */
  135. #pragma warn -ret
  136.         /* should use args to powl, but have no access */
  137.         return  __matherrl (_SI, "powl", NULL, NULL,
  138.                           (UNDERFLOW == _SI) ? 0.0 : _LHUGE_VAL);
  139. exp_OK:
  140.         __expld();
  141.         return;
  142. }
  143.  
  144. #pragma warn -rvl
  145. #pragma warn -powl
  146. #pragma warn -use
  147. long double _FARFUNC powl (long double x, long double y)
  148. {
  149.         long double temp;           /* also used as a 64-bit integer */
  150.         char        negate = 0;     /* boolean, negate after exp() ? */
  151.         volatile unsigned int Control;
  152.         volatile unsigned int Status;
  153.  
  154. asm     FLD     LONGDOUBLE (x)
  155. asm     mov     bx, 7FFFh                       /* mask just the exponent */
  156. asm     mov     ax, x [8]
  157. asm     and     ax, bx
  158. asm     jz      powl_ofZero
  159. asm     cmp     ax, bx
  160. asm     je      powl_ofInfinity
  161.  
  162. asm     FLD     LONGDOUBLE (y)
  163. asm     mov     ax, y [8]
  164. asm     and     ax, bx
  165. asm     jnz     temp1
  166. asm     jmp     powl_toZero
  167. temp1:
  168. asm     cmp     ax, bx
  169. asm     je      powl_toInfinity
  170. asm     jmp     powl_normal
  171.  
  172.  
  173. /***            Special cases           ***/
  174. /*
  175.   Raising any number to infinity is treated as a range error.
  176. */
  177. powl_toInfinity:
  178. asm     FSTP    LONGDOUBLE (temp)               /* propagate Y thru to result */
  179. asm     jmp     short   powl_discard
  180.  
  181. /*
  182.   Powers of infinity are range errors.
  183. */
  184. powl_ofInfinity:
  185. asm     FSTP    LONGDOUBLE (temp)               /* propagate X thru to result */
  186. asm     mov     ax, y[8]
  187. asm     or      ax, ax
  188. asm     jge     powl_overflow                   /* jump if exponent nonnegative */
  189. asm     mov     si, UNDERFLOW
  190.         temp = 0.0;
  191. asm     jmp     short   powl_complain
  192.  
  193. /*
  194.   Arrive here if Y is zero.  The zero'th power of any number is 1.
  195. */
  196. powl_toZero:
  197. asm     FSTP    ST(0)                           /* discard Y    */
  198. asm     FSTP    ST(0)                           /* discard X    */
  199. asm     FLD1
  200.         return; /* 1.0 */
  201.  
  202. powl_discard:
  203. asm     FSTP    ST(0)                           /* discard X    */
  204.  
  205. powl_overflow:
  206. asm     mov     si, OVERFLOW
  207. asm     jmp     short   powl_complain
  208.  
  209. /*
  210.   Powers of 0 are (EDOM, 1, 0) as Y ranges over (negative, zero, positive).
  211. */
  212. powl_ofZero:
  213. asm     FSTP    ST(0)                           /* discard  X   */
  214. asm     mov     ax, y [8]
  215. asm     or      ax, ax                          /* was Y positive ?     */
  216. asm     jg      powl_zero
  217. asm     mov     si, DOMAIN
  218. asm     je      powl_zz
  219.         temp = _LHUGE_VAL;
  220. asm     jmp     short   powl_complain
  221.  
  222. powl_zz:
  223.         temp = 1.0;
  224.  
  225. powl_complain:
  226.         return  __matherrl (_SI, "powl", &x, &y, temp);
  227.  
  228. powl_zero:
  229. asm     FLDZ
  230.         return; /* 0.0 */
  231.  
  232. /***            End of Special Cases            ***/
  233.  
  234.  
  235. /*
  236.   If arrived here then both x and y seem to be ordinary numbers.
  237. */
  238. powl_normal:
  239. asm     FCLEX
  240. asm     FRNDINT
  241. asm     FSTSW   W0 (Status)             /* is Y an integer      */
  242. asm     FWAIT
  243. asm     test    BY0 (Status), 20h       /* precision error if not */
  244. asm     jz      powl_integral
  245. asm     FSTP    ST(0)                   /* discard Y */
  246.  
  247. /*
  248.   Arrive here if the exponent exceeds integer range or if it contains
  249.   a fractional part.  Calculate using Log and Exp functions.  Just
  250.   x is on 87-stack.
  251. */
  252.  
  253. powl_fractional:
  254. /* make sure that x > 0 */
  255. asm     FTST
  256. asm     FSTSW   W0 (Status)             /* is Y an integer      */
  257. asm     FWAIT
  258. asm     mov     ax, Status
  259. asm     sahf
  260. asm     jae     powl_log
  261. asm     FSTP    ST(0)
  262.  
  263.         temp = *((long double *) NANLOGL);
  264. asm     mov     si, DOMAIN
  265. asm     jmp     short   powl_complain
  266.  
  267. powl_log:
  268. /* arg is > 0, so log cannot fail */
  269.  
  270. #ifdef _Windows
  271.         _f87_Log();
  272. #else
  273. asm     _FAST_  (_FLOG_)
  274. #endif
  275.  
  276. /* Disable underflow and overflow exceptions temporarily.
  277.  */
  278. asm     fstcw   Control                 /* save old control word */
  279. asm     fwait
  280. asm     mov     ax, Control
  281. asm     mov     cx, ax                  /* save copy in cx */
  282. asm     or      ax, 18h                 /* mask under/overflow */
  283. asm     mov     Control, ax
  284. asm     fldcw   Control
  285.  
  286. asm     FLD     LONGDOUBLE (y)
  287. asm     FMUL
  288.  
  289. /*  Check for exceptions.
  290.  */
  291. asm     fstsw   Status                  /* save the status */
  292. asm     fclex                           /* clear exceptions */
  293. asm     mov     Control, cx
  294. asm     fldcw   Control                 /* reload control word */
  295. asm     test    word ptr Status, 18h    /* did under/overflow occur? */
  296. asm     jz      powl_noexc              /* if so, discard result */
  297. asm     jmp     powl_discard
  298.  
  299. powl_noexc:
  300. asm     sub     sp, 10
  301. asm     mov     bx, sp
  302. asm     FSTP    LONGDOUBLE (SS_ [bx])
  303.         EXTPROC1 (__expl)
  304.  
  305.         if (negate)
  306.         {
  307.                 asm     FCHS
  308.         }
  309.         return;
  310.  
  311. /*
  312.   If arrived here then Y is some integer of up to 64 bits and has
  313.   been copied to temp. Y is ST(0), X is ST(1), AX is exponent of Y.
  314. */
  315. powl_integral:
  316. asm     sub     ax, 3FFFh                       /* remove the bias      */
  317. asm     cmp     ax, 63                          /* AX = n, the exponent */
  318. asm     jb      powl_trueIntegral
  319. asm     FSTP    ST(0)                           /* discard Y */
  320. powl_fracjmp:
  321. asm     jmp     short   powl_fractional
  322.  
  323. /*
  324.   The shift-and-add method is not accurate for extreme powers since
  325.   round off errors are magnified.  However, we cannot simply call for
  326.   evaluation like fractional powers because X may be negative and
  327.   fractional negative powers are treated as exceptions.
  328. */
  329. powl_trueIntegral:
  330. asm     cmp     al, 12
  331. asm     jb      powl_shiftAndAdd
  332. powl_unsafeRange:
  333. asm     FISTP   qword ptr (temp)                /* store an integer copy of Y */
  334. asm     test    BY0 (x [9]), 80h                /* X less than 0 ?      */
  335. asm     jz      powl_fracjmp                    /* X not signed, so no worry */
  336. asm     FCHS                                    /* make X absolute      */
  337. asm     test    BY0 (temp), 01h                 /* odd or even ?        */
  338. asm     jz      powl_fracjmp                    /* even powers are positive */
  339. /*
  340.   If we arrive here then X was negative and Y was odd.  Calculate with
  341.   abs(X) and then negate result.
  342. */
  343.         negate = 1;
  344. asm     jmp     short   powl_fracjmp
  345.  
  346.  
  347. /*
  348.   Arrive here for modest integral powers of any number.  We must also
  349.   check for overflow, by making a worst-case check on log (X^Y).  If
  350.   it has a potential to overflow, then we use the exp(log()) method.
  351. */
  352. powl_shiftAndAdd:
  353. asm     mov     bx, x [8]
  354. asm     shl     bx, 1
  355. asm     sub     bx, 7FFEh                       /* BX estimates log2 (X) */
  356. asm     mov     dx, bx
  357. asm     xchg    cx, ax
  358. asm     inc     cx                              /* 2^CL is max possible Y */
  359. asm     shl     bx, cl                          /* multiply BX by max Y */
  360. asm     sar     bx, cl
  361. asm     dec     cx
  362. asm     xchg    ax, cx
  363. asm     cmp     bx, dx                          /* did BX lose any bits ? */
  364. asm     jne     powl_unsafeRange
  365.  
  366. asm     FLD     ST (1)                          /* Z = X        */
  367. asm     mov     dx, y [6]
  368. asm     shl     dx,1
  369.  
  370. powl_iWhileBit:
  371. asm     dec     al
  372. asm     jl      powl_maybeInverse
  373.  
  374. asm     FMUL    ST(0), ST(0)                    /* Z *= Z       */
  375.  
  376. asm     shl     dx, 1
  377. asm     jnc     powl_iWhileBit
  378.  
  379. asm     FMUL    ST(0), ST(2)                    /* Z *= X       */
  380.  
  381. asm     jmp     short   powl_iWhileBit
  382.  
  383. powl_maybeInverse:
  384. asm     FSTP    ST(1)                           /* overwrite Y  */
  385. asm     test    BY0 (y [9]), 80h                /* was Y a negative power ? */
  386. asm     FSTP    ST(1)                           /* overwrite X */
  387. asm     jz      powl_iDone
  388. asm     FLD1
  389. asm     FDIVRP  ST(1), ST(0)                    /*   if so, invert result. */
  390. powl_iDone:
  391.         return;
  392. }
  393.