home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c065 / 2.ddi / MATH.ZIP / POW10.CAS < prev    next >
Encoding:
Text File  |  1990-06-07  |  6.2 KB  |  201 lines
  1. /*------------------------------------------------------------------------
  2.  * filename - pow10.cas
  3.  *
  4.  * function(s)
  5.  *        pow10 - power function, 10^p
  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. #pragma inline
  19. #include <asmrules.h>
  20.  
  21. #include <_math.h>
  22. #include <math.h>
  23. #include <errno.h>
  24. #include <stddef.h>
  25.  
  26. typedef unsigned short int  extend [5];     /* 80-bit constants */
  27.  
  28. static  const   float   e0to7 [8] =
  29. {
  30.     1, 1.e1,  1.e2,  1.e3,  1.e4,  1.e5,  1.e6,  1.e7,
  31. };
  32.  
  33. /* Exponents > 4932 become infinities.  Exponents < -4932 become 0. */
  34.  
  35. /* These values have been calculated with extra precision to ensure  */
  36. /* that the last bit is rounded correctly.                           */
  37. static  const   float   e8    = 1.e8;
  38. static  const   double  e16   = 1.e16;
  39. static  const   extend  e32   = {0xB59E, 0x2B70, 0xADA8, 0x9DC5, 0x4069};
  40. static  const   extend  e64   = {0xA6D5, 0xFFCF, 0x1F49, 0xC278, 0x40D3};
  41. static  const   extend  e128  = {0x8CE0, 0x80E9, 0x47C9, 0x93BA, 0x41A8};
  42. static  const   extend  e256  = {0xDE8E, 0x9DF9, 0xEBFB, 0xAA7E, 0x4351};
  43. static  const   extend  e512  = {0x91C7, 0xA60E, 0xA0AE, 0xE319, 0x46A3};
  44. static  const   extend  e1024 = {0x0C17, 0x8175, 0x7586, 0xC976, 0x4D48};
  45. static  const   extend  e2048 = {0x5DE5, 0xC53D, 0x3B5D, 0x9E8B, 0x5A92};
  46. static  const   extend  e4096 = {0x979B, 0x8A20, 0x5202, 0xC460, 0x7525};
  47. static  const   float    eINF  = 1.0/0.0;
  48.  
  49.  
  50. /*--------------------------------------------------------------------------*
  51.  
  52. Name            pow10 - power function, 10^p
  53.  
  54. Usage           double  pow10(int  p);
  55.  
  56. Prototype in    math.h
  57.  
  58. Description     Calculate 10  raised to power.  A lookup table  is used for
  59.                 values  from  10  through  10^7,  then this is augmented by
  60.                 multiplying with  table entries for  10^8/16/32/64/128/256,
  61.                 512/1024/2048/4096 which allows any power up to the
  62.                 implementation limit of 4932.
  63.                 The usual range of double precision is  10e308 but pow10
  64.                 has a wider  range so that it can be used for long double
  65.                 calculations when the result is retrieved off the TOS
  66.                 explicitly by inline code doing a store of an 80 bit #.
  67.  
  68.                 Negative powers are provided by a final division.
  69.  
  70.                 All registers  are preserved except   AX ! This  is done to
  71.                 enable  use by  xcvt(), which  was designed  to assume  its
  72.                 registers will be undisturbed.
  73.  
  74. Return value    pow10 returns 10^p.
  75.  
  76. *---------------------------------------------------------------------------*/
  77. #pragma warn -rvl
  78. double  pow10  (int  p)
  79. {
  80. #define MAX_87_EXP      4932
  81.  
  82. #ifdef __HUGE__
  83. asm     mov     ax, seg e0to7
  84. asm     mov     DS, ax
  85. #endif
  86.  
  87. /*--------------------------------------------------------------------------
  88.         Take care of all the easy special cases up front.
  89. --------------------------------------------------------------------------*/
  90. asm     mov     ax, p
  91.         if ((int)_AX < -MAX_87_EXP)     /* Extremely small -> Zero      */
  92.                 {
  93. asm             FLDZ
  94. asm             jmp     p10_end
  95.                 }
  96.         if ((int)_AX >  MAX_87_EXP)     /* Extremely large -> Infinity  */
  97.                 {
  98. asm             FLD     FLOAT(eINF)
  99. asm             jmp     p10_end
  100.                 }
  101.         if ((int)_AX == 0)              /* 10^0 -> 1.0                  */
  102.                 {
  103. asm             FLD1
  104. asm             jmp     p10_end
  105.                 }
  106.  
  107. /*--------------------------------------------------------------------------
  108.                 The non-trivial cases require some calculation.
  109. --------------------------------------------------------------------------*/
  110. /*asm   mov     ax, p*/
  111. asm     or      ax, ax
  112. asm     jnl     p10_abs
  113. asm     neg     ax
  114.  
  115. p10_abs:
  116. asm     mov     si, 7
  117. asm     and     si, ax
  118. asm     shl     si, 1
  119. asm     shl     si, 1
  120. asm     FLD     FLOAT (e0to7 [si])
  121.  
  122. asm     shr     ax, 1
  123. asm     shr     ax, 1
  124. asm     shr     ax, 1
  125.  
  126. p10_maybe8:
  127. asm     shr     ax, 1
  128. asm     jnc     p10_maybe16
  129. asm     FMUL    FLOAT (e8)
  130.  
  131. p10_maybe16:
  132. asm     jnz     keep_going
  133. asm     jmp     p10_checkSign           /* optimization, skip if all done */
  134. keep_going:
  135. asm     shr     ax, 1
  136. asm     jnc     p10_maybe32
  137. asm     FMUL    DOUBLE (e16)
  138.  
  139. p10_maybe32:
  140. asm     shr     ax, 1
  141. asm     jnc     p10_maybe64
  142. asm     FLD     LONGDOUBLE (e32)
  143. asm     FMUL
  144.  
  145. p10_maybe64:
  146. asm     shr     ax, 1
  147. asm     jnc     p10_maybe128
  148. asm     FLD     LONGDOUBLE (e64)
  149. asm     FMUL
  150.  
  151. p10_maybe128:
  152. asm     shr     ax, 1
  153. asm     jnc     p10_maybe256
  154. asm     FLD     LONGDOUBLE (e128)
  155. asm     FMUL
  156.  
  157. p10_maybe256:
  158. asm     shr     ax, 1
  159. asm     jnc     p10_maybe512
  160. asm     FLD     LONGDOUBLE (e256)
  161. asm     FMUL
  162.  
  163. p10_maybe512:
  164. asm     shr     ax, 1
  165. asm     jnc     p10_maybe1024
  166. asm     FLD     LONGDOUBLE (e512)
  167. asm     FMUL
  168.  
  169. p10_maybe1024:
  170. asm     shr     ax, 1
  171. asm     jnc     p10_maybe2048
  172. asm     FLD     LONGDOUBLE (e1024)
  173. asm     FMUL
  174.  
  175. p10_maybe2048:
  176. asm     shr     ax, 1
  177. asm     jnc     p10_maybe4096
  178. asm     FLD     LONGDOUBLE (e2048)
  179. asm     FMUL
  180.  
  181. p10_maybe4096:
  182. asm     shr     ax, 1
  183. asm     jnc     p10_checkSign
  184. asm     FLD     LONGDOUBLE (e4096)
  185. asm     FMUL
  186.  
  187. p10_checkSign:
  188. asm     test    BY1 (p), 80h
  189. asm     jz      p10_end
  190.  
  191. /* 10^(-n) = 1 / 10^n, so we need the reciprocal of TOS. */
  192.  
  193. asm     FDIVR   FLOAT (e0to7)           /* TOS = 1.0 / TOS */
  194.  
  195. /* Now the value 10^p is on TOS. */
  196.  
  197. p10_end:
  198.     return;
  199. }
  200. #pragma warn .rvl
  201.