home *** CD-ROM | disk | FTP | other *** search
/ EuroCD 3 / EuroCD 3.iso / Programming / vbcc / machines / amiga68k / libsrc / math / math_040 / pow.s < prev    next >
Encoding:
Text File  |  1998-06-24  |  3.1 KB  |  117 lines
  1. *
  2. *   $VER: pow.s 33.1 (19.1.97)
  3. *
  4. *   calculates the result of a to the power of b.
  5. *
  6. *   Version history:
  7. *
  8. *   33.1    19.1.97 laukkanen
  9. *
  10. *           - created (using the fact that pow(a,b) == exp(b*log(a))
  11. *
  12. *   33.2    20.1.97 laukkanen
  13. *
  14. *           - added pow2 and pow10
  15. *
  16. *   33.3    25.1.97 laukkanen
  17. *
  18. *           - hopefully fixed all the special cases in pow()
  19. *           - further optimization is possible if the cases
  20. *             where y is integer are handled separately.
  21. *
  22. *
  23.  
  24.     machine 68040
  25.     fpu     1
  26.  
  27.     XDEF    _pow
  28.     XDEF    @pow
  29.     XDEF    _pow2
  30.     XDEF    @pow2
  31.     XDEF    _pow10
  32.     XDEF    @pow10
  33.  
  34.     XREF    @exp
  35.     XREF    @log
  36.  
  37.         cnop            0,4
  38. _pow
  39.         fmove.d         (4,sp),fp0
  40.         fmove.d         (12,sp),fp1
  41. @pow
  42.         fmove.s         fp0,-(sp)               ; Let's check if
  43.         move.w          (sp),d0                 ; x is negative or zero
  44.         bmi.s           .neg                    ; and branch accordingly
  45.         addq.l          #4,sp
  46.         and.w           #$7f80,d0
  47.         beq.s           .zero
  48. .pos
  49.         fmove.x         fp2,-(sp)               ; x**y is calculated
  50.         fmove.x         fp1,fp2                 ; with exp(y*ln(x))
  51.         jsr             @log
  52.         fmul.x          fp2,fp0
  53.         jsr             @exp
  54.         fmove.x         (sp)+,fp2
  55.         rts
  56. .zero                                           ; x is zero
  57.         fmove.s         fp1,-(sp)
  58.         move.w          (sp),d0                 ; if y is zero
  59.         addq.l          #4,sp                   ;   return 1;
  60.         and.w           #$7f80,d0               ; else return 0;
  61.         beq             .zero1                  ; My calculator returns
  62.         rts                                     ; error for 0**0 but
  63. .zero1                                          ; "cephes" math package
  64.         fmove.s         #1,fp0                  ; did so. OTOH SAS/C
  65.         rts                                     ; scm881.lib just returns
  66.                                                 ; undefined number.
  67. .neg
  68.         fabs.x          fp0
  69.         addq.l          #4,sp
  70.         fmove.x         fp2,-(sp)
  71.         moveq           #0,d1
  72.         fmove.x         fp1,fp2
  73.         fabs.x          fp1
  74.         fmove.l         fp1,-(sp)
  75.         move.l          (sp)+,d0
  76.         and.b           #$01,d0
  77.         beq             .even
  78.         moveq           #1,d1
  79. .even
  80.         move.l          d1,-(sp)
  81.         jsr             @log
  82.         fmul.x          fp2,fp0
  83.         jsr             @exp
  84.         move.l          (sp)+,d1
  85.         beq             .even2
  86.         fneg.x          fp0
  87. .even2
  88.         fmove.x         (sp)+,fp2
  89.         rts
  90.  
  91.         cnop            0,4
  92. LOGTABLE
  93.  
  94.  
  95.  
  96. _pow2
  97.         fmove.d         (4,sp),fp0
  98. @pow2
  99.         fmul.d          (.ln2,pc),fp0
  100.         jsr             @exp
  101.         rts
  102.  
  103.         cnop            0,4
  104.  
  105. .ln2    dc.d            0.69314718055995
  106.  
  107. _pow10
  108.         fmove.d         (4,sp),fp0
  109. @pow10
  110.         fmul.d          (.ln10,pc),fp0
  111.         jsr             @exp
  112.         rts
  113.  
  114.         cnop            0,4
  115.  
  116. .ln10   dc.d            2.302585092994
  117.