;(setq p '(1 2 3) q '(7 8 9 ))
(setq p '(7 -7 0 0 1) q '(-7 0 3))
;(setq p '(4 0 1))
;(setq p '(9 8 7 6 5 4 3 2 1))
;(setq p '(1 1 1 1 1 1))

;;;;Polynomial algorithms

;;;Polynomial functions
(defun poly-add (p q)
  (cond ((and (null p) (null q)) nil)
	(t (cond ((null p) (cons (car q) (poly-add nil (cdr q))))
		 ((null q) (cons (car p) (poly-add (cdr p) nil)))
		 (t (cons (+ (car p) (car q)) (poly-add (cdr p) (cdr q))))))))

(defun poly-sub (p q)
  (cond ((and (null p) (null q)) nil)
	(t (cond ((null p) (cons (- (car q)) (poly-sub nil (cdr q))))
		 ((null q) (cons (car p) (poly-sub (cdr p) nil)))
		 (t (cons (- (car p) (car q)) (poly-sub (cdr p) (cdr q))))))))

(defun poly-mul (p q)
  (cond ((null p) nil)
	(t (poly-add (mapcar #'(lambda (x) (* (car p) x)) q) (cons 0 (poly-mul (cdr p) q))))))

(defun lc (p)
  (car (last p)))

(defun mulx (n p)
  (dotimes (i n p)
    (push 0 p)))

(defun poly-div (p q &optional quo)
  (cond ((< (length p) (length q)) (values quo p))
	(t (push (/ (lc p) (lc q)) quo)
	   (poly-div 
	    (butlast
	     (poly-sub p (mulx (- (length p) (length q)) (poly-mul (list (car quo)) q)))) q quo))))

(defun poly-pdiv (p q)
  (let ((degp (1- (length p)))
	(degq (1- (length q)))
	(lcq (lc q)))
	(poly-div (poly-mul (list (expt lcq (+ (- degp degq) 1))) p) q)))

(defun poly-eval (p x)
  (if (null p)
      0
    (+ (car p) (* x (poly-eval (cdr p) x)))))

;more to do
(defun SKF (p)
  (let ((m (floor (/ (1- (length p)) 2)))
	(y) (f))
    (dotimes (i (1+ m) y)
      (push (poly-eval p (* (expt -1 i) (floor (/ (+ i 1) 2)))) y))
    (dolist (i y f)
      (push (factors i) f))))

(defun lagrange (xsupp ysupp)
  (let ((li-lst nil)
	(temp '(0))
	(len (length xsupp)))
    (dotimes (i len li-lst)
      (setq li-lst (cons (lagrange-i xsupp (- len i 1)) li-lst)))
    (dotimes (i len temp)
      (setq temp (poly-add (poly-mul (list (nth i ysupp)) (nth i li-lst)) temp)))))

(defun lagrange-i (supp i)
  (let* ((wi (wi supp i)) 
	(wixi (poly-eval wi (nth i supp))))
    (mapcar #'(lambda (x) (/ x wixi)) wi)))

(defun wi (supp i &optional (curri 0))
  (cond ((null supp) '(1))
	((= i curri) (wi (cdr supp) i (1+ curri)))
	(t (poly-mul (cons (- (car supp)) '(1)) (wi (cdr supp) i (1+ curri))))))

(defun poly-diff (p &optional (i 1))
  (cond ((= i 1) (setq p (cdr p)) (cons (* i (car p)) (diff (cdr p) (1+ i))))
	(t (cond ((null p) nil)
		 (t (cons (* i (car p)) (diff (cdr p) (1+ i))))))))

(defun poly-pp (p)
  (cond ((= (length p) 1) (format t "~A~%" (car p)))
	(t (format t "~Ax^~A+" (car (last p)) (1- (length p)))
	   (poly-pp (butlast p)))))

;;Integer functions
(defun factors-prime (n &optional f (i 3))
  (cond ((<= n 1) f)
	((evenp n) (factors-prime (/ n 2) (union '(2) f)))
	(t (if (= 0 (mod n i))
	       (factors-prime (/ n i) (union (list i) f i))
	     (factors-prime n f (+ i 2))))))

(defun factors (n &optional (i 2) (f '(1 -1)))
  (cond ((> i n) f)
	(t (cond ((= 0 (mod n i)) (push (- i) f) (push i f) (factors n (1+ i) f))
		 (t (factors n (1+ i) f))))))