'*************************************************
'  Program to demonstrate the inverse hyperbolic
'  functions subroutines
'-------------------------------------------------
'Reference: BASIC Scientific Subroutines, Vol. II
'  By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981.
'  [BIBLI 01].
'*************************************************
defint i-n
defdbl a-h,o-z
cls
dim A(15),W(15)
F1$="  ##.#" : F2$="  ###.##########"
ligne=1
print
print "   X      ARCSINH(X)      ARCCOSH(X)      ARCTANH(X) "
print " ----------------------------------------------------"
for x=-3# to 3.2# step 0.2#
  print using F1$; X;
'get arcsinh(x)
  gosub 1000
  print using F2$; Y;
'get arccosh(x)
  gosub 1500
  print using F2$; Y;
'get arctanh(x)
  gosub 2000
  if ABS(Y)<1e12 then
    print using F2$; Y
  else
    print "   "; Y
  end if
  if ligne=20 then
    ligne=0
    input "",R$
    cls
    print:print
  end if
  ligne=ligne+1
next X
print
end

500 'Modified cordic logarithm subroutine
' This subroutine takes an input value and returns Y=LN(X)
' X may be any positive real value
' Get coefficients
  gosub 900
' If X<=0 then an error exists, return
  if X<=0 then return
  K=0
' Save X
  X1=X
' Reduce the range of X
510 if X<E then goto 520
' Divide out a power of E
    K=K+1
    X=X/E
  goto 510
' Test if X>=1, if so go to next step
' Otherwise, bring X to >1
520 if X>=1# then goto 530
  K=K-1
  X=X*E
  goto 520
' Determine the weighting coefficients, W(I)
530 gosub 800
' Calculate residual factor based on Z
' Want LN(Z), where Z is near unity
  Z=Z-1
  Z=Z*(1#-(Z/2#)*(1#+(Z/3#)*(1#-Z/4#)))
' Assemble results
  A=0.5#
  for I=1 to N
    Z=Z+W(I)*A
    A=A/2#
  next I
' Z is now the mantissa, K the characteristic
  Y=K+Z
' Restore X
  X=X1
return

800 'Weight determination subroutine
  Z=X
  for I=1 to N
    W(I)=0#
    if Z>A(I) then W(I)=1#
    if W(I)=1# then Z=Z/A(I)
  next I
return

900 'Exponential coefficients subroutine
  N=15
  E=2.718281828459045
  A(1)=1.648721270700128
  A(2)=1.284025416687742
  A(3)=1.133148453066826
  A(4)=1.064494458917859
  A(5)=1.031743407499103
  A(6)=1.015747708586686
  A(7)=1.007843097206448
  A(8)=1.003913889338348
  A(9)=1.001955033591003
  A(10)=1.000977039492417
  A(11)=1.000488400478694
  A(12)=1.000244170429748
  A(13)=1.000122077763384
  A(14)=1.000061037018933
  A(15)=1.000030518043791
return

'----------------------------------------------
'    Inverse hyperbolic sine subroutine
' ---------------------------------------------
' This subroutine calculates the inverse of
' the hyperbolic sine using the modified cordic
' natural logarithm subroutine 500.
' Input: X - Output: Y = ARCSINH(X)
' Formula: ARCSINH(X)=LN(X+SQRT(X*X+1))
'----------------------------------------------
' Test for zero argument
1000 if X<>0# then goto 1100
  Y=0#
  return
1100 'Save X
  X2=X
  X=ABS(X)
  X=X+SQR(X*X+1#)
' Get LN(X)
  gosub 500
' Insert sign
  Y=(X2/ABS(X2))*Y
' Restore X
  X=X2
return

1500 'Inverse hyperbolic cosine subroutine
'     ARCCOSH(X)=LN(X+SQRT(X*X-1))
' Test for argument less than or equal to unity
  if X>1# then goto 1550
  Y=0
  return
1550 'Save X
  X2=X
  X=ABS(X)
  X=X+SQR(X-1#)*SQR(X+1#)
  gosub 500 'Y=LN(X)
' Restore X
  X=X2
return

2000 'Inverse hyperbolic tangent subroutine
'     ARCTANH(X)=1/2*LN(1+x/1-x)
' Test for X>= +/- 1
  if ABS(X)<=0.999999 then goto 2100
' Y is BIG! (here +/- 1E18)
  Y=(X/ABS(X))*1000000#*1000000#*1000000
  return
2100 'Test for zero argument
  if X<>0# then goto 2150
  Y=0
  return
2150 'Save X
  X2=X
  X=(1#+X)/(1#-X)
  gosub 500 'Y=LN(X)
' Restore X
  X=X2
return

' End of file invhyper.bas