;;;f(x) = a*x^2+b*x+c = 0 ;
;;;Input: the coefficients a, b, c are real numbers ;
;;;Output: when a /= 0,one or Two solutions,all of them;
;;; like this: (x y),means: x + y * i,if it's a ;
;;; real number,then y = 0. ;
;;; Otherwise , return a real number or nil ;
;;;Ref: http://en.wikipedia.org/wiki/Quadratic_equation;
;;;----------------------------------------------------;
(defun Math:Quadratic_Equation (a b c / d e g)
(if (zerop a)
(if (not (zerop b))
(list (list (/ (- c) (float b)) 0))
)
(progn
(setq a (float a))
(if (not (equal a 1 1e-14))
(setq b (/ b a)
c (/ c a)
)
)
(setq d (- (* b b) (* 4 c)))
(setq e (* b -0.5))
(cond
( (equal d 0 1e-14)
(list (list e 0) (list e 0))
)
( (> d 0)
(setq g (* (sqrt d) -0.5))
(list (list (- e g) 0) (list (+ e g) 0))
)
( (< d 0)
(setq g (* (sqrt (- d)) -0.5))
(list (list e (- g)) (list e g))
)
)
)
)
)
;;;----------------------------------------------------; ;
;;;f(x) = a*x^3+b*x^2+c*x+d = 0 ;
;;;Input: the coefficients a, b, c, d are real numbers ;
;;;Output: when a /= 0,Three solutions,all of them like;
;;; this: (x y),means: x + y * i,if it's a real ;
;;; number,then y = 0; ;
;;; otherwise goto "Math:Quadratic_Equation" ;
;;;Ref: http://en.wikipedia.org/wiki/Cubic_function ;
;;;----------------------------------------------------;
(defun Math:Cubic_Equation (a b c d / e f g h u w x y)
(cond
( (zerop a)
(Math:Quadratic_Equation b c d)
)
( (zerop d)
(cons '(0 0) (Math:Quadratic_Equation a b c))
)
(t
(setq b (/ b a 3.0))
(setq c (/ c a 6.0))
(setq d (/ d a 2.0))
(setq e (- (* b (- (+ c c c) (* b b))) d)) ;Alpha
(setq f (- (* b b) c c)) ;Beta
(setq g (- (* e e) (* f f f))) ;Delta,The nature of the roots
(cond
( (equal g 0 1e-14)
(setq u (MATH:CubeRoot e))
(list (list (- (+ u u) b) 0.0)
(list (- (+ b u)) 0.0)
(list (- (+ b u)) 0.0)
)
)
( (> g 0)
(setq h (sqrt g))
(setq u (MATH:CubeRoot (+ e h)))
(setq w (MATH:CubeRoot (- e h)))
(setq x (- (+ b (* (+ u w) 0.5))))
(setq y (* (sqrt 3) (- u w) 0.5))
(list (list (+ u w (- b)) 0)
(list x y)
(list x (- y))
)
)
( (< g 0)
(setq x (/ e f (sqrt f)))
(setq y (sqrt (abs (- 1 (* x x)))))
(setq u (/ (atan y x) 3))
(setq w (/ (+ pi pi) 3))
(setq h (* 2 (sqrt f)))
(list (list (- (* (cos u) h) b) 0)
(list (- (* (cos (+ u w)) h) b) 0)
(list (- (* (cos (- u w)) h) b) 0)
)
)
)
)
)
)
;;;----------------------------------------------------; ;
;;;f(x) = a*x^4+b*x^3+c*x^2+d*x+e= 0 ;
;;;Input: the coefficients a,b,c,d,e are real numbers. ;
;;;Output: when a /= 0,Three solutions,all of them like;
;;; this: (x y),means: x + y * i,if it's a real ;
;;; number,then y = 0; ;
;;; otherwise goto "Math:Quadratic_Equation" ;
;;;Ref: http://en.wikipedia.org/wiki/Quartic_function ;
;;;----------------------------------------------------;
(defun Math:Quartic_Equation (a b c d e / B2 B3 B4 EPS F G H P Q R V W X Y Z S)
(setq eps 1e-8)
(cond
( (equal a 0 eps)
(Math:Cubic_Equation b c d e)
)
( (equal e 0 eps)
(cons '(0 0) (Math:Cubic_Equation a b c d))
)
( (and (equal b 0 eps) (equal d 0 eps))
(foreach x (Math:Quadratic_Equation a c e)
(foreach y (Math:CSqrt x)
(setq s (cons y s))
)
)
)
( (setq a (float a))
(setq b (/ b a))
(setq c (/ c a))
(setq d (/ d a))
(setq e (/ e a))
(setq b2 (* b b))
(setq b3 (* b2 b))
(setq b4 (* b3 b))
(setq f (+ (* b2 -0.375) c)) ;alpha
(setq g (+ (* b3 0.125) (* b c -0.5) d)) ;beta
(setq h (+ (* -0.01171875 b4) (* 0.0625 c b2) (* -0.25 b d) e))
(if (equal g 0 eps)
(progn
(setq x (* b -0.25))
(mapcar
(function (lambda (e) (list (+ (car e) x) (cadr e))))
(Math:Quartic_Equation 1 0 f 0 h) ;Recursion
)
)
(progn
(setq p (- (* f f 2) h))
(setq q (- (* f f f 0.5) (* f h 0.5) (* g g 0.125)))
(setq r (Math:Cubic_Equation 1 (* 2.5 f) p q))
(while (not (equal (cadar r) 0 eps))
(setq r (cdr r))
)
(setq r (caar r))
(setq w (sqrt (abs (+ f r r))))
(foreach i (list + -)
(foreach j (list + -)
(setq x (i (* b -0.25) (* w 0.5)))
(setq y (- (+ f f f r r (i (/ g 0.5 w)))))
(if (< y 0)
(setq z (list x (j (* (sqrt (- y)) 0.5))))
(setq z (list (+ x (j (* (sqrt y) 0.5))) 0))
)
(setq S (cons z S))
)
)
)
)
)
)
)
HI,Ribarm!
The above code is written by highflybir!
hope this helps you!