(in-package :maxima)

(defun require-positive-integer (n pos fun)
  (if (or (not (integerp n)) (< n 1))
      (merror "The ~:M argument of the function ~:M must be a positive integer" pos fun)))

(defvar lettervars (make-array 26))
(dotimes(i 26) (setf (aref lettervars i)
		     (concat #\$ (code-char (+ i #.(char-code #\a))))))

(defun random-identifier (n) (setf n ($random n))
  (if (< n 26) (aref lettervars n)
    (concat #\$ (code-char (+ n #.(char-code #\a))))))

(defun random-rational ()
  (div (* (if (= 1 ($random 2)) -1 1) (+ 1 ($random 100))) (+ 1 ($random 100))))

(defun random-maxima-number ()
  (let ((k ($random 2)))
    (cond ((= k 0) (random-rational))
	  ((= k 1) ($random 1.0))
	  (t ($bfloat ($random 1.0))))))
				  
;; Return a monomial of the form q x1^n1 x2^n2 .. xp^np, where q is rational 
;; and 0 <= nk <= max-power and 1 <= p <= max-vars. Thus the monomial is
;; nonconstant.

(defun $random_monomial (&optional (max-power 10) (max-vars 1))
  (require-positive-integer max-power "$second" "$random_monomial")
  (require-positive-integer max-vars "$third" "$random_monomial")

  (let ((acc nil) (x) (n  (+ 1 ($random max-vars))))
    (dotimes (i n)
      (while (member (setq x (random-identifier max-vars)) acc))
      (push x acc))
    (mul (random-rational) (apply 'mul (mapcar #'(lambda (s) (power s (+ 1 ($random max-power)))) acc)))))
	
;; Return an expanded random polynomial with rational coefficients.  Specifically,
;; random_poly(maxterms, maxpower, maxvars) returns a polynomial with terms of the
;; form  q  x1^n1 x2^n2 .... xp^np, where q is rational, each nk satisfies
;; 0 <= nk <= maxpower, and p = maxvars. Since any nk can vanish, each term can
;; have fewer than maxvars variables.

(defun $random_expanded_poly (&optional (max-terms 10) (max-power 2) (max-vars 1))
  (require-positive-integer max-terms "$first" "$random_poly")
  (require-positive-integer max-power "$second" "$random_poly")
  (require-positive-integer max-vars "$third" "$random_poly")

  (let ((n (+ 1 ($random max-terms))) (acc (list (random-rational))))
    (dotimes (i n `((mplus) ,@acc))
      (push ($random_monomial max-power max-vars) acc))))

(defvar *nary-moperators* `(mplus mplus mplus mplus mplus mplus mplus mplus mtimes))

(defun random-nary ()
  (nth ($random (length  *nary-moperators*))  *nary-moperators*))

(defun $random_poly (&optional (num-term 10) (max-power 5) (max-vars 1) (non-cnst nil) (non-power nil))
  ;(require-positive-integer max-vars "$second" "$random_expr")
  ;(setq num-term ($random num-term))
  (let ((k))

    (cond ((or (< num-term 2) (< max-power 1))
	   (if (or non-cnst (< ($random 5) 3))  (random-identifier max-vars) (random-rational)))
	  
	  ((and (not non-power) (< ($random 10) 2))
	   `((mexpt) ,($random_poly num-term max-power max-vars t t) ,($random max-power)))
	  
	  (t
	   (setq k (+ 1 ($random (- num-term 1))))
	       `((,(random-nary)) 
		 ,($random_poly k max-power max-vars nil non-power) 
		 ,($random_poly (- num-term k) max-power max-vars t non-power))))))

(defun $random_rational_expr (&optional (nlen 10) (max-power 5) (max-vars 1) (non-cnst nil) (non-power nil))
  ;(require-positive-integer max-vars "$second" "$random_expr")
  ;(setq nlen ($random nlen))
  `((mtimes) 
    ,($random_poly nlen max-power max-vars non-cnst non-power)
    ((mexpt) ,($random_poly nlen max-power max-vars non-cnst non-power) -1)))