home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c082_122 / 2.ddi / MATHSRC.ZIP / POW.CAS < prev    next >
Encoding:
Text File  |  1992-06-10  |  12.3 KB  |  385 lines
  1. /*------------------------------------------------------------------------
  2.  * filename - pow.cas
  3.  *
  4.  * function(s)
  5.  *        pow - 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   NANLOG [4] = {0,0,0x0480, 0xFFF8};
  28.  
  29. /*--------------------------------------------------------------------------*
  30.  
  31. Name            pow - power function, x^y
  32.  
  33. Usage           double  pow(double x, 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 "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, pow returns the
  83.                 value HUGE_VAL.
  84.  
  85.                 If the argument x passed to pow  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. Similar to exp(), but with
  106.         long double argument
  107.         error if result outside double range
  108.         error reported as a pow() error
  109. */
  110. long double near pascal __expl(long double x)
  111. {
  112. asm     FLD     LONGDOUBLE (x)
  113.  
  114. asm     mov     ax, 7FFFh
  115. asm     and     ax, x [8]               /* select exponent */
  116. asm     cmp     ax, 4008h
  117. asm     jnb     exp_tooBig      /* exp (+-709) is the limit for double */
  118.  
  119. exp_justFits:
  120.  
  121. #ifdef _Windows
  122.         _f87_Exp();
  123. #else
  124. asm     _FAST_  (_FEXP_)
  125. #endif
  126.  
  127.         return;
  128.  
  129. exp_tooBig:
  130. asm     mov     ax, 0FFFFh      /* force extreme */
  131. asm     ja      exp_excess
  132. asm     mov     ax, x [6]
  133.  
  134. exp_excess:
  135. asm     test    BY0 (x [9]), 80h
  136. asm     jnz     exp_tooTiny
  137. asm     cmp     ax, 0B172h
  138. asm     jb      exp_justFits
  139. asm     mov     si, OVERFLOW
  140. asm     jmp     short   exp_err
  141.  
  142. exp_tooTiny:
  143. asm     cmp     ax, 0B172h
  144. asm     jb      exp_justFits
  145. asm     mov     si, UNDERFLOW
  146.  
  147. exp_err:
  148. asm     FSTP    ST(0)           /* discard ST */
  149. #pragma warn -ret
  150.         /* should use args to pow, but have no access */
  151.         return  _matherr (_SI, "pow", NULL, NULL,
  152.                           (UNDERFLOW == _SI) ? 0.0 : HUGE_VAL);
  153. }
  154.  
  155. #pragma warn -rvl
  156. #pragma warn -pow
  157. #pragma warn -use
  158.  
  159. double _FARFUNC pow( double x, double y )
  160. {
  161.         double      temp;           /* also used as a 64-bit integer */
  162.         unsigned    sword;          /* status-word from testing y   */
  163.         char        negate = 0;     /* boolean, negate after exp() ? */
  164.  
  165. asm     FLD     DOUBLE (x)
  166. asm     mov     bx, 7FF0h                       /* mask just the exponent */
  167. asm     mov     ax, x [6]
  168. asm     and     ax, bx
  169. asm     jz      pow_ofZero
  170. asm     cmp     ax, bx
  171. asm     je      pow_ofInfinity
  172.  
  173. asm     FLD     DOUBLE (y)
  174. asm     mov     ax, y [6]
  175. asm     and     ax, bx
  176. asm     jz      pow_toZero
  177. asm     cmp     ax, bx
  178. asm     je      pow_toInfinity
  179. asm     jmp     pow_normal
  180.  
  181.  
  182. /***            Special cases           ***/
  183. /*
  184.   Raising any number to infinity is treated as a range error.
  185. */
  186. pow_toInfinity:
  187. asm     FSTP    DOUBLE (temp)                   /* propagate Y thru to result */
  188. asm     jmp     short   pow_discard
  189.  
  190. /*
  191.   Powers of infinity are range errors.
  192. */
  193. pow_ofInfinity:
  194. asm     FSTP    DOUBLE (temp)                   /* propagate X thru to result */
  195. asm     mov     ax, y[6]
  196. asm     or      ax, ax
  197. asm     jge     pow_overflow                    /* jump if exponent nonnegative */
  198. asm     mov     si, UNDERFLOW
  199.         temp = 0.0;
  200. asm     jmp     short   pow_complain
  201.  
  202. /*
  203.   Arrive here if Y is zero.  The zero'th power of any number is 1.
  204. */
  205. pow_toZero:
  206. asm     FSTP    ST(0)                           /* discard Y    */
  207. asm     FSTP    ST(0)                           /* discard X    */
  208. asm     FLD1
  209.         return; /* 1.0 */
  210.  
  211. pow_discard:
  212. asm     FSTP    ST(0)                           /* discard X    */
  213.  
  214. pow_overflow:
  215. asm     mov     si, OVERFLOW
  216. asm     jmp     short   pow_complain
  217.  
  218. /*
  219.   Powers of 0 are (EDOM, 1, 0) as Y ranges over (negative, zero, positive).
  220. */
  221. pow_ofZero:
  222. asm     FSTP    ST(0)                           /* discard  X   */
  223. asm     mov     ax, y [6]
  224. asm     or      ax, ax                          /* was Y positive ?     */
  225. asm     jg      pow_zero
  226. asm     mov     si, DOMAIN
  227. asm     je      pow_zz
  228.         temp = HUGE_VAL;
  229. asm     jmp     short   pow_complain
  230.  
  231. pow_zz:
  232.         temp = 1.0;
  233.  
  234. pow_complain:
  235.         return  _matherr (_SI, "pow", &x, &y, temp);
  236.  
  237. pow_zero:
  238. asm     FLDZ
  239.         return; /* 0.0 */
  240.  
  241. /***            End of Special Cases            ***/
  242.  
  243.  
  244. /*
  245.   If arrived here then both x and y seem to be ordinary numbers.
  246. */
  247. pow_normal:
  248. asm     FCLEX
  249. asm     FRNDINT
  250. asm     FSTSW   W0 (sword)              /* is Y an integer      */
  251. asm     FWAIT
  252. asm     test    BY0 (sword), 20h                /* precision error if not */
  253. asm     jz      pow_integral
  254. asm     FSTP    ST(0)                   /* discard Y */
  255.  
  256. /*
  257.   Arrive here if the exponent exceeds integer range or if it contains
  258.   a fractional part.  Calculate using Log and Exp functions.  Just
  259.   x is on 87-stack.
  260. */
  261.  
  262. pow_fractional:
  263. /* make sure that x > 0 */
  264. asm     FTST
  265. asm     FSTSW   W0 (sword)              /* is Y an integer      */
  266. asm     FWAIT
  267. asm     mov     ax, sword
  268. asm     sahf
  269. asm     jae     pow_log
  270. asm     FSTP    ST(0)
  271.  
  272.         temp = *((double *) NANLOG);
  273. asm     mov     si, DOMAIN
  274. asm     jmp     short   pow_complain
  275.  
  276. pow_log:
  277. /* arg is > 0, so log cannot fail */
  278.  
  279. #ifdef _Windows
  280.         _f87_Log();
  281. #else
  282. asm     _FAST_  (_FLOG_)
  283. #endif
  284.  
  285. asm     FMUL    DOUBLE (y)
  286.  
  287. asm     sub     sp, 10
  288. asm     mov     bx, sp
  289. asm     FSTP    LONGDOUBLE (SS_ [bx])
  290.         EXTPROC1 (__expl)
  291.  
  292.         if (negate)
  293.         {
  294.                 asm     FCHS
  295.         }
  296.         return;
  297.  
  298. /*
  299.   If arrived here then Y is some integer of up to 64 bits and has
  300.   been copied to temp. Y is ST(0), X is ST(1), AX is exponent of Y.
  301. */
  302. pow_integral:
  303. asm     mov     cl, 4
  304. asm     ror     ax, cl
  305. asm     sub     ax, 3FFh                        /* remove the bias      */
  306. asm     cmp     ax, 63                          /* AX = n, the exponent */
  307. asm     jb      pow_trueIntegral
  308. asm     FSTP    ST(0)                           /* discard Y */
  309. asm     jmp     short   pow_fractional
  310.  
  311. /*
  312.   The shift-and-add method is not accurate for extreme powers since
  313.   round off errors are magnified.  However, we cannot simply call for
  314.   evaluation like fractional powers because X may be negative and
  315.   fractional negative powers are treated as exceptions.
  316. */
  317. pow_trueIntegral:
  318. asm     cmp     al, 12
  319. asm     jb      pow_shiftAndAdd
  320. pow_unsafeRange:
  321. asm     FISTP   qword ptr (temp)                /* store an integer copy of Y */
  322. asm     test    BY0 (x [7]), 80h                /* X less than 0 ?      */
  323. asm     jz      pow_fractional                  /* X not signed, so no worry */
  324. asm     FCHS                                    /* make X absolute      */
  325. asm     test    BY0 (temp), 01h                 /* odd or even ?        */
  326. asm     jz      pow_fractional                  /* even powers are positive */
  327. /*
  328.   If we arrive here then X was negative and Y was odd.  Calculate with
  329.   abs(X) and then negate result.
  330. */
  331.         negate = 1;
  332. asm     jmp     short   pow_fractional
  333.  
  334.  
  335. /*
  336.   Arrive here for modest integral powers of any number.  We must also
  337.   check for overflow, by making a worst-case check on log (X^Y).  If
  338.   it has a potential to overflow, then we use the exp(log()) method.
  339. */
  340. pow_shiftAndAdd:
  341. asm     mov     bx, x [6]
  342. asm     shl     bx, 1
  343. asm     sub     bx, 7FE0h                       /* BX estimates log2 (X) */
  344. asm     mov     dx, bx
  345. asm     xchg    cx, ax
  346. asm     inc     cx                              /* 2^CL is max possible Y */
  347. asm     shl     bx, cl                          /* multiply BX by max Y */
  348. asm     sar     bx, cl
  349. asm     dec     cx
  350. asm     xchg    ax, cx
  351. asm     cmp     bx, dx                          /* did BX lose any bits ? */
  352. asm     jne     pow_unsafeRange
  353.  
  354. asm     FLD     ST (1)                          /* Z = X        */
  355. asm     mov     dx, y [4]
  356. asm     mov     bl, y [6]
  357. asm     and     bl,  0Fh                        /* most significant nibble */
  358. asm     and     dl, 0F0h                        /* DX is the next 12 bits */
  359. asm     or      dl, bl
  360. asm     ror     dx, cl                          /* top 16 bits of fraction */
  361.  
  362. pow_iWhileBit:
  363. asm     dec     al
  364. asm     jl      pow_maybeInverse
  365.  
  366. asm     FMUL    ST(0), ST(0)                    /* Z *= Z       */
  367.  
  368. asm     shl     dx, 1
  369. asm     jnc     pow_iWhileBit
  370.  
  371. asm     FMUL    ST(0), ST(2)                    /* Z *= X       */
  372.  
  373. asm     jmp     short   pow_iWhileBit
  374.  
  375. pow_maybeInverse:
  376. asm     FSTP    ST(1)                           /* overwrite Y  */
  377. asm     test    BY0 (y [7]), 80h                /* was Y a negative power ? */
  378. asm     FSTP    ST(1)                           /* overwrite X */
  379. asm     jz      pow_iDone
  380. asm     FLD1
  381. asm     FDIVRP  ST(1), ST(0)                    /*   if so, invert result. */
  382. pow_iDone:
  383.         return;
  384. }
  385.