share/affine/dim4-help.lisp

   1 ;;; -*- Mode:Lisp; Package:CL-MAXIMA; Base:10 -*-
   2 (in-package :maxima)
   3
   4 (defun $complexity(v &optional sum)
   5   (cond ((eql v 0)(if sum -1 '$z))
   6         ((mbagp v)
   7          (cond (sum
   8                 (loop for u in  (cdr v)
   9                       sum ($complexity u sum)))
  10                (t (cons (car v) (mapcar '$complexity (cdr v))))))
  11         (t (gen-pcomplexity (st-rat v)))))
  12
  13 (defun complexity-difference1 (x y)
  14  (cond ((eq x '$z) (setq x -1)))
  15   (cond ((eq y '$z) (setq y -1)))
  16  (cond ((numberp x)
  17         (- x y))
  18        ((mbagp x)
  19         (loop for v in (cdr x) for w in (cdr y)
  20               sum (complexity-difference1 v w)))))
  21
  22
  23 (defun $complexity_difference(a-list v )
  24   (cond ((eql v 0) -1)
  25         ((let ((x ($complexity  ($psublis a-list v)))
  26                (y  ($complexity  v)))
  27          (complexity-difference1 y x)))))
  28
  29 (defun $pdegree(v va)
  30   (cond ((eql v 0) '$z)
  31         ((mbagp v)
  32          (cons (car v)
  33                (loop for w in (cdr v) collect ($pdegree w va))))
  34         (t (pdegree (st-rat v) (car (st-rat va))))))
  35                                                 ;
  36
  37 (defun $delete_matrix_row(mat i)
  38   (cons (car mat)
  39         (loop for v in (cdr mat) for j from 1
  40                unless (eql j i) collect v)))
  41
  42 (defun $complexity_less_p (x y) (< ($complexity x t) ($complexity y t)))
  43
  44 (defun $my_rat (x)
  45   (cond ((mbagp x)
  46          (cons (car x) (mapcar '$my_rat (cdr x))))
  47         (t (header-poly (st-rat x)))))
  48
  49 (defun $mat_entry(mat i j)
  50   (nth  j (nth  i mat)))
  51
  52 (defvar $pivot_list '((mlist)))
  53
  54 (defun $eliminate_pivot (mat i j)
  55   (setq mat ($my_rat mat))
  56   (let ((pivot ($mat_entry mat i j))
  57         (pivot-row (nth  i mat)))
  58     (format t "~%Using pivot:")
  59     (sh pivot)
  60    (setq $pivot_list ($append $pivot_list `((mlist) ,(header-poly pivot) ,i ,j)))
  61     (cons (car mat)
  62     (loop for v in (cdr mat) for ii from 1
  63           when (not (eql i ii))
  64           collect
  65           (let* ((this (st-rat (nth j v))))
  66             (cond ((eql 0 this) v)
  67                   (t (let ((quot (nred this pivot)))
  68                        (cons (car v)
  69                              (loop for u in (cdr v)
  70                                    for a in (cdr pivot-row)
  71                                    collect
  72                                    (header-poly
  73                                      (n- (n* (function-numerator quot) a)
  74                                          (n* (function-denominator quot) u)))))))))))))
  75
  76
  77 (defun $psublis (a-list poly &optional(denom 1) palist)
  78   "use psublis([y=x^2,v=u^3],denom,poly)"
  79   (or palist
  80       (setq palist
  81             (loop for (u v  repl)  in (cdr a-list) by  'cdr
  82                   do (check-arg u (eq (car u) 'mequal) "Type a=repl")
  83                   collecting
  84                   (cons (p-var (st-rat v))
  85                         (st-rat repl)))))
  86   (cond ((mbagp poly) (cons (car poly)
  87                             (loop for v in (cdr poly) collect
  88                                   ($psublis a-list v denom palist))))
  89         (t (header-poly(psublis palist (st-rat denom) (st-rat poly))))))
  90
  91
  92 (defun $eliminate_pivot_list_hack (matrix pivot-list)
  93   (loop for (piv i j) on (cdr pivot-list) by 'cdddr
  94         do
  95         (format t "~%Pivot in list was :") (displa  ($matrix_entry matrix i j))
  96         (setq matrix ($eliminate_pivot matrix i j))
  97         finally (return matrix)))
  98
  99
 100 (defun $get_relation_matrix(relations mons)
 101   (setq mons (st-rat mons))
 102   (cons '($matrix)
 103         (loop for v in (cdr relations)
 104               do (setq v (st-rat v))
 105               collect
 106               (cons '(mlist)
 107                     (loop for w in mons
 108                           collect (header-poly  (pcoeff v w)))))))
 109
 110 (defun $switch_rows(mat i j)
 111   (let ((nmat (copy-list mat)))
 112     (setf (nth i nmat) (nth j mat))
 113     (setf (nth j nmat) (nth i mat))
 114     nmat))
 115
 116 (defun number_zeros(a)
 117   (loop for v in a when (eql v '$z) count t))
 118
 119 (defun $row_less (a b) (> (number_zeros a) (number_zeros b)))
 120
 121 (defun $row_sort (mat pred)
 122   (cons '($matrix) (sort (copy-list (cdr mat)) pred)))
 123
 124 (defun $reorder_matrix(mat &aux rows)
 125   (setq rows
 126         (loop for u in (cdr mat)
 127         collect (cons (loop for v in (cdr u) sum (gen-pcomplexity (st-rat v))) u)))
 128   (setq rows (sort rows #'< :key #'car))
 129   (cons '($matrix) (mapcar 'cdr rows)))
 130
 131 (defun $best_row(mat &aux  tem at)
 132         (loop for u in (cdr mat) for i from 1
 133               minimize(setq tem (loop for v in (cdr u) sum (gen-pcomplexity (st-rat v)))) into the-min
 134               when (eql tem the-min)
 135               do (setq at i)
 136               finally (return at)))
 137
 138
 139 (defun $best_piv (row &aux tem at )
 140   (loop for v in (cdr row) for i from 1
 141         unless (eql (setq v (st-rat v)) 0)
 142         minimize (setq tem (gen-pcomplexity (st-rat v))) into the-min
 143         and
 144         when (eql tem the-min) do (setq at  i)
 145         finally (return at)))
 146
 147 (defvar $current nil)
 148
 149 ;;sort by column.
 150 (defun $invertible_pivots(mat g &aux pivs)
 151   (setq g (st-rat g))
 152   (loop for ro in (cdr ($transpose mat)) for i from 1
 153          append
 154          (loop for u in (cdr ro) for j from 1
 155                    do (setq u (st-rat u))
 156                    when (may-invertp u g)
 157                    collect (list '(mlist) (header-poly u) j i ($complexity ro t) (pcomplexity u)))
 158          into all
 159          finally (setq pivs (sort all '< :key #'(lambda (x) (nth 4 x))))
 160          (return (cons '(mlist) pivs))))
 161
 162 ;(defun $invertible_pivots(mat g &aux pivs)
 163 ;  (setq g (st-rat g))
 164 ;  (loop for ro in (cdr mat) for i from 1
 165 ;        append
 166 ;        (loop for u in (cdr ro) for j from 1
 167 ;                  do (setq u (st-rat u))
 168 ;                  when (may-invertp u g)
 169 ;                  collect (list '(mlist) (header-poly u) i j ($complexity ro t) (pcomplexity u)))
 170 ;        into all
 171 ;        finally (setq pivs (sort all '< :key '(lambda (x) (nth 4 x))))
 172 ;        (return (cons '(mlist) pivs))))
 173
 174 (defun $reduce_by_complexity (&optional mat )
 175   (setq $current (or mat $current))
 176   (let* ((i ($best_row mat))
 177          (j ($best_piv (nth i mat))))
 178     (print (list i j))
 179     (setq $current ($eliminate_pivot mat i j))))
 180
 181
 182 (defun $sort_complexity(list)
 183   (loop for v in (cdr list)
 184          collect (cons (pcomplexity (st-rat v)) v) into all
 185          finally (setq all (sort all '< :key 'car))
 186          (return (cons (car list) (mapcar 'cdr all)))))
 187
 188 (defvar $gcds_used '((mlist)))
 189
 190 (defun $p_projective (lis &optional agcd &aux (zero-row t))
 191   "Cancel a gcd of AGCD and The rest of the elements of LIS.  If LIS is a
 192 list of lists call this function on each of the lists independently."
 193   (assert (mbagp lis))
 194   (and agcd (setq agcd (st-rat agcd)))
 195   (cond ((mbagp (second lis))
 196          (cons (car lis)
 197                (loop for w in (cdr lis)
 198                      collect ($p_projective w agcd))))
 199         (t
 200          (loop for v in (cdr lis)
 201                do (setq v (st-rat v))
 202                (assert (affine-polynomialp v))
 203                (cond ((and zero-row (not (pzerop v)) (setq zero-row nil))))
 204                (setq agcd (cond (agcd (pgcd agcd v))
 205                                 (t v)))
 206                finally (return
 207                          (cond ((numberp agcd) lis)
 208                                (zero-row lis)
 209                                (t
 210
 211                                 (setq $gcds_used ($append $gcds_used `((mlist) ,(header-poly agcd))))
 212                                 (format t "~%Found a non trivial factor:") (sh agcd)
 213                                 (cons '(mlist)
 214                                       (loop for vv in (cdr lis)
 215                                             collect (header-poly (pquotient (st-rat vv) agcd)))))))))))
 216
 217
 218
 219 (defun $cancel_pivot(lis piv)
 220  (setq piv (st-rat piv))
 221            (assert (affine-polynomialp piv))
 222   (cond ((mbagp lis)
 223          (cons (car lis) (loop for w in (cdr lis) collect($cancel_pivot  w piv))))
 224         (t
 225            (setq lis (st-rat lis))
 226            (assert (affine-polynomialp lis))
 227            (let ((tem (pgcd lis piv)))
 228              (cond ((numberp tem) (header-poly lis))
 229                    (t (header-poly (pquotient lis tem))))))))
 230
 231 (defun $linearize_nc (x)
 232   (cond ((consp x)
 233          (cond ((and (consp (car x)) (eq (caar x) 'mnctimes))
 234                 (setf (car x) '(mtimes))
 235                 (loop for v on (cdr x)
 236                    for term in (cdr x)
 237                    for i from 1
 238                    do (setf (car v) (intern (format nil "~a~d" term i)))))
 239                (t (cons ($linearize_nc (car x)) ($linearize_nc (cdr x))))))
 240         (t x)))
 241
 242 (defun $linearize_nc_to_nc (x)
 243   (cond ((consp x)
 244          (cond ((and (consp (car x)) (eq (caar x) 'mnctimes))
 245                 (loop for v on (cdr x)
 246                    for term in (cdr x)
 247                    for i from 1
 248                    do (setf (car v) (intern (format nil "~a~d" term i)))))
 249                (t (cons ($linearize_nc_to_nc (car x)) ($linearize_nc_to_nc (cdr x)))))
 250          x)
 251         (t x)))
 252
 253 (defun times_4_n (n)
 254   (if (zerop n)
 255       1
 256       (* 4 n)))
 257
 258 (defun hilbert_4(n)
 259   (round ($binomial (+ n 3) 3)))
 260
 261 (push 'hilbert_4 *all-rank-functions*)
 262
 263 ;(loop for i below 10 collect (list i (hilbert_tem i)))
 264
 265 (defun $minors2_2 (mat cols)
 266   (let ((nmat (apply '$submatrix mat (cdr cols))))
 267      (cons '(mlist)
 268            (loop for ro in (list-tableaux 2 (1- (length nmat)))
 269             do (show ro)
 270             collect (cons '($matrix) (loop for i in ro
 271                                         collect (nth i nmat)))))))
 272
 273 (defun $minors (n mat cols)
 274   (let ((nmat (apply '$submatrix mat (cdr cols))))
 275      (cons '(mlist)
 276            (loop for ro in (list-tableaux n (1- (length nmat)))
 277             do (show ro)
 278             collect (cons '($matrix) (loop for i in ro
 279                                         collect (nth i nmat)))))))
 280
 281 (loop for i below 5 collect i)
 282
 283 ;(setq he '(1 2))
 284
 285
 286 (defun sublist (l pred)
 287   (loop for v in l when (apply pred v) collect v))
 288
 289 (defun lis (l)
 290   (cons '(mlist) (loop for w in l collect (cons '(mlist) w))))
 291
 292 ;(setq $lis1 (lis (sublist (cartesian-product he he he he) #'(lambda ( i j l m) (and (< i j) (<= m l))))))
 293
 294 ;(setq $lis2 (lis                (sublist (cartesian-product he he he he) #'(lambda ( i j l m) (and (<= i j) (< m l))))))
 295
 296 ;; POLYS is a list of polynomials and a parallel list DEGS of degrees.
 297 ;; We compute all monomials in POLYS of deg l
 298 (defun monoms (polys degs &optional (cross '(nil)))
 299   (cond ((null polys) '(1))
 300         (t
 301           (loop for v in (monoms (cdr polys) (cdr degs))
 302                  append
 303                  (loop for j to (car degs)
 304                         collect (n* (pexpt (car polys) j) v))))))
 305
 306 (defun $extract_c_equations(poly vars-to-exclude)
 307   (let ((monoms (st-rat vars-to-exclude))
 308         (vars1  (mapcar 'car (st-rat vars-to-exclude))))
 309     (assert ($listp poly))
 310   (let (ans all-monoms some-monoms)
 311     (setq ans
 312     (loop for v in (cdr poly)
 313            do (setq v (function-numerator (st-rat v)))
 314            (setq degs (loop for w in vars1
 315                              collect (pdegree v w)))
 316            (setq some-monoms (monoms monoms degs))
 317            (setq all-monoms (append all-monoms some-monoms))
 318            nconc
 319            (loop for w in some-monoms
 320                   collect (new-disrep (pcoeff v w vars1)))))
 321     (list '(mlist)
 322       (cons '(mlist) ans)
 323       (cons '(mlist) (mapcar 'new-disrep all-monoms))))))
 324
 325 (defun $pb ()
 326   (save-linenumbers :file "/tmp/lines"))