(in-package "MAXIMA")
(use-package "SLOOP")

(defvar $todd_coxeter_state nil)
(proclaim '(type (vector (t))  $todd_coxeter_state))
;;  To turn on debug printing set to T
(defvar *debug* nil)
;;  When *debug* is not nil, this holds the multiplications for
;; the current row.
(defvar *this-row* nil)

(eval-when (compile eval)
(defmacro nvars () '(aref $todd_coxeter_state 0))
(defmacro ncosets() '(aref $todd_coxeter_state 1))
(defmacro multiply-table () '(aref $todd_coxeter_state 2))
(defmacro relations () '(aref $todd_coxeter_state 3))
(defmacro subgroup-generators () '(aref $todd_coxeter_state 4))
(defmacro row1-relations () '(aref $todd_coxeter_state 5))
(defmacro with-multiply-table (&body body)
  `(let ((nvars (nvars))(multiply-table (multiply-table)))
    (declare (fixnum nvars) (type (vector (t)) multiply-table))
    ,@ body)) 

(defmacro  undef (s) `(eql 0 ,s))

;; Multiply coset K times variable R
(defmacro mult (k r) `(the coset (aref (table ,r) ,k)))

; Force  k . r = s and  k = s . r^-1
(defmacro define-mult (k r s)
  `(progn (setf (mult ,k ,r) ,s)
	  (setf (mult ,s (- ,r)) ,k)))

;; cosets M < N are to be made equal
(defmacro push-todo (m n)
  `(progn (vector-push-extend ,m *todo*)
	  (vector-push-extend ,n *todo*)))

(defmacro f+ (a b) `(the fixnum (+ (the fixnum ,a) (the fixnum ,b))))
(defmacro f- (a &optional b) `(the fixnum (- (the fixnum ,a)
					     ,@ (if b `((the fixnum ,b))))))
  
;; The multiplication table for variable I 
(defmacro table (i)
  `(the (vector (coset)) (aref  multiply-table (f+ ,i nvars))))

;; Some optional declarations of functions.
(proclaim '(ftype (function (fixnum) t) doing-row)) 
(proclaim '(ftype (function (t t t) t) set-up todd-coxeter )) 
(proclaim
    '(ftype (function nil t) fill-in-inverses dprint-state next-coset
            dcheck-tables replace-coset-in-multiply-table)) 
(proclaim
    '(ftype (function (t) t) $todd_coxeter coerce-rel has-repeat)) 
(proclaim '(ftype (function (t t) t) my-print))
(proclaim '(ftype (function (t *) t) $todd_coxeter))


) ;; end of the macros and proclamations.

(eval-when (compile eval load)
(deftype coset nil 'fixnum)
(proclaim '(type (vector (coset)) *todo*))
)

;; The data type we use to enumerate cosets.


(defvar *todo* (make-array 10 :element-type 'coset :fill-pointer 0 :adjustable t))

;; NVARS is the number of of variables.   It should be the maximum
;; of the absolute values of the entries in the relations RELS.
;; The format of the relations is variables X,Y,.. correspond to
;; numbers 1,2,.. and X^^-1, Y^^-1 .. are -1,-2,...   RELS is
;; a list of lists in these variables.
;; Thus rels = '((1 -2 -1) (2 2 3) ..)  (ie [x1.x2^^-1 . x1^^-1, x2.x2.x3,.. ))
;; SUBGP is also a list of lists.
;; Returns order of G/H, where G is Free/(rels), and H is
;; This is the main entry point at lisp level.
;; Example: (TODD-COXETER 2 '((1 1) (1 2 1 2 1 2) (2 2)))
;; returns 6.   In (multiply-table) we find the current
;; state of the action of the variables on the cosets G/H.
;; For computing the symmetric group using the relations
;; p(i,j) :=concat(x,i).concat(x,j);
;; symet(n):=create_list(if (j - i) = 1 then (p(i,j))^^3 else
;;     if (not i = j) then (p(i,j))^^2 else p(i,i) , j,1,n-1,i,1,j);
;; todd_coxeter(symet(n)) == n!
;; the running time of the first version of this code is observed to be quadratic
;; in the number of cosets.  On a rios it is approx 5*10^-5 * (ncosets)^2.

(defun todd-coxeter (nvars rels subgp &aux (i 1) (c 0 ))
  (declare (fixnum i c))
  (set-up nvars rels subgp)
  (sloop while (>= (ncosets) i)
     do  (incf c) ;; count how many row tries.. 
     (cond
      ;; row still being done
      ((doing-row i) (replace-coset-in-multiply-table))
      ;; row finished but there is work to do
      ((> (the fixnum (fill-pointer *todo*)) 0)
       (incf i) (replace-coset-in-multiply-table))
      ;; row finished -- no work
      (t (incf i))))
  (format t "~%Rows tried ~a" c) 
  (ncosets))

;; Store the data in $todd_coxeter_state, and build multiply-table.
(defun set-up (nvars rels subgp)
  (setf $todd_coxeter_state  (make-array 6))
  (setf (nvars) nvars)
  (setf (ncosets) 1)
  (setf (relations) rels)
  (setf (subgroup-generators) subgp)
  (setf (row1-relations) (append (subgroup-generators)  (relations)))
  (setf (fill-pointer *todo*) 0)
  (setf (multiply-table) (make-array (* 2 (1+ (nvars)))))
  (with-multiply-table
   (sloop for rel in (row1-relations) do
      (sloop for v in rel
	 do (or (<= 1 (abs v)  nvars)
		(error
		 "Vars must be integers with absolute value between 1 and ~a"
		 nvars))))
   (sloop for i from (f- nvars)  to nvars
      when (not (zerop i))
      do (setf (table i) (make-array 10 :adjustable t :element-type 'coset)))
   ))

;; Starts multiplying coset i times the relations.  Basic fact is i . rel = i.
;; This gives a condition on the multiplication table.  Once we have made it all
;; the way through the relations for a given coset i, and NOT had any
;; incosistency in our current multiplication table, then we go on the the next
;; coset.  The coset 1 denotes H.  so for generators h of H we we have 1 . h = 1.
;; So when we do row 1, we add to the relations the generators of H.

;; When we do find an inconsistency eg: 7 . y = 1 and 4 . y = 1 or 7 = 1 . y^^-1
;; and 4 . y = 1, then we would know that 4 and 7 represent the same coset, and
;; so we put 4 and 7 in the *todo* vector and return t so that
;; replace-coset-in-multiply-table will identify them.  While we are running
;; inside doing-row, the multiply-table is accurate, up to our current state of
;; knowledge.  Note that once we find such a nonpermutation action of y, we could
;; not maintain the consistency of (table i) and (table -i).  We exit doing-row
;; with value t, to indicate replacements should be done, and that we must
;; return to complete row i.  (Actually we return t even in the case we were
;; finished the row and found the duplicate in the last step).

(defun doing-row ( i &aux (j 0) (k 0) (r 0)(s 0) *this-row* relations)
  (declare (fixnum j  k i r s))
  (setf relations (if (eql i 1) (row1-relations) (relations)))
  (with-multiply-table
   (sloop for rel in relations
      for v on relations
      do
      (setq k i)
      (sloop 
	 do
	 (setq r (car rel))
	 (setq s (mult  k r))
	 (cond
	  ((undef s)
	   (cond ((cdr rel)
		  (setq s  (next-coset))
		  (define-mult k r s))
		 (t (setq s (mult i  (- r)))
		    (cond ((undef s) (define-mult k r i))
			  ((< k s) (push-todo k s)(return-from doing-row (cdr v)))
			  ((> k s) (push-todo s k)(return-from doing-row (cdr v))))
		    (loop-finish)))))
	 (cond ((setq rel (cdr rel))
		(when *debug* 
		  (push s *this-row*)
		  (my-print (reverse *this-row*) i))
		(setq k s)
		(incf j))
	       ((< i s) (push-todo i s) (return-from doing-row (cdr v)))
	       ((> i s) (push-todo s i) (return-from doing-row (cdr v)))
	       (t			;rel is exhausted and it matched 
		(loop-finish))))))
  (when *debug*  (dcheck-tables) 
	(my-print (reverse *this-row*) i))
  nil)

;; FILL-IN-INVERSES not only completes the (table i) for i < 0
;; but at the same time checks that (table i) for i > 0
;; does not have repeats.   eg if  5 . y = 3 and 7 . y = 3,
;; then this would show up when we go to build the inverse.
;; if it does we add 5 and 7 to the *todo* vector.

(defun fill-in-inverses (&aux (s 0) (sp 0))
  (declare (fixnum s sp))
  ;(format t "~%In fillThere are now ~a cosets" (ncosets))
 (with-multiply-table
  (sloop for i from 1 to nvars
     do (let ((ta1 (table i))
	      (ta2 (table (- i))))
	  (declare (type (vector (coset)) ta1 ta2))
	  (sloop for j from 1 to (ncosets) do (setf (aref ta2 j) 0))
	  (sloop for j from 1 to (ncosets)  do
	     (setf s (aref ta1 j))
	     when (not (eql 0 s))
	     do
	     (setf sp (aref ta2 s))
	     (cond ((eql 0 sp) (setf (aref ta2 s) j))
		   (t	       ;; there's a duplicate!
		    (push-todo sp j)
		    (return-from fill-in-inverses t))))))))
		

;; set n (vector-pop *todo*) , m (vector-pop *todo*)
;; and replace n by m in multiply-table and in *todo*.
;; The replacement is done carefully so as not to lose ANY
;; information from multiply-table, without recording it in
;; *todo*.   It finishes by doing FILL-IN-INVERSES which may
;; in turn cause entries to be added to *todo*.

(defun replace-coset-in-multiply-table (&aux (m 0) (n 0) (s 0) (s2 0) )
  (declare (Fixnum m n s s2 ))
  (with-multiply-table	   
   (tagbody
    AGAIN
    (setf n (vector-pop *todo*))
    (setf m (vector-pop *todo*))
    (unless (eql m n)
	    (dprint-state)
	    (if *debug* (format t "     ~a --> ~a " n m))
	   
	    (sloop for i from 1 to nvars
	       do
	       (let ((ta (table i)))
		 (declare (type  (vector (coset)) ta))
		 (setq s2 (mult n i))
		 (unless (undef s2)
			 (setq s (mult m i))
			 (cond
			  ((undef s) (setf (mult m i) s2))
			  ((< s s2) (push-todo s s2))
			  ((> s s2)(push-todo s2 s))))
		 (sloop for  j downfrom (f- n 1) to 1
		    do (setq s (aref ta j))
		    (cond ((>  s n) (setf (aref ta j) (f- s 1)))
			  ((eql s n) (setf (aref ta j) m) )))
		 (sloop for  j from n below (ncosets)
		    do (setq s (aref ta (f+ j 1)))
		    (cond ((>  s n) (setf (aref ta j) (f- s 1)))
			  ((eql s n) (setf (aref ta j) m) )
			  (t  (setf (aref ta j) s))))))
		
	    (sloop for i downfrom (f- (fill-pointer *todo*) 1) to 0
	       do (setf s (aref *todo* i))
	       (cond ((> s n) (setf (aref *todo* i) (f- s 1)))
		     ((eql s n)(setf (aref *todo* i)  m))))
	    (decf (ncosets))
	    (dprint-state)
	    (progn nil)
	    )

    (if (>  (the fixnum (fill-pointer *todo*)) 0) (go AGAIN))
    ;;(format t "~%There are now ~a cosets" (ncosets))
    ;; check for new duplicates introduced!!
    (if (fill-in-inverses) (go AGAIN))
    )))

;; Get the next coset number, making sure the multiply-table will
;; have room for it, and is appropriately cleaned up.
(defun next-coset ()
  (let* ((n (f+ 1 (ncosets)))
	 (m 0))
    (declare (fixnum n))
    (with-multiply-table    
     (let ((ta (table 1)))
       (unless (> (the fixnum (array-total-size ta)) (f+ n 1))
	       (setf m (f+ n  (ash n -1)))
	       (sloop for i from (f- 0 nvars) to nvars 
		  when (not (eql i 0))
		  do (setf ta (table i))
		  (setf (table i) (adjust-array   ta m ))))
       (sloop for i from 1  to nvars
	  do (setf (aref (table i) n) 0)
	  (setf (aref (table (f- i)) n) 0))))
    (setf (ncosets) n)))



;;  $todd_coxeter parses maxima args
;; todd_coxeter(rels, subgrp) computes the
;; order of G/H where G = Free(listofvars(rels))/subgp_generated(rels));
;; and H is generated by subgp.   Subgp defaults to [].
;; todd_coxeter([x^^3,y.x.y^^-1 . x^^-1],[])  gives 6  the order of the symmetric group
;; on 3 elements.

(defun $todd_coxeter (rels &optional (subgp '((mlist))))
  (let ((vars ($sort ($listofvars rels)))
	(neg 1))
    (declare (special neg vars))
    (todd-coxeter ($length vars) (mapcar 'coerce-rel (cdr rels))
			  (mapcar 'coerce-rel (cdr subgp))
			  )))

(defun coerce-rel (rel )
      (declare (special vars neg))
  (cond ((atom rel)(list (* neg (position rel vars))))
	(t (case (caar rel)
	     (mnctimes (apply 'append (mapcar 'coerce-rel (cdr rel))))
	     (mncexpt
	      (let* ((n  (meval* (third rel)))
		     (neg (signum n))
		     (v (coerce-rel (second rel))))
		(declare (special neg))
		(sloop for i below (abs (third rel))
		   append v)))
	     (otherwise (error "bad rel"))))))

;; The following functions are for debugging purposes, and
;; for displaying the rows as they are computed.   They are
;; not required for correct running.
;#+debug
(progn
(defvar *names* '(nil x y z))
(defun my-print (ro i  &aux relations)
  (when *debug*
  (fresh-line)
;  (print ro)
  (format t "Row ~a " i)
  (setq relations (if (eql i 1) (row1-relations) (relations)))
  (sloop for rel in relations
     do 
     (sloop for v on rel
	do (format t
		   (if (> (car v) 0) "~a" "~(~a~)")
		   (nth (abs (car v)) *names*))
	(cond ((null ro) (return-from my-print)))
	(if (cdr v) (princ (pop ro))
	  (format t "~a | ~a" i i))))))
     
(defun has-repeat (ar &aux (j (+ 1 (ncosets)))  ans tem)
  (sloop for k from 1 to (ncosets)
     do (setq tem (aref ar k))
     (cond ((and  (not (eql tem 0))
		  (find tem ar :start (+ k 1) :end j))
	    (pushnew tem ans))))
     ans)

(defun dcheck-tables ( &aux tem )
  (when *debug*
    (with-multiply-table
    (sloop for i from 1 to (nvars)
       do (if (setq tem (has-repeat (table i )) )
	      (format t "~%Table ~a has repeat ~a " i tem))))))
     
(defun dprint-state ()
  (when *debug*
    (with-multiply-table
	  (format t "~%Ncosets = ~a, *todo*=" (ncosets) *todo*)
	  (sloop for i from 1 to (nvars) do
	     (format t "~%~a:~a" (nth i *names*)
		     (subseq (table i ) 1  (+ 1 (ncosets)))))
	  (my-print (reverse *this-row*) 0)
	  ))
))