;;
;; A computation of the places at which flip bifurcation occurrs in the
;; quadratic family
;; Note: period 3 finished quite easily, period 4 runs forever
;; Two methods: grobner and resultant
;;


(setf order (elimination-order
	     2
	     :primary-order #'grevlex>
	     :secondary-order #'grevlex>
	     ))

(setf f (cdr (string-read-poly "[a-x^2+b*y,x,a,b]" '(x y a b) :order order)))
(setf id (cdr (string-read-poly "[x,y,a,b]" '(x y a b) :order order)))
(setf one (string-read-poly "1" '(x y a b) :order order))

(defun f-composition (n) (poly-dynamic-power f n order))

;; Flip bifurcations occur when derivative is -1 at some fixed point
;; g
(defun g (n) (subseq (mapcar #'(lambda (x y) (poly- x y order))
			     (f-composition n) id) 0 2))
(defun f-jacobian (n) (characteristic-polynomial (jacobi-matrix (f-composition n) 2 2)
						 order))
(defun flip-value (n) (poly-scalar-composition
		       (f-jacobian n)
		       (cdr (string-read-poly "[x,y,a,b,-1]" '(x y a b) :order order))
		       order))


(defun ideal (n) (cons (flip-value n) (g n)))

(defun print-ideal (n) (poly-print (cons '[ (ideal n)) '(x y a b)) (terpri))

(defun bifurcation (n)
  (mapcar #'(lambda (x) (poly-contract x 2))
	  (elimination-ideal (ideal n) 2
			     :order order
			     )))

(defun print-bifurcation (n)
  (poly-print (cons '[ (bifurcation n)) '(a b))
  (terpri))