;;; -*- Mode:Common-Lisp; Package:ma; Base:10 -*-
;;
;; Written by: Tak W. Yan
;; modified by Richard Fateman

;; Contains: procedures that convert between lisp prefix form 
;;           and representation of rat's and also the simplifier itself

;;; (c) Copyright 1990 Richard J. Fateman, Tak Yan
;;; (c) 2005 RJF, modified for generic arithmetic

;;; Basically three functions are provided.
;;;     into-rat:  take an expression parsed by parser, simplify it, and
;;;                represent it as a "rat"
;;;     outof-rat: take a rat and translate it back to Mathematica (tm)-like
;;;                language
;;;     ratsimp:      simplify an expression by converting it back and forth

;;(provide 'simp1)

(declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))

(eval-when (compile load) 
  (require "poly" "poly.fasl")
  (require "rat1" "rat1.fasl")
  (provide 'simprat))
  
(in-package :ma)

(defvar varcount 0)
(defvar vartab (make-hash-table :test #'equal)) ;; map from kernels to integers
(defvar revtab (make-hash-table :test #'eql)) ;; map from integers to kernels
;; reptab should use eq, if we believe that expressions will be
;; "unique" coming from the parser or other transformation programs.
;; (see what happens coming out of "outof-rat" for example.)
;;
(defvar reptab (make-hash-table :test #'eq)) ;; map from expressions to rats
(defvar disreptab (make-hash-table :test #'eq)) ;; map from rats to expressions

;; look-up-var: look up a variable in the hash-table; if exists,
;;    	        return the value of it; otherwise, increment varcount
;;              and add variable to hash-table

(defun look-up-var (x)
  (cond ((gethash x vartab))
	(t (setf (gethash (setq varcount (1+ varcount)) revtab) x)
	   (setf (gethash x vartab) varcount))))

;; special hack to make the variable x the ``MOST'' MAIN variable.
;; useful for integration variable, for example.
;; Do this before any other use of the symbol z.

(defun make-main-var (z &optional (place #.most-positive-fixnum))
  (cond ((gethash z vartab) 
	 (if (not (eq z (gethash place revtab)))
		  (format t "~% ~s already is a variable with a different order;
   existing expressions may be broken" z))))
  (cond ((gethash place revtab)
	 (format t "~% ~s was most-main, previously."
		 (gethash place revtab))
	 (make-main-var (gethash place revtab) (1- place))))
  (setf (gethash z vartab) place)
  (setf (gethash place revtab) z))
	

;;; into the rat representation

;; into-rat: represent as a rat

(defun into-rat (x)
  (typecase x 
	    (integer (coef2rat x))
	    (ratio  (make-rat :numerator `((,(cl::numerator x) . 1))
			      :denominator `((,(cl::denominator x) . 1))))
	    (rat x) ;already a rat form
;;	    (array  whatever)
	    (t
	     (cond((gethash x reptab))       ; see if remembered
		  ((setf (gethash x reptab) (into-rat-1 x)))))
))

;; into-rat-1: actually do the work

(defun into-rat-1 (x)
   (cond ((symbolp x) (var2rat x))
	  ((and (consp x)(atom (car x)))
	   (case (car x)
	     (*  (reduce-rat #'rat* (into-rat (cadr x)) (cddr x)))
	     (/  (reduce-rat #'rat/ (into-rat (cadr x)) (cddr x)))
	     (+  (reduce-rat #'rat+ (into-rat (cadr x)) (cddr x)))
	     (-  (reduce-rat #'rat- (into-rat (cadr x)) (cddr x)))
	     (expt (into-rat-2 x))
	      (t(var2rat (umapcar #'simp x)))))

	  ((consp x)
	   (var2rat (umapcar #'simp x)))
	  ;; we could check for interval or bigfloat or ?? but then what?
	  ;; we give up on what this is. just wrap it up as a variable
	  (t (var2rat x))))

;; into-rat-2: handle the case of powering

(defun into-rat-2 (x)
  (let* ((exp (into-rat (caddr x)))
	 (exp-n (rat-numerator exp))
	 (exp-d (rat-denominator exp)))
    (cond ((and (fpe-coef-p exp-n) (fpe-coef-p exp-d))
	   (cond ((fpe-coefone-p exp-d) (rat^ (into-rat (cadr x))
					      (fpe-expand exp-n)))
		 (t (rat^ (var2rat (ulist 'expt
					  (simp (cadr x))
					  (simp (ulist '* 1
						 (ulist 
						  'expt
						  (fpe-expand exp-d) -1)))))
			  (fpe-expand exp-n)))))
	  (t (var2rat (ulist 'expt (simp (cadr x)) (outof-rat exp)))))))

;; var2rat: convert a variable into a rat

(defun var2rat (x)
  (make-rat :numerator (make-fpe (vector (look-up-var x) 0 1) 1)
	    :denominator (make-fpe (coefone) 1)))

;; coef2rat: convert a coef into a rat

(defun coef2rat (x)
  (make-rat :numerator (make-fpe x 1)
	    :denominator (make-fpe (coefone) 1)))

;; reduce-rat: apply fn to r and the result of applying fn recursively to
;;             the terms in l

(defun reduce-rat (fn r l)
  (cond ((null l) r)
	(t (funcall fn r (reduce-rat fn (into-rat (car l)) (cdr l))))))

;;; out of the rat representation

;; outof-rat: translate back to Mathematica (tm)-like lists

(defun outof-rat (r)
  (cond ((gethash r disreptab))
	((setf (gethash r disreptab)
	   (let ((n (fintol (rat-numerator r)))
		 (d (fintol (rat-denominator r))))    
	     ;; here we remember the output simplified form
	     ;; so that if it is fed back in, you get exactly the same rat rep.
	     ;; In fact, by re-collecting terms in a different order, it might
	     ;; sometimes change. E.g. (x^2+2x) results from (x+1)^2-1,
	     ;; but if you type it in, you'll get x*(x+2)!
	     (let ((ans (outof-rat-1 n d)))
	       (setf (gethash ans reptab) r) ;; controversy here
	       ans))))))


;; outof-rat-1: take 2 fpe's, n and d, and translate 
;;              into list form; n is numerator and d is denominator;


(defun outof-rat-1 (n d)
  (cond ((equal 1 d) n) ;; e.g. 3
	((equal 1 n) 
	 (cond ((and (consp d)(eq (car d) 'expt))
		(cond ((numberp (caddr d));; e.g. y^(-4)
		       (ulist 'expt (cadr d)(- (caddr d))))
		      (t ;; e.g. y^(-x)
		       (ulist 'expt (cadr d)(ulist '* -1 (caddr d))))))
	       ((numberp d)(expt d -1)) ;; 1/3
	       (t (ulist 'expt d -1)))) ;; x^-1
	((and (numberp n)(numberp d)) (/ n d)) ;; e.g. 1/2
	((numberp d)(ulist '* (outof-rat-1 1 d) n)) ;; e.g. 1/30(x^2+1)
	(t (ulist '* n 
		  (outof-rat-1 1 d)))))

;; fintol: convert a fpe into list form; the fpe must be normal,
 
(defun fintol (f)
  (let ((c (caar f))); constant term
    (do ((j (cdr f)(cdr j))
	 (lis nil))
	((null j)
	 (cond ((null lis) c)
	       ((and (coefonep c)(null(cdr lis)))
		(car lis))
	       (t (cond ((coefonep c) (ucons '* (uniq(nreverse lis))))
			(t (ucons '*(ucons c (uniq (nreverse lis)))))))))
	;; the loop body:
	(setq lis (cons (fintol2 (caar j)(cdar j)) lis)))))



;; fintol2: break a pair in a fpe into list form

(defun fintol2 (p e)
  (cond ((eq e 1) (intol p))
	(t (ulist 'expt (intol p) e))))

;; intol: convert a polynomial in vector form into list form

(defun intol (p)
       (cond ((vectorp p)
	      ;; look up what the kernel is from revtab
	      (intol+ (intol2 p (gethash (svref p 0) revtab
					 ;;in case (svref p 0) is not 
					 ;; in hashtable, use it literally
					 (svref p 0)))))
	     (t p)))

;; intol2: help intol

(defun outof-rat-check(x)
  (cond ((rat-p x) (outof-rat x))
	(t (intol x))))

(defun intol2 (p var)
  (let (res term)
    (do ((i (1- (length p)) (1- i)))
	((= i 0) res)
	(setq term (intol* (outof-rat-check (svref p i)) (intol^ (1- i) var)) )
	(cond ((eq term 0))
	      (t (setq res (ucons term res)))))))

;; intol^: handle the case for powering

(defun intol^ (n var)
       (cond ((zerop n) 1)
	     ((equal n 1) var)
	     (t (ulist 'expt var n))))

;; intol+: handle +

(defun intol+ (p)
  (cond ((null (cdr p)) (car p))
	((IsPlus (car p)) (uappend (car p) (cdr p)))
	(t (ucons '+ p))))

;; intol*: handle *

(defun intol* (a b)
  (cond ((eql a 0) 0)
	((eq a 1) b)
	((eq b 1) a)
	(t (ucons '*
		  (uappend (intol*chk a) (intol*chk b))))))

;; into*chk: help into*

(defun intol*chk (a)
  (if (IsTimes a) (cdr a) (ulist a)))

;; IsPlus: check if a is led by a +

(defun IsPlus (a) (and (consp a) (eq (car a) '+)))

;; IsTimes: check if a is led by *

(defun IsTimes (a) (and (consp a) (eq (car a) '*)))

;;; simplify by converting back and forth

;; simp: simplify the expression x
;; assumes that all kernels that are non-identical are algebraically
;; independent. Clearly false for e.g. Sin[x], Cos[x], E^x, E^(2*x)

(defmethod ratsimp ((x t)) (outof-rat (into-rat x)))

(defmethod ratexpand((u t))
  (outof-rat
  (let* ( ;;(*expand* t) ;; global flag to rat program
	 (x (into-rat u)))
    (make-rat :numerator (make-fpe (fpe-expand (rat-numerator x))1)
	       :denominator
	       (make-fpe (fpe-expand (rat-denominator x)) 1)))))


;; This gives a list of kernels assumed to be independent.
;; If they are NOT, then the simplification may be incomplete.
;; In general, the search for the "simplest" set of kernels is
;; difficult, and leads to (for example) the Risch structure
;; theorem, Grobner basis decompositions, solution of equations
;; in closed form algebraically or otherwise.  Don't believe me?
;; what if the kernels are Rootof[x^2-3],Rootof[y^4-9], Log[x/2],
;; Log[x], Exp[x], Integral[Exp[x],x] ....

#+ignore
(defun Kernelsin(x)(cons 'List
			  (mapcar #'(lambda(x)(gethash x revtab))
			    (collectvars (into-rat x)))))


#+ignore  ;; one way of putting floats into the system
(defun into-rat (x)
  (typecase x 
	    (integer (coef2rat x))
	    (ratio  (make-rat :numerator `((,(numerator x) . 1))
			      :denominator `((,(denominator x) . 1))))
	    (rat x) ;already a rat form
;;	    (array  whatever)
	    (float (coef2rat x))
	    (t
	     (cond((gethash x reptab))       ; see if remembered
		  ((setf (gethash x reptab) (into-rat-1 x)))))
))