share/affine/modsimp.lisp

   1 ;;; -*- Mode: LISP; Package: cl-maxima; Syntax: Common-lisp -*-
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;                                                                    ;;;;;
   4 ;;;     Copyright (c) 1984 by William Schelter,University of Texas     ;;;;;
   5 ;;;     All rights reserved                                            ;;;;;
   6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   7 (in-package :maxima)
   8
   9 ;;then have a compiler optimizer which does the selection.
  10 (defun nn+ (&rest l)
  11    (cond ((null l) 0)
  12          ((cddr l)
  13           (n+ (car l) (apply 'nn+ (cdr l))))
  14          (t (n+ (car l)  (second l)))))
  15
  16 ;;a and b should be polynomials or rational-functions with the following
  17 ;;restrictions:  Any noncommutative variables or monomials occur with
  18 ;;higher pointergp than any scalars.  This will be the case with
  19 ;;add-newvar if you stick to using a,b,%zeta, etc for the scalars and
  20 ;;x,y,z,.. for the noncommutative variables.  All the nc parts thus will
  21 ;;come first so that a polynomial will be scalar if its main variable is
  22 ;;scalar.  General assumption is that nc polynomial will have all the nc
  23 ;;parts before any of the nc parts and no non commutative denominators.
  24 ;;We also assume that  a noncommutative polynomial looks like
  25 ;;   ncprod-or-var*cof+ncpolynomial  where cof is a scalar polynomial
  26 ;;probably the default should be that variables are scalar, but
  27 ;;at present the central variables must satisfy ($scalarp var) ==> t
  28
  29 (defun n. (a b)
  30    (cond ((numberp a)
  31           (n* a b))
  32          ((numberp b) (n* a b))
  33          ((atom b) (setq b (st-rat b))(n. a b))
  34          ((affine-polynomialp a)
  35           (cond ((affine-polynomialp b)
  36                  (let ((va (get (p-var a) 'disrep)) vb ncprod)
  37                    (cond (($scalarp va)
  38                           (ptimes a b))
  39                          (($scalarp (setq vb (get (p-var b) 'disrep)))
  40                           (ptimes a b))
  41                          (t (check-arg a (and (<=  (length a) 5)
  42                                               (<= (p-deg a) 1)) "an nc polynomial")
  43                             (check-arg b (and (<=  (length b) 5)
  44                                               (<= (p-deg b) 1)) "an nc polynomial")
  45                             (setq ncprod(add-newvar  (ncmul* va vb)))
  46                             (nn+ (ptimes (list ncprod 1 1) (ptimes (p-cof a) (p-cof b)))
  47                                  (cond ((cdddr a) (n. (fifth a) (subseq b 0 3)))
  48                                        (t 0))
  49                                  (cond ((cdddr b) (n.   (subseq a 0 3) (fifth b) ))
  50                                        (t 0))
  51                                  (cond ((and (cdddr a) (cdddr b))
  52                                         (n. (fifth a ) (fifth b)))
  53                                        (t 0)))))))
  54                 ((rational-functionp b) (nred (n. a (num b)) (denom b)))
  55                 (t (n. a (st-rat b)))))
  56          ((rational-functionp a)
  57           (nred (n. (num a) b) (denom a)))
  58          (t  (setq a (st-rat a))
  59              (n. a b)
  60              )))
  61
  62 (defun $find_sygy (degn &aux monoms general-term eqns sols simp)
  63   (setq monoms ($mono $current_variables (1- degn)))
  64   (setq monoms
  65         (cons '(mlist)
  66               (loop for v in (cdr monoms)
  67                  collecting (ncmul* v '$x)
  68                  collecting (ncmul* v '$y))))
  69   (mshow monoms)
  70   (setq general-term ($general_sum monoms $aaaa))
  71   (setq simp ($ratdisrep ($dotsimp general-term)))
  72   (mshow simp)
  73   (mshow ($list_nc_monomials simp))
  74   (setq simp ($ratsimp ($numerator simp)))
  75   (setq eqns ($extract_linear_equations `((mlist) ,simp)))
  76   (setq sols ($fast_Linsolve eqns ($list_variables eqns "aa")))
  77   (mshow eqns sols)
  78   ($separate_parameters ($sublis sols general-term)))
  79
  80
  81 (defun $find_sygy (degn matrix
  82                    &aux monoms  eqns sols pref general-image  general-from from-n to-n simp scal)
  83   " find_sygy(2,matrix([x],[y]))===> pq"
  84   (setq from-n ($length matrix))
  85   (setq to-n ($length (second matrix)))
  86   (setq monoms ($mono $current_variables degn))
  87   (setq pref (loop for i
  88                 from 1 to from-n
  89                 collecting ($concat '$aa i)))
  90   (setq general-from
  91         (cons '(mlist) (loop for u in pref      do
  92                             (setq scal (loop for i below ($length monoms)
  93                                           collecting ($concat u i)))
  94                             ($declare_scalar_list (setq scal(cons '(mlist) scal)))
  95                           collecting ($general_sum monoms scal))))
  96   (mshow general-from)
  97   (setq general-image (ncmul* general-from matrix))
  98   (mshow general-image)
  99   (break t)
 100   (setq simp ($ratdisrep ($dotsimp general-image)))
 101   (displa simp)
 102   (mshow simp)
 103   (mshow ($list_nc_monomials simp))
 104   (setq simp ($ratsimp ($numerator simp)))
 105   (cond (($listp simp) nil)
 106         (($matrixp simp)(setq simp ($list_matrix_entries simp)))
 107         (t (setq simp (list '(mlist) simp))))
 108   (setq eqns ($extract_linear_equations simp))
 109   (setq sols ($fast_Linsolve eqns ($list_variables eqns "aa")))
 110   (mshow eqns sols)
 111   ($separate_parameters ($sublis sols general-from)))
 112
 113 (defmacro matrix-row (mat n)`(nth ,n ,mat))
 114
 115 (defun matrix-column (mat n)
 116  (cons '(mlist) (loop for v in (cdr mat)
 117         collecting (nth n v))))
 118
 119 (defvar $module_simp_for_kernel nil)
 120 (defun $find_kernel (in-deg matrix &optional multiply-on-left
 121                    &aux first eqns sols general-image $general_from   simp  )
 122 "if matrix is m x n then it finds the kernel of the MAPL on the right
 123 multiplication B:A^m-->A^n     v |---> vB or   the kernel of the left
 124 multiplication A^m <--- A^n:B    Bv<----|v"
 125   (cond (multiply-on-left (setq first (matrix-row matrix 1)))
 126         (t
 127          (setq first (matrix-column matrix 1))))
 128
 129   ;;in-deg refers to degree of image in 1-1 position
 130   (setq $general_from (cons '(mlist)
 131                             (loop for u in (cdr first)
 132                                   for i from 1
 133                            collecting ($scalar_sum ($concat '$ff i)
 134                                                    ($mono $current_variables
 135                                                           (- in-deg ($nc_degree u)))))))
 136   (mshow $general_from)
 137   (cond (multiply-on-left
 138          (setq general-image (ncmul*   matrix ($transpose  $general_from) )))
 139         (t      (setq general-image (ncmul*  $general_from matrix))))
 140   (mshow general-image)
 141   (setq simp ($ratdisrep ($dotsimp general-image)))
 142   (cond ($module_simp_for_kernel (format t "~%Simplifying by ~A" $module_simp_for_kernel)
 143          (setq simp($totaldisrep (funcall $module_simp_for_kernel simp)))))
 144   (mshow simp)
 145   (setq simp ($ratsimp ($numerator simp)))
 146   (cond (($listp simp) nil)
 147         (($matrixp simp)(setq simp ($list_matrix_entries simp)))
 148         (t (setq simp (list '(mlist) simp))))
 149   (setq eqns ($extract_linear_equations simp))
 150   (setq sols ($fast_Linsolve eqns ($list_variables eqns "ff")))
 151   (mshow eqns sols)
 152   ($separate_parameters ($sublis sols $general_from) "par" "ff"))
 153
 154 (defun $find_basis (list-vectors &aux sum eqns sols rk)
 155   (setq sum ($scalar_sum '$aa list-vectors))
 156   (setq eqns ($ratsimp ($numerator ($ratdisrep ($dotsimp sum)))))
 157   (setq eqns ($extract_linear_equations eqns))
 158   (setq sols ($fast_Linsolve eqns ($list_variables eqns "aa")))
 159   (setq rk (sp-number-of-pivots (pv-the-sparse-matrix $poly_vector)))
 160   ($separate_parameters ($sublis sols sum) "par" "aa"))
 161
 162 (defun $scalar_sum (prefix monoms &aux scal )
 163 ;  (declare (values general-sum scalars))
 164   (setq scal (loop for i below ($length monoms)
 165                    collecting ($concat prefix i)))
 166   ($declare_scalar_list (setq scal(cons '(mlist) scal)))
 167   (values ($general_sum monoms scal) scal))
 168
 169 (defun replace-parameters-by (expr prefix &aux olds pars tem )
 170   (setq pars (cdr ($list_variables expr "par")))
 171   (setq olds ($list_variables expr prefix))
 172   (loop for i from 0 while pars
 173         when (not (member (setq tem ($concat prefix i)) olds :test #'equal))
 174           collecting (cons (car pars) tem) into repl
 175         and
 176         do (setq pars (cdr pars))
 177         finally (cond (repl
 178                        (loop for v in repl do (format t "~%Replacing ~A by ~A" (car v) (cdr v)))
 179                        (return (sublis repl expr)))
 180                       (t (return expr)))))
 181
 182 (defun $leftsocle (in-degree module-relations rad &aux h newh relats allrelats
 183                    monoms eqns vari solns  hh condition tem)
 184   "currently the modulue is a cyclic left module M=A/A.module-relations and
 185    we test to find solutions for h of v.h-A.module-relations for v in rad
 186   everything is supposed to be homogeneous"
 187   (setq h ($scalar_sum '$bbb ($mono $current_variables in-degree)))
 188   (setq relats
 189   (loop for v in (cdr module-relations)
 190         for i from 0
 191         do
 192         (setq tem  ($scalar_sum ($concat '$ccc i)
 193                                 ($mono $current_variables
 194                                        (1+ (- in-degree ($nc_degree v))))))
 195         collecting (ncmul* tem v)))
 196  (setq allrelats (meval* (cons '(mplus) relats)))
 197   (mshow allrelats)
 198   (setq monoms ($mono $current_variables in-degree))
 199   ;;impose condition x.h-general-sum module-relations
 200   (loop for v in (cdr rad)
 201         with all-rel = (num (new-rat  allrelats))
 202         do (setq v (st-rat v))
 203            (setq hh (n. v (num (new-rat h))))
 204            (setq condition (n- hh all-rel))
 205 ;          (setq condition (sub* (ncmul* v h)
 206 ;                                allrelats))
 207 ;          (setq conditon (num (new-rat
 208 ;         (setq condition ($numerator condition))
 209 ;          (setq condition ($ratdisrep ($numerator ($dotsimp condition))))
 210         (setq condition (new-rat-dotsimp (cons condition 1)))
 211         (displa condition)
 212         (setq condition ($totaldisrep ($numerator condition)))
 213         (setq eqns ($extract_linear_equations (list '(mlist) condition)))
 214         (setq vari ($list_variables eqns "bbb" "ccc" ))
 215         (setq solns ($fast_linsolve eqns vari))
 216         (setq newh ($ratsimp($sublis solns h)))
 217         (cond (($zerop newh)
 218                (return '((mlist)))))
 219         (setq h (replace-parameters-by newh '$bbb))
 220         (displa h)
 221         finally (return ($separate_parameters newh "bbb" "par"))))
 222
 223 (defun $try_bezout( f g  &optional (linears   #$[z-a*x,y-b*x]$) &aux (deg ($nc_degree f))  eqns cofs cof-eqns answ det tem)
 224   "solve f=f1*(z-a*x)+f2*(z-b*x),g=g1*(z-a*x)+g2*(z-b*x)"
 225   (setq cofs (loop for v in '(f1 f2 g1 g2)
 226                      do (setq tem
 227                               ($scalar_sum ($concat '$ccc v)
 228                                            ($mono $current_variables (- deg 1))))
 229                      collecting tem))
 230   (setq eqns
 231         (loop for (aa bb) on cofs by #'cddr
 232               for ff in (list f g)
 233               collecting
 234                 (loop  for lin in (cdr linears)
 235                   for acof in (list aa bb)
 236                   with eqn = ff
 237                   do (setq eqn (sub* eqn (ncmul* acof lin)))
 238                   finally (return eqn))))
 239   (setq eqns (loop for v in eqns collecting ($ratdisrep ($dotsimp v))))
 240   (setq cof-eqns ($extract_linear_equations (cons '(mlist) eqns) ($mono $current_variables deg)))
 241   (mshow cof-eqns)
 242   (setq answ ($fast_linsolve cof-eqns ($list_variables cof-eqns "cc")))
 243   (setq det (sp-determinant (pv-the-sparse-matrix $poly_vector)))
 244   (list '(mlist) answ (new-disrep det)))
 245
 246 (defun $intersect_principals_variety ($f1 $f2 in-deg &optional in-free-ring
 247                               &aux  $g1 $g2 eqn eqns )
 248   "Af1 intersect f2A yielding g1.f1=f2.g2, and returns a equations.."
 249   (declare (special $g1 $g2 $f1 $f2))
 250          (setq $g1 ($scalar_sum '$gg1 ($mono $current_variables (- in-deg ($nc_degree $f1)))))
 251          (setq $g2 ($scalar_sum '$gg2 ($mono $current_variables (- in-deg ($nc_degree $f2)))))
 252                  (setq eqn #$[g1.f1-f2.g2]$)
 253          (cond (in-free-ring nil)
 254                (t  (setq eqn ($totaldisrep ($dotsimp eqn)))))
 255         (setq eqns ($extract_linear_equations eqn)))
 256
 257 (defun $intersect_principals ($f1 $f2 in-deg &optional modulo-free-ring
 258                               &aux free-result result $g1 $g2 eqn eqns solns)
 259   "Af1 intersect f2A yielding g1.f1=f2.g2, and returns a matrix([g1,g2],[g1p,g2p],...)"
 260   (declare (special $g1 $g2 $f1 $f2 $eqns))
 261   (cond (modulo-free-ring
 262          (let ($dot_simplifications )
 263            (setq free-result  ($intersect_principals $f1 $f2 in-deg)))
 264          (setq free-result ($totaldisrep ($dotsimp free-result)))
 265          (setq result  ($intersect_principals $f1 $f2 in-deg))
 266          ($nc_matrix_quotient result free-result))
 267         (t
 268          (setq $g1 ($scalar_sum '$gg1 ($mono $current_variables (- in-deg ($nc_degree $f1)))))
 269          (setq $g2 ($scalar_sum '$gg2 ($mono $current_variables (- in-deg ($nc_degree $f2)))))
 270          (setq eqn #$totaldisrep(dotsimp([g1.f1-f2.g2]))$)
 271          (mshow eqn)
 272          (setq $eqns (setq eqns ($extract_linear_equations eqn)))
 273          (setq solns ($fast_linsolve eqns ($list_variables eqns "gg")))
 274          (setq result ($sublis solns #$[g1,g2]$))
 275          (cons '($matrix) (cdr ($separate_parameters result "par" "gg")))
 276          )))
 277
 278 ;(defun $find_basis_module (in-deg matrix-relations &optional right)
 279 ;  "given matrix([[f1,f2,..,fr],[g1,g2,..,gr],...]) compute a basis
 280 ;  for the free module A^r modulo the left(resp right) submodule generated by the
 281 ;  rows of matrix.  Each row must satisfy nc_degree(fi)= nc_degree(fj)."
 282
 283 (defvar $module_simplifications nil)
 284
 285 (defun module-monom-must-replacep (monom &aux repl)
 286   (cond ((atom monom)
 287          (setq repl (member monom (cdr $module_simplifications) :test #'eq))
 288          (cond (repl (list monom (second repl)))))
 289         (t
 290          (loop for v on (cdr $module_simplifications) by #'cddr
 291                with leng =  (length monom)
 292                with leng-first-v
 293                when
 294                  (or (and (listp (car v))
 295                           (>= leng (setq leng-first-v(length (car v))))
 296                           (equal (nthcdr (- leng leng-first-v) (cdr monom)) (cdar v)))
 297                      (eql (car v) (nth (1- (length monom)) monom)))
 298                  do (return (list (first v) (cond ((numberp (second v)) (second v))
 299                                                   (t (cdr (second v))))))))))
 300 ;(defun module-must-replacep (poly)
 301 ;  (*catch 'found-replace
 302 ;    (mod-must-replace1 poly)))
 303
 304 (defun module-must-replacep (poly &aux monom)
 305   (loop while poly
 306         do
 307     (cond ((atom poly) (return nil))
 308           ((module-monom-must-replacep (setq monom (get (car poly) 'disrep)))
 309            (return (cons (car poly)   (module-monom-must-replacep monom))))
 310           (t (setq poly (fifth poly))))))
 311
 312 (defun rat-simplifications (simps)
 313  (cons '(mlist) (loop for (mon rep) on (cdr simps ) by #'cddr
 314         collecting mon
 315         collecting (cons *fake-rat* (new-rat rep)))))
 316
 317 (defun module-simp (fun &aux did-replace num-fun den-fun
 318                     mon actual-repl repl den prev-mon lft symbo num-repl     repl-data )
 319   "argument satisfies rational-functionp"
 320   (cond ((affine-polynomialp fun)(setq num-fun fun)(setq den-fun 1))
 321         ((rational-functionp fun)(desetq (num-fun . den-fun) fun))
 322         (t (fsignal "takes poly or rat'l function for argument")))
 323 ;  (show $module_simplifications)
 324   (loop while (setq repl-data (module-must-replacep num-fun)) do
 325     (desetq (symbo mon repl) repl-data)
 326     (setq prev-mon (get symbo 'disrep))
 327     (setq did-replace t)
 328     (setq lft (dot-right-quotient prev-mon mon))
 329     (setq actual-repl (n. (st-rat lft) repl))
 330     (cond ((affine-polynomialp actual-repl) (setq num-repl actual-repl) (setq den 1))
 331           ((rational-functionp actual-repl)(desetq (num-repl .  den) actual-repl))
 332           (t (fsignal "bad repl")))
 333     (format t "~%Module simp replacing ~A" (get symbo 'disrep))
 334     (setq num-fun (psublis(list  (cons symbo num-repl)) den num-fun))
 335     (setq den-fun (n* den den-fun))
 336         finally (return (cond (did-replace (ratreduce num-fun den-fun))
 337                       (t fun)))))
 338
 339 (defun delete-from-simps (simps mon repl)
 340   (loop for v on (cdr simps) by #'cdr
 341         when (and (equal (car v) mon)
 342                   (equal (second v) repl))
 343           do (return (cons '(mlist)(nconc tem (cddr v))))
 344           else
 345           nconc (subseq v 0 2) into tem))
 346
 347
 348 ;;;needs work:note that we have not been using new-rat-dotsimp and we need to!
 349 ;;;maybe the best plan is to sort the simps in increasing order of degree,
 350 ;;;then to (let the modsimps be just one term and one replacement and check
 351 ;;;alll remaining ones for subs.  Any replacement will be further down the
 352 ;;;list.   Also if you do a replacement you could not get earlier in the list
 353 ;;;if  they are sorted by degree, and within degree by worst monom
 354 ;(defun simplify-module-simplifications ()
 355 ;  (loop named sue for (mon rep) on (cdr $dot_simplifications) by #'cddr
 356 ;          do
 357 ;         (let (($module_simplifications (delete-from-simps $module_simplifications mon repl)))
 358 ;             (loop for (mon1 rep1) on (cdr $dot_simplifications) by #'cddr
 359 ;                   when (or (module-monom-must-replacep mon1)
 360 ;                            (module-must-replacep rep1))
 361 ;                     do
 362 ;               (setq relat (function-numerator (module-simp (make-relation mon1 rep1))))
 363 ;                (setq relat (module-simp relat))
 364 ;               (setq changed t)
 365 ;               (return-from sue 'start-over))
 366 ;               (setq new (plain-add-to-module-simps relat)))))
 367
 368 (defun make-relation (monom repl)
 369 ;  (iassert (or (numberp repl) (eql (caar repl) 'mrat)))
 370   (cond ((numberp repl) nil)
 371         ((affine-polynomialp repl)(setq repl (cons repl 1)))
 372         ((rational-functionp repl) nil)
 373         ((eql (caar repl) 'mrat) (setq repl (cdr repl)))
 374         (t (fsignal "bad repl")))
 375   (function-numerator  (n- monom repl)))
 376
 377
 378 (defun sort-simps-by-degree (simps)
 379  (cons (car simps) (sort-grouped-list (cdr simps) 2 #'(lambda (u v)
 380                                                         (cond ((atom u))
 381                                                          ((atom v) nil)
 382                                                          ((< (length u) (length v)) )
 383                                                          ((> (length u) (length v)) nil)
 384                                                          (( funcall $order_function  v u)))))))
 385
 386
 387 (defun module-and-dot-simp (pol  &aux (changed  t))
 388   (loop
 389     do
 390     (cond ((module-must-replacep pol)
 391            (setq pol (module-simp pol)) (setq changed t)))
 392     (cond (($must_replacep pol)
 393           (setq pol (cdr (new-rat-dotsimp pol)))
 394           (setq changed t)))
 395     while changed
 396     finally (return pol)))
 397
 398 (defun gen-dot-check-associative (a b c &optional (simplifier 'module-and-dot-simp) &aux answ)
 399   "simplifier should return a polynomial or rational function"
 400  (setq answ (n- (n. (funcall simplifier (n. a b)) c)
 401       (n. a (funcall simplifier a b))))
 402  (cond ((pzerop answ) (format t "~%It was associative"))
 403        (t answ)))
 404
 405
 406 ;(defun find-right-overlap (left right &aux answ)
 407 ;  (setq answ
 408 ;       (cond   ((atom left)
 409 ;
 410 ;                (cond ((eql left right)(list nil left nil))
 411 ;                      ((equal left (car (last right)))
 412 ;                       (list  (butlast right) left nil))
 413 ;                      (t nil)))
 414 ;               ((atom right)
 415 ;                (cond ((equal right (car (last left)))
 416 ;                       (list (butlast left) right nil))
 417 ;                      (t nil)))
 418 ;               (t
 419 ;                (multiple-value-bind (lap right-excess)
 420 ;                    (find-lap1 (cdr left) (cdr right))
 421 ;                (cond (lap (list (firstn (- (length left)  1 (length lap)) (cdr left))   lap right-excess)))))))
 422 ;  (show answ)
 423 ;  (and answ
 424 ;
 425 ;
 426 ;       (loop for v on answ
 427 ;            for a in answ
 428 ;            do
 429 ;
 430 ;        finally (return answ))))
 431
 432 (defun coerce-nctimes (a)
 433   "nil-->1, '(a) -->a '(x y z) --> x.y.z , x.y--> x.y"
 434          (cond ((null a)
 435                    1)
 436                ((listp a)
 437                 (cond ((listp (car a))
 438                        (cond ((eql (caar a) 'mnctimes)
 439                               (cond ((null (cdr a)) 1)
 440                                     ((cddr a) a)
 441                                     (t (second a))))
 442                              (t (fsignal "what am i like?"))))
 443                       ((not (cdr a))
 444                        (first a))
 445                       (t  (cons '(mnctimes) a))))
 446                (t a)))
 447
 448 (defun $lastn_ncdegree (mon n)
 449   "Returns an final segment of mon of degree n if possible, otherwise nil"
 450   (cond ((zerop n) 1)
 451         ((atom mon)
 452          (and (equal ($nc_degree mon) n)
 453               mon))
 454         (t (loop
 455                  for i downfrom (1- (length mon)) to 1
 456                  for u = (nth i mon)
 457                  summing ($nc_degree u) into tot
 458                  when (>= tot n)
 459                    do
 460                      (cond ((eql tot n)
 461                             (return (coerce-nctimes (nthcdr i mon))))
 462                            (t (return nil)))))))
 463
 464 (defun $firstn_ncdegree (mon n)
 465   "Returns an initial segment of mon of degree n if possible, otherwise nil"
 466   (cond ((zerop n) 1)
 467         ((atom mon)
 468          (and (equal ($nc_degree mon) n)
 469               mon))
 470         (t (loop for u in (cdr mon)
 471                  for i from 2
 472                  summing ($nc_degree u) into tot
 473                  when (>= tot n)
 474                    do
 475                      (cond ((eql tot n)
 476                             (return (coerce-nctimes (subseq mon 0 i))))
 477                            (t (return nil)))))))
 478
 479
 480
 481 (defun $find_right_lap (left right tot-deg &aux size-lap a-lap )
 482   "tot-deg is the degree of the lapped left and right. It will not catch left a subset
 483    of right since it is presumed the right is in reduced form "
 484   (let ((ld ($nc_degree left))
 485         (rd ($nc_degree right)))
 486     (setq size-lap (- (+ ld rd) tot-deg))
 487     (cond ((<= size-lap 0) nil)
 488           (t (setq a-lap ($lastn_ncdegree left size-lap))
 489             (cond((and a-lap
 490                  (or
 491                    (equal a-lap right)
 492                    (cond ((listp right)
 493                           (cond ((listp a-lap)
 494                                  (initial-seg (cdr a-lap) (cdr right)))
 495                                 (t (equal (second right) a-lap)))))))
 496                    (list ($firstn_ncdegree left (- ld size-lap))
 497                          a-lap
 498                          ($lastn_ncdegree right (- rd size-lap)))))))))
 499
 500 (defun $find_right_lap (left right tot-deg &aux size-lap a-lap )
 501   "tot-deg is the degree of the lapped left and right. It will not catch left a subset
 502    of right since it is presumed the right is in reduced form "
 503   (let ((ld ($nc_degree left))
 504         (rd ($nc_degree right)))
 505     (setq size-lap (- (+ ld rd) tot-deg))
 506     (cond ((<= size-lap 0) nil)
 507           (t (setq a-lap ($lastn_ncdegree left size-lap))
 508            (cond ((and a-lap  ;;have match
 509                     (cond ((equal a-lap right))
 510                           ((and (listp a-lap) (listp right)
 511                               (initial-seg (cdr a-lap) (cdr right))))
 512                           ((and (atom a-lap) (listp right)
 513                                 (equal (second right) a-lap)))))
 514                 (list ($firstn_ncdegree left (- ld size-lap))
 515                       a-lap
 516                       ($lastn_ncdegree right (- rd size-lap)))))))))
 517
 518 ;
 519 ;(defun $find_right_lap (left right tot-deg &aux size-lap a-lap )
 520 ;  "tot-deg is the degree of the lapped left and right. It will not catch left a subset
 521 ;   of right since it is presumed the right is in reduced form "
 522 ;  (let ((ld ($nc_degree left))
 523 ;       (rd ($nc_degree right)))
 524 ;    (setq size-lap (- (+ ld rd) tot-deg ))
 525 ;    (cond ((<= size-lap 0) nil)
 526 ;         (t (setq a-lap ($lastn_ncdegree left size-lap))
 527 ;;          (cond( a-lap
 528 ;;          (setq b-lap ($firstn_ncdegree right size-lap))))
 529 ;            (cond
 530 ;              ((and a-lap
 531 ;                    (or (equal a-lap right)
 532 ;                        (and (listp a-lap) (listp right)
 533 ;                             (initial-seg (cdr a-lap) (cdr right)))
 534 ;                        (listp right) (equal (second right) a-lap))
 535 ;                             )
 536 ;;              (equal (cdr a-lap) (cdr b-lap)))))
 537 ;               (list ($firstn_ncdegree left (- ld size-lap))
 538 ;                     a-lap
 539 ;                     ($lastn_ncdegree right (- rd size-lap)))))))))
 540
 541 (defun initial-seg (seg lis)
 542   "'(1 2) '(1 2 3) --> values: t '(3), '(1 2) '(0 1 2 3) --> nil, '(1 2) '(1) --> nil"
 543   (prog ()
 544         a
 545         (or (eql (car seg) (car lis)) (return nil))
 546         (setq lis (cdr lis))
 547         (setq seg (cdr seg))
 548         (cond ((null seg) (return (values t lis))))
 549         (and (null lis)  (return nil))
 550         (go a)))
 551
 552
 553
 554
 555 (defun plain-add-to-module-simps (poly)
 556   (cond ($module_simplifications nil)
 557         (t (setq $module_simplifications '((mlist)))))
 558            (setq $module_simplifications
 559                  (append $module_simplifications  (make-simp (function-numerator poly)))))
 560
 561
 562 ;;must do two things
 563 ;;I reduce module-simps leading terms wrt themselves.
 564 ;;II If M is a leading term of f and DM is a dotsimp replacement,
 565 ;; and there is a dotsimp repl R and the tail of R overlaps the left side of M
 566 ;;then we have two ways of simplifying and need to add the difference.
 567 ;;add in the simp dotsimp(Df).  Note this is essentially checking for overlaps
 568 ;;between the module simps and the dotsimps.
 569 ;;Have two functions: I one which checks for type I and if finds, returns the (changed list of simps , t)
 570 ;;else returns simps,nil
 571 ;;II one which checks for the II type and returns the new list of simps ,t or just old list of simps.
 572
 573
 574 (defun make-simp (relat &aux answ)
 575   "inverse of make-relation, except that constants will be cancelled."
 576  (setq answ (cond ((numberp relat)
 577          (cond ((zerop relat) nil)
 578                (t (list 1 0))))
 579         (t (list (get (car relat) 'disrep)
 580            (cond ((fifth relat)
 581                    (cons *fake-rat* (ratreduce (pminus (fifth relat)) (third relat))))
 582                  (t 0))))))
 583
 584  answ)
 585
 586 (defun dot-right-quotient (a b &aux tem answ )
 587   "a.b^^-1 with error if b not final segment"
 588   (cond ((atom b)
 589          (cond ((atom a)
 590                 (cond ((eql a b) 1)
 591                       (t (fsignal "doesn't divide evenly"))))
 592                ((equal (car (last a)) b)
 593                 (setq answ (butlast a)))
 594                (t (fsignal "not tail"))))
 595         (t
 596          (loop for v in (cdr b)
 597                for w in (nthcdr (setq tem  (1+ (- (length a) (length b)))) a)
 598 ;              do ;(show v w tem)
 599                when (not (equal v w))
 600                  do (fsignal "b is not a tail of a"))
 601          (setq answ (subseq a 0 tem))))
 602   (cond ((cddr answ) answ)
 603         ((cdr answ) (second answ))
 604         (t 1)))
 605
 606 (defvar *module-overlaps-checked* nil)
 607
 608 (defun simplify-module-simps (to-deg &aux changed new-simps ( present-deg 0))
 609   (declare (special present-deg))
 610   (user-supply *module-overlaps-checked*)
 611   (loop do (setq changed nil)
 612        (loop for f in '(type-i-simp type-ii-simp) do
 613             (multiple-value-setq
 614                 (new-simps changed)
 615               (funcall f present-deg to-deg ))
 616             (setq $module_simplifications new-simps)
 617             (cond (changed (return 'try-again))))
 618      while changed
 619      finally (return $module_simplifications)))
 620
 621 (defun type-i-simp (&rest ignore &aux old-simps relat new-relat)
 622   (declare (ignore ignore))
 623   "eliminate simps such that one leading term is a left multiple of another, ie common right overlap"
 624   (setq old-simps (setq $module_simplifications (sort-simps-by-degree $module_simplifications)))
 625   (loop named sue for tail-simps on (cdr old-simps) by #'cddr
 626         for i from 3 by  2
 627         do
 628     (let (($module_simplifications (cons '(mlist) (subseq tail-simps 0 2))))
 629       (loop for (mon repl) on (cddr tail-simps) by #'cddr
 630             for jj from (+ i 2) by 2
 631             when (or  (module-monom-must-replacep mon)
 632                       (and (listp repl)
 633                            (module-must-replacep (function-numerator (cdr repl)))))
 634               do (setq relat  (make-relation mon repl))
 635                  (setq new-relat (module-simp relat))
 636                  (setq new-relat (new-rat-dotsimp new-relat))
 637                  (cond ((numberp new-relat)
 638                         (cond ((zerop new-relat) nil)))
 639                        (t (setq new-relat (cdr new-relat))
 640                           (iassert (rational-functionp new-relat))
 641                           (setq new-relat (num new-relat))))
 642                  (return-from sue (values (nconc (subseq old-simps 0 (- jj 2))
 643                                                  (make-simp new-relat)
 644                                                  (nthcdr jj old-simps))
 645                                           t))
 646             finally (return-from sue (values old-simps nil))))))
 647
 648 (defun dotsimp (ratl-or-poly)
 649 ;  (declare (values ratl-function))
 650     (new-rat  (new-rat-dotsimp ratl-or-poly)))
 651
 652 (defun ratl-function-id (number-or-rational-function &aux tem)
 653   (cond ((affine-polynomialp number-or-rational-function)
 654          number-or-rational-function)
 655         ((rational-functionp number-or-rational-function)
 656          number-or-rational-function)
 657         ((and (listp number-or-rational-function)
 658               (listp (setq tem (car number-or-rational-function)))
 659               (eql (car tem) 'mrat))
 660          (cdr number-or-rational-function))
 661         (t (fsignal "not right type"))))
 662
 663 (defun type-ii-simp (from-deg to-deg &aux old-simps new-repl tem tem1 a b c)
 664   "checks overlaps of type (a b c) where a.b is a dotsimp and b.c is a modsimp and deg(a.b.c)<=i
 665    then replacement(a.b).c is a module relation to be added."
 666   (declare (special present-deg))
 667   (setq old-simps (setq $module_simplifications (sort-simps-by-degree $module_simplifications)))
 668   (loop for i from from-deg to to-deg do
 669        (block sue
 670          (setq present-deg i)
 671          (loop for (mon repl) on (cdr $dot_simplifications) by #'cddr do
 672               (loop for (monmod repmod) on (cdr old-simps) by #'cddr
 673                  when (setq tem ($find_right_lap mon monmod i))
 674                  do
 675                  ;;(show (list mon monmod))
 676                  ;;should really use
 677                  ;;(setq tem (list mon monmod i)) ;but only maximal overlaps need to be checked and since we go up in degree
 678                    (setq tem1 (cons mon monmod))
 679                    (cond ((not (member tem1 *module-overlaps-checked* :test #'equal))
 680                           (push tem1 *module-overlaps-checked*)
 681                           (iassert (nc-equal (ncmul* (first tem) (second tem)) mon))
 682                           (iassert (nc-equal (ncmul* (second tem) (third tem)) monmod))
 683                           (iassert (nc-equal ($nc_degree  (ncmul*  mon (third tem))) i))
 684                           (iassert (not (numberp (second tem))))
 685                           ;;(setq mod-prod (n. (first tem)(function-numerator repmod)))
 686                           ;;(setq new-repl (n-  (n. (function-numerator repl) (third tem)) mod-prod))
 687                           ;; rep.c - a.repmod is in the ideal defining the module!
 688                           (desetq (a b c) tem)
 689                           (setq new-repl (n- (n. (ratl-function-id repl) c)
 690                                              (n. a (ratl-function-id repmod))))
 691                           (setq new-repl (dotsimp new-repl))
 692                           (setq new-repl (function-numerator new-repl))
 693                           (setq new-repl (function-numerator (module-simp new-repl)))
 694                           (mshow new-repl)
 695                           (mapcar #'displa (make-simp new-repl))
 696                           (cond ((pzerop new-repl) nil)
 697                                 (t (return-from sue (values (append old-simps (make-simp new-repl)) t)))))))))
 698      finally (return (values old-simps nil))))
 699
 700 ;(defun type-II-simp ( to-deg &aux old-simps new-repl repl-data tem)
 701 ;  (setq old-simps (setq $module_simplifications (sort-simps-by-degree $module_simplifications)))
 702 ;  (loop for  (mon repl) on (cdr $dot_simplifications) by #'cddr
 703 ;       when (setq repl-data (module-monom-must-replacep mon))
 704 ;         do (cond ((not (member (setq tem (cons mon (car repl-data))) *module-overlaps-checked*))
 705 ;                   (push tem *module-overlaps-checked*)
 706 ;           (setq new-repl (n. (dot-right-quotient mon (car repl-data)) (apply 'make-relation repl-data)))
 707 ;            (setq new-repl (cadr (new-rat-dotsimp new-repl)))
 708 ;            (return (values (append old-simps (make-simp new-repl))
 709 ;                            t))))
 710 ;       finally (return (values old-simps nil))))
 711
 712 ;(let (($module_simplifications)(*module-overlaps-checked*))
 713 ;  (loop for v in (st-rat #$[x,-(y.z+x.x)]$)
 714 ;           do (plain-add-to-module-simps v))
 715 ;  ( simplify-module-simps))
 716
 717 (defun fake-header (form)
 718   (cond
 719          ((affine-polynomialp form)(cons *fake-rat* (cons form 1)))
 720          ((rational-functionp form) (cons *fake-rat* form))
 721          (t (fsignal "not poly or rat'l fun"))))
 722 (defun $modsimp(form)
 723   (cond ((atom form)(fake-header (module-simp (new-rat form))))
 724         ((or (affine-polynomialp form)
 725          (rational-functionp form))
 726          (fake-header (module-simp form)))
 727         ((mbagp form)(cons (car form) (mapcar '$modsimp (cdr form))))
 728         (t (fake-header (module-simp (new-rat form))))))
 729
 730 (defun $cyclic_module_basis (degree &optional (variables $current_variables) )
 731   (cyclic-module-basis variables degree ($replacements) ($module_replacements)))
 732
 733 (defun $module_replacements( &optional (module-simps $module_simplifications))
 734   (cons '(mlist) (loop for (mon rep) on (cdr module-simps) by #'cddr
 735                        collecting mon)))
 736
 737 (defremember  cyclic-module-basis (variables deg  dot-replacements module-replacements &aux tem)
 738   (cond ((eql deg 0)
 739          (cond ((not (or (member 1 dot-replacements :test #'equal)
 740                          (member 1 module-replacements :test #'equal)))
 741                 '((mlist) 1))
 742                (t '((mlist)))))
 743         (t
 744          (cons '(mlist)
 745                (loop for v in (cdr variables)
 746                       appending
 747                       (loop for w in   (cdr (cyclic-module-basis variables (1- deg) dot-replacements module-replacements))
 748                              do (setq tem (ncmul* v w ))
 749                              unless (or (module-monom-must-replacep tem) ($must_replacep tem))
 750                              collecting tem))))))
 751
 752
 753 (defun $module_dimensions (to-n &optional (variables $current_variables) &aux tem)
 754   (loop for i from 0 to to-n
 755         collecting (setq tem ($length ($cyclic_module_basis i variables))) into dims
 756         summing tem into tot
 757         do (format t "~%The dimension of the cyclic module in deg ~A is ~A." i tem)
 758            finally (format t "~%The total dimension through degree ~A is ~A." to-n tot)
 759                    (return (cons '(mlist) dims ))))
 760
 761
 762 (defun $set_up_module_simplifications(relations &optional to-deg)
 763   (setq $module_simplifications nil)
 764   (setq *module-overlaps-checked* nil)
 765   (loop for v in (st-rat relations) do
 766        (plain-add-to-module-simps  v))
 767   (cond (to-deg
 768          (simplify-module-simps to-deg)))
 769   $module_simplifications)
 770
 771 (defun $cyclic_module_basis (deg &optional (variables $current_variables))
 772  (cons '(mlist) (loop for v in (cdr  ($mono variables deg))
 773     when (not (module-monom-must-replacep v))
 774       collecting v)))
 775
 776 (defun $shift_sequence (seq n &aux (tem (make-list (abs n) :initial-element 0)))
 777   "shifts a sequence of degrees over by n to use for adding to check sums"
 778   (cond ((zerop n) seq)
 779         ((> n 0)
 780          (append ($rest seq  n) tem))
 781         (t ($append (cons '(mlist) tem) ($rest seq  n)))))
 782
 783 (defun $euler_sum (deg-seq size-of-module map-degs)
 784   (loop for v in (reverse (cdr map-degs))
 785         for s in (cdr size-of-module)
 786         for i from 1
 787         with ans = deg-seq
 788                  summing v into tot
 789         do (mshow ans)
 790         when (evenp i)
 791           do (setq ans (add* ans (mul* s ($shift_sequence deg-seq (- tot)))))
 792         else do (setq ans (sub* ans (mul* s ($shift_sequence deg-seq (- tot)))))
 793           finally (return ans)))
 794
 795 (defun $solve_sp (try eqns &aux  simp syst solns ld)
 796  (setq eqns  ($sublis try eqns))
 797  (setq ld  (make-ldata :eqns (mapcar 'function-numerator(st-rat  eqns))))
 798  (des ld)
 799  (setq simp (simplify-ldata ld))
 800  (des simp)
 801  (setq syst (loop for u in simp collecting (cons '(mlist) (mapcar 'new-disrep (ldata-eqns u)))))
 802  (setq solns (loop for v in syst collecting ($append try  ($fast_linsolve v ($list_variables v)))))
 803  (cons '(mlist) solns))
 804
 805
 806 (defun $nc_linear_factor (pol &aux answ (deg ($nc_degree pol)) lin fac2 result eqns solns general)
 807   (setq lin ($scalar_sum '$bb ($mono $current_variables 1)))
 808   (setq fac2 ($scalar_sum '$cc ($mono $current_variables (- deg 1))))
 809   (setq result (sub* pol (ncmul* lin fac2)))
 810   (setq result ($dotsimp result))
 811   (setq eqns  ($extract_linear_equations (list '(mlist) result)))
 812   (mshow eqns)
 813   (setq answ
 814         (loop
 815      with try
 816      do
 817           (format t "~%Supply a macsyma list of possible solns")
 818           (setq try (mread-noprompt))
 819           (mshow try)
 820           (setq solns ($solve_sp try eqns))
 821           (mshow solns)
 822           (setq solns ($ratsimp solns))
 823           (cond ((null (y-or-n-p "Try again?"))
 824                  (setq general (list '(mlist) lin fac2))
 825                  (return (cons '(mlist) (loop for v in (cdr solns)
 826                                               appending (cdr ($separate_parameters ($sublis v general) "bb" "cc" "par")))))))))
 827   (cond((y-or-n-p "Verify factors:?")
 828         (loop for v in (cdr answ)
 829               for i from 1
 830               with simp
 831               do (setq simp ($dotsimp (sub* pol  (ncmul* (meval*(second v)) (meval* (third v))))))
 832                  (cond ((pzerop simp)(format t "~%soln ~A is correct:" i)
 833                         (displa (list '(mequal) pol (cons '(mnctimes) (cdr v) ))))
 834                        (t(format t "~%soln ~A is not correct ****************"))))))
 835   answ)
 836
 837 (defun $degless (u v &aux (u-deg ($nc_degree u)) (v-deg ($nc_degree v)))
 838   (cond ((< u-deg v-deg))
 839         ((eql  u-deg v-deg) (funcall $order_function u v))
 840         (t nil)))
 841
 842