PROGRAM series2dbl
!
! Author: Br. David Carlson
!
! Date: December 31, 1999
!
! Revised: May 2, 2004
! March 9, 2016 to use double precision reals
! March 11, 2018 to require 16 decimal digits of precision and to
! report information about the type of reals used.
!
! This program prints out partial sums for a series, stopping when no further
! change is found by adding on another term (or when a set number of terms
! have been added -- to prevent an infinite loop). The particular series
! used is the well-known series for sin(x). Double precision reals are used
! and give about 16 digits of precision. Note the request for KIND=DBL where
! DBL is a constant that specifies at least 16 digits of precision and an
! exponent range sufficient to handle 10 to the 37 power. Besides outputing
! the partial sums, this program also prints information about the type of
! real numbers used.
IMPLICIT NONE
INTEGER, PARAMETER::DBL = SELECTED_REAL_KIND(p=16, r=37)
REAL(KIND=DBL)::x, xsquared, term, sum, sumold
INTEGER::n, count
WRITE (*, 20)
20 FORMAT (1X, 'Program to show partial sums of the power series for sin(x):')
WRITE (*, *) 'Enter an angle in radians: '
READ (*, *) x
sum = 0.0D0 !*** The D0 means that we have exponent 0 (like E0), but also says that this is a double precision real.
!*** Use this D instead of E anytime you want a double precision real number.
count = 0
n = 1
term = x
xsquared = x * x
WRITE (*, 50)
50 FORMAT (' ', T5, 'count', T17, 'term', T50, 'partial sum')
DO
count = count + 1
sumold = sum
sum = sum + term
WRITE (*, 100) count, term, sum
IF ((sum == sumold) .OR. (count == 200)) THEN
EXIT
END IF
term = -term * xsquared / ((n + 1) * (n + 2))
n = n + 2
END DO
100 FORMAT (' ', T5, I3, T12, D24.16, T45, D24.16) !*** Note the special D format (similar to E format) for
!*** double precision reals. Display 16 decimal places so
!*** we can see all those digits that double precision gives.
WRITE (*, 200) SIN(x)
200 FORMAT (/, 'The built-in sine function gives: ', D25.16)
! We write out information about the double precision reals that we used here.
! This would not normally be left in the program. It is here to show you what we
! know about these numbers.
WRITE (*, 300) KIND(x), PRECISION(x), RANGE(x)
300 FORMAT (/, 'kind = ', I2, ', precision = ', I2, ', range = ', I5)
END PROGRAM