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