program TEST_QUAD
!     Three tests of the basic quadruple precision module.
!     1.       (A + B)^2          = (A - B)^2 + 4.A.B
!     2.  (A^2 - B^2) / (A - B)   = A + B
!     3.  SQRT(A^2 + 2.A.B + B^2) = A + B

use QUADRUPLE_PRECISION
implicit none

type (QUAD)                        :: A, B, LHS, RHS, DIFF
real (kind=DP), parameter          :: HALF = 0.5_DP, SMALL = epsilon(HALF)
integer, allocatable, dimension(:) :: SEED
integer                            :: K, I

!     Set the random number seed.

call random_seed(size=K)
allocate (SEED(K))
call random_seed(get=SEED)
write(unit=*, fmt=*) "Old random number seeds: ", SEED

write(unit=*, fmt="(a, i4, a)")  &
                  " Enter ", K, " integers as random number seeds: "
read(unit=*, fmt=*) 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
                             ! (A + B)^2 = (A - B)^2 + 4.A.B
  LHS = (A + B) * (A + B)
  RHS = A * B
  RHS = 4.0_DP * RHS
  RHS = (A - B) * (A - B) + RHS
  DIFF = LHS - RHS
  write(unit=*, fmt="(a, g13.5, a, g12.4)")  &
                    " lhs =", LHS%HI, "  Diff. =", DIFF%HI

                             ! (A^2 - B^2) / (A - B) = (A + B)
  LHS = (A*A - B*B) / (A - B)
  RHS = A + B
  DIFF = LHS - RHS
  write(unit=*, fmt="(a, g13.5, a, g12.4)")  &
                    " lhs =", LHS%HI, "  Diff. =", DIFF%HI

                             ! SQRT(A^2 + 2.A.B + B^2) = A + B
  LHS = A * B
  LHS = 2.0_DP * LHS
  LHS = sqrt(A*A + LHS + B*B)
  if (RHS%HI < 0.0_DP) then             ! Force +ve square root
    RHS = -RHS
  end if
  DIFF = LHS - RHS
  write(unit=*, fmt="(a, g13.5, a, g12.4)")  &
                    " lhs =", LHS%HI, "  Diff. =", DIFF%HI

end do

stop
end program TEST_QUAD