subroutine polyan ( c, cm, n ) c*********************************************************************72 c c polyan obtains information about the location c of the roots of a polynomial by using c bound, radius and hrwtzr. c c is an n element array containing the coefficients c normalized so that the leading coefficient (which c is not included in c) is +1.0. c cm is a working array the same size as c. c n = degree of polynomial. c dimension c(n), cm(n) logical hrwtzr c test for zero root if ( c(n) .eq. 0.0 ) go to 50 c coefficients for reciprocal polynomial are put in cm. cm(n) = 1. / c(n) nm1 = n - 1 do i = 1, nm1 ni = n - i cm(i) = cm(n) * c(ni) end do rout = radius ( c, n ) rin = 1.0 / radius ( cm, n ) write ( 6, 201 ) rin, rout 201 format ( 40h roots are in an annulus of inner radius, & e10.3, 17h and outer radius, e10.3 ) rpu = bound ( c, n ) if ( rpu .ne. 0.0 ) go to 10 write ( 6, 202 ) 202 format ( 33h there are no real positive roots ) go to 20 10 rpl = 1. / bound ( cm, n ) write ( 6, 203 ) rpl, rpu 203 format & ( 40h the positive roots(if any) are between, & e10.3, 4h and, e10.3 ) c coefficients for negative reciprocal are put in cm. 20 do i = 1, n, 2 cm(i) = - cm(i) end do rnu = bound ( cm, n ) if ( rnu .ne. 0.0 ) go to 30 write ( 6, 204 ) 204 format & ( 33h there are no negative real roots ) go to 40 c coefficients for negative roots are put in cm. 30 x = -1.0 do i = 1, n cm(i) = x * c(i) x = - x end do rnu = -1.0 / rnu rnl = - bound ( cm, n ) write ( 6, 205 ) rnu, rnl 205 format & ( 44h the real negative roots(if any) are between, & e10.3, 4h and, e10.3 ) 40 if ( hrwtzr ( c, n ) ) write ( 6, 206 ) 206 format & ( 44h there are no roots with positive real parts ) if ( hrwtzr ( cm, n ) ) write ( 6, 207 ) 207 format & ( 44h there are no roots with negative real parts ) return 50 write ( 6, 208 ) 208 format & ( 41h polynomial has a zero root-reduce degree ) return end function radius ( c, n ) c*********************************************************************72 c c radius returns an upper limit for the modulus c of the roots of an n degree polynomial. c dimension c(n) radius = abs ( c(1) ) do i = 2, n if ( abs ( c(i) ) .gt. radius ) radius = abs ( c(i) ) end do radius = 1.0 + radius return end function bound ( c, n ) c*********************************************************************72 c c bound returns an upper limit for the c positive real roots of an n degree polynomial. c dimension c(n) m = 0 bound = 0.0 do i = 1, n if ( m .le. 0 ) then if ( c(i) .lt. 0.0 ) m = i if ( c(i) .lt. bound ) bound = c(i) end if end do if ( m .eq. 0 ) return bound = 1.0 + ( -bound )** ( 1.0 / float ( m ) ) return end logical function hrwtzr ( c, n ) c*********************************************************************72 c c hrwtzr returns .true. if all the roots have c negative real parts, otherwise .false. is returned. c if a real part is zero, then .false. is returned. c dimension c(n) hrwtzr = .false. c%% c%% modifications as suggested by driessen and hunt, c%% cacm volume 16, number 9, page 579, september 1973. c%% if ( c(1) .le. 0.0 .or. c(n) .le. 0. ) return c1 = c(1) m = n - 1 do i = 2, m do k = i, m, 2 c(k) = c(k) - c(k+1) / c1 end do c1 = c(i) / c1 if ( c1 .le. 0.0 ) return end do hrwtzr = .true. return end