; ==============================================================================
; PROJET LISP
; 
; Trouver les racines d'un polynom donnee a l'aide de l'algo. de Newton.
;
; Michael Vorburger, 16.7.97 - 18.7.97
; pour le cours de Programmation IV, Prof. Rapin, EPFL
; ==============================================================================


; =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
; NOTE: A "user" should call only function EQUATION_ALGEBRIQUE because only this 
;       guarantees correct "error" checks in parameters. So DON'T be surprised if 
;       you call for example function POLYZERO or SOLVE-2ND-DEGREE directly with
;       "non-standard" coefficients and get unexptexted behaviour.
; -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=


; ------------------------------------------------------------------------------
; (equation_algebrique polycoeffs epsi)
;
; This is the main function to find all solutions of a polynomial function.
; Call it with a list of coefficients (ordered by growing powers) as first,
; and the max. rel. error as second argument and you get a list of all solutions.

(defun equation_algebrique (polycoeffs epsi)
   (let* (( adjusted-poly  (polyadjust polycoeffs) )
          ( first-x0       (find-initial-value adjusted-poly) )
         )
         (if (null adjusted-poly)
         ;then
              nil
         ;else
              (let* ((x       (correct (polyzero adjusted-poly first-x0 epsi)) )
                     (reduced_poly  (polyreduce adjusted-poly x))
                    )
                    (if (not (null reduced_poly))
                        (cons x (equation_algebrique reduced_poly epsi) )
                        x )
              )
         )
   )
)


; ------------------------------------------------------------------------------
; (polyadjust polycoeffs)
;
; Removes zeros at the end of the list polycoeffs. This is crucial, because eg.
; '(1 2 0) would be interpreted as polynominal of second degree (1+2x+0x^2)
; instead as the linear (1+2x)

(defun polyadjust (polycoeffs)
   (if (not (null polycoeffs))
       (if (not (zerop (first (last polycoeffs))))   ; ATTENTION: last returns a 
           polycoeffs                                ;   LIST (with one element)
           (polyadjust (revrest polycoeffs) )        ;        che stupido...
       )
   )
)

; (revrest list) returns all but the last element of a list.
;

(defun revrest (a-list)
   (reverse (rest (reverse a-list))) )
   
                            
; ------------------------------------------------------------------------------
; (find-initial-value polycoeffs)
;
; Determines an initial value. If the polynomial is at least of 2nd degree,
; use the smaller (seen as absolute values) of the two solution. If this
; is not complex, add a small imaginary part to allow convergence.
;
; If it is a linear or constant function, simply return 0 because function
; polyzero will take special care of all those special cases.

(defun find-initial-value (polycoeffs)
   (if (> (length polycoeffs) 2)
       (let* ((solution        
                     (solve-2nd-degree polycoeffs) )
              (smaller-solution
                     (if (< (abs (first solution)) (abs (second solution)))
                         (first solution)
                         (second solution) ) )
             )
             (if (realp smaller-solution)
                 (+ smaller-solution #C(0 1.0d-8) )
                 smaller-solution
             )   
       )
       0.0d0                      ; linear or constant function, return 0
   )
)


; ------------------------------------------------------------------------------
; (correct number)
;
; If number is a complex, and if it's imag. part is < 1e-15, discard it.

(defun correct (z)
   (if (complexp z)
      (if (< (abs (imagpart z)) 1.0e-15)
         (realpart z)
         z
         )
      z
      )
   )


; ------------------------------------------------------------------------------
; (solve-2nd-degree polycoeffs)
;
; This function interpretes the elements of the list polycoeffs as coefficients
; of a quadratic function (ordered by growing powers) and returns a LIST
; containing the two solutions. The might be identical, however. -- Does NOT
; use the standard formula, but one which more appropriate for numeric reasons.
;
; SPECIAL BEHAVIOUR: If the third element of polycoeffs (second degree) is zero,
;                    a "small number" is taken instead. This is only for
;                    find-initial-value and does not play any role when called
;                    from polyzero, because such "pseudo-quadratic" polynomials
;                    would have been adjusted to linear polynomials by 
;                    polyadjust.
;
; CROSSREFERENCE NOTE: This function is called from find-initial-value
;                      and polyzero.

(defun solve-2nd-degree (polycoeffs)
   (let* ((c  (first polycoeffs)  )
          (b  (second polycoeffs) )
          (a- (third polycoeffs)  )
          (a  (if (zerop a-)  1.0d-8  a-) )
          (dr (sqrt (- (* b b) (* 4 (* a c)))) )    ; dr = sqrt(b*b - 4*a*c)
          (e  (* 2 c ) )
         )
         (list (/ e (- (- b) dr) )                  ; (e / (-b) - dr)
               (/ e (+ (- b) dr) )   )              ; (e / (-b) + dr)
   ) 
)

      
; ------------------------------------------------------------------------------
; (polyzero polycoeffs x0 epsi)
;
; Intepretes the list polycoeffs as coefficients (ordered by growing
; powers) of a polynominal function and finds a solution for it; one only.
; This is done by applying Newton's algorithm with x0 as initial value and
; epsi as error which is allowed.
; Returns a number or NIL if no solution could be found.

(defun polyzero (polycoeffs x0 epsi)
   (cond ((= (first polycoeffs) 0)   0)                 ; polynomial with const 0
         ((= (length polycoeffs) 1)  NIL )              ; constant polynomial
         ((= (length polycoeffs) 2)                     ; linear polynomial 
             (- (/ (first polycoeffs) (second polycoeffs))) )      ; x=-b/a
         ((= (length polycoeffs) 3)                     ; quadratic polynomial 
             (first (solve-2nd-degree polycoeffs)) )

          
         (t (let* (( x1        (newton polycoeffs x0) ) ; next iteration x1
                   ( dx        (abs (- x0 x1)) )        ; diff between x0 and x1
                   ( maxabserr (abs (* x0 epsi)) )      ; maximal absolute error
                  )
                  (cond ((nanp x1)                      ; if x1 cannot be calc,
                            (polyzero polycoeffs (+ x0 epsi) epsi) ) ; correction
                        ((< dx maxabserr)                            ; if precise
                             x1 )                                    ; then done
                        (t                 
                             (polyzero polycoeffs x1 epsi) )         ; else next
                  )
             )
         )
   )
)


; (newton polycoeffs xn) calculates the next iteration for Newton's
; approximation. This function just implements Newton's formula for
; ONE step and does not iterate. It is called from function polyzero.
; Returns a number or NaN if the next value could not be calculated.

(defun newton (polycoeffs xn)
   (let* (( yyd (polyeval polycoeffs xn) )
          ( y   (first yyd) )
          ( yd  (second yyd) )    ; PS: (nth 1 yyd) = (elt yyd 1) = (second yyd)
         )
         (- xn (/ y yd) )         ; x[n+1] := x[n] + f(x[n]) / f'(x[n])
   ) 
)


; -------------------------------------------------------------------------------
; (polyreduce polycoeffs x0)
;
; Intepretes the list polycoeffs as coefficients (ordered by growing
; powers) of a polynominal function and returns the polynomial which
; results from dividing polycoeffs by (x-x0)
; 
; This is done by using Horner's algorithm:
;   if    f(x) = u[n]*x^n + u[n-1]*x^(n-1) + ... + u[1]*x + u[0]
;   then  g(x) = f(x) / (x-x0)
;   is         = v[n-1]*x^(n-1) + v[n-2]*x^(n-2) + ... + v[1]*x + v[1]
;   where v[k] = v[k+1]*x0 + u[k]


(defun polyreduce (polycoeffs x0)
   (reverse (horner 0 (reverse polycoeffs) x0))
)

(defun horner (pv polycoeffs x0)       ; pv = previous v (from previous iteration)
   (if (> (length polycoeffs) 1)
       (let (( nv  (+ (* pv x0) (first polycoeffs)) )) ; (first polycoeffs) = u[k]
            (cons nv
                  (horner nv (rest polycoeffs) x0)
            )
       )
   )
)


; -------------------------------------------------------------------------------
; (polyeval polycoeffs x)
;
; Intepretes the list polycoeffs as coefficients (ordered by growing
; powers) of a polynominal function which is evaluted for x. The function
; value and the value of the derrived function at point x are returned as list.
;
; Source: This is copied from Ex. 19.

(defun polyeval (polycoeffs x)
   (cond ((null polycoeffs)          
             nil)                         ; return nil for empty lists
         ((null (rest polycoeffs))   
             (cons (first polycoeffs) '(0)) )  ; return c and 0 for const poly.
         (t
             (let* ((evalrest (polyeval (rest polycoeffs) x) )  ; otherwise...
                    (p_x     (first evalrest)  )
                    (p_prime (second evalrest) )

                   )
                   (cons (+ (* p_x x) (first polycoeffs))
                         (cons (+ p_x (* p_prime x)) nil)
                   )  
             )
      )
   )
)