home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c082_122 / 2.ddi / MATHSRC.ZIP / EXPLD.CAS < prev    next >
Encoding:
Text File  |  1992-06-10  |  4.1 KB  |  115 lines
  1. /*------------------------------------------------------------------------
  2.  * filename - expld.cas
  3.  *
  4.  * function(s)
  5.  *        expld - long double exponential function helper function
  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. static unsigned short log2hi[5] =
  26.         { 0x0000,0xD1D0,0x17F7,0xB172,0x3FFE };         /* +0.693E+0 */
  27. static unsigned short log2hi_ov_2[5] =
  28.         { 0x0000,0xD1D0,0x17F7,0xB172,0x3FFD };         /* +0.34657E+0 */
  29. static unsigned short log2lo[5] =
  30.         { 0xFC0D,0x4C67,0x361C,0x8654,0xBFCE };         /* -0.1864E-14 */
  31.  
  32.  
  33. /*--------------------------------------------------------------------------*
  34.  
  35. Name            __expld - long double exponential function
  36.  
  37. Usage           void near __expld(void)
  38.  
  39. Prototype in    _math.h
  40.  
  41. Description     expl calculates the exponent of the argument on the top
  42.                 of the FPU stack (e^TOS).  The result is left on the
  43.                 top of the FPU stack.  No error checking on the argument
  44.                 is performed.  This function is a helper for expl(), sinh(),
  45.                 and coshl().
  46.  
  47. Return value    __expld does not return a value, but leaves the result
  48.                 e^TOS on the TOS.
  49.  
  50. *---------------------------------------------------------------------------*/
  51.  
  52. void pascal near __expld(void)
  53. {
  54.         unsigned short status;
  55.  
  56. asm     FLD     LONGDOUBLE (log2hi)     /* st(0)= log2hi st(1)= x */
  57. asm     FLD     st(1)           /* st(0)= x; st(1)= log2hi; st(2)= x */
  58. asm     FPREM                   /* st(0)=hi; st(1)=log2hi; st(2)=x */
  59. asm     FLD     LONGDOUBLE (log2hi_ov_2)        /* st(0)=log2hi/2 */
  60. asm     FCOMP   st(1)           /* log2hi/2 >= hi ? */
  61. asm     FSTSW   status
  62. asm     FWAIT
  63. asm     mov     ah, status[1]
  64. asm     sahf
  65. asm     jae     expl_test_neg
  66. asm     FSUB    st(0), st(1)    /* st(0)= hi; st(1)= log2hi; st(2)= x */
  67. asm     jmp     short expl_begin
  68.  
  69. expl_test_neg:
  70. asm     FLD     LONGDOUBLE (log2hi_ov_2)        /* st(0)=-log2hi/2 */
  71. asm     FCHS
  72. asm     FCOMP   st(1)           /* -log2hi/2 <= hi ? */
  73. asm     FSTSW   status
  74. asm     FWAIT
  75. asm     mov     ah, status[1]
  76. asm     sahf
  77. asm     jbe     expl_begin
  78. asm     FADD    st(0), st(1)    /* st(0)=hi; st(1)=log2hi; st(2)=x */
  79.  
  80. expl_begin:
  81. asm     FXCH    st(2)           /* st(0)= x; st(1)= log2hi; st(2)= hi */
  82. asm     FSUB    st(0), st(2)    /* st(0)= x-hi; st(1)= log2hi; st(2)= hi */
  83. asm     FDIVRP  st(1), st(0)    /* st(0)= (x-hi)/log2hi; st(1)= hi */
  84. asm     FRNDINT                 /* st(0)= i= rint((x-hi)/log2hi); st(1)= hi */
  85. asm     FXCH    st(1)           /* st(0)= hi; st(1)= i */
  86. asm     FLD     LONGDOUBLE (log2lo)     /* st(0)=log2lo; st(1)=hi; st(2)=i */
  87. asm     FMUL    st(0), st(2)    /* st(0)= i*log2lo; st(1)= hi; st(2)= i */
  88. asm     FSUBP   st(1), st(0)    /* st(0)= r= x-i*log2; st(1)= i */
  89. asm     FLDL2E                  /* st(0)= log2(e); st(1)= r; st(2)= i */
  90. asm     FMUL                    /* st(0)= f; st(1)= i */
  91.  
  92. asm     FTST                    /* f >= 0? */
  93. asm     FSTSW   status
  94. asm     FWAIT
  95. asm     mov     ah, status[1]
  96. asm     sahf
  97. asm     jb      expl_red_neg
  98. asm     F2XM1                   /* st(0)= 2^f-1; st(1)= i */
  99. asm     FLD1                    /* st(0)= 1; st(1)= 2^f-1; st(2)= i */
  100. asm     FADDP   st(1), st(0)    /* st(0)= 2^f; st(1)= i */
  101. asm     jmp     short expl_red_next
  102.  
  103. expl_red_neg:
  104. asm     FCHS                    /* st(0)= -f; st(1)= i */
  105. asm     F2XM1                   /* st(0)= 2^(-f)-1; st(1)= i */
  106. asm     FLD1                    /* st(0)= 1; st(1)= 2^(-f)-1; st(2)= i */
  107. asm     FADD    st(1), st(0)    /* st(0)= 1; st(1)= 2^(-f); st(2)= i */
  108. asm     FDIVRP  st(1), st(0)    /* st(0)= 2^f; st(1)= i */
  109.  
  110. expl_red_next:
  111. asm     FSCALE                  /* st(0)= (2^i)*(2^f) */
  112. asm     FSTP    st(1)
  113.         return;
  114. }
  115.