home *** CD-ROM | disk | FTP | other *** search
/ C!T ROM 2 / ctrom_ii_b.zip / ctrom_ii_b / PROGRAM / PASCAL / TPL60N19 / ARISOURC / F48FEXP.ASM < prev    next >
Assembly Source File  |  1993-01-25  |  11KB  |  248 lines
  1.  
  2. ; *******************************************************
  3. ; *                                                     *
  4. ; *     Turbo Pascal Runtime Library Version 6.0        *
  5. ; *     Real Exponentiation                             *
  6. ; *                                                     *
  7. ; *     Copyright (C) 1989-1992 Norbert Juffa           *
  8. ; *                                                     *
  9. ; *******************************************************
  10.  
  11.              TITLE   F48FEXP
  12.  
  13.              INCLUDE SE.ASM
  14.  
  15.  
  16. CODE         SEGMENT BYTE PUBLIC
  17.  
  18.              ASSUME  CS:CODE
  19.  
  20. ; Externals
  21.  
  22.              EXTRN   RealAdd:NEAR,RealSub:NEAR,RealMul:NEAR,RealFloat:NEAR
  23.              EXTRN   RealTrunc:NEAR,RealSqr:NEAR,RealDivRev:NEAR
  24.              EXTRN   RealMulNoChk:NEAR,RealReduceMP:NEAR,ROverflow:NEAR
  25.              EXTRN   RealMulF:NEAR,RealMulFNoChk:NEAR
  26.  
  27.  
  28. ; Publics
  29.  
  30.              PUBLIC  RExp
  31.  
  32.              IFDEF   EXTENSIONS
  33.              PUBLIC  RPower2,RPower10
  34.              ENDIF
  35.  
  36. ;-------------------------------------------------------------------------------
  37. ; Routines RExp, RPower2, and RPower10 exponentiate their argument x to base e,
  38. ; base 2, and base 10, respectively. All routines do the exponentiation of an
  39. ; argument x by calculating e^g * 2^n, where g is in -0.5*ln(2)..0.5*ln(2).
  40. ; While 2^n is a simple scaling operation, e^g is computed using a rational
  41. ; approximation of the Pade type. The computation of i and g differs for the
  42. ; three functions. For RPower2, n is computed as n = Round (x), and g as g =
  43. ; (x - n) * ln (2). For RExp, the computation proceeds as follows: n = Round
  44. ; (x/ln(2)), g = x - ln(2)*n, where ln(2) is represented by two constants c1
  45. ; and c2 that provide more than 40 bits of precision, so that the actual
  46. ; computation becomes x - c1*n - c2*n. This way, loss of accuracy due to the
  47. ; argument reduction is prevented. The computation of RPower10 is similar to
  48. ; that of the RExp function n = Round (x/log10(2)), g = x - log10(2) * n by
  49. ; the multiple precision process described above. Finally, g = g * ln(10).
  50. ; Actually, these different argument reduction schemes are all derived from
  51. ; the same process by elimination of unnecessary steps: n = Round (x/logB (2)),
  52. ; g = (x - n*logB(2)) * ln (B), where B is the base for exponentation.
  53. ;
  54. ; If the result of the exponentiation routines underflows the floating point
  55. ; format, zero is returned. If the computation results in overflow, error 205
  56. ; is returned through the error handler.
  57. ;
  58. ; GPZ := (2.001347076967e1*g^2+8.408086858319e2)*g
  59. ; QZ  := (g^2+1.801617224813e2)*g^2+1.681617371664e3
  60. ;
  61. ; e^g := 2 * (0.5 + GPZ / (QZ - GPZ))
  62. ;
  63. ; This approximation has a theoretical maximum relative error of 2.98e-14.
  64. ; Maximum observed error when evaluated in REAL arithmetic is 1.74e-12.
  65. ;
  66. ; INPUT:     DX:BX:AX  argument
  67. ;
  68. ; OUTPUT:    DX:BX:AX  2^argument, e^argument, 10^argument depending on function
  69. ;
  70. ; DESTROYS:  AX,BX,CX,DX,SI,DI,Flags
  71. ;-------------------------------------------------------------------------------
  72.  
  73.              IFDEF   EXTENSIONS
  74.  
  75. RPower10     PROC    FAR
  76.              PUSH    BP                ; save caller's framepointer
  77.              MOV     BP,OFFSET Pow10Cst; address table of constants for pow10
  78.              JMP     $power_10         ; compute power of ten
  79. Pow10Cst     DW      0007Fh, 09A84h, 09A20h ; log(2)-eps = 3.010299955495e-1
  80.              DW      0005Fh, 0F799h, 0FBCFh ;        eps = 1.145110089921e-10
  81.              DW      0CD82h, 0784Bh, 0549Ah ; 1/log(2)   = 3.321928094887e+0
  82.              DW      0AB82h, 08DDDh, 0135Dh ; ln (10)    = 2.302585092994e+0
  83. RPower10     ENDP
  84.  
  85.              ALIGN   4
  86.  
  87. RPower2      PROC    FAR
  88.              PUSH    BP                ; save caller's framepointer
  89.              MOV     BP, OFFSET Pow2Cst-18;
  90.              PUSH    AX                ; save
  91.              MOV     SI, BX            ;  x
  92.              MOV     DI, DX            ;   ditto
  93.              MOV     CH, 0FFh          ; signal rounding desired
  94.              CALL    RealTrunc         ; compute n := round (x)
  95.              POP     CX                ; removed saved word
  96.              CMP     CL, 88h           ; abs (x) < 128 ?
  97.              JNB     $argumnt_err      ; no, overflow or underflow
  98.              PUSH    AX                ; save n
  99.              PUSH    CX                ; save lsw of x
  100.              CALL    RealFloat         ; compute float(n)
  101.              POP     CX                ; restore x to DI:SI:CX
  102.              XOR     DH, 80h           ; -n
  103.              CALL    RealAdd           ; compute t := -n + x = n - x
  104.              JMP     $power_2          ; compute power of two
  105.  
  106. Pow2Cst      DW      0D280h, 017F7h, 03172h ; ln (2) = 6.931471805599e-1
  107.  
  108. RPower2      ENDP
  109.  
  110.              ENDIF
  111.  
  112.  
  113. RExp_Ext     PROC    FAR
  114. $argumnt_err:POP     BP                ; restore caller's frame pointer
  115.              OR      DI, DI            ; argument negative ?
  116.              JS      $exp_zero         ; yes, result = 0 to machine precision
  117.              JMP     $exp_overfl       ; no, overflow, signal error
  118. $exp_zero:   XOR     AX, AX            ; load
  119.              MOV     BX, AX            ;  a
  120.              CWD                       ;   zero
  121.              RET                       ; done
  122. RExp_Ext     ENDP
  123.  
  124.              ALIGN   4
  125.  
  126. RExp         PROC    FAR
  127.              PUSH    BP                ; save frame pointer
  128.              MOV     BP, OFFSET PowECst; pointer to table of constants
  129. $power_10:   PUSH    DX                ; save
  130.              PUSH    BX                ;  argument x
  131.              PUSH    AX                ;   on stack
  132.              MOV     CX, CS:[BP+12]    ; load
  133.              MOV     SI, CS:[BP+14]    ;  constant
  134.              MOV     DI, CS:[BP+16]    ;   1/lg(2)
  135.              CALL    RealMulFNoChk     ; compute y := x/lg(2) (DI:SI:CX <> 0)
  136.              POP     CX                ; get
  137.              POP     SI                ;  back
  138.              POP     DI                ;   x
  139.              CMP     AL, 88h           ; abs (y) < 128 ?
  140.              JNB     $argumnt_err      ; no, underflow or overflow
  141.              PUSH    CX                ; save lsw of x
  142.              MOV     CH, -1            ; signal rounding
  143.              CALL    RealTrunc         ; n = round (y)
  144.              POP     CX                ; get back lsw of x
  145.              PUSH    AX                ; save n
  146.              PUSH    CX                ; save lsw of x
  147.              CALL    RealFloat         ; convert n to floating-point
  148.              POP     CX                ; get back lsw of x
  149.              CALL    RealReduceMP      ; compute t = x - n*c1 - n*c2
  150.  
  151.              IFDEF   EXTENSIONS
  152.  
  153.              CMP     BP, OFFSET PowECst; function = e^x ?
  154.              JE      $approx           ; yes, compute approximation with g = t
  155. $power_2:    MOV     CX, CS:[BP+18]    ; load
  156.              MOV     SI, CS:[BP+20]    ;  constant
  157.              MOV     DI, CS:[BP+22]    ;   1/lg(2)
  158.              CALL    RealMulNoChk      ; compute g := t/lg(2) (DI:SI:CX <> 0)
  159.  
  160.              ENDIF
  161.  
  162. $approx:     PUSH    DX                ; save
  163.              PUSH    BX                ;  g on
  164.              PUSH    AX                ;   stack
  165.              CALL    RealSqr           ; compute z := g^2
  166.              POP     CX                ; get
  167.              POP     SI                ;  g from
  168.              POP     DI                ;   stack
  169.              PUSH    DX                ; save
  170.              PUSH    BX                ;  z on
  171.              PUSH    AX                ;   stack
  172.              PUSH    DI                ; save
  173.              PUSH    SI                ;  g on
  174.              PUSH    CX                ;   stack
  175.              MOV     CX, 01A85h        ; load real
  176.              MOV     SI, 09690h        ;  constant
  177.              MOV     DI, 0201Bh        ;   2.001347076967e1
  178.              CALL    RealMulFNoChk     ; multiply by z  (DI:SI:CX <> 0)
  179.              MOV     CX, 0388Ah        ; load
  180.              MOV     SI, 0C182h        ;  real constant
  181.              MOV     DI, 05233h        ;   8.408086858319e2
  182.              CALL    RealAdd           ; make sum
  183.              POP     CX                ; get
  184.              POP     SI                ;  g from
  185.              POP     DI                ;   stack
  186.              CALL    RealMulF          ; compute GPZ
  187.              POP     CX                ; get
  188.              POP     SI                ;  z from
  189.              POP     DI                ;   stack
  190.              PUSH    DX                ; save
  191.              PUSH    BX                ;  GPZ
  192.              PUSH    AX                ;   on stack
  193.              MOV     AX, 00088h        ; load
  194.              MOV     BX, 066A5h        ;  real constant
  195.              MOV     DX, 03429h        ;   1.801617224813e2
  196.              PUSH    DI                ; save
  197.              PUSH    SI                ;  z on
  198.              PUSH    CX                ;   stack
  199.              CALL    RealAdd           ; add z
  200.              POP     CX                ; get
  201.              POP     SI                ;  z from
  202.              POP     DI                ;   stack
  203.              CALL    RealMulF          ; multiply result by z
  204.              MOV     CX, 0388Bh        ; load
  205.              MOV     SI, 0C182h        ;  real constant
  206.              MOV     DI, 05233h        ;   1.681617371664e3
  207.              CALL    RealAdd           ; compute QZ
  208.              POP     CX                ; get
  209.              POP     SI                ;  GPZ from
  210.              POP     DI                ;   stack
  211.              PUSH    DI                ; save
  212.              PUSH    SI                ;  GPZ on
  213.              PUSH    CX                ;   stack
  214.              CALL    RealSub           ; compute QZ - GPZ
  215.              POP     CX                ; get
  216.              POP     SI                ;  GPZ from
  217.              POP     DI                ;   stack
  218.              CALL    RealDivRev        ; compute GPZ / (QZ-GPZ)
  219.              MOV     CX, 80h           ; load
  220.              XOR     SI, SI            ;  real constant
  221.              MOV     DI, SI            ;   0.5
  222.              CALL    RealAdd           ; compute 0.5 + GPZ / (QZ-GPZ)
  223.              POP     CX                ; get n
  224.              INC     CX                ; compute n+1
  225.              ADD     AL, CL            ; adjust exponent by adding n+1
  226.              POP     BP                ; restore caller's frame pointer
  227.              JZ      $exp_overfl       ; overflow error if exponent overflows
  228.              RET                       ; done
  229.  
  230.              IFDEF   NOOVERFLOW
  231. $exp_overfl: XOR     CH, CH            ; make result positive
  232.              JMP     ROverflow         ; error 205, floating point overflow
  233.              ELSE
  234. $exp_overfl: JMP     ROverflow         ; error 205, floating point overflow
  235.              ENDIF
  236.  
  237. PowECst      DW      00080h, 017F7h, 0B172h ; ln(2)-eps = 6.931471803691e-1
  238.              DW      00060h, 079ACh, 0D1CFh ;       eps = 1.908214929385e-10
  239.              DW      05C81h, 03B29h, 038AAh ; 1/ln(2)   = 1.442695040889e+0
  240.  
  241. RExp         ENDP
  242.  
  243.              ALIGN   4
  244.  
  245. CODE         ENDS
  246.  
  247.              END
  248.