Math-Polynomial-Solve

 view release on metacpan or  search on metacpan

lib/Math/Polynomial/Solve.pm  view on Meta::CPAN

	# With the coefficents in ascending order,
	# pretend it was always that way for the next
	# function calls.
	#
	my $temp_ascending_flag = $ascending_flag;
	$ascending_flag = 1;

	if ($option{hessenberg} or $#coefficients > 4)
	{
		#
		# QR iterations from the matrix.
		#
		@x = hqr_eigen_hessenberg(
			balance_matrix(build_companion(@coefficients))
			);
	}
	elsif ($#coefficients == 4)
	{
		@x = quartic_roots(@coefficients);
	}
	elsif ($#coefficients == 3)

lib/Math/Polynomial/Solve.pm  view on Meta::CPAN

					# Complex or twin pair.
					#
					push @roots, $x + $p - $y * i;
					push @roots, $x + $p + $y * i;
				}

				$n -= 2;
				next ROOT;
			}

			croak "Too many iterations ($its) at n=$n\n" if ($its >= $iteration{hessenberg});

			if ($its && $its % 10 == 0)
			{
				#
				# Form exceptional shift.
				#
				### Exceptional shift at: $its
				#

				$t += $x;

lib/Math/Polynomial/Solve.pm  view on Meta::CPAN

		else
		{
			#
			# We've divided the roots between two ranges. Do it
			# again until each range has a single root in it.
			#
			push @boundaries, sturm_bisection($chain_ref, $from, $mid);
			push @boundaries, sturm_bisection($chain_ref, $mid, $to);
			last ROOT;
		}
		croak "Too many iterations ($its) at mid=$mid\n" if ($its >= $iteration{sturm_bisection});
		$its++;
	}
	return @boundaries;
}


=head3 sturm_bisection_roots()

Return the I<real> roots counted by L</sturm_real_root_range_count()>.
Uses L</sturm_bisection()> to bracket the roots of the polynomial,

lib/Math/Polynomial/Solve.pm  view on Meta::CPAN

			#
			my $dx = $n/($g + $f);

			$x -= $dx;
			if (abs($dx) <= $tolerance{laguerre})
			{
				push @roots, $x;
				last ROOT;
			}

			croak "Too many iterations ($its) at dx=$dx\n" if ($its >= $iteration{laguerre});
			$its++;
		}

		### root found at iteration $its
		#### $x
	}

	return @roots;
}

lib/Math/Polynomial/Solve.pm  view on Meta::CPAN

				push @roots, $x;
				last ROOT;
			}

			#
			#### At Iteration: $its
			#### x: $x
			#### f(x): $y
			#### f'(x): $dy
			#
			croak "Too many iterations ($its) at dx=$dx\n" if ($its >= $iteration{newtonraphson});
			$its++;
		}

		### root found at iteration $its
		#### $x
	}

	return @roots;
}

=head3 poly_iteration()

Sets the limit to the number of iterations that a solving method may go
through before giving up trying to find a root. Each method of root-finding
used by L</poly_roots()>, L</sturm_bisection_roots()>, and L</laguerre()>
has its own iteration limit, which may be found, like L</poly_option()>,
simply by looking at the return value of poly_iteration().

    #
    # Get all of the current iteration limits.
    #
    my %its_limits = poly_iteration();

reference/qralg/qralg.f  view on Meta::CPAN

      integer cnt(n)
      integer rtn_code
Comments:
*     n       : degree of the monic polynomial.
*     c       : non-principal coefficients of the polynomial.
*               i-th degree coefficients is stored in c(i).
*     macheps : double precision machine epsilon.
*     zr,zi   : Re and Im part of output roots.
*     detil   : The accuracy hint.
*     a       : work matrix.
*     cnt     : work area for counting the qr-iterations.
*     rtn_code: return code from hqr_eigen_hessenberg.
*end comments;
*================
      integer i
      integer iter
      double precision afnorm
*
      if(n.le.0)then
        print *, 'qr_alg_solver: wrong arguments. (n<=0)'
        rtn_code=3

reference/qralg/qralg.f  view on Meta::CPAN

*
* Build the companion matrix A.
      call build_companion(n,a,c)
*
* Balancing the A itself.
      call balance_companion(n,a)
*
* Compute the Frobenius norm of the balanced companion matrix A.
      call frobenius_norm_companion(n,a,afnorm)
*
* QR iterations from A.
      call hqr_eigen_hessenberg(n,a,macheps,  zr,zi,cnt,rtn_code)
      if(rtn_code.ne.0)then
        print*, 'in calling hqr_eigen_hessenberg abnormal condition'
        print*, 'rtn_code=',rtn_code
        if(rtn_code.eq.1) print *, 'matrix is completely zero.'
        if(rtn_code.eq.2) print *, 'QR iteration does not converged.'
        if(rtn_code.gt.3) print *, 'arguments violate the condition.'
        return
      endif
*

reference/qralg/qralg.f  view on Meta::CPAN

*       Numer. Math. 14, 219-231(1970).
*
************************************************************************
      subroutine hqr_eigen_hessenberg(n0,h,macheps,  wr,wi,cnt,rtn_code)
      implicit none
*Comment: Finds the eigenvalues of a real upper Hessenberg matrix,
*         H, stored in the array h(1:n0,1:n0), and stores the real
*         parts in the array wr(1:n0) and the imaginary parts in the
*         array wi(1:n0). macheps is the relative machine precision.
*         The procedure fails if any eigenvalue takes more than 
*         30 iterations.
*
      integer n0
      double precision h(n0,n0)
      double precision macheps
*
      double precision wr(n0),wi(n0)
      integer cnt(n0)
      integer rtn_code
*===
      integer i,j,k,l,m,na,its