;; -*- Mode:Lisp; Syntax:ANSI-Common-LISP; Coding:us-ascii-unix; fill-column:158 -*-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;; @file      use-poly.lisp
;; @author    Mitch Richling <https://www.mitchr.me>
;; @brief     Polynomials over complex, real, rational, and integers.@EOL
;; @std       Common Lisp
;; @see       tst-poly.lisp
;; @copyright
;;  @parblock
;;  Copyright (c) 1994,1997,1998,2004,2008,2012,2013,2014,2015, Mitchell Jay Richling <https://www.mitchr.me> All rights reserved.
;;
;;  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
;;
;;  1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
;;
;;  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
;;     and/or other materials provided with the distribution.
;;
;;  3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
;;     without specific prior written permission.
;;
;;  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
;;  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
;;  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
;;  OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
;;  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
;;  DAMAGE.
;;  @endparblock
;; @todo      mjr_poly_root-solve-search-deflate: Want generic way to specifiy search, refine, and xform methods.@EOL@EOL
;; @todo      Make some kind of convention for the type (int, gausian, float, complex, etc...) of polynomials a function works with.@EOL@EOL
;; @todo      Better way to find solutions that are in radical extension fields -- i.e. $a,b,c\in\mathbb{Q}$ but $a+b\sqrt{c}\not\in\mathbb{Q}$.@EOL@EOL
;; @todo      Better control/seporation/deliniation of Z[x], Q[x], & R[x] across functions.  Make this uniform.@EOL@EOL
;; @todo      Update unit tests for new functionality.@EOL@EOL
;; @todo      mjr_poly_decompose.@EOL@EOL
;; @todo      mjr_poly_factor-irreducible.@EOL@EOL
;; @todo      mjr_poly_root-solve-jenkins-traub.@EOL@EOL
;; @todo      mjr_poly_root-solve-continued-fraction.@EOL@EOL
;; @todo      mjr_poly_root-solve-descartes-bisection.@EOL@EOL
;; @todo      mjr_poly_root-separate-real.@EOL@EOL
;; @warning   Much of this code has become quite experimental!  Use the old version if you want to be carefull.@EOL@EOL
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defpackage :MJR_POLY
  (:USE :COMMON-LISP
        :MJR_NLEQ
        :MJR_NUMU
        :MJR_CMP
        :MJR_CHK
        :MJR_PRIME
        :MJR_VEC
        :MJR_PRNG
        :MJR_COMBE
        :MJR_INTRP
        :MJR_COMBC
        :MJR_VVEC
        :MJR_GPOLY
        :MJR_EPS)
  (:DOCUMENTATION "Brief: Univariate Polynomials over R or C.;")
  (:EXPORT #:mjr_poly_help
           ;; Autogenerated from #:mjr_gpoly
           #:mjr_poly_coeff #:mjr_poly_diff #:mjr_poly_- #:mjr_poly_+ #:mjr_poly_* #:mjr_poly_iexpt #:mjr_poly_scale
           #:mjr_poly_truncate #:mjr_poly_rem #:mjr_poly_mod #:mjr_poly_divides?  #:mjr_poly_gcd #:mjr_poly_degree #:mjr_poly_leading-coeff
           #:mjr_poly_constant-coeff #:mjr_poly_zerop #:mjr_poly_onep #:mjr_poly_constantp #:mjr_poly_simplify #:mjr_poly_eval #:mjr_poly_subst
           #:mjr_poly_density #:mjr_poly_index
           ;; Implemented here..

           #:mjr_poly_imul

           #:mjr_poly_count-sign-changes

           #:mjr_poly_print #:mjr_poly_code

           #:mjr_poly_intg

           #:mjr_poly_gcd-primitive

           ;; Polynomial metrics
           #:mjr_poly_height #:mjr_poly_length #:mjr_poly_mahler-measure
           #:mjr_poly_scount-descartes #:mjr_poly_scount-sturm #:mjr_poly_scount-fourier #:mjr_poly_scount-budan

           #:mjr_poly_root-structure #:mjr_poly_root-structure-print

           #:mjr_poly_root-count-distinct-real #:mjr_poly_root-count-distinct-interval #:mjr_poly_root-count-distinct-negative-zero-positive

           #:mjr_poly_root-bound-all #:mjr_poly_root-bound-positive #:mjr_poly_root-bound-all-cauchy

           #:mjr_poly_eval-poly-and-first-n-derivatives
           #:mjr_poly_2func
           #:mjr_poly_seq-eval

           #:mjr_poly_root-separate-real

           #:mjr_poly_root-solve-rational #:mjr_poly_root-solve-search-deflate #:mjr_poly_root-solve-low-degree

           #:mjr_poly_2square-free
           #:mjr_poly_2primitive
           #:mjr_poly_2monic

           #:mjr_poly_root-search-bsect #:mjr_poly_root-search-newton #:mjr_poly_root-search-laguerre #:mjr_poly_root-search-largest-real

           #:mjr_poly_root-separate-largest-real

           #:mjr_poly_make-from-roots

           ;; Make special polynomials
           #:mjr_poly_make-lagrange #:mjr_poly_make-chebyshev #:mjr_poly_make-legendre #:mjr_poly_make-bernstein #:mjr_poly_make-hermite
           #:mjr_poly_make-laguerre

           ;; Make sequences of polynomials
           #:mjr_poly_seq-make-lagrange #:mjr_poly_seq-make-chebyshev #:mjr_poly_seq-make-legendre
           #:mjr_poly_seq-make-bernstein #:mjr_poly_seq-make-hermite #:mjr_poly_seq-make-laguerre
           #:mjr_poly_seq-make-fourier #:mjr_poly_seq-make-sturm-canonical

           #:mjr_poly_deflate #:mjr_poly_shift #:mjr_poly_reflect #:mjr_poly_zap-zero-roots

           ;; EXPERIMENTAL!!!!
           #:mjr_poly_test-property
           #:mjr_poly_find-integer-factor
           #:mjr_poly_factor-over-integers

           ;; MORE EXPERIMENTAL!!!!
           #:mjr_poly_cubic-depress
           #:mjr_poly_tschirnhaus-3-2
           #:mjr_poly_resultant
           #:mjr_poly_discriminant-low-degree
           #:mjr_poly_discriminant-high-degree
           #:mjr_poly_discriminant
           #:mjr_poly_factor-square-free
           #:mjr_poly_factor-irreducible
           #:mjr_poly_content
           #:mjr_poly_primitive-part
           #:mjr_poly_rationalize
           ))

(in-package :MJR_POLY)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_help ()
  "Univariate, Dense Polynomial Library

Polynomials are represented as vectors containing the coefficients -- the position in the array implies the power of the monomial in the term to which the
coefficient belongs.  For example, x^2+4 would be represented as #(1 0 4).

The first coefficient should never be zero unless the polynomial is the zero polynomial.  No function in this package should ever produce a non-constant
polynomial with a zero in the first position; however, all functions will gladly accept such polynomials and silently fix them.

Note that while this vector representation is quite adequate for dense polynomials (i.e. polynomials that have mostly non-zero coefficients), this data
structure is inefficient in both space and time for sparse polynomials.

This library is a work in progress, but it has enough functionality to be useful.  Still, some care should be exercised when using it.

Tags: (G) means the function is provided by GPOLY and a (!) means the function is not yet implemented.  Function names have been abbreviated in this list by
replacing 'mjr_poly_' with 'MP_'.

  * Strings & Printing
  ** MP_print MP_code
  * Evaluation
  ** MP_eval(G) MP_eval-poly-and-first-n-derivatives MP_seq-eval
  * Arithmetic
  ** MP_+(G) MP_-(G) MP_*(G) MP_iexpt(G) MP_rem(G) MP_mod(G) MP_truncate(G) MP_gcd(G) MP_gcd-primitive MP_imul MP_divides?(G)
  * Coefficient access
  ** MP_coeff(G) MP_leading-coeff(G) MP_constant-coeff(G)
  * Recognizing special polynomials
  ** MP_onep(G) MP_constantp(G) MP_zerop(G)
  * Polynomial factorization
  ** MP_decompose(!) MP_factor-square-free MP_factor-irreducible(!)
  * Things we do with roots:
  **  MP_root-structure
  ** Solve polynomial (find all roots):
  *** MP_root-solve-integer MP_root-solve-rational MP_root-solve-low-degree MP_root-solve-search-deflate
  *** MP_root-solve-jenkins-traub(!) MP_root-solve-continued-fraction(!) MP_root-solve-descartes-bisection(!)
  ** Finding roots:
  *** MP_root-search-bsect MP_root-search-laguerre MP_root-search-newton MP_root-search-largest-real
  *** MP_root-search-max-magnitude
  ** Bounding roots:
  *** MP_root-bound-all MP_root-bound-positive
  ** Separating roots (a root is 'separated' by finding an interval that contains it but no other roots):
  *** MP_root-separate-largest-real MP_root-separate-real(!)
  ** Counting distinct roots
  *** MP_root-count-distinct-real MP_root-count-distinct-interval MP_root-count-distinct-positive
  * Construct sequences of polynomials
  ** MP_seq-make-fourier MP_seq-make-sturm-canonical MP_seq-make-chebyshev MP_seq-make-legendre MP_seq-make-lagrange
     MP_seq-make-bernstein MP_seq-make-hermite MP_seq-make-laguerre
  * Construct polynomials
  ** MP_make-chebyshev MP_make-legendre MP_make-lagrange MP_make-bernstein MP_make-hermite MP_make-laguerre
  * Polynomial metrics
  ** Sign change based
  *** MP_scount-descartes MP_scount-sturm MP_scount-fourier MP_scount-budan
  ** Other
  *** MP_height MP_length MP_density(G) MP_index(G) MP_degree(G) MP_discriminant-low-degree
      MP_discriminant-high-degree MP_discriminant MP_content MP_mahler-measure
  * Other Stuff
  ** mjr_poly_resultant
  * Polynomial xforms
  ** MP_deflate MP_shift MP_reflect MP_diff(G) MP_intg MP_subst(G) MP_zap-zero-roots MP_2square-free MP_2primitive
     MP_simplify(G) MP_scale(G) MP_rationalize
  ** MP_cubic-depress mjr_poly_tschirnhaus-3-2"
  (documentation 'mjr_poly_help 'function))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(mjr_gpoly_make-coeff          "")
(mjr_gpoly_make-simplify       "")
(mjr_gpoly_make-scale          "")
(mjr_gpoly_make-diff           "")
(mjr_gpoly_make-eval           "")
(mjr_gpoly_make-+              "")
(mjr_gpoly_make--              "")
(mjr_gpoly_make-*              "")
(mjr_gpoly_make-iexpt          "")
(mjr_gpoly_make-subst          "")
(mjr_gpoly_make-leading-coeff  "")
(mjr_gpoly_make-degree         "")
(mjr_gpoly_make-density        "")
(mjr_gpoly_make-index          "")
(mjr_gpoly_make-constant-coeff "")
(mjr_gpoly_make-truncate       "")
(mjr_gpoly_make-rem            "")
(mjr_gpoly_make-mod            "")
(mjr_gpoly_make-onep           "")
(mjr_gpoly_make-zerop          "")
(mjr_gpoly_make-constantp      "")
(mjr_gpoly_make-divides?       "")
(mjr_gpoly_make-gcd            "")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_print (poly &key (var "x") (add-op "+") (add-apad " ") (add-bpad " ") (expt-op "^") (expt-apad "")
                       (expt-bpad "") (mul-op "*") (mul-apad "") (mul-bpad "") (cof-fmt "~s") (cof-bpad "") (cof-apad "")
                       (pow-fmt "~d") (pow-bpad "") (pow-apad "") (ply-apad "") (ply-bpad "") (suppress-zero-terms 't)
                       (suppress-unit-expt 't) (suppress-unit-prod 't) (suppress-zero-powers 't) (return-string nil))
  "Print a polynomial.  Return POLY.
Arguments:
  - return-string .......... If non-NIL, then do not print but return a string rep of POLY with no following newline.
                             If NIL, then print POLY with a following newline, and return POLY
  - var .................... variable to use (a string)
  - add-op ................. string to print for + signs
  - add-bpad ............... string to print before + signs
  - add-apad ............... string to print after + signs
  - cof-fmt ................ format string for coefficients
  - cof-bpad ............... string to print before coefficients
  - cof-apad ............... string to print after coefficients
  - pow-fmt ................ format string powers
  - pow-bpad ............... string to print before powers
  - pow-apad ............... string to print after powers
  - mul-op ................. string to print for * signs
  - mul-bpad ............... string to print before * signs
  - mul-apad ............... string to print after * signs
  - expt-op ................ string to print for ^ signs
  - expt-bpad .............. string to print before ^ signs
  - expt-apad .............. string to print after ^ signs
  - suppress-zero-terms .... do not print 0*x^n
  - suppress-unit-expt ..... do not print ^1
  - suppress-unit-prod ..... do not print 1*
  - suppress-zero-powers ... do not print x^0"
  (let ((da-str (with-output-to-string (str-out)
                  (if (vectorp poly)
                      (let ((plen (length poly)))
                        (if (< 0 plen)
                            (loop with not-first-print = nil
                                  for power from (1- plen) downto 0
                                  for term across poly
                                  initially (format str-out ply-bpad)
                                  finally   (format str-out ply-apad)
                                  do (if (or (not suppress-zero-terms) (mjr_cmp_!=0 term) (and (= power 0) (not not-first-print)))
                                         (progn (if not-first-print (format str-out "~a~a~a" add-bpad add-op add-apad))
                                                (if (or (not suppress-unit-prod) (mjr_cmp_!= 1 term) (zerop power))
                                                    (format str-out (concatenate 'string cof-bpad cof-fmt cof-apad) term))
                                                (if (and (or (not suppress-unit-prod) (mjr_cmp_!= 1 term))
                                                         (or (not suppress-zero-powers) (not (zerop power))))
                                                    (format str-out "~a~a~a" mul-bpad mul-op mul-apad))
                                                (if (or (not suppress-zero-powers) (not (zerop power)))
                                                    (progn (format str-out "~a" var)
                                                           (if (or (not suppress-unit-expt) (not (= 1 power)))
                                                               (progn (format str-out "~a~a~a" expt-bpad expt-op expt-apad)
                                                                      (format str-out  (concatenate 'string pow-bpad pow-fmt pow-apad) power)))))
                                                (setq not-first-print 't))))))))))
    (if return-string
        (if (vectorp poly)
            da-str)
        (progn (if (vectorp poly)
                   (format 't "~a~%" da-str))
               poly))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_code (poly &key (lang :lang-matlab))
  "Return a string using the syntax of the selected programming language or computational environment."
  (if (not (vectorp poly))
      (error "mjr_poly_code: Argument must be a polynomial (a vector)!"))
  (with-output-to-string (str-out)
    (flet ((gpp (&key (var "x") (add-op "+") (add-apad "") (add-bpad "") (expt-op "^") (expt-apad "")
                      (expt-bpad "") (mul-op "*") (mul-apad "") (mul-bpad "") (cof-bpad "") (cof-apad "")
                      (pow-bpad "") (pow-apad "") (ply-apad "") (ply-bpad ""))
             (loop with plen = (length poly)
                   for not-first-print = nil then 't
                   for power from (1- plen) downto 0
                   for term across poly
                   initially (format str-out ply-bpad)
                   finally   (format str-out ply-apad)
                   when not-first-print
                   do (format str-out "~a~a~a" add-bpad add-op add-apad)
                   do (progn (format str-out "~a~a~a~a" cof-bpad (mjr_numu_code term :lang lang) cof-apad cof-bpad)
                             (if (not (zerop power))
                                 (format str-out "~a~a~a~a~a~a~a~a~a~a" mul-bpad mul-op mul-apad var pow-bpad expt-bpad expt-op expt-apad (mjr_numu_code power :lang lang) pow-apad))))))
      (case lang
        ((:lang-povray
          :lang-matlab
          :lang-octave
          :lang-idl
          :lang-r
          :lang-hp48
          :lang-mathematica
          :lang-maple
          :lang-maxima)      (gpp))
        ((:lang-latex
          :lang-pdflatex
          :lang-amstex
          :lang-tex)         (gpp :ply-bpad "$$" :ply-apad "$$" :mul-op " " :pow-bpad "{" :pow-apad "}"))
        ((:lang-lisp)        (progn (format str-out "(+")
                                    (loop for power from (1- (length poly)) downto 0
                                          for term across poly
                                          do (if (mjr_cmp_!=0 term)
                                                 (let ((tstr (if (or (mjr_cmp_!= 1 term) (zerop power))
                                                                 (format nil " ~a" term)))
                                                       (vstr (if (not (zerop power))
                                                                 (if (= 1 power)
                                                                     (format nil " x")
                                                                     (format nil " (expt x ~a)" power)))))
                                                   (if (and tstr vstr)
                                                       (format str-out " (*~a~a)" tstr vstr)
                                                       (format str-out "~a" (or tstr vstr))))))
                                    (format str-out ")")))
        ('t                  (error "mjr_poly_code: Language unsupported!"))))
    str-out))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_intg (poly a b)
  "Integrate over the interval $[a,b]$."
  (let* ((poly  (mjr_poly_simplify poly))
         (plen  (length poly))
         (ipoly (make-array (1+ plen) :initial-element 0)))
    (loop for i from 0 upto (1- plen)
          do (setf (aref ipoly i) (/ (aref poly i) (- plen i))))
    (- (mjr_poly_eval ipoly b) (mjr_poly_eval ipoly a))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_zap-zero-roots (poly &optional eps)
  "Return two values: 1) a new polynomial, $\\text{POLY}/x^n$, where $n$ is the number of times $0$ is a root of POLY, and 2) $n$.
EPS is used with mjr_cmp_!=0 to detect non-zero coefficients."
  (let* ((plen (length poly))
         (zcnt (loop for i from (1- plen) downto 0
                     for ip = (aref poly i)
                     until (mjr_cmp_!=0 ip eps)
                     count 1))
         (slim (- plen zcnt)))
    (values (if (zerop slim)
                #(1)
                (subseq poly 0 slim))
            zcnt)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_deflate (poly const)
  "Return the quotient and the remainder from POLY / (X - CONST).
Note: The remainder is equal to POLY(CONST).
While the result is equivalent to (MJR_POLY_TRUNCATE POLY (VECTOR 1 (- CONST))), this function is faster."
  (cond ((not (numberp const)) (error "mjr_poly_eval: CONST must be a number!")))
  (let* ((poly    (mjr_poly_simplify poly))
         (plen    (length poly))
         (polyval 0)
         (newpoly (mjr_vec_make-const (max 1 (1- plen)))))
    (loop for i from 0 upto (1- plen)
          for ht = (+ (* const polyval) (aref poly i))
          do (setq polyval ht)
          when (not (= i (1- plen)))
          do (setf (aref newpoly i) ht)
          finally (return (values newpoly polyval)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_shift (b poly)
  "Shift POLY right B units (or left if B is negative).
References:
;; MJR TODO NOTE <2011-11-23 13:07:23 CST> mjr_poly_shift: Add a reference here.  The algorithm is the 'fast taylor shift'.  The pink book or sure.  Perhaps Prasolov too..."
  (let ((new-poly (mjr_poly_simplify (copy-seq poly)))
        (b        (- b)))
    (mjr_poly_simplify (dotimes (i (length new-poly) new-poly)
                         (loop for k from i downto 1
                               do (incf (aref new-poly k) (* b (aref new-poly (1- k)))))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_reflect (poly)
  "Reflect POLY across the Y-axis (i.e. subst -x)."
  (let ((new-poly (mjr_poly_simplify (copy-seq poly))))
    (loop for cur-coef across poly
          for odd-pow = (not (oddp (length poly))) then (not odd-pow)
          for i from 0
          when odd-pow
          do (setf (aref new-poly i) (- cur-coef))
          finally (return (mjr_poly_simplify new-poly)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_make-from-roots (&rest the-roots)
  "Return the the monic polynomial with the given roots.
The roots may be given as individual arguments, or a single list may be provided."
  (let ((the-roots (if (and (car the-roots) (listp (car the-roots))) (car the-roots) the-roots)))
    (let* ((n      (length the-roots))
           (a-lst  (make-array (1+ n) :initial-element 0))
           (a-nxt  (make-array (1+ n) :initial-element 0)))
      (loop initially (setf (aref a-lst n) (- (car the-roots)))
            for k from 1 upto (1- n)
            for zk in (cdr the-roots)
            do (loop initially (setf (aref a-nxt (- n k)) (- (aref a-lst (- n (1- k))) zk))
                     for j from 1 upto (1- k)
                     do (setf (aref a-nxt (- n j)) (- (aref a-lst (1+ (- n j))) (* zk (aref a-lst (- n j)))))
                     finally (setf (aref a-nxt n) (* (- zk) (aref a-lst n))))
            do (rotatef a-lst a-nxt)
            finally (progn (setf (aref a-lst 0) 1)))
      a-lst)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_eval-nth-derivative (poly x &optional (order 1))
  "Evaluate polynomial derivative."
  (mjr_poly_eval (mjr_poly_diff poly order) x))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_eval-poly-and-first-n-derivatives (poly x &optional (order 1))
  "Return value of poly and the values of the first ORDER derivatives at x.

Algorithm due to Pankiewicz:
  W. Pankiewicz (1968); Algorithm 337: calculation of a polynomial and its derivative values by Horner scheme; Communications of the ACM; Vol 11, Issue 9, pp 633"
  (let* ((len-1  (1- (length poly)))
         (odr-1  order)
         (pd     (make-array (1+ order) :initial-element 0)))
    (loop for i from (1- len-1) downto 0
          for nnd = (min odr-1 (- len-1 i))
          initially (setf (aref pd 0) (aref poly 0))
          do (loop for j from nnd downto 1
                   do (setf (aref pd j) (+ (* (aref pd j) x) (aref pd (1- j)))))
          do (setf (aref pd 0) (+ (* (aref pd 0) x) (aref poly (- len-1 i)))))
    (loop with cnst = 1
          for i from 2 upto odr-1
          do (setf cnst (* cnst i))
          do (setf (aref pd i) (* (aref pd i) cnst)))
    (values-list (concatenate 'list pd))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_2func (poly &optional num-derivatives)
  "Return a function that evaluates the POLY and its first NUM-DERIVATIVES."
  (if (and num-derivatives (> num-derivatives 0))
      (eval `(lambda (x) (mjr_poly_eval-poly-and-first-n-derivatives ,poly x ,num-derivatives)))
      (eval `(lambda (x) (mjr_poly_eval ,poly x)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-eval (seq x &optional result-type)
  "Evaluate the polynomials in the sequence (list or vector) at $x$, and return the results in a sequence.

If result-type is NIL, then the returned type will be the same as the type of seq."
  (map (or result-type (if (vectorp seq) 'vector 'list))
       (lambda (p) (mjr_poly_eval p x))
       seq))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_height (poly)
  "This is the infinity-norm (maximum absolute value of coefficients)"
  (mjr_vec_norm-infinity poly))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_length (poly)
  "This is the one-norm (sum of absolute values of coefficients)"
  (mjr_vec_norm-one poly))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_2square-free (p)
  "Return a new polynomial with the same roots as poly, but with no multiple roots."
  (mjr_poly_truncate p (mjr_poly_gcd p (mjr_poly_diff p))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_2monic (poly)
  "Transform POLY, via multiplication, into a monic polynomial that shares the roots of the original."
  (mjr_vec_/ poly (aref poly 0)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_rationalize (poly)
  "Transform POLY so that all coefficients are rational -- imaginary parts of coefficients too."
  (map 'vector
       (lambda (x) (if (complexp x)
                       (complex (rationalize (realpart x))
                                (rationalize (imagpart x)))
                       (rationalize x)))
       poly))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_content (poly)
  "Compute the content of an integer or rational polynomial.  POLY is rationalized if necessary.

If $p\\in\\mathbb{Z}[x]$, then the content is the GCD of the coefficients.  If $p\\in\\mathbb{Q}[x], then the content is the GCD of the numerators divided by
the LCM of the denominators.  i.e. the content is the unique $q\\in\\mathbb{Q}$ such that $p/q$ is a primitive polynomial in $\\mathbb{Z}[x]$ -- all the
coefficients are integers have have GCD of $1$."
  (let* ((poly (mjr_poly_simplify poly)))
    (cond ((every #'integerp   poly) (reduce #'gcd poly))
          ((every #'rationalp  poly) (/ (apply #'gcd (map 'list #'numerator   poly))
                                        (apply #'lcm (map 'list #'denominator poly))))
          ((some  #'complexp poly)   (error "mjr_poly_content: Polynomial must not be complex!"))
          ((every #'numberp  poly)   (mjr_poly_content (mjr_poly_rationalize poly)))
          ('t                        (error "mjr_poly_content: POLY is not a polynomial!")))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_primitive-part (poly)
  "Returns the primitive part of the polynomial and the content.

Note: poly = cont(poly)*primitive-part(poly)"
  (let* ((poly (mjr_poly_simplify poly))
         (c    (mjr_poly_content poly)))
    (values (mjr_vec_/ poly c)
            c)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_2primitive (poly)
  "Transform POLY, via multiplication by a rational number, into an primitive polynomial that shares the roots of the original.

The rational number used to transform POLY is the second return value.

Results are exact if the polynomial has only integer and/or rational coefficients.

A simple generalization of the idea of primitive is used when the polynomial has complex coefficients -- The real and imaginary parts will both be transformed
to integers and the entire collection of integers (real parts and imaginary parts together) will not share a common factor other than one.

Floating point numbers in the coefficients will be rationalized. The finite precision of floating point arithmetic implies that polynomials made up of floats,
doubles, and long doubles are but a simple scalar product away from their integer doppelganger.  Unfortunately the vulgarities of floating point arithmetic
can destroy this almost-truth.  Still, we do the best we can.

NOTE: A polynomial $p(x)=\sum_{j=0}^na_jx^j\in\mathbb{Z}[x]$ is primitive if $\mathrm{GCD}(a_0,\cdots,a_n)=1$ -- i.e. the coefficients are relatively
prime. Gauss's lemma: If $p$ and $q$ are primitive, then $p\cdot q$ primitive, and if a $p\in\mathbb{Z}[x]$ is irreducible over the $\mathbb{Z}$, then the
same polynomial considered in $\mathbb{Q}[x]$ is also irreducible over $\mathbb{Q}$.

References:
  Carl Friedrich Gauss(1801); Disquisitiones Arithmeticae; Article 42"
  (let* ((poly (mjr_poly_simplify poly))
         (p1   (mjr_poly_rationalize poly))
         (L    (apply #'lcm
                      (mapcar #'denominator (loop for coef across p1
                                                  collect (realpart coef)
                                                  when (complexp coef)
                                                  collect (imagpart coef)))))
         (p2   (mjr_vec_* p1 L))
         (G    (apply #'gcd (loop for coef across p2
                                  collect (realpart coef)
                                  when (complexp coef)
                                  collect (imagpart coef))))
         (p3   (mjr_vec_/ p2 G)))
    (values p3 (/ L G))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-sturm-canonical (p)
  "Return a list of polynomials representing the 'Canonical Sturm sequence'.

The canonical Sturm sequence of a polynomial $p$ is the intermediate results of Euclid's algorithm applied to to $p$ and $p'$:
  $$
  \\begin{array}{lccclcl}
     p_0(x)  &   &   & = & p(x)                          &   &                                          \\\\
     p_1(x)  &   &   & = & p'(x)                         &   &                                          \\\\
     p_2(x)  &   &   & = & -\\text{rem}(p_0,     p_1)    & = & p_1(x)    \\cdot q_0(x)     - p_0(x)     \\\\
     p_3(x)  &   &   & = & -\\text{rem}(p_1,     p_2)    & = & p_2(x)    \\cdot q_1(x)     - p_1(x)     \\\\
     \\vdots &   &   &   &                               &   &                                          \\\\
     p_m     &   &   & = & -\\text{rem}(p_{m-2}, p{m-1}) & = & p{m-1}(x) \\cdot q_{m-2}(x) - p_{m-2}(x) \\\\
     p_{m+1} & = & 0 & = & -\\text{rem}(p_{m-1}, p_m)    &   &                                          \\\\
  \\end{array}
  $$

where $\\text{rem}(p_i,p_j)$ and $q_i$ are the remainder and the quotient of the polynomial long division of $p_i$ by $p_j$, and where $m\\le\\text{deg}{p}$
is the minimal number of polynomial divisions needed to obtain a zero remainder.

Note: The sequence is well defined even when $p$ is not square-free; however, it may not be a Sturm sequence in this case.  Confusingly enough, this sequence
is always called the canonical Sturm sequence even when it is not really a Sturm sequence."
  (let ((pd (mjr_poly_diff p)))
    (nconc (list p pd)
           (loop for p-2 = p  then p-1
                 for p-1 = pd then p-0
                 for p-0 = (mjr_poly_- (mjr_poly_rem p-2 p-1))
                 collect p-0
                 until (zerop (mjr_poly_degree p-0))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-fourier (p)
  "Return a list of polynomials representing the 'Fourier sequence'.

Let $p$ be a real polynomial of degree $n>0$, then the Fourier sequence of $p$ is:
   $$F_\\text{seq}(x)=\\big\\{ p(x), p^{(1)}(x),\\ldots,p^{(n)}(x)\\big\\}$$"
  (loop for i from 0 upto (mjr_poly_degree p)
        for p-0 = p then (mjr_poly_diff p-0)
        collect p-0))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_count-sign-changes (seq)
  "Return the number of sign changes in seq (vector or list) -- zeros are ignored"
  (if (vectorp seq)
      (loop for cur-elt across seq
            for lst-sgn = nil then (if (zerop cur-sgn) lst-sgn cur-sgn)
            for (cur-sgn cur-amb) = (multiple-value-list (mjr_cmp_signum cur-elt))
            count (and lst-sgn (not (zerop cur-sgn)) (not (= cur-sgn lst-sgn))) into count-sgn
            count cur-amb                                                       into count-amb
            finally (return (values count-sgn (if (not (zerop count-amb)) count-amb))))
      (loop for cur-elt in seq
            for lst-sgn = nil then (if (zerop cur-sgn) lst-sgn cur-sgn)
            for (cur-sgn cur-amb) = (multiple-value-list (mjr_cmp_signum cur-elt))
            count (and lst-sgn (not (zerop cur-sgn)) (not (= cur-sgn lst-sgn))) into count-sgn
            count cur-amb                                                       into count-amb
            finally (return (values count-sgn (if (not (zerop count-amb)) count-amb))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_scount-sturm (p l r &optional sturm-sequence)
  "Compute the difference in the sign changes for $S(x+l)$ and $S(x+h)$ where S is the given Sturm sequence or the canonical one.

The return will be the number of distinct real roots the polynomial has in the interval $(l,h]$ when:
  1) sturm-sequence is NIL or the canonical Sturm sequence
  2) when sturm-sequence is a valid Sturm sequence and $p$ is square-free.

A Sturm sequence is a set of polynomials $p_0, p_1, \\dots, p_m\\in\\mathbb{R}[x]$ such that:
  \\begin{itemize}
    \\item $p_i$ are all polynomials
    \\item $m$ is a finite, non-negative integer
    \\item if $0\\ge i< j\\le m$ then $\\text{deg}(p_i)<\\text{deg}(p_j)$
    \\item $p_0$ is square free (no repeated roots)
    \\item if $p(x_0)=0$, then $\\text{sign}(p_1(x_0))= \\text{sign}(p'(x_0))$
    \\item if $p_i(\\xi)=0$ for $0<i<m$ then $\\text{sign}(p_{i-1}(\\xi))= -\\text{sign}(p_{i+1}(\\xi))$
    \\item $p_m$ does not change its sign.
  \\end{itemize}

Sturm's theorem\\\\
  Let $p$ be a real, square-free polynomial of degree $n>0$ and $l,r\\in \\mathbb{R}$ with $l<r$.  Let $S_\\text{seq}(x)=\\big\\{ p_0(x), p_1(x), ...,
  p_m(x)\\big\\}$ be a Sturm sequence with $p_0=p$.  Let $v_l$ and $v_r$ be the sign variations in the sequences $S_\\text{seq}(l)$ and $S_\\text{seq}(r)$.
  Let $\\rho$ be the number of the real roots of $p$ in $(l,r]$. Let $v_\\delta=v_l - v_r$.  Then we have $\\rho = v_\\delta$

Sturm's Canonical Sequence theorem\\\\
  Let $p$ be a real polynomial of degree $n>0$ and $l,r\\in \\mathbb{R}$ with $l<r$ and neither $l$ or $r$ is a multiple root of $p$.  Let
  $S_\\text{seq}(x)=\\big\\{ p_0(x), p_1(x), ..., p_m(x)\\big\\}$ be the canonical Sturm sequence of $p$.  Let $v_l$ and $v_r$ be the sign variations in the
  sequences $S_\\text{seq}(l)$ and $S_\\text{seq}(r)$.  Let $\\rho$ be the number of distinct real roots of $p$ in $(l,r]$. Let $v_\\delta=v_l - v_r$.  Then
  we have $\\rho = v_\\delta$

References:
  Jacques Charles Francois Sturm(1829); Memoire sur la resolution des equations numeriques; Bulletin des Sciences de Ferussac 11; 419-425"
  ;; MJR TODO NOTE <2013-01-01 19:01:18 CST> mjr_poly_scount-sturm: Add refs (pink and blue polynomial books)
  (let ((ss (or sturm-sequence (mjr_poly_seq-make-sturm-canonical p))))
    (- (mjr_poly_count-sign-changes (mjr_poly_seq-eval ss l))
       (mjr_poly_count-sign-changes (mjr_poly_seq-eval ss r)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_scount-budan (p l h)
  "Compute the difference in the sign changes for $p(x+l)$ and $p(x+h)$

Budan's theorem:\\\\
  Let $p\\in\\mathbb{R}[x]$ with degree $n>0$ and $l,r\\in \\mathbb{R}$ with $l<r$ and $p(r)\\ne0$. Let $v_l$ and $v_r$ be the sign variations in $p(x+l)$ and
  $p(x+r)$ respectively. Let $\\rho$ be the number of the real roots of $p$ in $(l,r)$. Let
  $v_\\delta=v_l - v_r$.
    \\begin{enumerate}
      \\item $v_l \\ge v_r$
      \\item $\\rho \\le v_\\delta$
      \\item if $\\rho < v_\\delta$, then $\\rho = v_\\delta -2\\lambda$ where $\\lambda \\in \\mathbb{Z}_+$
    \\end{enumerate}

References:
  Francois D. Budan (1807); Nouvelle methode pour la resolution des equations numeriques; Paris: Courcier"
  (if (not (mjr_chk_!=0 (mjr_poly_eval p h)))
      (warn "budan-count: Polynomial should not be zero on h if Budan's theorem is being used!"))
  (- (mjr_poly_count-sign-changes (mjr_poly_shift (- l) p))
     (mjr_poly_count-sign-changes (mjr_poly_shift (- h) p))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_scount-fourier (p l h)
  "Compute the difference in the sign changes for $F(x+l)$ and $F(x+h)$ where F is the 'Fourier Sequence' for $p$

Fourier's theorem:\\\\
  Let $p$ be a real polynomial of degree $n>0$ and $l,r\\in \\mathbb{R}$ with $l<r$ and $p(r)\\ne0$. Let $F_\\text{seq}(x)$ be the Fourier sequence of $p$.
  Let $v_l$ and $v_r$ be the sign variations in the sequences $F_\\text{seq}(l)$ and $F_\\text{seq}(r)$.  Let $\\rho$ be the number of the real roots of $p$
  in $(l,r)$. Let $v_\\delta=v_l - v_r$.
    \\begin{enumerate}
      \\item $v_l \\ge v_r$
      \\item $\\rho \\le v_\\delta$
      \\item if $\\rho < v_\\delta$, then $\\rho = v_\\delta -2\\lambda$ where $\\lambda \\in \\mathbb{Z}_+$
    \\end{enumerate}

References:
  Jean Baptiste Joseph Fourier (1820); Sur l'usage du theoreme de Descartes dans la recherche des limites des racines; Bulletin des Sciences, par la Societe Philomatique de Paris: 156-165
  E.J. Barbeau (2003); Polynomials (Problem Books in Mathematics); pp 174
  Victor V. Prasolov (2004); Polynomials; ISBN: 3540407146; pp 27"
  (if (not (mjr_chk_!=0 (mjr_poly_eval p h)))
      (warn "fourier-count: Polynomial should not be zero on h if Budan's theorem is being used!"))
  (let ((d  (mjr_poly_degree p)))
    (- (mjr_poly_count-sign-changes (multiple-value-list (mjr_poly_eval-poly-and-first-n-derivatives p l d)))
       (mjr_poly_count-sign-changes (multiple-value-list (mjr_poly_eval-poly-and-first-n-derivatives p h d))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_scount-descartes (poly &optional (c 0))
  "Apply Descartes rule of signs to the shifted polynomial.

Four values returned:
 * Max number roots ($P$) to the right of $C$   -- Actual number will be $P-2\\cdot m$ where $m$ is a non-negative integer.
 * Exact number of roots ($Z$) at $C$           -- Computed from a direct analysis of the polynomial.
 * Max number of roots ($N$) to the left of $C$ -- Actual number will be $N-2\\cdot n$ where $n$ is a non-negative integer.
 * Number of ambiguous zeros, or nil            -- This occurs when a value is assumed zero because it was small, but not zero

Theorem (Descartes):
   Let
   $$p(x)=\\sum\\limits_{j=0}^n a_jx^j\\in\\mathbb{R}[x]$$
   Let $P=\\sigma(p)$ and $N=\\sigma(p(-x))$ where the $\\sigma$ function is defined as the number of sign changes in the polynomial coefficients ignoring
   zero coefficients.  Let $Z$ be the number of trailing zeros -- i.e. number of zero coefficients at the end of the polynomial.  Then the polynomial $p$ has
   $P-2m$ positive roots, $Z$ zero roots, and $N-2n$ negative roots where $m$ and $n$ non-negative integers.

References:
  Rene Descartes (1637); La Geometrie (an appendix to Discours de la methode -- Discourse on Method)
  Michael Mahoney (1979); Geometry (a translation of La Geometrie)
  Grabiner (1999); Descartes Rule of Signs: Another Construction; Amer. Math. Monthly 106, 854-855
  Victor V. Prasolov (2004); Polynomials; ISBN: 3540407146; pp 28"
  (let* ((poly (if (zerop c) (mjr_poly_simplify poly) (mjr_poly_shift c poly)))
         (plen (length poly)))
    (cond ((< plen 2) (error "mjr_poly_scount-descartes: Polynomial must be at least degree 1!")))
    (loop for cur-coef across poly
          for odd-pow = (not (oddp plen)) then (not odd-pow)
          for count=0 = 0 then (if (mjr_cmp_=0 cur-coef) (1+ count=0) 0)
          for lst-sign   = nil then (if (zerop cur-sign)   lst-sign cur-sign)
          for lst-sign-x = nil then (if (zerop cur-sign-x) lst-sign-x cur-sign-x)
          for (cur-sign cur-ambiguous) = (multiple-value-list (mjr_cmp_signum cur-coef))
          for cur-sign-x  = (if odd-pow (- cur-sign) cur-sign)
          count (and lst-sign   (not (zerop cur-sign))   (not (= cur-sign   lst-sign)))   into count+x
          count (and lst-sign-x (not (zerop cur-sign-x)) (not (= cur-sign-x lst-sign-x))) into count-x
          count cur-ambiguous                                                             into count-ambiguous-zeros
          finally (return (values count+x count=0 count-x (if (not (zerop count-ambiguous-zeros)) count-ambiguous-zeros))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-bound-positive (poly &key (algorithm :rb-lagrange))
  "Return an upper bound (perhaps strict) on the positive real roots of POLY.  Returns NIL if anything goes wrong."
  (let* ((poly (mjr_poly_simplify poly))
         ;; MJR TODO NOTE <2013-01-01 18:59:12 CST> mjr_poly_root-bound-positive: Add references, formula, and notes that it is also a global bound...
         ;; MJR TODO NOTE <2013-01-01 18:59:47 CST> mjr_poly_root-bound-positive: Move into mjr_poly_root-bound-all????
         (pdeg (mjr_poly_degree   poly))
         (pn   (aref poly 0)))
    (if (mjr_chk_!=0 pn)
        (case algorithm
          (:rb-lagrange    (and (> pdeg 0)
                                (loop with max1 = nil
                                      with max2 = nil
                                      for j from 0 upto (1- pdeg)
                                      for pjr = (/ (aref poly (- pdeg j)) pn)
                                      when (< pjr 0)
                                      do (let ((pjv (expt (abs pjr) (/ (- pdeg j)))))
                                           (if (null max1)
                                               (setf max1 pjv)
                                               (if (> pjv max1)
                                                   (rotatef max2 max1
                                                            max1 pjv)
                                                   (if (or (null max2) (> pjv max2))
                                                       (setf max2 pjv)))))
                                      finally (return (and max1 max2 (+ max1 max2))))))
          (otherwise       (error "mjr_poly_root-bound-positive: Unknown algorithm!"))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-bound-all (poly &key (algorithm :rb-fujiwara))
  "Return an upper bound (perhaps strict) on the modulus of the roots of POLY using the given algorithm

If the test can not be performed because the polynomial violates a constraint of the algorithm, then NIL is returned.

Note: $$P(z)=\\sum\\limits_{i=0}^n a_ix^i$$

The algorithm argument may be one of, or a list of, the following symbols:
  * :rb-fujiwara    - 8.1.11 - Fujiwara (1916)     - $$2\\max_{0\\le j\\le n-1}\\left|\\frac{a_j}{a_n}\\right|$$
  * :rb-cauchy      - 8.1.9  - Cauchy (1829, p122) - $$1+\\max_{0\\le j\\le n-1}\\left|\\frac{a_j}{a_n}\\right|$$
                             - M,M,R (1994, p224)  - $$1+\\frac{1}{\\vert a_n\\vert}\\left(\\max\\limits_{ 0\\leq i\\leq n-1}\\vert a_i\\vert\\right)$$
  * :rb-kojima      - 8.1.17 - Kojima (1917)       - $$\\max\\left\\{ \\left|\\frac{a_0}{a_1}\\right|, 2\\left|\\frac{a_1}{a_2}\\right|, \\cdots, 2\\left|\\frac{a_{n-1}}{a_n}\\right| \\right\\}$$
  * :rb-hirst-macey - N/A    - Hirst,Macey (1997)  - $$ ? $$
  * :rb-lagrange    - N/A    - Lagrange (????)     - $$ ? $$

Not yet supported, but will be in a future version
  * :rb-carmichal   - 8.1.12 - Carmichal (1918)    - $$\\sum_{j=0}^{n-1}\\left(\\left|\\frac{a_j}{a_n}\\right|^\\frac{1}{n-j}\\right)$$
  * :rb-westerfield - N/A    - Westerfield (1931)  - $$\\hat{M}+\\check{M}$$ where $\\hat{M}$ and $\\check{M}$ are the two largest elements of $$\\left\\{ \\left|\\frac{a_j}{a_n}\\right|^\\frac{1}{n-j} : j=0,...,{n-1}\\right\\}$$

References:
  Fujiwara (1916); Uber die obere Schranke des absoluten Betrages der Wurzeln einer algebraischen Gleichung; Tohoku Math J 10; p167-171
  Kojima (1917); On the limits of the roots of an algebraic equation; Tohoku Math J 11; p119-127
  Cauchy (1829); Exercises de mathematique. Oeuvres 2 (9)
  Milovanovic, Mitrinovic, Rassias (1994); Topics in polynomials: extremal problems, inequalities, zeros; ISBN: 981020499X
  Hirst & Macey (1997); Bounding the roots of polynomials; Coll Math J 28 (4); pp292"
  ;; MJR TODO NOTE <2013-01-01 18:57:44 CST> mjr_poly_root-bound-all: IMPLEMENT :rb-carmichal and :rb-westerfield.
  ;; MJR TODO NOTE <2013-01-01 18:58:19 CST> mjr_poly_root-bound-all: Add refs for Carmichal and Westerfield.
  ;; MJR TODO NOTE <2013-01-01 18:58:37 CST> mjr_poly_root-bound-all: Add ref for blue polynomial book and/or blue computer algebra book.
  (let* ((poly (mjr_poly_simplify poly))
         (pdeg (mjr_poly_degree   poly))
         (pn   (aref poly 0)))
    (if (mjr_chk_!=0 pn)
        (flet ((cmpbnd (alg)
                 (case alg
                   (:rb-fujiwara    (* 2 (max (expt (abs (/ (aref poly pdeg) pn)) (/ pdeg))
                                              (loop for j from 1 upto (1- pdeg)
                                                    maximize (expt (abs (/ (aref poly (- pdeg j)) pn)) (/ (- pdeg j)))))))
                   (:rb-kojima      (if (mjr_chk_!=0 (aref poly (1- pdeg)))
                                        (* 2 (max (* 1/2 (abs (/ (aref poly pdeg) (aref poly (1- pdeg)))))
                                                  (loop for j from 1 upto (1- pdeg)
                                                        for pj+1 = (aref poly (- pdeg (1+ j)))
                                                        maximize (if (mjr_chk_!=0 pj+1)
                                                                     (abs (/ (aref poly (- pdeg j)) pj+1))
                                                                     (return-from cmpbnd nil)))))))
                   (:rb-cauchy      (1+ (/ (loop for i from 1 upto pdeg
                                                 maximize (abs (aref poly i)))
                                           (abs pn))))
                   (:rb-lagrange    (mjr_poly_root-bound-positive poly :algorithm :rb-lagrange))
                   (:rb-hirst-macey (max 1
                                         (loop for j from 0 upto (1- pdeg)
                                               for pj = (aref poly (- pdeg j))
                                               maximize (if (mjr_chk_!=0 pj)
                                                            (abs (/ pj pn))
                                                            (return-from cmpbnd nil)))))
                   (otherwise       (error "mjr_poly_root-bound-all: Unknown algorithm(~a)!" alg)))))
          (if (listp algorithm)
              (mapcar #'cmpbnd algorithm)
              (cmpbnd algorithm))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-bound-all-cauchy (poly &key (xeps 0.0001) (yeps 0.0001) (max-itr 1000))
  "Compute an approximation of the cauchy bound on the roots of the polynomial."
  (loop with soly = (mjr_poly_zap-zero-roots (mjr_poly_simplify poly))
        with coly = (apply #'vector (loop for ai  across soly
                                          for aai = (- (abs ai)) then (abs ai)
                                          collect aai))
        with ubnd = (or (mjr_poly_root-bound-positive coly :algorithm :rb-lagrange)
                        (mjr_poly_root-bound-all      coly :algorithm :rb-fujiwara))
        with x0   = 0
        with x1   = ubnd
        with y0   = (mjr_poly_eval coly x0)
        with y1   = (mjr_poly_eval coly x1)
        for x-mid = (/ (+ x0 x1) 2)
        for y-mid = (mjr_poly_eval coly x-mid)
        for i-cur from 1 upto max-itr
        until (or (mjr_eps_= 0  y-mid yeps)
                  (mjr_eps_= x0 x1    xeps)
                  (mjr_eps_= x0 x-mid xeps)
                  (mjr_eps_= x1 x-mid xeps)
                  (mjr_eps_= y0 y1    yeps)
                  (mjr_eps_= y0 y-mid yeps)
                  (mjr_eps_= y1 y-mid yeps))
        while (cond ((mjr_cmp_< y0 y1) (if (mjr_cmp_< y-mid 0)
                                           (setq x0 x-mid y0 y-mid)
                                           (setq x1 x-mid y1 y-mid)))
                    ((mjr_cmp_> y0 y1) (if (mjr_cmp_< y-mid 0)
                                           (setq x1 x-mid y1 y-mid)
                                           (setq x0 x-mid y0 y-mid))))
        finally (return x1)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-count-distinct-interval (p a b &optional eps)
  "Use mjr_poly_scount-sturm to compute the total number of distinct, real roots of the polynomial in $(a,b)$

If polynomial is zero at a or b, then those roots will be removed via deflation.  If floating point comparisons are required to determine zeros at a and/or b,
then eps will be used.  Polynomial rationalization via mjr_poly_2primitive is recommended if floating point round-off is a problem."
  (if (mjr_chk_!=0 (mjr_poly_eval p a) eps)
      (if (mjr_chk_!=0 (mjr_poly_eval p b) eps)
          (mjr_poly_scount-sturm p a b)
          (mjr_poly_root-count-distinct-interval (mjr_poly_deflate p b) a b eps))
      (mjr_poly_root-count-distinct-interval (mjr_poly_deflate p a) a b eps)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-count-distinct-negative-zero-positive (p &optional eps)
  "Use mjr_poly_scount-sturm to compute the total number of distinct, real roots of the polynomial in $(a,b)$

mjr_poly_zap-zero-roots and EPS are used to count zero roots (mjr_poly_2primitive is recommended)."
  (multiple-value-bind (poly nzr) (mjr_poly_zap-zero-roots p eps)
    (let ((u  (ceiling (mjr_poly_root-bound-all poly)))
          (ss (mjr_poly_seq-make-sturm-canonical poly)))
      (values (mjr_poly_scount-sturm p (- -1 u) 0        ss)
              nzr
              (mjr_poly_scount-sturm p 0        (+ 1 u)  ss)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-count-distinct-real (p)
  "Use mjr_poly_scount-sturm and mjr_poly_root-bound-all to compute the total number of distinct, real roots of the polynomial."
  (let ((u (ceiling (mjr_poly_root-bound-all p))))
    (mjr_poly_scount-sturm p (- -1 u) (+ 1 u))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-structure (poly &optional use-sturm-sequence)
  "Return a list of vectors with each vector containing possible counts for negative, zero, positive, and complex roots for POLY."
  (let* ((poly (mjr_poly_simplify poly))
         (num-roots  (1- (length poly))))
    (multiple-value-bind (max-pos-roots num-zero-roots max-neg-roots) (mjr_poly_scount-descartes poly)
      (multiple-value-bind (min-neg-roots num-zero-roots-s min-pos-roots) (if use-sturm-sequence
                                                                              (mjr_poly_root-count-distinct-negative-zero-positive poly)
                                                                              (values 0 nil 0))
        (if (and num-zero-roots-s (not (= num-zero-roots-s num-zero-roots)))
            (error "mjr_poly_root-count-distinct-negative-zero-positive: Inconsistent zero root counts!"))
        ;;(format 't "~10d ~10d ~10d~%" max-pos-roots num-zero-roots max-neg-roots)
        ;;(format 't "~10d ~10d ~10d~%" min-pos-roots num-zero-roots-s min-neg-roots)
        (loop for num-not0-roots from (- num-roots num-zero-roots) downto (+ min-pos-roots min-neg-roots) by 2
              append (loop for num-neg-roots from max-neg-roots downto min-neg-roots by 2
                           append (loop for num-pos-roots from max-pos-roots downto min-pos-roots by 2
                                        for num-cmplx-roots = (- num-roots num-not0-roots num-zero-roots)
                                        when (and (<= 0 num-cmplx-roots)
                                                  (= num-not0-roots (+ num-pos-roots num-neg-roots)))
                                        collect (vector num-neg-roots
                                                        num-zero-roots
                                                        num-pos-roots
                                                        (- num-roots num-not0-roots num-zero-roots)))))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-structure-print (rstruct)
  "Print out all the list returned from mjr_poly_root-structure in a nice humanly readable table."
  (format 't "  Negative       Zero   Positive    Complex~%")
  (loop for cnts in rstruct
        do (format 't "~10d ~10d ~10d ~10d~%" (aref cnts 0) (aref cnts 1) (aref cnts 2) (aref cnts 3))
        finally (return rstruct)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-solve-rational (poly &key only-integer-roots)
  "Find rational roots of a polynomial.

Morally speaking this only works on a polynomials with integer coefficients; however, this function uses MJR_POLY_2PRIMITIVE to transform POLY into one with
integer coefficients if required.

Returns:
  * A list of found roots
  * The polynomial actually solved (we call this IPOLY) -- the first return of MJR_POLY_2PRIMITIVE on POLY
  * The logical factor required to transform POLY into IPOLY -- the second return of MJR_POLY_2PRIMITIVE on POLY
  * the result of deflating IPOLY with each root found"
  (cond ((not (vectorp poly))       (error "mjr_poly_root-solve-rational: POLY must be a vector"))
        ((not (every #'realp poly)) (error "mjr_poly_root-solve-rational: POLY must have real coefficients.")))
  (multiple-value-bind (ipoly ipolyf) (mjr_poly_2primitive poly)
    (let ((plen (length poly)))
      (if (zerop (aref ipoly (1- plen))) ;; Here we take care of the case with zero constant coefficient....
          (let ((last-zero-idx (loop for i from (1- plen) downto 0
                                     until (not (zerop (aref ipoly i)))
                                     minimize i)))
            (if (zerop last-zero-idx)
                (error "mjr_poly_root-solve-rational: POLY is zero -- infinitely many roots!!!")
                (multiple-value-bind (roots ipoly-throw-away ipolyf-throw-away wpoly)
                    (mjr_poly_root-solve-rational (subseq ipoly 0 last-zero-idx))
                  (declare (ignore ipoly-throw-away ipolyf-throw-away)) ;; Just keep compiler from squawking
                  (values (sort (append (make-list (- plen last-zero-idx) :initial-element 0) roots) #'<)
                          ipoly ;; The original -- not ipoly-throw-away..
                          ipolyf;; The original -- not ipolyf-throw-away..
                          wpoly))))
          ;; constant coefficient not zero.
          (let* ((wpoly            (copy-seq ipoly))
                 (first-coef       (abs (aref ipoly 0)))
                 (last-coef        (abs (aref ipoly (1- plen))))
                 (first-factors    (mjr_prime_all-factors first-coef))
                 (last-factors     (mjr_prime_all-factors last-coef))
                 (proot-candidates (remove-duplicates (loop for first-factor in first-factors
                                                            append (loop for last-factor in last-factors
                                                                         for c = (/ last-factor first-factor)
                                                                         when (or (not only-integer-roots) (integerp c))
                                                                         collect c))))
                 (root-candidates  (sort (append (mapcar #'- proot-candidates) proot-candidates) #'<)))
            (values
             (loop for x in root-candidates
                   until (> 2 (length wpoly))
                   append (loop until (not (multiple-value-bind (dfpoly pval) (mjr_poly_deflate wpoly x)
                                             (and (= 0 pval) (setq wpoly dfpoly))))
                                collect x))
             ipoly
             ipolyf
             wpoly))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-solve-integer (poly)
  "Find integer roots of a polynomial.

Morally speaking this only works on a polynomials with integer coefficients; however, this function uses MJR_POLY_2PRIMITIVE to transform POLY into one with
integer coefficients if required.

Returns:
  * A list of found roots
  * The polynomial actually solved (we call this IPOLY) -- the first return of MJR_POLY_2PRIMITIVE on POLY
  * The logical factor required to transform POLY into IPOLY -- the second return of MJR_POLY_2PRIMITIVE on POLY
  * the result of deflating IPOLY with each root found"
  (mjr_poly_root-solve-rational poly :only-integer-roots 't))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-search-bsect (poly x0 x1 &rest kw-args &key xeps yeps max-itr show-progress)
  "Use bsect's method to find a root of POLY between X0 and X2.
See the documentation for MJR_NLEQ_ROOT-BSECT"
  (and xeps yeps max-itr show-progress)
  (flet ((fp  (x) (mjr_poly_eval poly x)))
    (apply #'mjr_nleq_root-bsect #'fp x0 x1 kw-args)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-search-newton (poly x0 &rest kw-args &key xeps yeps max-itr show-progress)
  "Use newton's method to find a root of POLY near (hopefully anyhow) X0.
See the documentation for MJR_NLEQ_ROOT-NEWTON"
  (and xeps yeps max-itr show-progress) ;; Variables not directly used
  (flet ((fdf (x) (mjr_poly_eval-poly-and-first-n-derivatives poly x 1)))
    (apply #'mjr_nleq_root-newton #'fdf x0 kw-args)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-search-laguerre (poly x0 &rest kw-args &key xeps yeps max-itr show-progress)
  "Use laguerre's method to find a root of POLY near (hopefully anyhow) X0.
See the documentation for MJR_NLEQ_ROOT-LAGUERRE"
  (and xeps yeps max-itr show-progress) ;; Variables not directly used
  (flet ((fdfddf (x) (mjr_poly_eval-poly-and-first-n-derivatives poly x 2)))
    (apply #'mjr_nleq_root-laguerre #'fdfddf (mjr_poly_degree poly) x0 kw-args)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-separate-largest-real (poly &key (xeps 1) max-itr yeps show-progress)
  "Use the canonical sturm-sequence to find an interval containing the largest real root and no other roots.

An open interval $(a,b)$ containing the largest real root (and no other roots) is found.  This interval will be no wider than :XEPS, and will contain no roots
other than the largest real one.  The return of this function will two values corresponding to the left and right endpoints of the interval.  In the case of a
polynomial with no real roots, then (values nil nil) will be returned. If the polynomial is rational, then the interval endpoints will be as well.

When removed from the vulgarities of real world computing, this algorithm is guaranteed to work; however, some things can go wrong with real world
implementations.  When the inputs are rational, and enough memory is available for the necessary precision, the algorithm will always work unless the :MAX-ITR
count is violated (the error is 'Maximum iteration count exceeded!').  When using floating point arithmetic various things can go wrong, but when the error
'Could not find a non-root' occurs a lower :YEPS (the only use of this argument is to detect when the polynomial is zero at prospective new interval
endpoints) might help."
  (let ((d (mjr_poly_degree poly)))
    (cond ((some #'complexp poly) (error "mjr_poly_root-separate-largest-real: POLY must be a non-complex coefficients!"))
          ((< d 2)                (error "mjr_poly_root-separate-largest-real: POLY must be of at least degree two!")))
    (let* ((rb (max 1 (ceiling (mjr_poly_root-bound-all poly))))
           (ss (mjr_poly_seq-make-sturm-canonical poly))
           (b  (* 11/10 rb))
           (a  (- b))
           (va (mjr_poly_count-sign-changes (mjr_poly_seq-eval ss a)))
           (vb (mjr_poly_count-sign-changes (mjr_poly_seq-eval ss b))))
      (loop for itr from 1
            for n = (- va vb)
            do (if show-progress (format 't "~5d ~3d ~40a ~40a ~40a ~%" itr n a "" b))
            do (if (and max-itr (> itr max-itr)) (error "mjr_poly_root-separate-largest-real: Maximum iteration count exceeded!"))
            do (if (> n 0)
                   (loop with delta = (/ (- b a) (* 2 (+ d 1)))
                         for i from 0 upto d
                         for c = (/ (+ a b) 2) then (- c delta)
                         when (mjr_cmp_not-zerop (mjr_poly_eval poly c) yeps)
                         do (return (let ((vc (mjr_poly_count-sign-changes (mjr_poly_seq-eval ss c))))
                                      (if (> (- vc vb) 0)
                                          (setq a c
                                                va vc)
                                          (setq b c
                                                vb vc))))
                         do (if show-progress (format 't "~5d ~3a ~40a ~40a ~40a ~%" "" "" "" c "")) ;; Only show after first try
                         finally (error "mjr_poly_root-separate-largest-real: Could not find a non-root!"))
                   (return-from mjr_poly_root-separate-largest-real (values nil nil)))
            until (and (= n 1)
                       (< (- b a) xeps)))
      (values a b))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-search-largest-real (poly &rest rest)
  "Return the largest real root or NIL if the polynomial has no real roots.

When successful this function returns the center of the interval found by mjr_poly_root-separate-largest-real.  All arguments are passed directly to
mjr_poly_root-separate-largest-real, so consult that function for usage information."
  (multiple-value-bind (a b) (apply #'mjr_poly_root-separate-largest-real poly rest)
    (if (and a b)
        (/ (+ a b) 2))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-solve-low-degree (poly &optional eps)
  "Return list of root(s) to the small degree (<=3) polynomial.
Multiple roots are listed multiple times.
The argument EPS is used to avoid division by zero."
  (if (<= (mjr_poly_degree poly) 3)
      (let* ((p3   (mjr_poly_coeff poly 3))
             (p2   (mjr_poly_coeff poly 2))
             (p1   (mjr_poly_coeff poly 1))
             (p0   (mjr_poly_coeff poly 0)))
        (if (mjr_chk_!=0 p3 eps)
            (if (every #'realp (list p3 p2 p1 p0)) ;; ............................................................... Cubic(deg=3)
                (let* ((m  (+ (* 2 p2 p2 p2) (* -9 p3 p2 p1) (* 27 p3 p3 p0))) ;; ..................................... Real Poly
                       (n  (* 4 (expt (- (* p2 p2) (* 3 p3 p1)) 3)))
                       (r  (mjr_numu_sqrt (- (* m m) n)))
                       (p  (mjr_numu_cubert (/ (+ m r) 2)))
                       (q  (mjr_numu_cubert (/ (- m r) 2)))
                       (i1 (/ (+ 1 (complex 0 (sqrt 3))) 2))
                       (i2 (/ (- 1 (complex 0 (sqrt 3))) 2)))
                  (list (/ (+ p2 p q) (* -3 p3))
                        (/ (+ (- p2) (* i1 p) (* i2 q)) (* 3 p3))
                        (/ (+ (- p2) (* i2 p) (* i1 q)) (* 3 p3))))
                (let* ((a (/ p2 p3)) ;; ............................................................................... Complex
                       (b (/ p1 p3))
                       (c (/ p0 p3))
                       (q  (/ (- (* a a) (* 3 b)) 9))
                       (r  (/ (+ (* 2 a a a) (* -9 a b) (* 27 c)) 54))
                       (r2 (* r r))
                       (q3 (expt q 3)))
                  (if (and (not (complexp r)) (not (complexp q)) (< r2 q3))
                      (let* ((an  (acos (/ r (mjr_numu_sqrt q3)))))
                        (list (- (* -2 (mjr_numu_sqrt q) (cos (/ an 3)))              (/ a 3))
                              (- (* -2 (mjr_numu_sqrt q) (cos (/ (+ an (* 2 pi)) 3))) (/ a 3))
                              (- (* -2 (mjr_numu_sqrt q) (cos (/ (- an (* 2 pi)) 3))) (/ a 3))))
                      (let* ((ca   (if (>= (realpart (* (conjugate r) (mjr_numu_sqrt (- r2 q3)))) 0)
                                       (mjr_numu_cubert (- (+ r (mjr_numu_sqrt (- r2 q3)))))
                                       (mjr_numu_cubert (- (- r (mjr_numu_sqrt (- r2 q3)))))))
                             (cb (if (mjr_chk_!=0 ca)
                                     (/ q ca)
                                     0)))
                        (list (- (+ ca cb) (/ a 3))
                              (+ (- (* -1/2 (+ ca cb)) (/ a 3)) (* #C(0 1) (/ (sqrt 3) 2) (- ca cb)))
                              (- (- (* -1/2 (+ ca cb)) (/ a 3)) (* #C(0 1) (/ (sqrt 3) 2) (- ca cb))))))))
            (if (mjr_chk_!=0 p2 eps)
                (let* ((dsr (mjr_numu_sqrt (- (* p1 p1) (* 4 p2 p0)))) ;; ........................................... Quadratic(deg=2)
                       (tmp (/ (mjr_cmp_abs-max (+ p1 dsr) (- p1 dsr)) -2)))
                  (if (mjr_chk_!=0 tmp eps)
                      (list (/ tmp p2) (/ p0 tmp))
                      (list 0 0)))
                (if (mjr_chk_!=0 p1 eps)
                    (list (/ (- p0) p1)) ;; ......................................................................... Linear(deg=1)
                    (if (mjr_chk_!=0 p0 eps)
                        nil
                        (error "mjr_poly_root-solve-low-degree: Infinitely many solutions!"))))))
      (error "mjr_poly_root-solve-low-degree: Polynomial order is greater than 3!")))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_make-lagrange (points 1point &optional (pval 1))
  "Return the unique polynomial that is PVAL (default 1) on the 1point'th element (zero indexed) of POINTS, and zero all the
  other elements of POINTS"
  (let ((poly 1)
        (divs 1)
        (1p   (elt points 1point)))
    (if (vectorp points)
        (loop for i from 0
              for p across points
              when (not (= i 1point))
              do (setq poly (mjr_poly_* poly (vector 1 (- p)))
                       divs (* divs (- 1p p))))
        (loop for i from 0
              for p in points
              when (not (= i 1point))
              do (setq poly (mjr_poly_* poly (vector 1 (- p)))
                       divs (* divs (- 1p p)))))
    (mjr_poly_* poly (/ pval divs))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-lagrange (points &optional (pval 1))
  "Return a list of polynomials such that the i'th one has value PVAL on the i'th element of POINTS and zero on all other elements"
  (loop for i from 0 upto (1- (length points))
        collect (mjr_poly_make-lagrange points i pval)))
;; MJR TODO NOTE <2012-11-12 22:30:11 CST> mjr_poly_seq-make-lagrange: Optimize this function!

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-chebyshev (num &optional (kind 1))
  "Compute first NUM Chebyshev Polynomials (as a vector of vectors).
kind==1 corresponds to Chebyshev Polynomials of the first kind $T_n$.
kind==2 corresponds to Chebyshev Polynomials of the second kind $U_n$."
  (if (= num 1)
      (make-array 1 :initial-contents '(#(1)))
      (loop with p = (make-array num)
            for i from 2 upto (1- num)
            initially (setf (aref p 0) #(1)
                            (aref p 1) (if (= kind 1) #(1 0) #(2 0)))
            finally (return p)
            do (setf (aref p i) (mjr_poly_- (mjr_poly_* 2 #(1 0) (aref p (- i 1))) (aref p (- i 2)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_make-chebyshev (power &optional (kind 1))
  "Compute Chebyshev Polynomial of the given power.
kind==1 corresponds to Chebyshev Polynomials of the first kind $T_n$.
kind==2 corresponds to Chebyshev Polynomials of the second kind $U_n$."
  (let* ((p0 #(1))
         (p1 (if (= kind 1) #(1 0) #(2 0))))
    (case power
      (0         p0)
      (1         p1)
      (otherwise (loop for i from 2 upto (1+ power)
                       for pa = p0 then pb
                       for pb = p1 then pn
                       for pn = (mjr_poly_- (mjr_poly_* 2 #(1 0) pb) pa)
                       when (= power i) do (return pn))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_make-legendre (n)
  "Compute the Nth legendre polynomial using $$\\frac{1}{2^n n!}\\cdot\\frac{d^n}{dx^n}(x^2-1)^n$$"
  (mjr_poly_* (mjr_poly_diff (mjr_poly_iexpt #(1 0 -1) n) n) (/ (* (expt 2 n) (mjr_combe_! n)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-legendre (n)
  "Return a list of the first n legendre polynomials (starting with the 0'th one)"
  (loop for i from 0 upto n
        collect (mjr_poly_make-legendre i)))
;; MJR TODO NOTE <2012-11-12 22:40:48 CST> mjr_poly_seq-make-legendre: Optimize this function!

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_make-bernstein (nu n)
  "The nu'th Bernstein basis polynomial of degree n."
  (let ((poly (make-array (1+ n) :initial-element 0))
        (tmp  (mjr_poly_* (mjr_combe_comb n nu)
                          (mjr_poly_iexpt #(-1 1) (- n nu)))))
    (dotimes (i (length tmp) poly)
      (setf (aref poly i) (aref tmp i)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-bernstein (n &key show-progress)
  "The Bernstein basis polynomials of degree n (all n of them in a list)."
  (reverse (loop with xn = (make-array (1+ n) :initial-element 0)
                 initially (setf (aref xn 0) 1)
                 for nu from n downto 0
                 for c = (mjr_combe_comb n nu)
                 for p1 = (subseq xn 0 (1+ nu))
                 for p2 = #(1) then (mjr_poly_* #(-1 1) p2)
                 do (if show-progress (format 't "~4d: ~10d ~4d ~80a ~%" nu c (mjr_poly_degree p1) p2))
                 collect (mjr_poly_* c p1 p2))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_imul (poly n)
  "Multiply POLY times the integer N.

This function is provided for consistency with other packages that implement arithmetic for various rings -- after all integer
multiplication is not terribly interesting for polynomials over Z, Q, R, or C."
  (mjr_poly_* n poly))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_root-solve-search-deflate (poly &key
                                                  (xeps 1.0d-5) (yeps 1.0d-5)
                                                  (xepsr 1.0d-10) (yepsr 1.0d-10) (refine-steps 10) (only-refined nil)
                                                  (yepsp nil)
                                                  (retry-factor 4) (float-type 0.0d0)
                                                  (show-warnings nil) (show-progress nil))
  "Find floating point approximations for all of the roots of a polynomial.
Three values are returned: list of roots found, list of poly values on roots, and the part of the poly that couldn't be solved.

             _____THIS FUNCTION SHOULD BE CONSIDERED EXPERIMENTAL_____

This function works most of the time, but the roots are generally only accurate to a few decimal points.  Retrying the function after failure can lead to
success (different random seeds).  Lowering the :XEPS and :YEPS can help in initial root location.  If :ONLY-REFINED is used, then increasing :RETRY-FACTOR
may be required.  On LISPs that support it, setting :FLOAT-TYPE to 0.0L0 may increase accuracy and help with convergence problems due to round off errors.  As
a last resort, the :YEPSP argument may be set to a value larger than :YEPS in cases where numerical instability leads to deflated polynomials that have roots
which are not within :YEPS of being a root of the original polynomial but are within :YEPSP of being a root (i.e.  the absolute value of the original
polynomial at the proposed root is less than :YEPSP).  This last option can lead not-quite-roots that are quite far off, but still useful.

The algorithm:
   0) Set w<-p
   1) Use laguerre to find a root xb of w (:xeps & :yeps)
   2) Use laguerre to refine xb (:refine-steps, :xepsr, & :xepsr)
   3) If refinement leads to a non-root of w and :only-refined is non-NIL, then NEXT
   4) If p is real and xb is complex, check to see if re(xb) is a root of p and w.  If so xb<-re(xb)
   5) If xb is a root of p and w, then we deflate:
      a) If p is real and xb is complex, divide out xb and conj(xb) and update w
      b) otherwise just divide out xb
   6) Jump to 1) if w is deg 1 or higher and we haven't gone beyond (retry-factor)*deg(poly) tries"
  (let ((rpoly (every #'realp poly))
        (wpoly (copy-seq poly))
        (roots nil)
        (pvals nil)
        (yepsp (or yepsp yeps))
        (trys  (* retry-factor (length poly))))
    (loop with refined = nil
          for guess = (mjr_prng_float-co (float -10 float-type) (float 10 float-type))
          for itr from 1 upto trys
          while (< 1 (length wpoly))
          do (multiple-value-bind (x-bst w-bst ex-why)
                 (progn (if show-progress (format 't "Root search: ~f~%" guess))
                        (mjr_poly_root-search-laguerre wpoly guess :yeps yeps :xeps xeps :show-progress show-progress))
               (if show-progress (format 't "Root candidate: ~a~%" x-bst))
               ;; Refine the root.
               (setq refined nil)
               (if refine-steps
                   (multiple-value-bind (x-bst2 p-bst2 ex-why2)
                       (progn (if show-progress (format 't "Root refine:~%"))
                              (mjr_poly_root-search-laguerre poly x-bst :yeps yepsr :xeps xepsr :max-itr refine-steps :show-progress show-progress))
                     (if show-progress (format 't "Refinement candidate: ~a~%" x-bst2))
                     (if (mjr_eps_=0 p-bst2 yepsp)
                         (let ((w-bst2 (mjr_poly_eval wpoly x-bst2)))
                           (if (mjr_eps_=0 w-bst2 yeps)
                               (progn (if show-progress (format 't "Refinement success~%"))
                                      (setf x-bst   x-bst2
                                            w-bst   w-bst2
                                            ex-why  ex-why2
                                            refined 't))
                               (if show-warnings (warn "mjr_poly_root-solve-search-deflate: laguerre refinement failed to find new root: ~a" x-bst))))
                         (if show-warnings (warn "mjr_poly_root-solve-search-deflate: laguerre refinement failed to converge: ~a" x-bst)))))
               (if (or refined (null only-refined))
                   (progn
                     ;; Snip off the complex bit if the real bit is a good root
                     (if (and rpoly (complexp x-bst))
                         (let* ((x-bst2 (realpart x-bst))
                                (w-bst2 (mjr_poly_eval wpoly x-bst2)))
                           (if (and (mjr_eps_=0 w-bst2 yeps)
                                    (mjr_eps_=0 (mjr_poly_eval poly x-bst2) yepsp))
                               (progn (if show-progress (format 't "Root real-ifyed~%"))
                                      (setf x-bst  x-bst2
                                            w-bst  w-bst2)))))
                     ;; Snip off the fractional bit if this results in a good root
                     (let* ((x-bst2 (complex (round (realpart x-bst)) (round (imagpart x-bst))))
                            (w-bst2 (mjr_poly_eval wpoly x-bst2)))
                           (if (and (mjr_eps_=0 w-bst2 yeps)
                                    (mjr_eps_=0 (mjr_poly_eval poly x-bst2) yepsp))
                               (progn (if show-progress (format 't "Root int-ifyed~%"))
                                      (setf x-bst  x-bst2
                                            w-bst  w-bst2))))
                     ;; Check the "roots", and deflate if things look good
                     (if (mjr_eps_=0 w-bst yeps)
                         (let ((p-bst (mjr_poly_eval poly x-bst)))
                           (if (mjr_eps_=0 p-bst yepsp)
                               (let ((wpoly2 (mjr_poly_deflate wpoly x-bst)))
                                  (setf roots (cons x-bst roots)
                                        pvals (cons p-bst pvals))
                                  (if (and rpoly (complexp x-bst))
                                      (let* ((x-bst-c (conjugate x-bst)))
                                        (if (mjr_eps_=0 (mjr_poly_eval wpoly2 x-bst-c) yeps)
                                            (setf roots (cons x-bst-c roots)
                                                  pvals (cons (mjr_poly_eval poly x-bst-c) pvals)
                                                  wpoly (mjr_poly_truncate wpoly
                                                                           (vector 1
                                                                                   (* -2 (realpart x-bst))
                                                                                   (mjr_numu_abssqr x-bst))))
                                            (progn
                                              (if show-warnings (warn "mjr_poly_root-solve-search-deflate: Numerical instability detected: real poly but conjugate not a root!: ~a" x-bst))
                                              (setf wpoly wpoly2))))
                                      (setf wpoly wpoly2)))
                               (if show-warnings (warn "mjr_poly_root-solve-search-deflate: Numerical instability detected: root of deflated poly not a root of original poly!"))))
                         (if show-warnings (warn "mjr_poly_root-solve-search-deflate: laguerre failed to converge: ~a ~a ~a" x-bst w-bst ex-why))))
                   (if show-warnings (warn "mjr_poly_root-solve-search-deflate: Ignoring unrefined root: ~a" x-bst))))
          do (if (and (= itr (truncate trys 2)) show-warnings) (warn "mjr_poly_root-solve-search-deflate: laguerre having difficulty converging"))
          do (if show-progress (format 't "ROOT:  ~a~%" roots)))
    (values roots pvals wpoly)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_find-integer-factor (poly degree-or-sample-points &key show-progress)
  "Return an integer polynomial factor of the integer polynomial specified by POLY of the specified size or NIL if none exist.

If a factor was found, then a second return value returned is POLY divided by the factor that was found.

If degree-or-sample-points is a vector, then the points will be used in Kronecker's method.  As the value of the polynomial on each of these points must be
factored, the polynomial must not be zero on any of the points in the vector.

The size of the factors may be specified as an integer indicating the degree of the desired factors or a vector of sample points to be used by Kronecker's
method (in which case the degree of the factors will be one less than the length of the list)

References:
  Victor V. Prasolov (2004); Polynomials; ISBN: 3540407146; pp 49-50"
  (let* ((degree   (if (vectorp degree-or-sample-points)
                       (1- (length degree-or-sample-points))
                       degree-or-sample-points))
         (x-points (if (vectorp degree-or-sample-points)
                       degree-or-sample-points
                       (make-array (1+ degree))))
         (y-points (if (vectorp degree-or-sample-points)
                       (map 'vector (lambda (x) (mjr_poly_eval poly x)) x-points)
                       (apply #'vector (loop with npts = 0 ;; find degree+1 points where the poly is not zero
                                             for x = 0 then (if (not (plusp x)) (1+ (abs x)) (- x))
                                             for y = (mjr_poly_eval poly x)
                                             do (if show-progress (format 't "testing sample point ~d => ~d~%" x y))
                                             when (not (zerop y))
                                             collect (progn (setf (aref x-points npts) x)
                                                            (incf npts)
                                                            y)
                                             until (= npts (1+ degree)))))))
    (if show-progress (format 't "Using x sample points: ~a~%" x-points))
    (if show-progress (format 't "Using y sample points: ~a~%" y-points))
    (mjr_combc_gen-all-cross-product (loop for i from 0
                                           for x across x-points
                                           for y across y-points
                                           for pf = (mjr_prime_all-factors (abs y))
                                           do (if show-progress (format 't "factor ~d from ~d into ~d factors~%" y x (length pf)))
                                           collect (concatenate 'vector pf (if (not (zerop i)) (mapcar #'- pf))))
                                     :func (lambda (y-vals) (let ((p (mjr_intrp_poly x-points y-vals)))
                                                              (if (and (= (1+ degree) (length p)) (every #'integerp p))
                                                                  (let ((quo (mjr_poly_divides? p poly :require-integer-quotient 't)))
                                                                    (if quo
                                                                        (return-from mjr_poly_find-integer-factor (values p quo))))))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_factor-over-integers (poly &key show-progress)
  "Return a list of integer polynomial factors of the integer polynomial specified by POLY.

Algorithm:
  1) The polynomial is transformed into a primitive, integer polynomial if required (see: MJR_POLY_2PRIMITIVE)
     The multiplicative factor required for this transformation is included in the returned factor list as a constant polynomial
     NOTE: This will rationalize the polynomial if it contains floating point coefficients (real or complex)
  2) All linear factors are found and removed
  3) Kronecker's method is applied to the remainder
     NOTE: Kronecker's method can be quite slow when the polynomial sample values are large integers with many factors.

References:
  Victor V. Prasolov (2004); Polynomials; ISBN: 3540407146; pp 49-50
  Joel S. Cohen (2003); Computer Algebra and Symbolic Computation: Mathematical Methods; ISBN: 1568811594; pp362-367"
  (multiple-value-bind (int-roots ipoly ipolyf pleft) (mjr_poly_root-solve-rational poly)
    (declare (ignore ipoly))
    (setq pleft (mjr_poly_* (/ (reduce #'* (mapcar #'denominator int-roots))) pleft))
    (loop with factors = (append (mapcar (lambda (r) (vector (denominator r) (- (numerator r)))) int-roots)
                                 (if (not (= 1 ipolyf)) (list (vector (/ ipolyf)))))
          for degree = 2 then (if gotfac 2 (1+ degree))
          for gotfac = nil
          for sample-points = (mjr_vvec_to-list (list :vvec-type :vvt-aseq :start (floor (- (truncate degree 2))) :step 1 :len (+ degree 1)))
          for pleftlen = (length pleft)
          finally (return (append factors (if (or (< 1 pleftlen) (not (= (aref pleft 0) 1))) (list pleft))))
          do (if show-progress (format 't "factoring(~d) ~a~%" degree pleft))
          do (if show-progress (format 't "factors so far ~a ~%" factors))
          while (and (<= degree (floor (/ pleftlen 2))) (< 1 (length pleft)) (or (< 1 pleftlen) (not (= (aref pleft 0) 1))))
          do (mjr_combc_gen-all-cross-product (loop for i from 0
                                                    for x in sample-points
                                                    for y = (mjr_poly_eval pleft x)
                                                    for pf = (mjr_prime_all-factors (abs y))
                                                    do (if show-progress (format 't "factor ~d from ~d into ~d factors~%" y x (length pf)))
                                                    collect (concatenate 'vector pf (if (not (zerop i)) (mapcar #'- pf))))
                                              :func (lambda (y-vals) (mjr_intrp_poly sample-points y-vals))
                                              :exit-if (lambda (p) (if (and (= (1+ degree) (length p)) (every #'integerp p))
                                                                       (loop for quo = (mjr_poly_divides? p pleft :require-integer-quotient 't)
                                                                             for j from 1
                                                                             finally (return gotfac)
                                                                             while quo
                                                                             do (progn (push p factors)
                                                                                       (setq pleft quo)
                                                                                       (setq gotfac 't)))))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_gcd-primitive (poly1 poly2)
  "Return primitive GCD of poly1 and poly2. See: MJR_POLY_GCD"
  (mjr_poly_2primitive (mjr_poly_gcd poly1 poly2)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_test-property (poly &rest pp-list)
  "Return non-nil if the given polynomial has all of the specified polynomial properties.

This is not a fast function, but it sure is handy.  A few property tests exist as separate functions because of performance.

  Symbol              Description                          elt_comp
  :PP-PRIMITIVE       Integer, primitive polynomial        exact
  :PP-SQUARE-FREE     No double factors                    fizzy
  :PP-IRREDUCIBLE     No factors                           exact
  :PP-MONIC           Leading coefficient is 1             fizzy
  :PP-NUMBER          All elements are numbers             N/A
  :PP-REAL            All elements are real                N/A
  :PP-INTEGER         All elements are integers            N/A
  :PP-RATIONAL        All elements are rational            N/A
  :PP-COMPLEX         All elements are complex             N/A
  :PP-POSITIVE        All elements numeric and positive    cmpzy
  :PP-EVEN            The polynomial is an even function   cmpzy
  :PP-ODD             The polynomial is an odd function    cmpzy
  :PP-ZERO            The zero polynomial                  cmpzy
  :PP-IDENTITY        The constant 1 polynomial            cmpzy
  :PP-CONSTANT        A constant polynomial                cmpzy"
  (let ((poly (mjr_poly_simplify poly))
        (deg  (mjr_poly_degree poly)))
    (if (cdr pp-list)
        (every (lambda (pp) (mjr_poly_test-property poly pp)) pp-list)
        (case (car pp-list)
          (:pp-complex      (every #'complexp   poly))
          (:pp-real         (every (lambda (x) (and (numberp x) (not (complexp x)))) poly))
          (:pp-rational     (every #'rationalp  poly))
          (:pp-integer      (every #'integerp   poly))
          (:pp-number       (every #'numberp    poly))
          (:pp-positive     (every #'plusp      poly))
          (:pp-monic        (mjr_cmp_= (mjr_poly_leading-coeff poly) 1))
          (:pp-odd          (and (oddp deg)
                                 (loop for c across poly
                                       for i from 0 upto deg
                                       finally (return 't)
                                       when (and (oddp i) (mjr_cmp_!=0 c))
                                       do (return nil))))
          (:pp-even         (and (evenp deg)
                                 (loop for c across poly
                                       for i from 0 upto deg
                                       finally (return 't)
                                       do (format 't "~a~%" c)
                                       when (and (oddp i) (mjr_cmp_!=0 c))
                                       do (return nil))))
          (:pp-irreducible  (if (< deg 2)
                                't
                                (if (or (not (mjr_poly_test-property poly :pp-rational))
                                        (mjr_poly_test-property poly :pp-complex))
                                    nil
                                    (error "mjr_poly_test-property: Irreducibly tests not implemented for integer or rational polynomials yet!"))))
          (:pp-primitive    (and (mjr_poly_test-property poly :pp-integer)
                                 (= 1 (reduce #'gcd poly))))
          (:pp-square-free  (mjr_poly_test-property (mjr_poly_gcd (mjr_poly_diff poly) poly) :pp-identity))
          (:pp-zero         (mjr_eps_= #(0) poly))
          (:pp-identity     (mjr_eps_= #(1) poly))
          (:pp-constant     (= 0 deg))
          (otherwise        (error "mjr_poly_test-property: Unknown property"))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-hermite (num)
  "Compute first NUM Hermite Polynomials (as a vector of vectors)."
  (loop with p = (make-array num)
        for i from 2 upto (1- num)
        initially (progn (setf (aref p 0) #(1))
                         (if (> num 1)
                             (setf (aref p 1) #(1 0))))
        finally (return p)
        do (setf (aref p i) (mjr_poly_- (mjr_poly_* #(1 0) (aref p (- i 1)))
                                        (mjr_poly_diff (aref p (- i 1)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_make-hermite (power)
  "Compute the Hermite Polynomial of the given power."
  (let* ((p0 #(1))
         (p1 #(1 0)))
    (case power
      (0         p0)
      (1         p1)
      (otherwise (loop for i from 2 upto (1+ power)
                       for pa = p0 then pb
                       for pb = p1 then pn
                       for pn = (mjr_poly_- (mjr_poly_* #(1 0) pb) (mjr_poly_diff pb))
                       when (= power i) do (return pn))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_make-laguerre (power)
  "Compute the Laguerre Polynomial of the given power."
  (let* ((p0 #(1))
         (p1 #(-1 1)))
    (case power
      (0         p0)
      (1         p1)
      (otherwise (loop for i from 2 upto (1+ power)
                       for pa = p0 then pb
                       for pb = p1 then pn
                       for pn = (mjr_poly_* (/ i)
                                            (mjr_poly_- (mjr_poly_* (mjr_poly_+ (* 2 (1- i)) #(-1 1)) pb)
                                                        (mjr_poly_* (1- i) pa)))
                       when (= power i) do (return pn))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_seq-make-laguerre (num)
  "Compute first NUM Hermite Polynomials (as a vector of vectors)."
  (loop with p = (make-array num)
        for i from 2 upto (1- num)
        initially (progn (setf (aref p 0) #(1))
                         (if (> num 1)
                             (setf (aref p 1) #(-1 1))))
        finally (return p)
        do (setf (aref p i) (mjr_poly_* (/ i)
                                        (mjr_poly_- (mjr_poly_* (mjr_poly_+ (* 2 (1- i)) #(-1 1)) (aref p (1- i)))
                                                    (mjr_poly_* (1- i) (aref p (- i 2))))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_cubic-depress (poly)
  "Return the depressed form of the given cubic polynomial.

If $p(x)=ax^3+bx^2+cx+d$, then divide $p$ by $a$ and substitute $t=\frac{b}{-3a}$."
  (let* ((cpoly (mjr_poly_simplify poly))
         (deg   (1- (length cpoly))))
    (if (= 3 deg)
        (let* ((a (aref cpoly 0))
               (a2 (* a a))
               (b (aref cpoly 1))
               (b2 (* b b))
               (c (aref cpoly 2))
               (d (aref cpoly 3)))
          (vector 1
                  0
                  (/ (- (* 3 a c) b2) (* 3 a2))
                  (/ (+ (* 2 b2 b) (* -9 a b c) (* 27 a2 d)) (* 27 a2 a))))
        (error "mjr_poly_cubic-depress: Polynomial must be a cubic!"))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_tschirnhaus-3-2 (poly)
  "Tschirnhaus transformation (wipe out the power 2 term of a degree 3 polynomial)

If $p(x)=ax^3+bx^2+cx+d$, then substitute $t=\frac{b}{3a}$."
  (mjr_poly_simplify (mjr_poly_subst (vector 1 (/ (mjr_poly_coeff poly 2) (* -3 (mjr_poly_coeff poly 3)))) poly)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_discriminant-low-degree (poly)
  "Return discriminant of a non-linear, small degree (<=5) polynomial via a direct formula"
  ;; Maxima Code to compute discriminant
  ;;     degree(f,x) := hipow(f,x)$
  ;;     discriminant(f,x) := block([n:degree(f,x),fp,d,an],
  ;;      fp : diff(f,x),
  ;;      an : coeff(f,x,n),
  ;;      w : mod(n,4),
  ;;      w : w*(w-1)/2,
  ;;      d : (-1)^w/an,
  ;;      d : d*resultant(f,fp,x)
  ;;     )$
  (let* ((cpoly (mjr_poly_simplify poly))
         (deg   (1- (length poly)))
         (f     (mjr_poly_coeff cpoly 5))
         (e     (mjr_poly_coeff cpoly 4))
         (d     (mjr_poly_coeff cpoly 3))
         (c     (mjr_poly_coeff cpoly 2))
         (b     (mjr_poly_coeff cpoly 1))
         (a     (mjr_poly_coeff cpoly 0)))
    (case deg
      (5          (let* ((a2 (* a  a)) (a3 (* a2 a))
                         (b2 (* b  b)) (b3 (* b2 b)) (b4 (* b3 b))
                         (c2 (* c  c)) (c3 (* c2 c)) (c4 (* c3 c))
                         (d2 (* d  d)) (d3 (* d2 d)) (d4 (* d3 d))
                         (e2 (* e  e)) (e3 (* e2 e)) (e4 (* e3 e))
                         (f2 (* f  f)) (f3 (* f2 f)))
                    ;; 3125*a^4*f^4-2500*a^3*b*e*f^3-3750*a^3*c*d*f^3+2000*a^2*b^2*d*f^3+2250*a^2*b*c^2*f^3-1600*a*b^3*c*f^3+256*b^5*
                    ;; f^3+2000*a^3*c*e^2*f^2-50*a^2*b^2*e^2*f^2+2250*a^3*d^2*e*f^2-2050*a^2*b*c*d*e*f^2+160*a*b^3*d*e*f^2-900*a^2*c^3*
                    ;; e*f^2+1020*a*b^2*c^2*e*f^2-192*b^4*c*e*f^2-900*a^2*b*d^3*f^2+825*a^2*c^2*d^2*f^2+560*a*b^2*c*d^2*f^2-128*b^4*d^2*
                    ;; f^2-630*a*b*c^3*d*f^2+144*b^3*c^2*d*f^2+108*a*c^5*f^2-27*b^2*c^4*f^2-1600*a^3*d*e^3*f+160*a^2*b*c*e^3*f-36*a*b^3*
                    ;; e^3*f+1020*a^2*b*d^2*e^2*f+560*a^2*c^2*d*e^2*f-746*a*b^2*c*d*e^2*f+144*b^4*d*e^2*f+24*a*b*c^3*e^2*f-6*b^3*c^2*e^2*
                    ;; f-630*a^2*c*d^3*e*f+24*a*b^2*d^3*e*f+356*a*b*c^2*d^2*e*f-80*b^3*c*d^2*e*f-72*a*c^4*d*e*f+18*b^2*c^3*d*e*f+108*a^2*
                    ;; d^5*f-72*a*b*c*d^4*f+16*b^3*d^4*f+16*a*c^3*d^3*f-4*b^2*c^2*d^3*f+256*a^3*e^5-192*a^2*b*d*e^4-128*a^2*c^2*e^4+144*a*
                    ;; b^2*c*e^4-27*b^4*e^4+144*a^2*c*d^2*e^3-6*a*b^2*d^2*e^3-80*a*b*c^2*d*e^3+18*b^3*c*d*e^3+16*a*c^4*e^3-4*b^2*c^3*e^3-
                    ;; 27*a^2*d^4*e^2+18*a*b*c*d^3*e^2-4*b^3*d^3*e^2-4*a*c^3*d^2*e^2+b^2*c^2*d^2*e^2
                    (+ (* -3750 c d a3 f3)   (* -2500 b e a3 f3)   (* -2050 b c d e a2 f2) (* -1600 a c b3 f3)  (* -1600 d f a3 e3)
                       (* -900 b a2 d3 f2)   (* -900 e a2 c3 f2)   (* -746 a c d f b2 e2)  (* -630 a b d c3 f2) (* -630 c e f a2 d3)
                       (* -192 b d a2 e4)    (* -192 c e b4 f2)    (* -128 a2 c2 e4)       (* -128 b4 d2 f2)    (* -80 a b d c2 e3)
                       (* -80 c e f b3 d2)   (* -72 a b c f d4)    (* -72 a d e f c4)      (* -50 a2 b2 e2 f2)  (* -36 a f b3 e3)
                       (* -27 a2 d4 e2)      (* -27 b2 c4 f2)      (* -27 b4 e4)           (* -6 a b2 d2 e3)    (* -6 f b3 c2 e2)
                       (* -4 a c3 d2 e2)     (* -4 f b2 c2 d3)     (* -4 b2 c3 e3)         (* -4 b3 d3 e2)      (* 16 a f c3 d3)
                       (* 16 a c4 e3)        (* 16 f b3 d4)        (* 18 a b c d3 e2)      (* 18 c d b3 e3)     (* 18 d e f b2 c3)
                       (* 24 a b f c3 e2)    (* 24 a e f b2 d3)    (* 108 a c4 c f2)       (* 108 f a2 d4 d)    (* 144 a c b2 e4)
                       (* 144 c a2 d2 e3)    (* 144 d f b4 e2)     (* 144 d b3 c2 f2)      (* 160 a d e b3 f2)  (* 160 b c f a2 e3)
                       (* 256 a3 e4 e)       (* 256 b4 b f3)       (* 356 a b e f c2 d2)   (* 560 a c b2 d2 f2) (* 560 d f a2 c2 e2)
                       (* 825 a2 c2 d2 f2)   (* 1020 a e b2 c2 f2) (* 1020 b f a2 d2 e2)   (* 2000 c a3 e2 f2)  (* 2000 d a2 b2 f3)
                       (* 2250 b a2 c2 f3)   (* 2250 e a3 d2 f2)   (* 3125 a3 a f3 f)      (* b2 c2 d2 e2))))
      (4          (let* ((a2 (* a  a))
                         (b2 (* b  b)) (b3 (* b2 b))
                         (c2 (* c  c)) (c3 (* c2 c))
                         (d2 (* d  d)) (d3 (* d2 d))
                         (e2 (* e  e)))
                    ;; 256*a^3*e^3-192*a^2*b*d*e^2-128*a^2*c^2*e^2+144*a*b^2*c*e^2-27*b^4*e^2+144*a^2*c*d^2*e-6*a*b^2*d^2*e-80*a*b*c^2*d*e+
                    ;; 18*b^3*c*d*e+16*a*c^4*e-4*b^2*c^3*e-27*a^2*d^4+18*a*b*c*d^3-4*b^3*d^3-4*a*c^3*d^2+b^2*c^2*d^2
                    (+ (* -192 b d a2 e2)    (* -128 a2 c2 e2)     (* -80 a b d e c2)      (* -27 a2 d3 d)      (* -27 b3 b e2)
                       (* -6 a e b2 d2)      (* -4 a c3 d2)        (* -4 e b2 c3)          (* -4 b3 d3)         (* 16 a e c3 c)
                       (* 18 a b c d3)       (* 18 c d e b3)       (* 144 a c b2 e2)       (* 144 c e a2 d2)    (* 256 a2 a e2 e)
                       (* b2 c2 d2))))
      (3          (let* ((b2 (* b  b))
                         (c2 (* c  c)))
                    ;; -27*a^2*d^2-(4*b^3-18*a*b*c)*d-4*a*c^3+b^2*c^2
                    (+ (* -27 a a d d)       (* -4 a c2 c)         (* 18 a b c d)          (* -4 b2 b d)        (* b2 c2))))
      (2          (let ((b2 (* b b)))
                    ;;  b^2-4*a*c
                    (- b2 (* 4 a c))))
      (otherwise  (error "mjr_poly_discriminant-low-degree: Only polynomials of degree 2, 3, 4, and 5 are supported!")))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_resultant (poly1 poly2)
  "Return resultant the two polynomials via a recursive GCD-like computation."
  (labels ((resultant (poly1 poly2)
             (let ((deg1 (mjr_poly_degree poly1))
                   (deg2 (mjr_poly_degree poly2)))
               (if (> deg1 deg2)
                   (mjr_poly_* (if (evenp (* deg1 deg2)) 1 -1) (mjr_poly_resultant poly2 poly1))
                   (let ((lc (mjr_poly_leading-coeff poly1)))
                     (if (zerop deg1)
                         (vector (expt lc deg2))
                         (let ((r (mjr_poly_rem poly2 poly1)))
                           (if (mjr_poly_zerop r)
                               (vector 0)
                               (mjr_poly_* (expt lc (- deg2 (mjr_poly_degree r))) (mjr_poly_resultant poly1 r))))))))))
    (mjr_poly_leading-coeff (resultant poly1 poly2))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_discriminant-high-degree (poly)
  "Return discriminant of a non-linear polynomial via a recursive resultant computation."
  (let ((pdeg (mjr_poly_degree   poly)))
    (* (if (evenp (/ (* pdeg (1- pdeg)) 2)) 1 -1)
       (/ (mjr_poly_leading-coeff poly))
       (mjr_poly_resultant poly (mjr_poly_diff poly)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_discriminant (poly)
  "Return discriminant of a non-linear polynomial using the most efficient method available."
  (if (> (mjr_poly_degree poly) 4)
      (mjr_poly_discriminant-high-degree poly)
      (mjr_poly_discriminant-low-degree poly)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_factor-square-free (poly)
  "Return the square-free factorization of POLY as an assoc array with (poly . power).

Any factors equal to $1$ are not in the returned list.

If the input polynomial is $p\\in\\mathbb{Q}[x]$, then a square-free factorization is

     $$p=c\\prod_{j=1}^nu_j^j$$

Where $u_j\\in\\mathbb{Q}[x]$, $u^j \\vert p$, and $u^{(j+1)} \\nmid p$.

References:
  Antonio Machi (2012); Algebra for Symbolic Computation; ISBN: 978-88-470-2396-3; pp 108-109"
  ;; MJR TODO NOTE mjr_poly_factor-square-free: Add reference to Wikipedia and Cohen.
  ;; MJR TODO NOTE mjr_poly_factor-square-free: Replace with Yun's algorithm, use this for regression tests.  See: http://en.wikipedia.org/wiki/Square-free_polynomial#Yun.27s_algorithm
  ;; MJR TODO NOTE mjr_poly_factor-square-free: Implement finite field version of this for polygfp.
  (loop with d = (mjr_poly_gcd poly (mjr_poly_diff poly))
        with v = (mjr_poly_truncate poly d)
        with w = (mjr_poly_diff v)
        with z = (mjr_poly_truncate (mjr_poly_diff poly) d)
        for j from 1 upto (mjr_poly_degree poly)
        for uj = (mjr_poly_gcd v (mjr_poly_- z (mjr_poly_* j w)))
        when (not (mjr_poly_onep uj))
        collect (cons uj j)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_factor-irreducible (poly &key show-progress)
  "Return the irreducible factorization of POLY as an assoc array with (poly . power).

VERY SLOW.  Uses kronecker's method."
  ;; MJR TODO NOTE mjr_poly_factor-irreducible: Implement a modern method, and use this for regression tests.
  (loop with poly-left = poly
        with d         = 1
        with pdeg      = (mjr_poly_degree poly)
        for (f pl) = (multiple-value-list (mjr_poly_find-integer-factor poly-left d))
        if f
        collect (let ((daPow (1+ (loop initially (setf poly-left pl)
                                       for v = (mjr_poly_divides? f poly-left)
                                       while v
                                       count (setf poly-left v)))))
                  (if show-progress (format 't "left: ~5d deg: ~5d pow: ~5d factor: ~30a~%" (mjr_poly_degree poly-left) d daPow f ))
                  (cons f daPow))
        else
        do (if show-progress (format 't "left: ~5d deg: ~5d pow:     0 factor: DONE~%" (mjr_poly_degree poly-left) d))
        and
        do (incf d)
        until (mjr_poly_onep poly-left)
        until (>= d pdeg)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_poly_mahler-measure (poly &optional (method-class :solve) (method-function #'mjr_poly_root-solve-search-deflate) method-args)
  "Compute the Mahler measure of the polynomial POLY

The optional arguments to determine how the computation is preformed:

   method-class   method-function                          method-args
   :solve         #'mjr_poly_root-solve-search-deflate     nil                  <= DEFAULT
   :integrate     #'mjr_intg_simple-gauss-kronrod          '(:order 21)
   :integrate     #'mjr_intg_simple-gauss-kronrod          '(:order 31)
   :integrate     #'mjr_intg_simple-gauss-legendre         '(:order 20)
   :integrate     #'mjr_intg_glb-adp-composite-romberg     '(:the-err 1.0D-5)
   :integrate     #'mjr_intg_glb-adp-composite-trapezoidal '(:the-err 1.0D-5)

Definition

  Let $p\\in\\mathbb{C}[z]$ with $\\mathrm{deg}(p)=d$.  For definiteness, we express  $p$ like so:

  $$p(z)=\\sum_{j=0}^d a_jx^j = a_d\\prod_{j=1}^d (z-\\alpha_j)$$

  We then define the Mahler measure to be:

  $$M(p) =
    \\vert a_d\\vert\\prod_{\\vert\\alpha_j\\vert\\geq1} \\vert\\alpha_j\\vert =
    \\exp\\left(\\frac{1}{2\\pi}\\int_0^{2\\pi}\\ln(\\vert p(e^{i\\theta})\\vert)\\,\\mathrm{d}\\theta\\right)$$

Examples:

  * (mahler-measure #(1 1 0 -1 -1 -1 -1 -1 0 1 1) :integrate #'mjr_intg_simple-gauss-kronrod          '(:order 21      )) => 1.4751609391661360d0
  * (mahler-measure #(1 1 0 -1 -1 -1 -1 -1 0 1 1) :integrate #'mjr_intg_simple-gauss-kronrod          '(:order 31      )) => 0.7012142844290665d0
  * (mahler-measure #(1 1 0 -1 -1 -1 -1 -1 0 1 1) :integrate #'mjr_intg_simple-gauss-legendre         '(:order 20      )) => 1.0601232234056748d0
  * (mahler-measure #(1 1 0 -1 -1 -1 -1 -1 0 1 1) :integrate #'mjr_intg_glb-adp-composite-romberg     '(:the-err 1.0D-5)) => 1.1762811031170224d0
  * (mahler-measure #(1 1 0 -1 -1 -1 -1 -1 0 1 1) :integrate #'mjr_intg_glb-adp-composite-trapezoidal '(:the-err 1.0D-5)) => 1.1762806954644658d0
  * (mahler-measure #(1 1 0 -1 -1 -1 -1 -1 0 1 1) :solve     #'mjr_poly_root-solve-search-deflate     '(               )) => 1.1762808182599178d0"
  (case method-class
    (:solve     (let* ((roots (mapcar (lambda (x) (max 1 (abs x))) (apply method-function poly method-args)))
                       (numr  (length roots)))
                  (if (= numr (mjr_poly_degree poly))
                      (reduce #'* roots)
                      (error "mahler-measure: Could not find all roots.  Found ~d" numr))))
    (:integrate (exp (/ (apply method-function
                                 (lambda (x) (log (abs (mjr_poly_eval poly (complex (cos x) (sin x))))))
                                 :start 0.0D0
                                 :end (* 2.0d0 pi)
                                 method-args)
                        (* 2 pi))))
    (otherwise  (error "mahler-measure: method-class unsupported"))))