Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagesmc
 Path: blob/master/src/ext/maxima/willis/topoly.lisp
8815 views
1
;; Written by Barton Willis, [email protected]

2
;; and emailed to maxima users group email list

3
;; on Sep 2006

4
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

5
(defun suppress-multiple-zeros (q)

6
  (let ((acc 1) ($factorflag nil))

7
    (setq q ($factor q))

8
    (setq q (if (mtimesp q) (margs q) (list q)))

9
    (dolist (qi q acc)

10
      (setq acc (mul acc (cond ((mnump qi) (if (eq t (meqp qi 0)) 0 1))

11
			       ((mexptp qi) (nth 1 qi))

12
			       (t qi)))))))

13


14
;; When the second argument is true, convert constant expressions to polynomial form.

15
;; The default for the second argument is false. Examples:

16
;;

17
;;(%i1) topoly(x-sqrt(2),false);

18
;;(%o1) x-sqrt(2)=0

19
;;(%i2) topoly(x-sqrt(2),true);

20
;;(%o2) x^2-2=0

21
;;(%i7) topoly(sqrt(2)+sqrt(3)-x,true);

22
;;(%o7) x^4-10*x^2+1=0

23
;;(%i8) solve(%,x);

24
;;(%o8) [x=-sqrt(2*sqrt(6)+5),x=sqrt(2*sqrt(6)+5),x=-sqrt(5-2*sqrt(6)),x=sqrt(5-2*sqrt(6))]

25
;;(%i9) map(lambda([e],topoly(e,true)),%);

26
;;(%o9) [x^4-10*x^2+1=0,x^4-10*x^2+1=0,x^4-10*x^2+1=0,x^4-10*x^2+1=0]

27


28
(defun $topoly (p &optional (cnst nil))

29


30
  (if (not (or (eq t cnst) (eq nil cnst)))

31
      (merror "The second argument to 'topoly' must be either 'true' or 'false'"))

32


33
  (let ((subs) (q) (nv `((mlist)))) ;; new variables

34
    (setq p (meqhk p))

35
    (setq q ($ratdenom p))

36
    (if (not ($constantp q)) (mtell "Assuming that ~:M " `((mnotequal) ,q 0)))

37
    (setq p ($ratdisrep ($ratnumer p)))

38


39
    (setq p (to-polynomial p nil cnst))

40
    (setq subs (second p))

41
    (setq p (first p))

42
    (dolist (sk subs)

43
      (setq nv ($append nv ($listofvars ($lhs sk)))))

44
    (setq p (if (null subs) p ($first (mfuncall '$eliminate `((mlist) ,p ,@subs) nv))))

45
    `((mequal) ,(suppress-multiple-zeros p) 0)))

46


47
(defun to-polynomial (p subs cnst)

48
  (cond (($mapatom p) (list p subs))

49
	((and (not cnst) ($constantp p)) (list p subs))

50


51
	((mexptp p)

52
	 (let ((n (nth 2 p)) (b (nth 1 p)) (nv)(l))

53
	   (cond ((integerp n)

54
		  (setq b (to-polynomial b nil cnst))

55
		  (setq subs (append (second b) subs))

56
		  (setq b (first b))

57
		  (if (> n 0) (list (power b n) subs) (merror "Unable to convert to a polynomial equation")))

58


59
		 (($ratnump n)

60


61
		  (setq b (to-polynomial b nil cnst))

62
		  (setq subs (append (second b) subs))

63
		  (setq b (first b))

64
		  (setq nv (gensym))

65
		  (setq subs (cons `((mequal) ,(power nv ($denom n)) ,(power b ($num n))) subs))

66
		  (list nv subs))

67
		 (t (merror "Nonalgebraic argument given to 'topoly'")))))

68


69
	((op-equalp p 'mabs)

70
	 (setq b (to-polynomial (first (margs p)) nil cnst))

71
	 (setq subs (append (second b) subs))

72
	 (setq b (first b))

73
	 (setq nv (gensym))

74
	 (list nv (cons `((mequal) ,(power nv 2) ,(power b 2)) subs)))

75


76
	((mtimesp p)

77
	 (let ((z 1) (acc nil))

78
	   (setq p (mapcar #'(lambda (s) (to-polynomial s nil cnst)) (margs p)))

79
	   (dolist (pk p)

80
	     (setq z (mul z (first pk)))

81
	     (setq acc (append acc (second pk))))

82
	   (list z acc)))

83


84
	 ((mplusp p)

85
	  (let ((z 0) (acc nil))

86
	    (setq p (mapcar #'(lambda (s) (to-polynomial s nil cnst)) (margs p)))

87
	    (dolist (pk p)

88
	      (setq z (add z (first pk)))

89
	      (setq acc (append acc (second pk))))

90
	    (list z acc)))

91


92
	 (t (merror "Nonalgebraic argument given to 'topoly'"))))

93


94


95


96


97


98


99