PROGRAM test_inv_trig
!  1. a + b = asin [ sin(a).cos(b) + cos(a).sin(b) ]
!  2. a + b = acos [ cos(a).cos(b) - sin(a).sin(b) ]
!  3. a + b = atan [ tan(a) + tan(b) ] / [ 1 - tan(a).tan(b) ]

USE quadruple_precision
IMPLICIT NONE

TYPE (quad)          :: a, b, lhs, rhs, diff, sina, sinb, cosa, cosb,  &
                        tana, tanb, qone
REAL (r10)           :: half = 0.5_r10, small = EPSILON(half)
INTEGER, ALLOCATABLE :: seed(:)
INTEGER              :: k, i

qone = quad(1.0_r10, 0.0_r10)

!     Set the random number seed.

CALL RANDOM_SEED(size=k)
ALLOCATE (seed(k))
CALL RANDOM_SEED(get=seed)
WRITE(*, *)'Old random number seeds: ', seed

WRITE(*, '(1x, a, i4, a)') 'Enter ', k, ' integers as random number seeds: '
READ(*, *) seed
CALL RANDOM_SEED(put=seed)

DO i = 1, 10
  CALL RANDOM_NUMBER(a%hi)
  a%hi = (a%hi - half) / a%hi
  a%lo = a%hi * small
  CALL RANDOM_NUMBER(b%hi)
  b%hi = (b%hi - half) / b%hi
  b%lo = b%hi * small
                             ! asin [ sin(a).cos(b) + cos(a).sin(b) ]
                             ! Range (-pi/2, +pi/2)
  diff = (a + b)
  CALL longmodr(diff, pi, k, lhs)
  IF (k /= 2*(k/2)) lhs = -lhs
  sina = SIN(a)
  sinb = SIN(b)
  cosa = COS(a)
  cosb = COS(b)
  rhs = longasin( sina*cosb + cosa*sinb )
  diff = lhs - rhs
  WRITE(*, '(" lhs =", g13.5, "  Diff. =", g12.4)') lhs%hi, diff%hi
                             ! acos [ cos(a).cos(b) - sin(a).sin(b) ]
                             ! Range (0, pi)
  diff = (a + b)
  CALL longmodr(diff, pi, k, lhs)
  IF (k /= 2*(k/2)) THEN
    IF (lhs%hi > 0._r10) THEN
      lhs = pi - lhs
    ELSE
      lhs = pi + lhs
    END IF
  END IF
  IF (lhs%hi < 0._r10) lhs = -lhs
  rhs = longacos( cosa*cosb - sina*sinb )
  diff = lhs - rhs
  WRITE(*, '(" lhs =", g13.5, "  Diff. =", g12.4)') lhs%hi, diff%hi
                             ! atan [ tan(a) + tan(b) ] / [ 1 - tan(a).tan(b) ]
                             ! Range (-pi/2, +pi/2)
  diff = (a + b)
  CALL longmodr(diff, pi, k, lhs)
  tana = TAN(a)
  tanb = TAN(b)
  rhs = longatan( (tana + tanb) / (qone - tana*tanb) )
  diff = lhs - rhs
  WRITE(*, '(" lhs =", g13.5, "  Diff. =", g12.4)') lhs%hi, diff%hi
END DO

STOP
END PROGRAM test_inv_trig