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