real a, b, c, x1, x2, r1, r2 real*8 aa, bb, cc, xx1, xx2, rr1, rr2 integer i c c for real c write(6,*) ' *** using real ***' write(6,'(3x,a,11x,a,14x,a)') 'exact solution','calc. solution', $ 'relative error' a = 1.0 r1 = 1.0 do i = 1 , 10 r2 = r1/10.0**i b = - ( r1 + r2 ) c = r1 * r2 c x1 = ( -b + sqrt( b*b - 4.0*a*c ))/ ( 2.0*a ) x2 = ( -b - sqrt( b*b - 4.0*a*c ))/ ( 2.0*a ) c write(6,'(2e10.3,4e15.8)') $ r1, r2, x1, x2, (x1-r1)/r1, (x2-r2)/r2 end do c write(6,*) write(6,*) ' *** another algorihm ***' write(6,'(3x,a,11x,a,14x,a)') 'exact solution','calc. solution', $ 'relative error' a = 1.0 r1 = 1.0 do i = 1 , 10 r2 = r1/10.0**i b = - ( r1 + r2 ) c = r1 * r2 c if( b .ge. 0.0 ) then x1 = ( -b - sqrt( b*b - 4.0*a*c ))/ ( 2.0*a ) x2 = c / ( a*x1 ) else x1 = ( -b + sqrt( b*b - 4.0*a*c ))/ ( 2.0*a ) x2 = c / ( a*x1 ) end if c write(6,'(2e10.3,4e15.8)') $ r1, r2, x1, x2, (x1-r1)/r1, (x2-r2)/r2 end do c c for double precision c write(6,*) write(6,*) ' *** using double ***' write(6,'(3x,a,11x,a,14x,a)') 'exact solution','calc. solution', $ 'relative error' aa = 1.0d0 rr1 = 1.0d0 do i = 1 , 10 rr2 = rr1/10.0d0**i bb = - ( rr1 + rr2 ) cc = rr1 * rr2 c xx1 = ( -bb + sqrt( bb*bb - 4.0d0*aa*cc ))/ ( 2.0d0*aa ) xx2 = ( -bb - sqrt( bb*bb - 4.0d0*aa*cc ))/ ( 2.0d0*aa ) c write(6,'(2e10.3,4e15.8)') $ rr1, rr2, xx1, xx2, (xx1-rr1)/rr1, (xx2-rr2)/rr2 end do c end