home *** CD-ROM | disk | FTP | other *** search
/ EuroCD 3 / EuroCD 3.iso / Programming / vbcc / machines / amiga68k / libsrc / math / math_040 / sinh.s < prev    next >
Encoding:
Text File  |  1998-06-24  |  5.3 KB  |  120 lines
  1. *
  2. *   $VER: sinh.s 33.1 (22.1.97)
  3. *
  4. *   hyperbolic sine
  5. *
  6. *   Version history:
  7. *
  8. *   33.1    22.1.97 (c) Motorola
  9. *
  10. *           - snipped from M68060SP sources
  11. *
  12.  
  13.     machine 68040
  14.     fpu     1
  15.  
  16.     XDEF    _sinh
  17.     XDEF    @sinh
  18.  
  19.     XREF    @exp
  20.     XREF    @expm1
  21.  
  22. *************************************************************************
  23. * sin():  computes the hyperbolic sine of a normalized input            *
  24. *                                                                       *
  25. * INPUT *************************************************************** *
  26. *       fp0 = extended precision input                                  *
  27. *                                                                       *
  28. * OUTPUT ************************************************************** *
  29. *       fp0 = sinh(X)                                                   *
  30. *                                                                       *
  31. * ACCURACY and MONOTONICITY ******************************************* *
  32. *       The returned result is within 3 ulps in 64 significant bit,     *
  33. *       i.e. within 0.5001 ulp to 53 bits if the result is subsequently *
  34. *       rounded to double precision. The result is provably monotonic   *
  35. *       in double precision.                                            *
  36. *                                                                       *
  37. * ALGORITHM *********************************************************** *
  38. *                                                                       *
  39. *       SINH                                                            *
  40. *       1. If |X| > 16380 log2, go to 3.                                *
  41. *                                                                       *
  42. *       2. (|X| <= 16380 log2) Sinh(X) is obtained by the formula       *
  43. *               y = |X|, sgn = sign(X), and z = expm1(Y),               *
  44. *               sinh(X) = sgn*(1/2)*( z + z/(1+z) ).                    *
  45. *          Exit.                                                        *
  46. *                                                                       *
  47. *       3. If |X| > 16480 log2, go to 5.                                *
  48. *                                                                       *
  49. *       4. (16380 log2 < |X| <= 16480 log2)                             *
  50. *               sinh(X) = sign(X) * exp(|X|)/2.                         *
  51. *          However, invoking exp(|X|) may cause premature overflow.     *
  52. *          Thus, we calculate sinh(X) as follows:                       *
  53. *             Y       := |X|                                            *
  54. *             sgn     := sign(X)                                        *
  55. *             sgnFact := sgn * 2**(16380)                               *
  56. *             Y'      := Y - 16381 log2                                 *
  57. *             sinh(X) := sgnFact * exp(Y').                             *
  58. *          Exit.                                                        *
  59. *                                                                       *
  60. *       5. (|X| > 16480 log2) sinh(X) must overflow. Return             *
  61. *          sign(X)*Huge*Huge to generate overflow and an infinity with  *
  62. *          the appropriate sign. Huge is the largest finite number in   *
  63. *          extended format. Exit.                                       *
  64. *                                                                       *
  65. *************************************************************************
  66.  
  67. T1      dc.l            $40C62D38,$D3D64634     ; 16381 LOG2 LEAD
  68. T2      dc.l            $3D6F90AE,$B1E75CC7     ; 16381 LOG2 TRAIL
  69.  
  70. _sinh
  71.         fmove.d         (4,sp),fp0
  72. @sinh
  73.         fmove.x         fp0,-(sp)
  74.         move.l          (sp),d1
  75.         move.w          (4,sp),d1
  76.         move.l          d1,a1                   ; save (compacted) operand
  77.         lea             (12,sp),sp
  78.         and.l           #$7FFFFFFF,d1
  79.         cmp.l           #$400CB167,d1
  80.         bgt.b           .SINHBIG
  81.  
  82. ;--THIS IS THE USUAL CASE, |X| < 16380 LOG2
  83. ;--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) )
  84.  
  85.         fabs.x          fp0
  86.  
  87.         move.l          a1,-(sp)                ; {a1}
  88.         jsr             @expm1                  ; FP0 IS Z = EXPM1(Y)
  89.         move.l          (sp)+,a1                ; {a1}
  90.  
  91.         fmove.x         fp0,fp1
  92.         fadd.s          #$3F800000,fp1          ; 1+Z
  93.         fmove.x         fp0,-(sp)
  94.         fdiv.x          fp1,fp0                 ; Z/(1+Z)
  95.         move.l          a1,d1
  96.         and.l           #$80000000,d1
  97.         or.l            #$3F000000,d1
  98.         fadd.x          (sp)+,fp0
  99.         move.l          d1,-(sp)
  100.         fmul.s          (sp)+,fp0               ; last fp inst - possible exceptions set
  101.         rts
  102.  
  103. .SINHBIG
  104.         cmp.l           #$400CB2B3,d1
  105.         bgt             .t_ovfl
  106.  
  107.         fabs.x          fp0
  108.         fsub.d          (T1,pc),fp0             ; (|X|-16381LOG2_LEAD)
  109.         clr.l           -(sp)
  110.         move.l          #$80000000,-(sp)
  111.         move.l          a1,d1
  112.         and.l           #$80000000,d1
  113.         or.l            #$7FFB0000,d1
  114.         move.l          d1,-(sp)                ; EXTENDED FMT
  115.         fsub.d          (T2,pc),fp0             ; |X| - 16381 LOG2, ACCURATE
  116.         jsr             @exp
  117.         fmul.x          (sp)+,fp0               ; possible exception
  118. .t_ovfl
  119.         rts
  120.