From b6690c482697f1e5d759e00eecaee6ffc9e82a4e Mon Sep 17 00:00:00 2001
From: Stephen Heumann <stephenheumann@gmail.com>
Date: Tue, 30 Nov 2021 19:15:21 -0600
Subject: [PATCH] Implement acosh().

This is basically the implementation recommended in Apple Numerics Manual Ch. 9, except that there is an added case for large values that would otherwise cause an overflow.
---
 math2.asm    | 108 +++++++++++++++++++++++++++++++++++++++++++++++++++
 math2.macros |  37 ++++++++++++++++++
 2 files changed, 145 insertions(+)

diff --git a/math2.asm b/math2.asm
index 0c45612..2fd2ff7 100644
--- a/math2.asm
+++ b/math2.asm
@@ -220,6 +220,114 @@ lb1      sta   mask
          creturn 2:mask
          end
 
+****************************************************************
+*
+*  double acosh(double x);
+*
+*  Returns the inverse hyperbolic cosine of x.
+*
+****************************************************************
+*
+acosh    start
+acoshf   entry
+acoshl   entry
+         using MathCommon2
+         
+         csubroutine (10:x),0
+         
+         phb
+         phk
+         plb
+         
+         pha                            save env & set to default
+         tsc
+         inc   a
+         pea   0
+         pha
+         FPROCENTRY
+         
+         lda   x                        y = sqrt(x-1)
+         sta   y
+         lda   x+2
+         sta   y+2
+         lda   x+4
+         sta   y+4
+         lda   x+6
+         sta   y+6
+         lda   x+8
+         sta   y+8
+         ph4   #one
+         ph4   #y
+         FSUBI
+         ph4   #y
+         FSQRTX
+
+         lda   x                        t1 = sqrt(x+1)
+         sta   t1
+         lda   x+2
+         sta   t1+2
+         lda   x+4
+         sta   t1+4
+         lda   x+6
+         sta   t1+6
+         lda   x+8
+         sta   t1+8
+         ph4   #one
+         ph4   #t1
+         FADDI
+         ph4   #t1
+         FSQRTX
+
+         ph4   #y                       t1 = ln(1+y*(y+t1))
+         ph4   #t1
+         FADDX
+         ph4   #y
+         ph4   #t1
+         FMULX
+         ph4   #t1
+         FLN1X
+         
+         lda   t1+8                     if t1 = +inf
+         cmp   #32767
+         bne   ret
+         lda   t1+6
+         asl   a
+         ora   t1+4
+         ora   t1+2
+         ora   t1
+         bne   ret
+         
+         pea   0                          clear exceptions
+         FSETENV
+         lda   x                          t1 = ln(x) + ln(2)
+         sta   t1
+         lda   x+2
+         sta   t1+2
+         lda   x+4
+         sta   t1+4
+         lda   x+6
+         sta   t1+6
+         lda   x+8
+         sta   t1+8
+         ph4   #t1
+         FLNX
+         ph4   #ln2
+         ph4   #t1
+         FADDX
+
+ret      FPROCEXIT                      restore env & raise any new exceptions
+         plb
+         lda   #t1                      return t1
+         sta   x
+         lda   #^t1
+         sta   x+2
+         creturn 4:x
+
+y        ds    10                       temporary variable
+one      dc    i'1'                     constants
+ln2      dc    e'0.69314718055994530942'
+         end
+
 ****************************************************************
 *
 *  double cbrt(double x);
diff --git a/math2.macros b/math2.macros
index 8c11b43..c666b3f 100644
--- a/math2.macros
+++ b/math2.macros
@@ -580,3 +580,40 @@
  LDX #$090A
  JSL $E10000
  MEND
+ MACRO
+&LAB FADDX
+&LAB PEA $0000
+ LDX #$090A
+ JSL $E10000
+ MEND
+ MACRO
+&LAB FADDI
+&LAB PEA $0400
+ LDX #$090A
+ JSL $E10000
+ MEND
+ MACRO
+&LAB FSUBI
+&LAB PEA $0402
+ LDX #$090A
+ JSL $E10000
+ MEND
+ MACRO
+&LAB FMULX
+&LAB PEA $0004
+ LDX #$090A
+ JSL $E10000
+ MEND
+ MACRO
+&LAB FSQRTX
+&LAB PEA $0012
+ LDX #$090A
+ JSL $E10000
+ MEND
+ MACRO
+&LAB FLNX
+&LAB PEA $0000
+ LDX #$0B0A
+ JSL $E10000
+ MEND
+