?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Description of maxima_unicode_display in reference manual,
[maxima.git] / share / affine / sparsemat.lisp
blob187c37774f276a1a458c3d09fcf8ec8b8e3f46b4
1 ;;; -*- mode: lisp; package: cl-maxima; syntax: common-lisp ;base: 10; -*-
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; ;;;;;
4 ;;; Copyright (c) 1984 by William Schelter,University of Texas ;;;;;
5 ;;; All rights reserved ;;;;;
6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
7
8 (in-package :maxima)
9
10 (defstruct (sparse-matrix :named (:conc-name sp-))
11 rows
12 number-of-rows
13 type-of-entries
14 current-row
15 current-row-number
16 current-row-length
17 pivot-row
18 pivot-row-number
19 characteristic
20 inverse-array
21 special-inverse
22 special-plus
23 (minimum-size-to-grow 1)
24 (allow-reorder t)
25 last-good-row
26 columns-used-to-pivot ;;key is the column ==>row number
27 column-used-in-row ;;(aref column-used-in-row row-num)
28 ;;==>column for row row-num
29 list-of-all-columns-occurring
30 pivot-entry
31 (pivot-test-list '(unit min))
32 row-number-before-swap
33 special-multiply
34 transpose
35 (reduced nil)
36 sort-pivot
37 (sign-of-row-permutation 1)
38 current-column-above-pivot-row-number
39 columns-with-no-pivot ;;list of them
40 rows-with-no-pivot ;;has a 1 in slot i --> no pivot in row i
41 solutions
42 constants-column-number
43 special-solution ;; a sparse-matrix row which dots with rows
44 (constants-column nil)) ;to give the constants-column.
45
46 (defstruct (solution-row :named (:conc-name solution-row-))
47 number
48 data)
49
50
51 (defmacro row-length (arow)
52 `(array-dimension-n 1 ,arow))
53
54 ;(defmacro set-keywords (listt key-list &aux body)
55 ; (let ((pack package))
56 ; (loop for v in key-list
57 ; do(show pack)
58 ; (show v)
59 ; collecting `(, v (setq ,(intern-local v pack) (car (setq key (cdr key)))))
60 ; into tem
61 ; finally (setq body
62 ; (append tem (list '(otherwise
63 ; (merror "unrecognized key word"))))))
64 ; `(do ((key (cdr ,listt) (cdr key)))
65 ; ((null key))
66 ; (case (car key)
67 ; ,@ body))))
68 ;
69 ;(defun foo ()(set-keywords '(:hi 4 :bye 5) (:hiii :hi :bye :extra :byye :there)))
70 ;(DO ((KEY (CDR '(:HI 4)) (CDR KEY)))
71 ; ((NULL KEY))
72 ; (CASE (CAR KEY)
73 ; (:HI (SETQ HI (CAR (SETQ KEY (CDR KEY)))))
74 ; (:BYYE (SETQ :BYYE (CAR (SETQ KEY (CDR KEY)))))
75 ; (:THERE (SETQ THERE (CAR (SETQ KEY (CDR KEY)))))
76 ; (OTHERWISE (FERROR "unrecognized key word"))))
77
78 (deff sp-rational-quotient #'/)
79
80 (defmacro without-double-evaluation (list-of-forms &body body &aux repl)
81 (setq repl (loop for v in list-of-forms
82 when (atom v) nconc nil
83 else collecting (list (gensym) v)))
84 (cond ((> (length body) 1)
85 (setq body `(progn ,@body)))
86 (t (setq body (car body))))
87 (cond (repl
88 `(let ,repl
89 ,(sublis (loop for u in repl collecting (cons (second u) (car u)))
90 body :test 'equal)))
91 (t `,body)))
92
93 ;;this is really like vector-push-extend except it may speed up later references.
94 (defmacro array-push-extend-replace (aarray data &key amount-to-grow replace &aux forms)
95 (declare (ignore replace))
96 `(without-double-evaluation (,data)
97 (cond ((vector-push ,data ,aarray))
98 (t (let ((.new. (adjust-array ,aarray (+ ,(or amount-to-grow 100) (array-total-size ,aarray))
99 :fill-pointer (fill-pointer ,aarray))))
100 (vector-push ,data .new.),@ forms)))))
101
102 ;(defun sp-add* (&rest llist)
103 ; (let ((varlist))
104 ; (apply 'add* llist)))
105 ;
106 ;(defun sp-minus* (a)
107 ; (let ((varlist))
108 ; (simplifya (list '(mminus) a) nil)))
109 ;
110 ;(defun sp-sub* (a b )
111 ; (sp-add* a (sp-minus* b)))
112 ;
113 ;(defun sp-mul* (&rest a-list)
114 ; (let ((varlist ))
115 ; (cond ((< (length a-list) 3)(apply 'mul* a-list))
116 ; (t (sp-mul* (car a-list) (apply 'sp-mul* (cdr a-list)))))))
117 ;
118 ;(defun sp-div* (a b)
119 ; (let ((varlist))
120 ; (simplifya (list '(mquotient) a b) nil)))
121 ;;see polyc
122
123 (defmacro sp-number-of-rows (self)
124 `(fill-pointer (sp-rows ,self)))
125
126 ;(defun make-one-dimensional (aarray &aux ans)
127 ; (cond ((eql (array-#-dims aarray) 1) aarray)
128 ; (t (setq ans (make-array (apply '* (array-dimensions aarray)) :fill-pointer 0 :adjustable t))
129 ; (loop for i below (row-length aarray)
130 ; when (aref aarray i 0)
131 ; do
132 ; (array-push-extend-replace ans (aref aarray i 0))
133 ; (array-push-extend-replace ans (aref aarray i 1)))
134 ; ans)))
135
136 (defmacro with-characteristic (&body body &aux body1 body2 body3)
137 (setq body1 (sublis '((special-times . *) (special-plus . +)) body))
138 (setq body2 (sublis '((special-times . finite-characteristic-times)
139 (special-plus . finite-characteristic-plus)) body))
140 (setq body3 (sublis '((special-times . sp-mul*)
141 (special-plus . sp-add*)) body))
142 `(cond ((equal (sp-type-of-entries sp-mat) ':any-macsyma) ,@ body3)
143 ((zerop (sp-characteristic sp-mat)) ,@ body1)
144 (t (let ((.characteristic. (sp-characteristic sp-mat)))
145 ,@ body2))))
146
147 (defmacro with-once-only (variables &body body1 &aux ll var reset)
148 "This macro is to save re-evaluation of instance variables again and
149 again. It unwind protects the resetting of variables in case of throws etc.
150 Variables may be generalized"
151 ; (setq body1 (cdr body))
152 ; (setq variables (car body))
153 (loop for u in variables
154 when (appears-in body1 u)
155 do
156 (setq ll (cons (setq var (list (gensym) u)) ll))
157 (setq body1 (subst (car var) (cadr var) body1))
158 (setq reset (append reset (list `(setf ,@(reverse var))))))
159 (cond (ll `(let (,@ ll)
160 (unwind-protect (progn ,@ body1) ,@ reset)))
161 (t `(progn ,@ body1))))
162
163 ;;(with-once-only ((sp-rows sp-mat) hi)
164 ;; (setq hi 4) (+ 1 2) 3)
165 ;;((LAMBDA (#:G2972)
166 ;; (UNWIND-PROTECT (PROGN (SETQ #:G2972 4)
167 ;; (+ 1 2)
168 ;; 3)
169 ;; (SETF HI #:G2972)))
170 ;; HI)
171
172 (defmacro finite-characteristic-times (&rest l)
173 `(mod (* ,@l) .characteristic.))
174
175 (defmacro finite-characteristic-plus (&rest l)
176 `(mod (+ ,@l) .characteristic.))
177
178 (defmacro special-minus (a) `(special-times -1 ,a))
179
180 (defmacro special-times (&rest l)
181 `(cond ((equal (sp-type-of-entries sp-mat) ':any-macsyma)(sp-mul* ,@ l))
182 ((zerop (sp-characteristic sp-mat))(* ,@ l))
183 (t (mod (* ,@ l) (sp-characteristic sp-mat)))))
184
185 (defmacro special-plus (&rest l)
186 `(cond ((eq (sp-type-of-entries sp-mat) ':any-macsyma) (sp-add* ,@ l))
187 ((zerop (sp-characteristic sp-mat)) (+ ,@ l))
188 (t (mod (+ ,@ l) (sp-characteristic sp-mat)))))
189
190 (defmacro non-negative-remainder (x y)
191 `(progn (setf ,x (rem ,x ,y))
192 (cond ((minusp ,x) (+ ,x ,y))
193 (t ,x))))
194 ;
195 ;(defsubst row-entry ( arow j )
196 ; (let ((.this-row. arow)(.j. j) ind)
197 ; (catch 'entry
198 ; (loop for ii below (array-active-length .this-row.) by 2
199 ; when (and (setf ind (aref .this-row. ii )) (eq ind .j.))
200 ; do (throw 'entry (aref .this-row. (1+ ii )))))))
201
202 (defun row-entry ( arow j &aux col)
203 (catch 'entry
204 (loop for ii below (length (the cl:array arow)) by 2
205 when (and (setf col (aref arow ii )) (eql j col))
206 do (throw 'entry (aref arow (1+ ii ))))))
207
208 ;(defmacro row-dot (.this-row. brow ) "dot product of rows arow and brow"
209 ; `(let ((.val.)
210 ; (.char. characteristic) .prod.)
211 ; (loop for .iii. below (row-length ,arow)
212 ; when (and (aref ,arow .iii. 0) (setf .val. (row-entry ,brow (aref ,arow .iii. 0))))
213 ; summing (special-times
214 ; .val.
215 ; (aref ,arow .iii. 1)))))
216
217 ;(defun row-dot (arow brow ) "dot product of rows arow and brow"
218 ; (let ((val)
219 ; prod)
220 ; (loop for iii below (row-length arow)
221 ; when (and (aref arow iii 0) (setq val (row-entry brow (aref arow iii 0))))
222 ; summing (special-times
223 ; val
224 ; (aref arow iii 1)))))
225
226
227 ;(defun make-number (n)
228 ; (cond ((null n) 0)
229 ; (t n)))
230
231 (defun sp-row-dot (sp-mat arow brow &aux aindex (ans 0))
232 (with-characteristic
233 (let (val)
234 (loop for iii below (length (the cl:array arow)) by 2
235 when (and (setf aindex (aref arow iii ))
236 (setq val (row-entry brow aindex)))
237 do (setq ans
238 (special-plus ans
239 (special-times
240 val
241 (aref arow (1+ iii) )))))))
242 ans)
243
244
245 (defun sp-row-dot-row-numbers (sp-mat i j)
246 (sp-row-dot sp-mat (aref (sp-rows sp-mat) i) (aref (sp-rows sp-mat) j)))
247
248 (defmacro below-fill (a)
249 `(max (- (length (the cl:array ,a)) 2) 0))
250 ;(defsubst set-fill-pointer (aarray n)
251 ; (store-array-leader n aarray 0))
252 ;
253 ;(defsubst maybe-move-back-fill-pointer (arow)
254 ; (cond ((equal (array-active-length arow) 0) nil)
255 ; (t
256 ; (let ((this-row arow))
257 ; (loop for i downfrom (below-fill this-row) to 0 by 2
258 ; when (aref this-row i)
259 ; do (setf (fill-pointer this-row) (+ i 2)) (return 'done)
260 ; finally (setf (fill-pointer this-row) 0))))))
261
262 (defun maybe-move-back-fill-pointer (arow)
263 (cond ((equal (length (the cl:array arow)) 0) nil)
264 (t
265 (loop for i downfrom (below-fill arow) to 0 by 2
266 when (aref arow i)
267 do (setf (fill-pointer arow) (+ i 2)) (return 'done)
268 finally (setf (fill-pointer arow) 0)))))
269
270 (defun sp-reduce-elements-for-type (sp-mat )
271 (cond ((fixnump (sp-type-of-entries sp-mat))
272 (setf (sp-characteristic sp-mat) (sp-type-of-entries sp-mat))
273 (loop for iii below (sp-number-of-rows sp-mat)
274 do (let ((this-row (aref (sp-rows sp-mat) iii)))
275 (loop for ii below (length (the cl:array this-row)) by 2
276 when (aref this-row ii )
277 do (let ((value (aref this-row (1+ ii))))
278 (cond (value
279 (setq value (mod value (sp-type-of-entries sp-mat)))
280 (cond ((eql value 0)(setf (aref this-row ii ) nil))
281 (t (setf (aref this-row (1+ ii)) value)))))))
282 ;(setf;(aref this-row ii 1) value)))))))
283 (maybe-move-back-fill-pointer this-row))))))
284
285 (defun sp-reset-list-of-all-columns-occurring (sp-mat &aux ind)
286 (setf (sp-list-of-all-columns-occurring sp-mat)
287 (loop for ii below (sp-number-of-rows sp-mat)
288 for this-row = (aref (sp-rows sp-mat) ii)
289 appending
290 (loop for jj below (length (the cl:array this-row)) by 2
291 when (and (setq ind (aref this-row jj ))
292 (not (member ind temp :test #'eq)))
293 collecting ind)
294 into temp
295 finally (return temp))))
296
297 (defun sp-set-rows (sp-mat the-rows &optional list-of-columns)
298 (setf (sp-rows sp-mat)
299 (if (array-has-fill-pointer-p the-rows)
300 the-rows
301 (make-array (array-dimension the-rows 0) :adjustable t
302 :displaced-to the-rows :fill-pointer t)))
303 (setf (sp-number-of-rows sp-mat) (fill-pointer (sp-rows sp-mat)))
304 ; (cond ((zerop (fill-pointer the-rows)) (break 'why-no-rows)))
305 ; (show (sp-rows sp-mat))
306 (setf (sp-last-good-row sp-mat) (1- (sp-number-of-rows sp-mat)))
307 (cond ((null list-of-columns) ( sp-reset-list-of-all-columns-occurring sp-mat))
308 (t (setf (sp-list-of-all-columns-occurring sp-mat) list-of-columns)))
309 (setf (sp-column-used-in-row sp-mat)
310 (MAKE-ARRAY (fill-pointer (sp-rows sp-mat)) :adjustable t))
311 (setf (sp-rows-with-no-pivot sp-mat)
312 (make-ARRAY (fill-pointer (sp-rows sp-mat)) :element-type '(mod 2)))
313 (sp-check-set-up sp-mat)
314 (cond ((and (sp-columns-used-to-pivot sp-mat)
315 (hash-table-p (sp-columns-used-to-pivot sp-mat) ))
316 (clrhash (sp-columns-used-to-pivot sp-mat) ))
317 (t
318 (setf (sp-columns-used-to-pivot sp-mat)
319 (make-hash-table :size (sp-number-of-rows sp-mat))))))
320
321 (defun ml-sort (lis &optional (predicate #'alphalessp))
322 (sort lis predicate))
323
324 (defun sp-check-set-up (sp-mat &aux tem)
325 (setf (sp-current-row sp-mat) nil)
326 (setf (sp-reduced sp-mat ) nil)
327 (setf (sp-last-good-row sp-mat) (sp-number-of-rows sp-mat))
328 (setf ( sp-current-column-above-pivot-row-number sp-mat) nil)
329 (cond ((setq tem (sp-columns-used-to-pivot sp-mat))
330 (clrhash tem ))
331 (t (setf (sp-columns-used-to-pivot sp-mat) (make-hash-table))))
332 (setf (sp-columns-with-no-pivot sp-mat) nil)
333 ; (setf (sp-column-used-in-row sp-mat) nil)
334 )
335
336
337 ;(defun compare-structures (a b &aux slotb slota)
338 ; (let ((slots (fourth (get (ml-typep a) 'si:defstruct-description))))
339 ; (loop for v in slots
340 ; for i from 0
341 ; when (not (equalp (setq slota (funcall (car (last v)) a))
342 ; (setq slotb (funcall (car (last v)) b))))
343 ; do (show (car (last v)) slota slotb))))
344
345 (defun sp-set-current-row (sp-mat i)
346 (setf (sp-current-row-number sp-mat) i)
347 (setf (sp-current-row sp-mat) (aref (sp-rows sp-mat) i))
348 (setf (sp-current-row-length sp-mat)
349 (length (the cl:array (sp-current-row sp-mat))))
350 (sp-current-row sp-mat))
351
352 (eval-when
353 (:load-toplevel :compile-toplevel :execute)
354 (defmacro sp-row (sp-mat i)
355 `(aref (sp-rows ,sp-mat) ,i)))
356
357 (defun sp-set-pivot-row (sp-mat i)
358 (setf (sp-pivot-row-number sp-mat) i)
359 (setf (sp-pivot-row sp-mat) (aref (sp-rows sp-mat) i))
360
361 (sp-pivot-row sp-mat))
362
363 (defun sp-entry (sp-mat i j &aux (this-row (aref (sp-rows sp-mat) i)))
364 (catch 'entry
365 (loop for ii below (length (the cl:array this-row)) by 2
366 when (eql j (aref this-row ii ))
367 do (throw 'entry (aref this-row (1+ ii))))))
368
369 (defun fix-even (n)
370 (* 2 (round (quotient n 2.0))))
371
372 (defun sp-grow-current-row (sp-mat &optional (ratio 1.3))
373 (let ((new-length
374 (max (fix-even (* ratio (row-length
375 (aref (sp-rows sp-mat)
376 (sp-current-row-number sp-mat)) )))
377 (sp-minimum-size-to-grow sp-mat))))
378 (setf (aref
379 (sp-rows sp-mat) (sp-current-row-number sp-mat))
380 (array-grow (sp-current-row sp-mat) new-length ))
381 ( sp-set-current-row sp-mat (sp-current-row-number sp-mat))))
382
383 (defun sp-grow-row (sp-mat row-number &optional (ratio 1.3))
384 (let ((new-length (max (fix-even (* ratio (row-length (aref (sp-rows sp-mat) row-number) )))
385 (sp-minimum-size-to-grow sp-mat))))
386 (setf (aref
387 (sp-rows sp-mat) row-number)
388 (array-grow (aref (sp-rows sp-mat) row-number) new-length ))))
389
390 ;(defmacro current-rowset ( element index sslot)
391 ; `(progn (cond ((equal ,element 0) (aset nil (sp-current-row sp-mat) ,sslot )
392 ; (maybe-move-back-fill-pointer (sp-current-row sp-mat)))
393 ; (t (cond (< ,sslot
394 ; (aset ,element (sp-current-row sp-mat) (1+ ,sslot) )
395 ; (aset ,index (sp-current-row sp-mat) , sslot )))))
396
397 ;(defmacro this-row-set ( element index sslot)
398 ; `(progn (cond ((equal ,element 0) (aset nil this-row ,sslot )
399 ; (maybe-move-back-fill-pointer this-row))
400 ; (t
401 ; (aset ,element this-row (1+ ,sslot) )
402 ; (aset ,index this-row , sslot )))))
403 ;
404 ;(defsubst set-entry ( value arow index)
405 ; (catch 'entry-is-set
406 ; (let* ((.this-row. arow)
407 ; (.ind. index) j (.val. value)
408 ; first-empty-slot
409 ; (active-length (array-active-length .this-row.)))
410 ; (setf first-empty-slot
411 ; (loop for ii below (array-active-length .this-row.) by 2
412 ; do
413 ; (cond ((null (setq j (aref .this-row. ii))) (return ii))
414 ; ((eql j .ind.)
415 ; (cond (($zerop .val.) (aset nil .this-row. .ind.)
416 ; (if (eql ii (- active-length 2))
417 ; (maybe-move-back-fill-pointer .this-row.))
418 ; )
419 ; (t
420 ; (aset .val. .this-row. (1+ ii ))
421 ; (aset .ind. .this-row. ii)))
422 ; (throw 'entry-is-set t)))))
423 ; (cond (first-empty-slot
424 ; (loop for ii from first-empty-slot below active-length by 2
425 ; when (eql .ind. (aref .this-row. ii))
426 ; do
427 ; (cond (($zerop .val.) (aset nil .this-row. .ind.)
428 ; (if (eql ii (- active-length 2))
429 ; (maybe-move-back-fill-pointer .this-row.))
430 ; )
431 ; (t
432 ; (aset .val. .this-row. (1+ ii ))
433 ; (aset .ind. .this-row. ii)))
434 ; (throw 'entry-is-set t)
435 ; finally
436 ; (cond (($zerop .val.) nil)
437 ; (t (aset .val. .this-row. (1+ first-empty-slot))
438 ; (aset .ind. .this-row. first-empty-slot)))))
439 ; ((not ($zerop .val.))
440 ; (array-push-extend .this-row. .ind.)
441 ; (array-push-extend .this-row. .val.))))))
442
443 (defun set-entry (value arow index &aux (aarow arow) ;; so can declare it an array-register
444 j first-empty-slot (active-length (fill-pointer arow)))
445 (catch 'entry-is-set
446 (setf first-empty-slot
447 (loop for ii below (fill-pointer aarow) by 2
448 do
449 (cond ((null (setq j (aref aarow ii))) (return ii))
450 ((eql j index)
451 (cond (($zerop value)(setf (aref aarow index) nil)
452 (if (eql ii (- active-length 2))
453 (maybe-move-back-fill-pointer aarow))
454 )
455 (t
456 (setf (aref aarow (1+ ii )) value)
457 (setf (aref aarow ii) index)))
458 (throw 'entry-is-set t)))))
459 (cond (first-empty-slot
460 (loop for ii from first-empty-slot below active-length by 2
461 when (eql index (aref aarow ii))
462 do
463 (cond (($zerop value)(setf (aref aarow index) nil)
464 (if (eql ii (- active-length 2))
465 (maybe-move-back-fill-pointer aarow))
466 )
467 (t
468 (setf (aref aarow (1+ ii )) value)
469 (setf (aref aarow ii) index)))
470 (throw 'entry-is-set t)
471 finally
472 (cond (($zerop value) nil)
473 (t(setf (aref aarow (1+ first-empty-slot)) value)
474 (setf (aref aarow first-empty-slot) index)))))
475 ((not ($zerop value))
476 (vector-push-extend index aarow)
477 (vector-push-extend value aarow)))))
478
479 ;;; did not seem to be calling this
480 ;(defsubst set-entry-without-moving-back-fill-pointer ( value arow index)
481 ;
482 ;(defun sp-rem-current-row-entry (sp-mat index)
483 ; (aset nil (sp-current-row sp-mat) index 0))
484
485 (defun sp-set-current-row-entry (sp-mat value ind)
486 (set-entry value (sp-current-row sp-mat) ind))
487
488 (defun sp-set-entry (sp-mat value i j)
489 (set-entry value (aref (sp-rows sp-mat) i) j))
490
491 ;(defun row-set-entry (value this-row j)
492 ; (let ((first-empty-slot nil))
493 ; (catch 'finished
494 ; (setq first-empty-slot
495 ; (catch 'first
496 ; (loop for ii below (row-length this-row)
497 ; do (cond ((eql (aref this-row ii 0) j) (this-row-set value j ii)
498 ; (throw 'finished t))
499 ; ((null (aref this-row ii 0)) (throw 'first ii))
500 ; (t nil)))))
501 ; (cond ( first-empty-slot
502 ; (loop for ii from first-empty-slot below (row-length this-row)
503 ; do (cond ((eql (aref this-row ii 0) j) (this-row-set value j ii)
504 ; (throw 'finished t))))
505 ; (this-row-set value j first-empty-slot))
506 ; (t
507 ; (let ((last-spot (row-length this-row)))
508 ; (array-grow this-row (list (round (* 1.3 (row-length this-row))) 2))
509 ; (this-row-set value j last-spot)))))))
510
511 ;(defmacro pivot-row-ref (index)
512 ; `(catch 'pivot-row-ref
513 ; (loop for .ii. below (array-active-length pivot-row) by 2
514 ; do (cond ((eql ,index (aref pivot-row .ii.))
515 ; (throw 'pivot-row-ref (aref pivot-row (1+ .ii. ))))))))
516
517 ;(defmacro current-row-ref (index)
518 ; `(catch 'current-row-ref
519 ; (loop for .ii. below (array-active-length (sp-current-row sp-mat)) by 2
520 ; do (cond ((eql ,index (aref (sp-current-row sp-mat) .ii. ))
521 ; (throw 'current-row-ref (aref (sp-current-row sp-mat) (1+ .ii.) )))))))
522
523 ;(defmacro current-row-slot (index)
524 ; `(catch 'current-row-slot
525 ; (loop for .ii. below (array-active-length (sp-current-row sp-mat)) by 2
526 ; do (cond ((eql ,index (aref (sp-current-row sp-mat) .ii. ))
527 ; (throw 'current-row-slot .ii.))))))
528
529
530
531 ;(defmacro matrix-entry ( i j) "This finds the entry of a sparse-matrix
532 ; the ith row and jth column"
533 ; `(catch 'matrix-entry
534 ; (let ((.row. (aref (sp-rows sp-mat) ,i)))
535 ; (loop for .ii. below (row-length .row.)
536 ; do (cond ((eql ,j (aref .row. .ii. 0))
537 ; (throw 'matrix-entry (aref .row. .ii. 1))))))))
538
539 ;(defmacro set-matrix-entry (value i j) "This sets the entry of a sparse-matrix
540 ; the ith row and jth column"
541 ; `(catch 'set-matrix-entry
542 ; (setq first-empty-slot (catch 'where
543 ; (let ((.row. (aref (sp-rows sp-mat) ,i)))
544 ; (loop for .ii. below (row-length .row.)
545 ; do (cond ((null (aref .row. .ii. 0) (throw 'where .ii.))))))))))
546 ; ((eql ,j (aref .row. .ii. 0))
547 ; (aset value .row. .ii. 1)
548 ; (throw 'set-matrix-entry t)))))) nil)
549
550 (defmacro row-slot (row index)
551 "Returns the slot that the INDEX (column) appears in ROW"
552 `(catch 'row-slot
553 (let ((.row. ,row))
554 (loop for .ii. below (length (the cl:array .row.)) by 2
555 do (cond ((eql ,index (aref .row. .ii. ))
556 (throw 'row-slot .ii.)))))))
557
558 (defmacro set-current-row-entry-in-new-column (value row index)
559 "Puts new entry in a column not occurring in row. Would work
560 for (sp-rows sp-mat) other than current-row except for wanting the row number
561 to be able to replace it in ROWS if ROW is grown."
562 `(let ((.val. ,value)
563 (.ind. ,index))
564 (cond (($zerop .val.) nil)
565 (t
566 (loop for ii below (length (the cl:array ,row)) by 2
567 when (null (aref ,row ii))
568 do(setf (aref ,row ii) .ind.)
569 (setf (aref ,row (1+ ii)) .val.)
570 (return 'entry-is-set)
571 finally (array-push-extend-replace ,row .ind. :replace ((aref (sp-rows sp-mat) (sp-current-row-number sp-mat))))
572 (array-push-extend-replace ,row .val. :replace ((aref (sp-rows sp-mat) (sp-current-row-number sp-mat)))))))))
573
574 ;;Note: the with-once-only makes local variables for the pivot-row
575 ;;and the current-row and it is these that are declared
576 (defun sp-row-operation (sp-mat factor &aux current-row-slot new-value
577 piv-col temp)
578 "this replaces the current-row by current-row + factor pivot-row"
579 (with-characteristic
580 (with-once-only
581 ((sp-pivot-row sp-mat) (sp-current-row sp-mat))
582 (cond
583 ((equal factor 0) nil)
584 (t
585 (loop for ii below (length (the cl:array (sp-pivot-row sp-mat))) by 2
586 when (setq piv-col (aref (sp-pivot-row sp-mat) ii ))
587 do (cond ((setq current-row-slot
588 (row-slot (sp-current-row sp-mat) piv-col))
589 (setq new-value
590 (special-plus
591 (special-times
592 factor (aref (sp-pivot-row sp-mat) (1+ ii )))
593 (aref (sp-current-row sp-mat)
594 (1+ current-row-slot))))
595 (cond (($zerop new-value )
596 (setf (aref (sp-current-row sp-mat) current-row-slot ) nil))
597 (t(setf (aref
598 (sp-current-row sp-mat)
599 (1+ current-row-slot )) new-value))))
600 (t (set-current-row-entry-in-new-column
601 (special-times factor
602 (aref (sp-pivot-row sp-mat) (1+ ii)))
603 (sp-current-row sp-mat) piv-col))))
604 (maybe-move-back-fill-pointer (sp-current-row sp-mat)))))
605 (cond ((and (sp-constants-column sp-mat)
606 (not ($zerop (setq temp (aref (sp-constants-column sp-mat)
607 (sp-pivot-row-number sp-mat))))))
608 (setf (aref
609 (sp-constants-column sp-mat) (sp-current-row-number sp-mat)) (special-plus (aref (sp-constants-column sp-mat)
610 (sp-current-row-number sp-mat))
611 (special-times factor temp)))))))
612
613 ;(with-characteristic
614 ; (special-times 3 4))
615
616 (defun get-rows-from-array ( matrix)
617 (let* ((dims (array-dimensions matrix))
618 (the-rows (MAKE-ARRAY (car dims) :adjustable t :fill-pointer (car dims)))
619 arow)
620 (loop for i below (car dims)
621 do (setq arow (MAKE-ARRAY (* (second dims) 2) :fill-pointer 0 :adjustable t))
622 (setf (aref the-rows i) arow)
623 (loop for j below (second dims)
624 do (let ((entry (aref matrix i j)))
625 (cond ((not (equal entry 0))
626 (vector-push j arow)
627 (vector-push entry arow)))))
628 finally (return the-rows))))
629
630 (defun random-matrix (m n &optional (random-size 5))
631 (let ((mat (MAKE-ARRAY (list m n) :adjustable t)) )
632 (loop for i below m
633 do (loop for j below n
634 do(setf (aref mat i j) (random random-size)))
635 finally (return mat))))
636
637 (defun sp-gcd-row (sp-mat i)
638 (let ((ans ( sp-first-element-of-row sp-mat i))
639 (this-row (aref (sp-rows sp-mat) i)))
640 (loop for ii below (length (the cl:array this-row)) by 2
641 when (aref this-row ii )
642 do (setq ans (gcd ans (aref this-row (1+ ii ))))
643 finally (return ans))))
644 (defun sp-first-element-of-row (sp-mat ii)
645 (catch 'first-element
646 (let ((this-row (aref (sp-rows sp-mat) ii)))
647 (loop for i below (length (the cl:array this-row)) by 2
648 when (aref this-row i )
649 do (throw 'first-element (aref this-row (1+ i )))))))
650
651
652 (defmacro show-row (arow)
653 `(let ((this-row ,arow) ind)
654 (loop for jj below (length (the cl:array this-row)) by 2
655 when (setq ind (aref this-row jj))
656 do (format t "~%In slot ~D ~D-->~D "
657 (truncate jj 2) ind (aref this-row (1+ jj))))))
658
659 (defun sp-show-row (sp-mat i)
660 (let ((this-row (aref (sp-rows sp-mat) i)))
661 (format t "~%Row ~D is ~A." i this-row)
662 (show-row this-row)))
663
664 (defun sp-show-current-and-pivot (sp-mat )
665 (let ((crn (sp-current-row-number sp-mat)))
666 (format t "~%Current row is ~D which is ~A." crn (sp-current-row sp-mat))
667 (show-row (sp-current-row sp-mat))
668 (format t "~%Pivot Row is ~D which is ~A." (sp-pivot-row-number sp-mat)
669 (sp-pivot-row sp-mat))
670 (show-row (sp-pivot-row sp-mat))))
671
672 ;(defmacro special-inverse (x)
673 ; `(cond ((equal (sp-type-of-entries sp-mat)
674 ; ':any-macsyma)(sp-div* 1 ,x))
675 ; ((numberp (sp-type-of-entries sp-mat))
676 ; (aref (sp-inverse-array sp-mat)
677 ; (mod ,x (sp-type-of-entries sp-mat))))
678 ; (t (sp-rational-quotient 1 ,x))))
679
680 (defmacro special-inverse (x)
681 `(cond ((equal (sp-type-of-entries sp-mat)
682 ':any-macsyma)(sp-div* 1 ,x))
683 ((numberp (sp-type-of-entries sp-mat))
684 (let ((modulus (sp-type-of-entries sp-mat)))
685 (crecip ,x)))
686 (t (sp-rational-quotient 1 ,x))))
687
688 (defun sp-rat (x)
689 (cond ((or (affine-polynomialp x)(rational-functionp x)) x)
690 (($bfloatp x) x)
691 (t (new-rat x))))
692
693
694 (defun sp-choose-type-of-entries (sp-mat &aux rational float tem1
695 (rows (sp-rows sp-mat)))
696 (loop named sue for i below (length (the cl:array rows))
697 do
698 (let ((a-row (aref rows i)))
699 (loop for ii below (length (the cl:array a-row)) by 2
700 when (aref a-row ii)
701 do
702 (cond ((numberp (setq tem1 (aref a-row (1+ ii))))
703 (cond ((or float (floatp tem1))(setq float t))
704 ((rationalp tem1) (setq rational t))))
705 (t (sp-set-type-of-entries sp-mat ':any-macsyma)
706 (return-from sue 'done)))))
707 finally
708 ;;if any floating use float else if any rational use rational
709 ;;else use integer
710 (cond (float (sp-set-type-of-entries sp-mat ':float))
711 (rational (sp-set-type-of-entries sp-mat :rational))
712 (t (sp-set-type-of-entries sp-mat :integer)))))
713
714 (defun sp-set-type-of-entries (sp-mat y)
715 (setf (sp-type-of-entries sp-mat) y)
716 (cond ((equal y ':integer)
717 (setf (sp-characteristic sp-mat) (or modulus 0))
718 (setf (sp-pivot-test-list sp-mat) '(unit gcd-column )))
719 ;;any is bad do gcd of col etc.
720 ((equal y ':rational)
721 (cond (modulus (setf (sp-characteristic sp-mat) modulus)
722 (setf (sp-type-of-entries sp-mat)modulus))
723 (t
724 (setf (sp-pivot-test-list sp-mat) '(unit any)) ;;removed gcd
725 (setf (sp-characteristic sp-mat) 0))))
726 ((equal y ':float)
727 (setf (sp-pivot-test-list sp-mat) '(any))
728 (setf (sp-characteristic sp-mat) (or modulus 0)))
729 ((equal y ':any-macsyma)
730 (setf (sp-pivot-test-list sp-mat) '(:any-macsyma))
731 (setf (sp-characteristic sp-mat) (or modulus 0))
732 (loop for i below (array-total-size (sp-rows sp-mat))
733 do
734 (let ((this-row (aref (sp-rows sp-mat) i)))
735 (loop for j below (array-total-size this-row) by 2
736 when (aref this-row j)
737 do(setf (aref
738 this-row (1+ j)) (sp-rat (aref this-row (1+ j))))))))
739 ((and (fixnump y) (> y 0))
740 (setf (sp-characteristic sp-mat) y))
741 (t (merror "At present only positive integers (ie finite char)
742 ,integer, rational, or':any-macsyma are valid types"))))
743
744
745 (defun sp-enter-pivot-data (sp-mat col)
746 (setf (gethash
747 col
748 (sp-columns-used-to-pivot sp-mat))
749 (sp-pivot-row-number sp-mat))
750 (setf (aref
751 (sp-column-used-in-row sp-mat) (sp-pivot-row-number sp-mat)) col)
752 (fmat t "Row ~D has pivot ~D in column ~D." (sp-pivot-row-number sp-mat)
753 (sp-pivot-entry sp-mat)
754 col))
755
756 (defun sp-remove-zero-entries-from-row (sp-mat i)
757 (let ((this-row (aref (sp-rows sp-mat) i)))
758 (loop for ii below (length (the cl:array this-row)) by 2
759 when (aref this-row ii )
760 do (if (null (aref this-row (1+ ii) ))(setf (aref this-row ii 0) nil)))
761 (maybe-move-back-fill-pointer this-row)))
762
763 (defun show-hash (tabl)
764 (loop for x being the hash-keys in tabl
765 using (hash-value y)
766 do (format t "~%~A -->~A" x y)))
767
768 (defun sp-show-rows(sp-mat &rest l)
769 (cond ((null l) (setq l (loop for i below (sp-number-of-rows sp-mat)
770 collecting i))))
771 (loop for i in l do ( sp-show-row sp-mat i))
772 (format t "~%Current row is row ~D." (sp-current-row-number sp-mat))
773 (format t "~%Pivot row is row ~D." (sp-pivot-row-number sp-mat)))
774
775 (defun sp-list-pivots (sp-mat &optional pred &aux entry col)
776 (cond ((null pred) (setq pred (function (lambda (ign)ign t))))
777 ((equal pred :non-unit)
778 (setq pred (function (lambda (x) (not (eql (abs x) 1)))))))
779 (loop for i below (sp-number-of-rows sp-mat)
780 when (and (setq col (aref (sp-column-used-in-row sp-mat) i))
781 (funcall pred
782 (setq entry ( sp-entry sp-mat i col))))
783 collecting entry into tem
784 ; (format t "~%Row ~D has pivot ~D in column ~D." i entry col)
785 finally (return tem)))
786
787 (defun sp-show-pivots (sp-mat &optional pred &aux entry col)
788 (cond ((null pred) (setq pred (function (lambda (ignor)ignor t))))
789 ((equal pred :non-unit)
790 (setq pred (function (lambda (x) (not (eql (abs x) 1)))))))
791 (loop for i below (sp-number-of-rows sp-mat)
792 when (and (setq col (aref (sp-column-used-in-row sp-mat) i))
793 (funcall pred (setq entry ( sp-entry sp-mat i col))))
794 do
795 (format t "~%Row ~D has pivot ~D in column ~D." i entry col)))
796
797 (defmacro null-to-zero (x)
798 `(cond ((null ,x) 0)
799 (t ,x)))
800
801 (defun sp-show-matrix(sp-mat )
802 (cond ( (sp-list-of-all-columns-occurring sp-mat) nil)
803 (t
804 (setf (sp-list-of-all-columns-occurring sp-mat)
805 (sort (sp-list-of-all-columns-occurring sp-mat) '<))))
806 (format t "~% " )
807 (loop for i in (sp-list-of-all-columns-occurring sp-mat) do
808 (format t "C~D " i))
809 (cond ((sp-constants-column sp-mat)
810 (format t "Const~D" (sp-constants-column-number sp-mat))))
811 (loop for i below (sp-number-of-rows sp-mat)
812 do (format t "~%R~2D " i)
813 (loop for j in (sp-list-of-all-columns-occurring sp-mat)
814 do (format t "~5D"
815 (null-to-zero ( sp-entry sp-mat i j))))
816 (cond ((sp-constants-column sp-mat)
817 (format t "~5D" (aref (sp-constants-column sp-mat) i)))))
818 (cond ((and (sp-special-solution sp-mat) )
819 (format t "~%SpSol")
820 (loop for j in (sp-list-of-all-columns-occurring sp-mat)
821 do (format t "~5D"
822 (null-to-zero (row-entry
823 (sp-special-solution sp-mat) j)))))))
824
825 ;(defun sp-reduce (sp-mat &optional type-of-entry &aux current-row-entry)
826 ; (cond (type-of-entry ( sp-set-type-of-entries sp-mat type-of-entry)))
827 ; (cond ((typep (sp-type-of-entries sp-mat) :fixnum)
828 ; ( sp-reduce-elements-for-type sp-mat)))
829 ; ( sp-clear-pivots sp-mat)
830 ; (with-characteristic
831 ; (loop for i below (sp-number-of-rows sp-mat)
832 ; do
833 ; ( sp-set-pivot-row sp-mat i)
834 ; ( sp-best-pivot sp-mat)
835 ; (if (sp-pivot-entry sp-mat)
836 ; (let ((pivoting-column
837 ; (aref (sp-column-used-in-row sp-mat)
838 ; (sp-pivot-row-number sp-mat)))
839 ; (minus-inverse-pivot-entry
840 ; (special-minus (special-inverse (sp-pivot-entry sp-mat)))))
841 ; (loop for ii from (+ 1 i) below (sp-number-of-rows sp-mat)
842 ; do
843 ; ( sp-set-current-row sp-mat ii)
844 ; (cond ((setq current-row-entry
845 ; ( sp-entry sp-mat ii
846 ; pivoting-column ))
847 ; ( sp-row-operation sp-mat
848 ; (special-times
849 ; current-row-entry
850 ; minus-inverse-pivot-entry)
851 ; )))
852 ;;;this gets rid of the denoms but you should gcd and.. what about constants
853 ; (cond ((and (eql (sp-type-of-entries sp-mat) :integer)
854 ; (not (typep minus-inverse-pivot-entry :fixnum)))
855 ;
856 ; (sp-multiply-row sp-mat ii (sp-pivot-entry sp-mat))))
857 ;;; (sp-show-matrix sp-mat)
858 ;
859 ; )))))
860 ;
861 ; (Setf (sp-reduced sp-mat) t))
862
863 (defun sp-reduce (sp-mat &optional type-of-entry &aux current-row-entry
864 (original-number-of-rows (sp-number-of-rows sp-mat)))
865 (cond (type-of-entry ( sp-set-type-of-entries sp-mat type-of-entry)))
866 (cond ((fixnump (sp-type-of-entries sp-mat))
867 ( sp-reduce-elements-for-type sp-mat)))
868 (sp-clear-pivots sp-mat)
869 (with-characteristic
870 (loop for i from 0 while (< i (sp-number-of-rows sp-mat))
871 do
872 ( sp-set-pivot-row sp-mat i)
873 ( sp-best-pivot sp-mat)
874 (if (sp-pivot-entry sp-mat)
875 (let ((pivoting-column
876 (aref (sp-column-used-in-row sp-mat)
877 (sp-pivot-row-number sp-mat)))
878 (minus-inverse-pivot-entry
879 (special-minus (special-inverse
880 (sp-pivot-entry sp-mat)))))
881 (loop for ii from (+ 1 i) below (sp-number-of-rows sp-mat)
882 do
883 ( sp-set-current-row sp-mat ii)
884 (cond ((setq current-row-entry
885 (sp-entry sp-mat ii
886 pivoting-column ))
887 (sp-row-operation sp-mat
888 (special-times
889 current-row-entry
890 minus-inverse-pivot-entry)
891 ))))
892 ))))
893 (sp-kill-extra-rows sp-mat original-number-of-rows)
894 (setf (sp-reduced sp-mat) t))
895
896 (defun sp-kill-extra-rows (sp-mat original-number-of-rows &aux up-bound)
897 (let ((rows (sp-rows sp-mat)))
898 (loop for i from original-number-of-rows below (sp-number-of-rows sp-mat)
899 when (not (zerop (fill-pointer (aref rows i)) ))
900 do (setq up-bound (1+ i))
901 finally(cond(up-bound nil) (t (setq up-bound original-number-of-rows)))
902 (setf (fill-pointer rows) up-bound)
903 (setf (sp-number-of-rows sp-mat) up-bound))))
904
905
906 (defun sp-multiply-row (sp-mat row-number factor)
907 (let ((this-row (aref (sp-rows sp-mat) row-number)))
908 (cond ((equal 0 (special-times 1 factor))
909 (loop for i below (fill-pointer this-row)
910 do (setf (aref this-row i ) nil)
911 (setf (fill-pointer this-row) 0)))
912 ((equal 1 factor) nil)
913 (t
914 (with-characteristic
915 (loop for i below (length (the cl:array this-row)) by 2
916 when (aref this-row i )
917 do(setf (aref this-row (1+ i) ) (special-times
918 (aref this-row (1+ i ))
919 factor ))))))))
920
921
922 (defmacro dsend (object &rest body)
923 `(progn (send ,object .,body)
924 (send ,object :show-matrix)))
925
926 (defun sp-clear-pivots(sp-mat )
927 (clrhash (sp-columns-used-to-pivot sp-mat))
928 (fillarray (sp-column-used-in-row sp-mat) nil)
929 (fillarray (sp-rows-with-no-pivot sp-mat) '(0))
930 (loop for i below (sp-number-of-rows sp-mat)
931 do(setf (aref (sp-column-used-in-row sp-mat) i) nil)))
932
933 (defun sp-get-rows-from-array (sp-mat aarray &optional (type :integer))
934 ( sp-set-rows sp-mat (get-rows-from-array aarray))
935 ( sp-set-type-of-entries sp-mat type)
936 (sp-check-set-up sp-mat))
937
938
939 ;(defun h (i j) ($concat '$dd i j))
940 ;(defun z
941 ; (loop for u in a-list do (apply 'aset 0 aarray u)))
942 ;(zero-out aa '((0 1) (0 3) (0 5) (0 7)(1 0)(1 2)(1 4) (1 6)(3 1) (3 3)(3 5)(3 7)
943 ; (4 0)(4 2)(4 4)(4 6)))
944
945 (defun genmatrix (m n fun)
946 (let ((ans (MAKE-ARRAY (list m n) :adjustable t)))
947 (loop for i below m
948 do (loop for j below n do
949 (setf (aref ans i j) (funcall fun i j ))))
950 ans))
951
952 (defun macsyma-array-to-row ( aarray &optional (grow-ratio 1))
953 (let* ((aarray (symbol-array (mget aarray 'hashar)))
954 (actual-size (aref aarray 1))
955 (index 0)
956 ans val )
957 (setq ans (MAKE-ARRAY (round (* actual-size 2 grow-ratio)) :fill-pointer 0 :adjustable t))
958 (loop for i from 3 below (car (array-dimensions aarray))
959 do (cond ((setq val (aref aarray i))
960 (loop for u in val do
961 (array-push-extend-replace ans (caar u))
962 (array-push-extend-replace ans (cdr u))))))
963 (cond ((not (equal (setq index (sp-rational-quotient (fill-pointer ans) 2) )
964 actual-size))
965 (format t "We only found ~D elements while
966 there should have been ~D elements."
967 index actual-size)))
968 ans))
969
970 (defun list-of-macsyma-arrays-to-rows (llist)
971 (let* ((size (length (cdr llist)))
972 (rows (MAKE-ARRAY size :fill-pointer 0 :adjustable t)))
973 (loop for u in (cdr llist)
974 do (vector-push-extend (funcall 'macsyma-array-to-row u) rows))
975 rows))
976
977 ;(defun get-array (&quote macsyma-array)
978 ; (fsymeval (mget macsyma-array 'hashar)))
979 ;
980 ;(defmacro mydump ( name-of-object-to-dump &optional filename file-atribute-list )
981 ; (cond ((null filename) (setq filename (string name-of-object-to-dump))))
982 ; `(sys:dump-forms-to-file ,filename
983 ; (list (list 'setq ',name-of-object-to-dump
984 ; (list 'quote ,name-of-object-to-dump)))
985 ; ,file-atribute-list))
986 ;
987 ;(defun te (n c &aux a aa)
988 ; (setq a c) ;;96 ms. this is best.93 ms. if put the aa let outside the do loop
989 ; (loop for i below n do (let(( aa (* i 2))) (cond ((zerop a) aa)
990 ; (t (rem aa a))))))
991 ;
992 ;(defun te (n c &aux d) ;;130ms.
993 ; (local-declare ((special d)) (setq d c)
994 ; (loop for i below n do ((lambda(a b)(cond ((eql d 0) (* a b))
995 ; (t (rem (* a b) d)))) i 2 ))))
996 ;
997 ;(defun te (n c &aux d)
998 ; (let ((d c)) ;;108 ms. slower than if you set up variable for (* i 2)
999 ; (loop for i below n do (cond ((zerop d) (* i 2))
1000 ; (t (rem (* i 2) d))))))
1001 ;(defun te (n c &aux d) ;;faster than without the setq d! approx 75ms.
1002 ; (loop for i below n do (setq d (aref c 500)) ))
1003 ;
1004 ;(defun te (n c &aux sp)
1005 ; (loop for i below n do (setq c ( i 2))))
1006 ;
1007 ;(defmacro with-speed (&body body)
1008 ; `(let ((.char. characteristic) .prod. .sum.)
1009 ; ,@body))
1010 ;
1011 ;
1012 ;(defun te (n c &aux )
1013 ; (let ((d c)) ;;120 ms.
1014 ; (loop for i below n do ((lambda(a b d)(cond ((eql d 0) (* a b))
1015 ; (t (rem (* a b) d)))) i 2 d ))))
1016 ;
1017 ;(progn 'compile (setf (function sp-ti) (function (lambda(a b d)(cond ((eql d 0) (* a b))
1018 ; (t (rem (* a b) d)))) )))
1019 ;(defun sp-test (sp-mat n)
1020 ; (let (((sp-type-of-entries sp-mat) type-of-entries)) ;;fairly-fast
1021 ; (loop for i below n do (hh i 2 type-of-entries))))
1022 ;
1023 ;(defun sp-test (sp-mat n &aux a)
1024 ; (let ((type-of-entries type-of-entries))
1025 ; (loop for i below n do
1026 ; (setq a (* i 2)) (cond ((zerop type-of-entries) a)
1027 ; (t (rem type-of-entries a))))))
1028 ;
1029 ;(defun sp-test (sp-mat n)
1030 ; (with-characteristic ; 77ms. in characteristic=0 this is best!!
1031 ; (loop for i below n do (special-times i 2))))
1032 ;
1033 ;
1034 ;(defun-method h sparse-matrix (a b) ;;slow
1035 ; (cond ((eql type-of-entries 0) (* a b))
1036 ; (t (rem (* a b) type-of-entries))))
1037
1038 (defun sp-number-of-pivots(sp-mat &optional (below-row (sp-number-of-rows sp-mat)))
1039 (let ((count 0))
1040 (loop for i below below-row
1041 when (aref (sp-column-used-in-row sp-mat) i) do (setf count (1+ count))
1042 finally (return count))))
1043
1044 ;(defun te (&rest l &aux with bar five)
1045 ; (keyword-extract l vv ((:with bar) five ))
1046 ; (print (list bar five)))
1047 ;(for element of row with index in slot do body)
1048 ;(defmacro for (u &rest l &aux element row index slot)
1049 ; (keyword-extract l vv ( (:of row) (:with index) (:in slot) (:do body)))
1050 ; `(let ((,row (aref rows ,i)) ,element ,index )
1051 ; (loop for ,slot below (row-length ,row)
1052 ; when (aref ,row ,slot 0)
1053 ; do (setq ,element (aref ,row ,slot 1))
1054 ; (setq ,index (aref ,row ,slot 0))
1055 ; . ,body)))
1056
1057 (defun sp-solve (sp-mat &key reduce reset-list-of-columns-occurring)
1058 (cond (reset-list-of-columns-occurring ( sp-reset-list-of-all-columns-occurring sp-mat)))
1059 (cond ((or reduce (null (sp-reduced sp-mat)))
1060 ( sp-clear-pivots sp-mat)
1061 ( sp-reduce sp-mat)))
1062 (setf (sp-columns-with-no-pivot sp-mat)
1063 (loop for u in (sp-list-of-all-columns-occurring sp-mat)
1064 when (null (gethash u (sp-columns-used-to-pivot sp-mat)))
1065 collecting u))
1066 (let* ((number-of-columns (length (sp-list-of-all-columns-occurring sp-mat)))
1067 (number-of-pivots ( sp-number-of-pivots sp-mat))
1068 (number-of-solutions (- number-of-columns number-of-pivots))
1069 a-solution a-new-entry a-special-solution
1070 (solution-rows
1071 (MAKE-ARRAY (length (sp-columns-with-no-pivot sp-mat))
1072 :fill-pointer 0 :adjustable t)))
1073 (cond ((not (equal (length (sp-columns-with-no-pivot sp-mat))
1074 number-of-solutions))
1075 (format t "~%There are ~D columns and ~D pivots and ~D columns with no pivot,
1076 something is wrong" (length (sp-list-of-all-columns-occurring sp-mat)) number-of-pivots (length (sp-columns-with-no-pivot sp-mat)))))
1077
1078 (cond
1079 ((> number-of-solutions 0)
1080 (with-characteristic
1081
1082 (loop for u in (sp-columns-with-no-pivot sp-mat)
1083 do (setf a-solution (MAKE-ARRAY (* (1+ number-of-pivots) 2)
1084 :fill-pointer 0 :adjustable t))
1085
1086 (vector-push a-solution solution-rows)
1087 (fmat t "~%solution-rows ~A" (listarray solution-rows))
1088 (set-entry 1 a-solution u)
1089 (loop for nn downfrom (1- (fill-pointer (sp-rows sp-mat))) to 0
1090 when (aref (sp-column-used-in-row sp-mat) nn)
1091 do
1092 (setf (sp-pivot-entry sp-mat)
1093 ( sp-entry sp-mat nn
1094 (aref (sp-column-used-in-row sp-mat) nn)))
1095 (setq a-new-entry
1096 (special-times
1097 -1
1098 (special-inverse (sp-pivot-entry sp-mat))
1099 ( sp-row-dot sp-mat a-solution (aref (sp-rows sp-mat) nn))))
1100
1101 (cond ((not ($zerop a-new-entry))
1102 (vector-push (aref (sp-column-used-in-row sp-mat) nn) a-solution)
1103 (vector-push a-new-entry a-solution))))))))
1104 (cond ((sp-constants-column sp-mat)
1105 (loop for i below (array-total-size (sp-column-used-in-row sp-mat))
1106 when (null (aref (sp-column-used-in-row sp-mat) i))
1107 do
1108 (cond ((not ($zerop (aref (sp-constants-column sp-mat) i)))
1109 (format t "~%~% ******************************~%")
1110 (format t "~%Row ~D has no pivot but the (sp-constants-column sp-mat) has entry ~D.
1111 ~%Warning: The equations were INCONSISTENT.
1112 ~%The special-solution is NOT VALID.
1113 ~% ******************************" i (aref (sp-constants-column sp-mat) i)))))
1114 (show (sp-type-of-entries sp-mat))
1115 (with-characteristic
1116 (setq a-special-solution (MAKE-ARRAY (* (1+ number-of-pivots) 2)
1117 :fill-pointer 0 :adjustable t))
1118 (loop for nn downfrom (1- (fill-pointer (sp-rows sp-mat))) to 0
1119 when (aref (sp-column-used-in-row sp-mat) nn)
1120 do
1121 (setf (sp-pivot-entry sp-mat)
1122 ( sp-entry sp-mat nn (aref
1123 (sp-column-used-in-row sp-mat) nn)))
1124 (setq a-new-entry
1125 (special-times -1
1126 (special-inverse (sp-pivot-entry sp-mat))
1127 (special-plus
1128 ( sp-row-dot sp-mat a-special-solution
1129 (aref (sp-rows sp-mat) nn))
1130 (aref (sp-constants-column sp-mat) nn))))
1131 (cond ((not ($zerop a-new-entry))
1132 (vector-push
1133 (aref (sp-column-used-in-row sp-mat) nn) a-special-solution)
1134 (vector-push a-new-entry a-special-solution)))))))
1135 (cond ((and (sp-solutions sp-mat) (eq (type-of (sp-solutions sp-mat)) 'sparse-matrix)))
1136 (t (setf (sp-solutions sp-mat) (make-sparse-matrix ))))
1137 (sp-set-rows (sp-solutions sp-mat) solution-rows)
1138 (sp-set-type-of-entries (sp-solutions sp-mat)
1139 (sp-type-of-entries sp-mat))
1140 (setf (sp-special-solution (sp-solutions sp-mat)) a-special-solution)))
1141
1142 ;
1143 ;(defun sp-solve (sp-mat &rest l &aux reset-list-of-columns-flag reduce-flag)
1144 ; (keyword-extract l vv nil ((:reduce reduce-flag)
1145 ; (:reset-list-of-columns-occurring reset-list-of-columns-flag)))
1146 ; (cond (reset-list-of-columns-flag ( sp-reset-list-of-all-columns-occurring sp-mat)))
1147 ; (cond ((or reduce-flag (null (sp-reduced sp-mat)))
1148 ; ( sp-clear-pivots sp-mat)
1149 ; ( sp-reduce sp-mat)))
1150 ;
1151 ; (setf (sp-columns-with-no-pivot sp-mat)
1152 ; (loop for u in (sp-list-of-all-columns-occurring sp-mat)
1153 ; when (null (gethash u (sp-columns-used-to-pivot sp-mat))) collecting u))
1154 ; (let* ((number-of-columns (length (sp-list-of-all-columns-occurring sp-mat)))
1155 ; (number-of-pivots ( sp-number-of-pivots sp-mat))
1156 ; (number-of-solutions (- number-of-columns number-of-pivots))
1157 ; a-solution a-new-entry a-special-solution
1158 ; (solution-rows (make-array (length (sp-columns-with-no-pivot sp-mat)) :fill-pointer 0 :adjustable t)))
1159 ; (cond ((not (equal (length (sp-columns-with-no-pivot sp-mat)) number-of-solutions))
1160 ; (format t "~%There are ~D columns and ~D pivots and ~D columns with no pivot,
1161 ;something is wrong" (length (sp-list-of-all-columns-occurring sp-mat)) number-of-pivots (length (sp-columns-with-no-pivot sp-mat)))))
1162 ;
1163 ; (cond
1164 ; ((> number-of-solutions 0)
1165 ; (with-characteristic
1166 ; (loop for u in (sp-columns-with-no-pivot sp-mat)
1167 ; do (setf a-solution (make-array (* (1+ number-of-pivots) 2)
1168 ; :leader-list (list 0 u)))
1169 ;
1170 ; (array-push solution-rows a-solution)
1171 ; (fmat t "~%solution-rows ~A" (listarray solution-rows))
1172 ; (set-entry 1 a-solution u)
1173 ; (loop for nn downfrom (1- (row-length (sp-rows sp-mat))) downto 0
1174 ; when (aref (sp-column-used-in-row sp-mat) nn)
1175 ; do
1176 ; (setf (sp-pivot-entry sp-mat) ( sp-entry sp-mat nn (aref (sp-column-used-in-row sp-mat) nn)))
1177 ; (setq a-new-entry (special-times -1
1178 ; (special-inverse (sp-pivot-entry sp-mat))
1179 ; ( sp-row-dot sp-mat a-solution (aref (sp-rows sp-mat) nn))))
1180 ;; (+ ( sp-row-dot sp-mat a-solution (aref (sp-rows sp-mat) nn))
1181 ;; (cond ((null (sp-constants-column sp-mat)) 0)
1182 ;; (t (aref (sp-constants-column sp-mat) nn))))))
1183 ; (cond ((not ($zerop a-new-entry))
1184 ; (array-push a-solution (aref (sp-column-used-in-row sp-mat) nn))
1185 ; (array-push a-solution a-new-entry))))))))
1186 ; (cond ((sp-constants-column sp-mat)
1187 ; (loop for i below (array-total-size (sp-column-used-in-row sp-mat))
1188 ; when (null (aref (sp-column-used-in-row sp-mat) i))
1189 ; do
1190 ; (cond ((not ($zerop (aref (sp-constants-column sp-mat) i)))
1191 ; (format t "~%~% ******************************~%")
1192 ; (format t "~%Row ~D has no pivot but the (sp-constants-column sp-mat) has entry ~D.
1193 ; ~%WARNING. The equations were INCONSISTENT.
1194 ; ~%The special-solution is NOT VALID.
1195 ; ~% ******************************" i (aref (sp-constants-column sp-mat) i)))))
1196 ; (with-characteristic
1197 ; (setq a-special-solution (make-array (* (1+ number-of-pivots) 2)
1198 ; :leader-list (list 0 )))
1199 ; (loop for nn downfrom (1- (row-length (sp-rows sp-mat))) downto 0
1200 ; when (aref (sp-column-used-in-row sp-mat) nn)
1201 ; do
1202 ; (setf (sp-pivot-entry sp-mat) ( sp-entry sp-mat nn (aref (sp-column-used-in-row sp-mat) nn)))
1203 ; (setq a-new-entry
1204 ; (special-times -1
1205 ; (special-inverse (sp-pivot-entry sp-mat))
1206 ; (special-plus ( sp-row-dot sp-mat a-special-solution (aref (sp-rows sp-mat) nn))
1207 ; (aref (sp-constants-column sp-mat) nn))))
1208 ; (cond ((not ($zerop a-new-entry))
1209 ; (array-push a-special-solution (aref (sp-column-used-in-row sp-mat) nn))
1210 ; (array-push a-special-solution a-new-entry)))))))
1211 ; (cond ((and (sp-solutions sp-mat)(typep (sp-solutions sp-mat) 'sparse-matrix))
1212 ;;;why do they not work
1213 ;; ( sp-set-rows sp-mat solution-rows )
1214 ;; ( sp-set-type-of-entries sp-mat (sp-type-of-entries sp-mat)))
1215 ; (sp-set-rows (sp-solutions sp-mat) solution-rows)
1216 ; (sp-set-type-of-entries (sp-solutions sp-mat) (sp-type-of-entries sp-mat))
1217 ;; (setf (sp-rows (sp-solutions sp-mat))solution-rows)
1218 ;; (setf (sp-type-of-entries (sp-solutions sp-mat) ) (sp-type-of-entries sp-mat))
1219 ; )
1220 ; (t (setf (sp-solutions sp-mat) (make-sparse-matrix rows solution-rows
1221 ; type-of-entries (sp-type-of-entries sp-mat)))))
1222 ; (setf (sp-special-solution (sp-solutions sp-mat)) a-special-solution)))
1223
1224 (defun bring-to-left-side (possible-eqn)
1225 "Converts a=b into a-b, or a into a also doing nothing if b is 0 "
1226 (cond ((atom possible-eqn) possible-eqn)
1227 ((eq (caar possible-eqn) 'mequal)
1228 (cond ((eql (third possible-eqn) 0)(second possible-eqn))
1229 (t (sub (second possible-eqn) (third possible-eqn)))))
1230 (t possible-eqn)))
1231
1232
1233
1234 (defun sp-show-solutions(sp-mat )
1235 ( sp-show-matrix sp-mat))
1236
1237 (defun sp-verify-solutions (sp-mat &aux tem temm spec)
1238 (loop for i below (sp-number-of-rows sp-mat)
1239 do
1240 (loop for j below (sp-number-of-rows (sp-solutions sp-mat))
1241 do
1242 (format t "~%The dot product of row ~D and solution ~D is ~D"
1243 i j ( sp-row-dot sp-mat
1244 ( sp-row sp-mat i)
1245 (sp-row (sp-solutions sp-mat) j)))))
1246 (cond ((setq spec (sp-special-solution (sp-solutions sp-mat)))
1247 (let ((const (sp-constants-column sp-mat)))
1248 (loop for i below (sp-number-of-rows sp-mat)
1249 do
1250 (setq tem ( sp-row-dot sp-mat
1251 ( sp-row sp-mat i)
1252 spec))
1253 (format t "~%The dot product of the special solution and row ~D is ~A while the constant is ~A"
1254 i tem (aref const i))
1255 (iassert (or ($zerop(setq temm (special-plus (aref const i) tem )))
1256 (and (numberp temm)
1257 (< (abs temm) .00001))))
1258 finally (format t "~%The special solution is ok"))))))
1259 (defun nil-lessp (x y)
1260 (cond ((null y) (cond ((null x) nil)
1261 (t t)))
1262 ((null x) nil)
1263 (t (< x y))))
1264
1265
1266
1267 (defun sp-sort-row (sp-mat i)
1268 (sort-grouped-array (aref (sp-rows sp-mat) i) 2 'nil-lessp))
1269
1270
1271
1272 (defmacro may-have-pivot (row-number)
1273 `(not (eql 1 (aref (sp-rows-with-no-pivot sp-mat) ,row-number))))
1274
1275 (defun desirable-row (x) x nil) ;
1276 ;;?? fix array-leader stuff..
1277 ;(defun desirable-row (x) (eq (array-leader-length x) 2))
1278
1279 (defun sp-best-pivot (sp-mat &aux value temp-col)
1280 (setf (sp-pivot-entry sp-mat) nil)
1281 (setf (sp-row-number-before-swap sp-mat) nil)
1282 (fmat t "~%What is best pivot in Row ~D? " (sp-pivot-row-number sp-mat))
1283 (cond
1284 ((eql 1 (aref (sp-rows-with-no-pivot sp-mat) (sp-pivot-row-number sp-mat)))
1285 (fmat t "Row ~D was previously found to have no possible pivot."
1286 (sp-pivot-row-number sp-mat)))
1287 (t
1288 (catch 'done
1289 (loop
1290 for test-name in (sp-pivot-test-list sp-mat)
1291 when (may-have-pivot (sp-pivot-row-number sp-mat))
1292 do
1293 (setf (sp-last-good-row sp-mat) (1- (sp-number-of-rows sp-mat)))
1294 (let ((fresh-row t)
1295 best no-gcd test col-to-pivot)
1296 (loop
1297 while fresh-row
1298 do
1299 (cond
1300 ((and (desirable-row (sp-pivot-row sp-mat))
1301 (equal (sp-type-of-entries sp-mat) :integer))
1302 (setf test-name 'gcd-column)))
1303 (case test-name
1304 (any (setq test (function (lambda (ignor ignor1)ignor
1305 ignor1 t)))
1306 (setq best 1))
1307 (:any-macsyma ( sp-best-pivot-macsyma sp-mat )
1308 (throw 'done 'here-maybe-pivot))
1309 (gcd (multiple-value (best no-gcd)
1310 (catch-error ( sp-gcd-row sp-mat
1311 (sp-pivot-row-number sp-mat)) nil))
1312 (setq test (function equal)))
1313 (min (setq best ( sp-smallest-possible-pivot sp-mat))
1314 (setq test (function equal)))
1315 (gcd-column (setq best ( sp-smallest-possible-pivot sp-mat))
1316 (setq test (function equal)))
1317 (unit (cond ((numberp (sp-type-of-entries sp-mat)) (setq best 1)
1318 (setq test (function (lambda (ignor ignor1)ignor
1319 ignor1 t))))
1320 (t (setq best 1) (setq test (function equal)))))
1321 (otherwise (merror "~A is not a possible test." test-name)))
1322 (setq value nil)
1323 (cond
1324 (no-gcd nil)
1325 ((or (null best) ($zerop best))
1326 (fmat t "~%Row ~D was found to be zero while looking for ~A."
1327 (sp-pivot-row-number sp-mat)
1328 test-name)
1329 ( sp-put-back-row-with-no-pivot sp-mat))
1330 ((and (eq test-name 'gcd) no-gcd) nil)
1331 (t
1332 (with-once-only
1333 ((sp-pivot-row sp-mat))
1334 (loop
1335 for ii below (length (the cl:array (sp-pivot-row sp-mat))) by 2
1336 when (aref (sp-pivot-row sp-mat) ii)
1337 ; (and (aref (sp-pivot-row sp-mat) ii)
1338 ; (null (gethash ;;these entries should be zero!!
1339 ; (aref (sp-pivot-row sp-mat) ii )
1340 ; (sp-columns-used-to-pivot sp-mat))))
1341 do (setq value (abs (aref (sp-pivot-row sp-mat) (1+ ii ))))
1342
1343 (cond ((funcall test value best)
1344 (setq col-to-pivot (aref (sp-pivot-row sp-mat) ii))
1345 (cond ((equal test-name 'gcd-column)
1346 (cond
1347 ((and
1348 (desirable-row (sp-pivot-row sp-mat))
1349 (setq temp-col
1350 (
1351 sp-find-good-column-to-pivot sp-mat)))
1352 (setq col-to-pivot temp-col)))
1353
1354 ( sp-force-pivot-row-to-contain-gcd sp-mat
1355 col-to-pivot)
1356 (setf (sp-pivot-entry sp-mat)
1357 (sp-entry sp-mat (sp-pivot-row-number sp-mat)
1358 col-to-pivot))
1359 (setf (sp-pivot-row sp-mat)
1360 (aref (sp-rows sp-mat)
1361 (sp-pivot-row-number sp-mat))))
1362 (t (setf (sp-pivot-entry sp-mat)
1363 (aref (sp-pivot-row sp-mat) (1+ ii) ))))
1364 ( sp-enter-pivot-data sp-mat col-to-pivot )
1365 (throw 'done 'here-is-pivot)))))))
1366 (cond ( no-gcd nil)
1367 ((null value)
1368 ( sp-put-back-row-with-no-pivot sp-mat)
1369 (fmat t "It is zero.")))
1370 ; (fmat t "~%Row ~D has no possible pivot." (sp-pivot-row-number sp-mat))
1371 ; (aset 1 (sp-rows-with-no-pivot sp-mat) (sp-pivot-row-number sp-mat))
1372 ; (throw 'done 'no-pivot)))
1373 ; (cond ((member test-name '(min any gcd-column) :test #'eq)
1374 ; (fmat t "There is no possible pivot in row ~D." (sp-pivot-row-number sp-mat))
1375 ; (aset 1 (sp-rows-with-no-pivot sp-mat) (sp-pivot-row-number sp-mat))))
1376 (cond (no-gcd (format t "~%Row ~D has no gcd." (sp-pivot-row-number sp-mat)))
1377 ((eq test-name 'gcd)
1378 (fmat t "~%Row ~D has gcd ~D but no pivot of that size."
1379 (sp-pivot-row-number sp-mat) best ))
1380 ((eq test-name 'unit)
1381 (fmat t "~%Row ~D does not have a unit for a pivot."
1382 (sp-pivot-row-number sp-mat))))
1383 (setf fresh-row ( sp-swap-pivot-row-with-later-one sp-mat)))))))))
1384
1385 (defmacro swap-rows-and-constants (m n)
1386 `(prog (tem) (rotatef (aref (sp-rows sp-mat) ,m) (aref (sp-rows sp-mat) ,n))
1387 (setf (sp-sign-of-row-permutation sp-mat) (- (sp-sign-of-row-permutation sp-mat)))
1388 (cond ((setq tem (sp-current-column-above-pivot-row-number sp-mat))
1389 (rotatef (aref tem ,m) (aref tem ,n))))
1390 ;(cond (-boundp sign)(setq sign (- sign))) ;;keep track of sign for det
1391 (cond ((sp-constants-column sp-mat)
1392 (rotatef (aref (sp-constants-column sp-mat) ,m) (aref (sp-constants-column sp-mat) ,n))))))
1393
1394 (defun sp-put-back-row-with-no-pivot(sp-mat )
1395 (cond ((sp-row-number-before-swap sp-mat)
1396 (swap-rows-and-constants (sp-pivot-row-number sp-mat) (sp-row-number-before-swap sp-mat))
1397 ;(rotatef (aref (sp-rows sp-mat) (sp-pivot-row-number sp-mat)) (aref (sp-rows sp-mat) (sp-row-number-before-swap sp-mat) ))
1398 (fmat t "There is no pivot in row ~D so exchanging it back"
1399 (sp-row-number-before-swap sp-mat) )
1400 (setf (sp-pivot-row sp-mat) (aref (sp-rows sp-mat) (sp-pivot-row-number sp-mat)))
1401 (setf (aref (sp-rows-with-no-pivot sp-mat) (sp-row-number-before-swap sp-mat) ) 1)
1402 (setf (sp-row-number-before-swap sp-mat) nil))))
1403
1404 (defun sp-swap-pivot-row-with-later-one(sp-mat )
1405 "Returns t if it did nil if it did not"
1406 (catch 'exchanged
1407 (loop for j downfrom (sp-last-good-row sp-mat)
1408 until (<= j (sp-pivot-row-number sp-mat))
1409 when ($zerop (aref (sp-rows-with-no-pivot sp-mat) j))
1410 do (swap-rows-and-constants (sp-pivot-row-number sp-mat) j)
1411 ;(rotatef (aref (sp-rows sp-mat) (sp-pivot-row-number sp-mat)) (aref (sp-rows sp-mat) j))
1412 (setf (sp-pivot-row sp-mat) (aref (sp-rows sp-mat) (sp-pivot-row-number sp-mat)))
1413
1414 (setf (sp-last-good-row sp-mat) (1- j) )
1415 (setf (sp-row-number-before-swap sp-mat) j)
1416 (fmat t " Exchanging Row ~D for Row ~D." (sp-pivot-row-number sp-mat) j)
1417 (throw 'exchanged t))
1418 nil))
1419 ;
1420 ;(defun sp-rat (x)($vrat x))
1421 ;;the new generic functions n* n+ etc. do not require rat form entries, they convert to
1422 ;;them.
1423 ;
1424 ;(defun sp-best-pivot-macsyma(sp-mat )
1425 ; (loop for i below (array-active-length (sp-pivot-row sp-mat)) by 2
1426 ; when (aref (sp-pivot-row sp-mat) i)
1427 ; do
1428 ; (setf (sp-pivot-entry sp-mat) (aref (sp-pivot-row sp-mat) (1+ i)))
1429 ; (cond ((numberp (sp-pivot-entry sp-mat)) nil)
1430 ; (t (cond (($numberp (sp-pivot-entry sp-mat))
1431 ; (setf (sp-pivot-entry sp-mat)
1432 ; (cond ((atom (sp-pivot-entry sp-mat))(sp-pivot-entry sp-mat))
1433 ; ((or (affine-polynomialp (sp-pivot-entry sp-mat))
1434 ; (rational-functionp (sp-pivot-entry sp-mat)))
1435 ; (sp-pivot-entry sp-mat))
1436 ; (t
1437 ; ($ratsimp (sp-pivot-entry sp-mat))))))
1438 ; (t
1439 ; (setq (sp-pivot-entry sp-mat) (sp-rat (sp-pivot-entry sp-mat)))
1440 ; ))
1441 ; (aset (sp-pivot-entry sp-mat) (sp-pivot-row sp-mat) (1+ i) )))
1442 ; (send self :enter-pivot-data (aref (sp-pivot-row sp-mat) i))
1443 ; (return 'done)
1444 ; finally (cond (( sp-swap-pivot-row-with-later-one sp-mat)
1445 ; (send self :best-pivot-macsyma)))))
1446 ;
1447 ;(defun sp-best-pivot-macsyma(sp-mat )
1448 ; (cond ((loop for i below (array-active-length (sp-pivot-row sp-mat)) by 2
1449 ; when (and (aref (sp-pivot-row sp-mat) i)
1450 ; ($numberp
1451 ; (aref (sp-pivot-row sp-mat) (1+ i))))
1452 ; do
1453 ; (setq (sp-pivot-entry sp-mat) (aref (sp-pivot-row sp-mat) (1+ i)))
1454 ; (send self :enter-pivot-data (aref (sp-pivot-row sp-mat) i))
1455 ; (return 'done)))
1456 ; (t
1457 ; (loop for i below (array-active-length (sp-pivot-row sp-mat)) by 2
1458 ; when (aref (sp-pivot-row sp-mat) i)
1459 ; do
1460 ; (setq (sp-pivot-entry sp-mat) (aref (sp-pivot-row sp-mat) (1+ i)))
1461 ; (send self :enter-pivot-data (aref (sp-pivot-row sp-mat) i))
1462 ; (return 'done)
1463 ; finally (cond ((send self :swap-pivot-row-with-later-one)
1464 ; (send self :best-pivot-macsyma)))))))
1465
1466
1467 (defun sp-best-pivot-macsyma(sp-mat )
1468 (loop named sue for test in (list #'(lambda (x) (or (eql x -1)(eql x 1)))
1469 #'(lambda (x) (and (numberp x) (< (abs x) 100)))
1470 #'(lambda (x) (numberp x))
1471 #'(lambda (x) ($numberp x))
1472 #'(lambda (ignor) ignor t))
1473 do
1474 (loop for i below (length (the cl:array (sp-pivot-row sp-mat))) by 2
1475 when (and (aref (sp-pivot-row sp-mat) i)
1476 (funcall test
1477 (aref (sp-pivot-row sp-mat) (1+ i))))
1478 do
1479 (setf (sp-pivot-entry sp-mat) (aref (sp-pivot-row sp-mat) (1+ i)))
1480 ( sp-enter-pivot-data sp-mat (aref (sp-pivot-row sp-mat) i))
1481 (return-from sue 'done))
1482 finally (cond (( sp-swap-pivot-row-with-later-one sp-mat)
1483 ( sp-best-pivot-macsyma sp-mat)))))
1484
1485
1486
1487 (defun sp-smallest-possible-pivot(sp-mat &aux
1488 (.pivot-row. (sp-pivot-row sp-mat)))
1489 (loop for ii below (length (the cl:array .pivot-row.)) by 2
1490 when (aref .pivot-row. ii)
1491 minimize (abs (aref .pivot-row. (1+ ii )))))
1492 ;
1493 ;(defun test (m n &key mat
1494 ; (type :rational) solve
1495 ; constants constants-column sparse-matrix &aux const sp)
1496 ; (cond (sparse-matrix (setq sp sparse-matrix))
1497 ; (t (setq sp (make-sparse-matrix))))
1498 ; (cond (mat nil)
1499 ; (t
1500 ; (setq mat (random-matrix m n))))
1501 ; (cond (constants
1502 ; (setq
1503 ; constants-column
1504 ; (make-array m
1505 ; :fill-pointer m))
1506 ; (loop for i below m
1507 ; do (setf (aref constants-column i) (random 6)))))
1508 ;
1509 ; (sp-get-rows-from-array sp mat type)
1510 ; (setf (sp-constants-column sp) constants-column)
1511 ; (sp-show-matrix sp)
1512 ; (cond (solve
1513 ; (cond ((eq type :float) (sp-float sp)))
1514 ; (if constants (setq const (copy-atomic-structure constants-column 4)))
1515 ; (sp-solve sp :reduce t)
1516 ; (sp-get-rows-from-array sp mat type)
1517 ; (cond ((eq type :float) (sp-float sp)))
1518 ; (setf (sp-constants-column sp) const)
1519 ; (format t "~%Verifying the solutions with the original matrix and constants:")
1520 ; (sp-verify-solutions sp))
1521 ; (t
1522 ; (sp-reduce sp)))
1523 ; (sp-show-matrix sp)
1524 ;
1525 ; mat )
1526 (defun sp-verify-solutions-for-original-array (sp-mat aarray)
1527 ( sp-get-rows-from-array sp-mat aarray (sp-type-of-entries sp-mat))
1528 ( sp-verify-solutions sp-mat))
1529
1530
1531 (defun sp-init (sp-mat plist)
1532 (let ((initial-data (get plist ':rows))
1533 (atype-of-entries (get plist ':type-of-entries))
1534 (asolutions (get plist :solutions))
1535 (acolumn-used-in-row (get plist :column-used-in-row))
1536 (acolumns-used-to-pivot (get plist :columns-used-to-pivot)))
1537 (if initial-data ( sp-set-rows sp-mat initial-data))
1538 (if atype-of-entries ( sp-set-type-of-entries sp-mat atype-of-entries))
1539 (if asolutions (setf (sp-solutions sp-mat) asolutions))
1540 (if acolumns-used-to-pivot
1541 (setf (sp-columns-used-to-pivot sp-mat) acolumns-used-to-pivot))
1542 (if acolumn-used-in-row
1543 (setf (sp-column-used-in-row sp-mat) acolumn-used-in-row))))
1544
1545 ;(once (piv div ) (setq div 4)(setq piv 3));;watch out if control leaves in middle of body the
1546 ;values of the variables may not be restored!!
1547
1548 ;;The t in the instantiate-flavor tells it to send the :init message to
1549 ;;the newly created flavor. It uses the values present in the options-present
1550 ;;list. The catch-error is for those instance variables which may be unbound.
1551 ;;We do not want them on the options-present list.
1552
1553 ;(defun sp-fasd-form (sp-mat )
1554 ; (let ((options-present
1555 ; (loop for u in (list :solutions :column-used-in-row
1556 ; :columns-used-to-pivot :rows :type-of-entries)
1557 ; when (catch-error ( u) nil)
1558 ; appending (list u (send self u)))))
1559 ; (setf options-present (cons nil options-present))
1560 ; `(instantiate-flavor 'sparse-matrix ',options-present t)))
1561
1562 (defun sp-make-transpose (sp-mat )
1563 "Creates the transpose of the sparse-matrix. The result is stored in the transpose
1564 instance-variable. The original column is stored in slot 1 of the array-leader of the row"
1565 (let ((transpose-rows (MAKE-ARRAY
1566 (length (sp-list-of-all-columns-occurring sp-mat)) :adjustable t))
1567 rowj val)
1568 (setf (sp-list-of-all-columns-occurring sp-mat)
1569 (sort (sp-list-of-all-columns-occurring sp-mat) '<))
1570 (loop for j in (sp-list-of-all-columns-occurring sp-mat)
1571 for i below (length (sp-list-of-all-columns-occurring sp-mat))
1572 do
1573 (setq rowj (make-solution-row))
1574 (setf (solution-row-data rowj ) (MAKE-ARRAY 100 :fill-pointer 0 :adjustable t))
1575 (setf (solution-row-number rowj) j)
1576 (loop for ii below (sp-number-of-rows sp-mat)
1577 when (setq val ( sp-entry sp-mat ii j))
1578 do
1579 (array-push-extend-replace (solution-row-data rowj) ii)
1580 (array-push-extend-replace (solution-row-data rowj) val))
1581 (setf (aref transpose-rows i) (solution-row-data rowj)))
1582
1583 (setf (sp-transpose sp-mat) (make-sparse-matrix
1584 :rows transpose-rows
1585 :type-of-entries (sp-type-of-entries sp-mat)))))
1586 ;conflicts never called
1587 ;(defmacro make-sparse-matrix (aarray name)
1588 ; ` (progn (setf ,name (make-instance 'sparse-matrix))
1589 ; (send ,name :get-rows-from-array ,aarray)))
1590
1591
1592 ;; To handle integer matrices properly we will have to find the pivot
1593 ;;in a given column. To do this we will find the gcd of that column.
1594 ;;Then find a row where it occurs or create a linear combination of rows
1595 ;;where it occurs. WE add that to the set of rows. Then use that row to
1596 ;;clean the given column. The span over Z of the set of rows is the
1597 ;;same as before (since adding a row did not hurt and cleaning out with
1598 ;;respect to a given row is a reversible operation, since all entries in
1599 ;;that column are integer multiples of the given entry).
1600
1601 ;; We repeat the process applying it to the "set of rows occurring
1602 ;;after our row and to the remaining columns" Note that the previous
1603 ;;pivot column does not occur anyway. Eventually we end up with a set
1604 ;;of rows which span the same lattice as the original, but which are
1605 ;;clearly independent over the rationals. Therefore the number of them
1606 ;;must be correct and they must be a basis.
1607
1608 ;; Alternateley we could take our row and choose some column say j.
1609 ;;Then we could perform row' -q*row + remainder. We would then replace
1610 ;;row' by remainder. Here row' is another row that has an entry in
1611 ;;column j which is not a multiple of the the entry in row column j.
1612 ;;The remainder will have an entry smaller than that of row. Then we swap
1613 ;; the new row' and row. If row now has entry in column j
1614 ;;equal to the gcd of the column, fine we use the remainder as our
1615 ;;pivot. Otherwise we repeat with remainder taking the place of row and
1616 ;;we find some row' whose entry in column j is not a multiple of the
1617 ;;entry in row column j.
1618
1619 (defun sp-set-current-column-above-pivot-row-number (sp-mat j)
1620 ; (setq current-column-number j)
1621 (cond ((sp-current-column-above-pivot-row-number sp-mat)
1622 (setf (sp-current-column-above-pivot-row-number sp-mat)
1623 (adjust-array (sp-current-column-above-pivot-row-number sp-mat)
1624 (sp-number-of-rows sp-mat)
1625 :fill-pointer
1626 (fill-pointer
1627 (sp-current-column-above-pivot-row-number
1628 sp-mat))
1629 ))
1630 (fillarray (sp-current-column-above-pivot-row-number sp-mat) nil))
1631 (t (setf (sp-current-column-above-pivot-row-number sp-mat)
1632 (MAKE-ARRAY (sp-number-of-rows sp-mat) :adjustable t ))))
1633 (loop for i from (sp-pivot-row-number sp-mat) below (sp-number-of-rows sp-mat)
1634 do(setf (aref
1635 (sp-current-column-above-pivot-row-number sp-mat) i) ( sp-entry sp-mat i j))))
1636
1637 (defun sp-set-constants-column (sp-mat j &aux this-row)
1638 (setf (sp-constants-column-number sp-mat) j)
1639 (cond ((sp-constants-column sp-mat)
1640 (if (< (array-total-size (sp-constants-column sp-mat))
1641 (array-total-size (sp-rows sp-mat)))
1642 (setf (sp-constants-column sp-mat)
1643 (adjust-array (sp-constants-column sp-mat)
1644 (array-total-size (sp-rows sp-mat))
1645 :fill-pointer
1646 (fill-pointer (sp-constants-column sp-mat))
1647 ))))
1648 (t (setf (sp-constants-column sp-mat)
1649 (MAKE-ARRAY (array-total-size (sp-rows sp-mat))
1650 :fill-pointer 0 :adjustable t))))
1651 (loop for i below (array-total-size (sp-constants-column sp-mat))
1652 do (setf (aref (sp-constants-column sp-mat) i) 0))
1653 (setf (sp-list-of-all-columns-occurring sp-mat)
1654 (delete (sp-constants-column-number sp-mat)
1655 (sp-list-of-all-columns-occurring sp-mat) :test #'equal))
1656 (loop for i below (sp-number-of-rows sp-mat)
1657 do
1658 (setf this-row (aref (sp-rows sp-mat) i))
1659 (loop for ii below (length (the cl:array this-row)) by 2
1660 when (eql (aref this-row ii) j)
1661 do
1662 (setf (aref this-row ii) nil)
1663 (setf (aref (sp-constants-column sp-mat) i) (aref this-row (1+ ii))))))
1664
1665 ;(defun where-the-min (an-array &aux temp the-min where-the-min)
1666 ; "Where an-array has its minimum"
1667 ; (loop for i below (array-total-size an-array) do
1668 ; (cond ((setq the-min (aref an-array i)) (setq where-the-min i) (return 'done))))
1669 ; (cond ((null where-the-min) nil)
1670 ; (t (setq the-min (abs the-min))
1671 ; (loop for i below (array-total-size an-array)
1672 ; when (and (setq temp (aref an-array i)) (< (abs temp) the-min))
1673 ; do (setq where-the-min i))))
1674 ; where-the-min)
1675
1676 (defun where-the-min (an-array &aux null-vector tem the-min where-the-min)
1677 (loop for i below (array-total-size an-array)
1678 when (aref an-array i)
1679 do (setq the-min (abs (aref an-array i)))
1680 (setq where-the-min i)
1681 (return 'done)
1682 finally (setq null-vector t))
1683 (cond ((Not null-vector)
1684 (loop for i below (array-total-size an-array)
1685 when (setq tem (aref an-array i))
1686 do
1687 (cond ((< (abs tem ) the-min)(setq where-the-min i)
1688 (setq the-min (abs tem))))
1689 finally (return where-the-min)))
1690 (t nil)))
1691
1692 (defun sp-force-pivot-row-to-contain-gcd (sp-mat j &aux
1693 temp
1694 piv-entry
1695 tem constant crow
1696 gcdfacts new-row
1697 the-gcd where-the-min the-min-gcd where-the-min-gcd)
1698 "J is column number. Elementary row operations and swaps of rows
1699 are performed with rows occurring after (sp-pivot-row-number sp-mat) to ensure that
1700 the pivot-row contains the gcd of the entries in rows greater than
1701 the (sp-pivot-row-number sp-mat) and in column j."
1702 (declare (special *verbose*))
1703 ( sp-set-current-column-above-pivot-row-number sp-mat j)
1704 (let* ((colj (sp-current-column-above-pivot-row-number sp-mat))
1705 (first-entry (aref colj (sp-pivot-row-number sp-mat))))
1706
1707 (cond ((not (eql colj (sp-current-column-above-pivot-row-number sp-mat))) (merror 'wow)))
1708 (setq temp first-entry)
1709 (setq the-gcd temp)
1710 (loop for i below (row-length
1711 (sp-current-column-above-pivot-row-number sp-mat))
1712 when (setq temp (aref colj i))
1713 do (setq the-gcd (gcd temp the-gcd)))
1714
1715 (loop
1716 do
1717 (setq where-the-min (where-the-min colj))
1718
1719 (cond ((not (eql where-the-min (sp-pivot-row-number sp-mat)))
1720 (swap-rows-and-constants (sp-pivot-row-number sp-mat)
1721 where-the-min)
1722 ; (rotatef (aref (sp-rows sp-mat) (sp-pivot-row-number sp-mat)) (aref (sp-rows sp-mat) where-the-min))
1723 ;(rotatef (aref colj (sp-pivot-row-number sp-mat)) (aref colj where-the-min))
1724 (fmat t "~%Exchanging rows ~D and ~D to help force gcd ~D into pivot-row."
1725 (sp-pivot-row-number sp-mat) where-the-min the-gcd)
1726 ( sp-set-pivot-row sp-mat (sp-pivot-row-number sp-mat))))
1727 (cond
1728 ((equal the-gcd (abs (null-to-zero (aref colj (sp-pivot-row-number sp-mat) ))))
1729 (fmat t "~%Row ~D has ~D in column ~D which is the gcd of the column above row ~D."
1730 (sp-pivot-row-number sp-mat)
1731 (aref colj (sp-pivot-row-number sp-mat)) j (sp-pivot-row-number sp-mat))
1732 (return the-gcd))
1733 (t (setq piv-entry (aref colj (sp-pivot-row-number sp-mat)))
1734 (setq the-min-gcd (abs piv-entry))
1735
1736 (loop for i below (array-total-size colj)
1737 when (aref colj i)
1738 do
1739
1740 (cond((eql (setq tem (abs (gcd piv-entry (aref colj i)))) the-gcd)
1741 (setq where-the-min-gcd i)(return i)))
1742 (cond ( (< tem the-min-gcd)
1743 (setq the-min-gcd tem)
1744 (setq where-the-min-gcd i))))
1745
1746 (setq gcdfacts (general-gcd piv-entry (aref colj where-the-min-gcd)))
1747 ;;add new row (first gcdfacts)*pivot-row +(second gcdfacts)*(sp-row self (+ pivot-row-number where-the-min-gcd)
1748 ;;and fix constants. redo force pivot row;
1749 ; (show (listarray colj))
1750 ; (sp-show-matrix sp-mat)
1751 (sp-set-current-row sp-mat where-the-min-gcd)
1752 (setq crow (sp-current-row sp-mat))
1753 (setq new-row
1754 (MAKE-ARRAY (fill-pointer crow) :adjustable t :fill-pointer (fill-pointer crow)))
1755 (fillarray new-row crow)
1756 (cond ((setq tem (sp-constants-column sp-mat))
1757 (setq constant (aref tem (sp-current-row-number sp-mat)))))
1758
1759 (sp-add-new-row sp-mat new-row :constant constant)
1760 (sp-multiply-row sp-mat (sp-current-row-number sp-mat) (second gcdfacts))
1761 (sp-row-operation sp-mat (first gcdfacts) )
1762 (cond (*verbose* (format t "~%The gcd is ~A Adding a new row with entry in col ~A of ~A" the-gcd j
1763 (sp-entry sp-mat (sp-current-row-number sp-mat) j ) )))
1764 ; (sp-show-matrix sp-mat)
1765 ; (show (listarray colj))
1766 (sp-force-pivot-row-to-contain-gcd sp-mat j)(return 'start-over)))))
1767
1768 'done)
1769
1770 (defun sp-add-new-row (self a-row &key constant &aux ok)
1771 (let ((dim (1+ (sp-number-of-rows self))))
1772 (loop for i below (1- dim)
1773 do
1774 (maybe-move-back-fill-pointer (sp-row self i))
1775 when (eql (fill-pointer (sp-row self i)) 0)
1776 do
1777 (cond (constant (cond ((eql 0 (null-to-zero (aref (sp-constants-column self)
1778 i)))
1779 (setf (aref (sp-constants-column self) i) constant)
1780 (setq ok t))))
1781 (t (setq ok t)))
1782 (cond (ok
1783 (setf (aref (sp-rows-with-no-pivot self) i) 0)
1784 (setf (sp-row self i) a-row))))
1785 (cond ((not ok)
1786 (setf (sp-current-column-above-pivot-row-number self) nil)
1787 (setf (sp-rows-with-no-pivot self)
1788 (array-grow (sp-rows-with-no-pivot self) dim))
1789 (setf (sp-column-used-in-row self)
1790 (array-grow (sp-column-used-in-row self) dim))
1791 (vector-push-extend a-row (sp-rows self))
1792 (setf (sp-number-of-rows self) dim)
1793 (cond (constant (vector-push-extend constant (sp-constants-column self))))))))
1794
1795 ;
1796 ;(defun sp-force-pivot-row-to-contain-gcd (sp-mat j)
1797 ; "J is column number. Elementary row operations and swaps of rows
1798 ; are performed with rows occurring after (sp-pivot-row-number sp-mat) to ensure that
1799 ; the pivot-row contains the gcd of the entries in rows greater than
1800 ; the (sp-pivot-row-number sp-mat) and in column j."
1801 ; ( sp-set-current-column-above-pivot-row-number sp-mat j)
1802 ; (let* ((colj (sp-current-column-above-pivot-row-number sp-mat))
1803 ; (first-entry (aref colj (sp-pivot-row-number sp-mat))) temp
1804 ; rem factor where-the-smallest-remainder
1805 ; the-gcd where-the-min smallest-remainder )
1806 ; (setq temp first-entry)
1807 ; (setq the-gcd temp)
1808 ; (loop for i below (row-length (sp-current-column-above-pivot-row-number sp-mat))
1809 ; when (setq temp (aref colj i))
1810 ; do (setq the-gcd (gcd temp the-gcd)))
1811 ; (loop
1812 ; do
1813 ; (setq where-the-min (where-the-min colj))
1814 ; (cond ((not (eql where-the-min (sp-pivot-row-number sp-mat)))
1815 ; (swap-rows-and-constants (sp-pivot-row-number sp-mat) where-the-min)
1816 ; ; (rotatef (aref (sp-rows sp-mat) (sp-pivot-row-number sp-mat)) (aref (sp-rows sp-mat) where-the-min))
1817 ; ;(rotatef (aref colj (sp-pivot-row-number sp-mat)) (aref colj where-the-min))
1818 ; (fmat t "~%Exchanging rows ~D and ~D to help force gcd ~D into pivot-row."
1819 ; (sp-pivot-row-number sp-mat) where-the-min the-gcd)
1820 ; ( sp-set-pivot-row sp-mat (sp-pivot-row-number sp-mat))))
1821 ; (cond
1822 ; ((equal the-gcd (abs (aref colj (sp-pivot-row-number sp-mat) )))
1823 ; (fmat t "~%Row ~D has ~D in column ~D which is the gcd of the column above row ~D."
1824 ; (sp-pivot-row-number sp-mat) (aref colj (sp-pivot-row-number sp-mat)) j (sp-pivot-row-number sp-mat))
1825 ; (return the-gcd))
1826 ; (t
1827 ; (setq smallest-remainder (abs (aref colj (sp-pivot-row-number sp-mat))))
1828 ; (loop for i below (array-total-size colj)
1829 ; when (and
1830 ; (aref colj i)
1831 ; (not ($zerop (setq rem (mod (aref colj i)
1832 ; (aref colj (sp-pivot-row-number sp-mat))))))
1833 ; (< rem smallest-remainder))
1834 ; do (setq where-the-smallest-remainder i) (setq smallest-remainder rem))
1835 ; (setq factor (special-times -1
1836 ; (aref colj where-the-smallest-remainder)
1837 ; (special-inverse (aref colj (sp-pivot-row-number sp-mat)))))
1838 ; ( sp-set-current-row sp-mat where-the-smallest-remainder)
1839 ; ( sp-row-operation sp-mat factor)
1840 ; (aset (special-plus (special-times (aref colj (sp-pivot-row-number sp-mat))
1841 ; factor)
1842 ; (aref colj (sp-current-row-number sp-mat)))
1843 ; colj (sp-current-row-number sp-mat)))))
1844 ; (aref colj (sp-pivot-row-number sp-mat))))
1845
1846 (defun gcd-array (an-array &aux where the-gcd temp)
1847 "Gcd of all non-nil elements occurring in an array"
1848 (loop for ii below (array-total-size an-array)
1849 until (setq the-gcd (aref an-array ii)))
1850 (loop for ii below (array-total-size an-array)
1851 when (setq temp (aref an-array ii))
1852 do (setq the-gcd (gcd the-gcd temp))
1853 (cond ((equal the-gcd (abs temp)) (setq where ii)
1854 (cond ((eql the-gcd 1) (return 'done))))))
1855 (values the-gcd where))
1856
1857 (defun rational-inverse (x)
1858 (sp-rational-quotient 1 x))
1859
1860 (defun divides (a b)
1861 ($zerop (mod b a)))
1862
1863 (defun sp-find-good-column-to-pivot (sp-mat &aux the-gcd)
1864 "Tries to find a column where the gcd is equal to the entry of the pivot-row
1865 without changing the pivot row"
1866 (loop for ii below (length (the cl:array (sp-pivot-row sp-mat))) by 2
1867 when (aref (sp-pivot-row sp-mat) ii)
1868 do ( sp-set-current-column-above-pivot-row-number sp-mat
1869 (aref (sp-pivot-row sp-mat) ii))
1870 (setq the-gcd (gcd-array (sp-current-column-above-pivot-row-number sp-mat)))
1871 (cond ((eql (abs (aref (sp-pivot-row sp-mat) (1+ ii))) the-gcd)(return (aref (sp-pivot-row sp-mat) ii))))))
1872
1873 (defun sp-gcd-column (sp-mat j )
1874 (setf (sp-pivot-row-number sp-mat) 0)
1875 ( sp-set-current-column-above-pivot-row-number sp-mat j)
1876 (gcd-array (sp-current-column-above-pivot-row-number sp-mat)))
1877
1878 (defun show-gcd-of-columns (a-sparse-matrix)
1879 (let (the-gcd where (all-col ( sp-list-of-all-columns-occurring a-sparse-matrix)))
1880 (loop for i in all-col do
1881 (multiple-value (the-gcd where) ( sp-gcd-column a-sparse-matrix i))
1882 (format t "~%Column ~D has gcd ~D at it occurs in row ~A. "i the-gcd where))))
1883
1884 (defun sp-show-row-array-leaders(sp-mat )
1885 (loop for i below (array-total-size (sp-rows sp-mat))
1886 do
1887 (format t "~% Row ~D has array-leader ~A." i (list-array-leader
1888 ( sp-row sp-mat i)))))
1889
1890 (defun sp-make-rows-desirable (sp-mat &rest a-list &aux temp this-row)
1891 (loop for i in a-list
1892 do
1893 (setq this-row (aref (sp-rows sp-mat) i))
1894 ; (setq temp (make-array (array-total-size (aref (sp-rows sp-mat) i))
1895 ; :leader-list (list (fill-pointer this-row ) nil)))
1896 (setq temp (make-sparse-solution :row (MAKE-ARRAY (array-total-size (aref (sp-rows sp-mat) i)) :adjustable t
1897 :fill-pointer (fill-pointer this-row ) )))
1898 (fillarray (sparse-solution-row temp) this-row)
1899 (setf (aref (sp-rows sp-mat) i) temp)))
1900
1901 ;si:
1902 ;(defun-hash-table sp-get-key (sp-mat value)
1903 ; (catch 'key
1904 ; (send self :map-hash `(lambda (u v) (cond ((eq v ,value) (throw 'key u)))))))
1905
1906 (defvar *verbose* nil)
1907 (defun fmat (&rest l)
1908 (cond (*verbose* (apply 'format l))
1909 (t nil)))
1910
1911
1912 (defun sp-show-column (sp-mat j)
1913 (loop for i below (sp-number-of-rows sp-mat) do
1914 (format t "~% entry in row ~D col ~D is ~A" i j ( sp-entry sp-mat i j))))
1915 (defun sp-verify-columns-above-pivots-are-zero(sp-mat &aux (ok t) col )
1916 "Verifies the entries above the pivot are indeed 0"
1917 (loop for i below (sp-number-of-rows sp-mat)
1918 when (setq col (aref (sp-column-used-in-row sp-mat) i))
1919 do
1920 (setq ok t)
1921 (loop for ii from (1+ i) below (sp-number-of-rows sp-mat)
1922 when ( sp-entry sp-mat ii col)
1923 do (format t "~%Pivot in row ~A has nonzero entry above it in column ~A row ~D."
1924 i col ii)
1925 (setq ok nil))
1926 (cond (ok
1927 (format t "~%For Row ~D column ~D is ok above the pivot" i col)))))
1928
1929 (defun sp-reduce-row-with-respect-to-rows (sp-mat a-row
1930 &aux (not-done t)
1931 piv-row col factor)
1932 "Cleans out all the columns in a-row such that rows has a pivot in that column."
1933 (setf (sp-current-row sp-mat) a-row)
1934 (setf (sp-current-row-number sp-mat) "Fake row number")
1935 (loop while not-done
1936 do
1937 (loop for ii below (length (the cl:array a-row)) by 2
1938 when (and (setq col (aref a-row ii))(setq piv-row (gethash col (sp-columns-used-to-pivot sp-mat)
1939 )))
1940 do
1941 ( sp-set-pivot-row sp-mat piv-row)
1942 (setf (sp-pivot-entry sp-mat) ( sp-entry sp-mat piv-row col))
1943 (setq factor
1944 (special-times -1
1945 (special-inverse (sp-pivot-entry sp-mat))
1946 (aref a-row (1+ ii))))
1947 ( sp-row-operation sp-mat factor)
1948 (return 'try-again)
1949 finally (setq not-done nil)))
1950 ( sp-set-current-row sp-mat 0)
1951 a-row)
1952
1953
1954
1955 ;;((0 0 2 1 1 0) (2 6 6 4 0 6) (4 0 5 0 0 5) (4 1 2 0 2 3) (5 6 5 4 5 6) (5 1 0 0 6 6))
1956 ;;need to fix this, or the reduce function: If given an
1957 ;;integer matrix it may return a matrix with more rows than
1958 ;;it started with, Actually I guess the reduce function should
1959 ;;be fixed for example with the above matrix it returns to big a thing..
1960 (defun sp-determinant(sp-mat &aux answer col entry (correct 0) tem)
1961 (setf (sp-sign-of-row-permutation sp-mat) 1)
1962 (cond ((eq (sp-type-of-entries sp-mat) :integer)
1963 (format t "Integer type is not good for determinant")
1964 (iassert (null (sp-reduced sp-mat)))
1965 (sp-set-type-of-entries sp-mat :any-macsyma)))
1966 ( sp-reduce sp-mat)
1967 (loop for i below (sp-number-of-rows sp-mat)
1968 do (show i)
1969 when (eql 0 (setq tem (aref (sp-column-used-in-row sp-mat) i)))
1970 do (return 'done)
1971 when (null tem) do (return (setq answer 0))
1972 finally (setq correct -1))
1973 (cond ((null answer)
1974 (setq answer (sp-sign-of-row-permutation sp-mat))
1975 (loop for i below (sp-number-of-rows sp-mat)
1976 when (setq col (aref (sp-column-used-in-row sp-mat) i))
1977 do
1978 (setq entry ( sp-entry sp-mat i col))
1979 (setq answer (sp-mul* answer entry))
1980 and
1981 collecting (+ col correct) into tem
1982 else
1983 do (format t "~%There is no pivot in row ~A" i)(return (setq answer 0 ))
1984 finally
1985 (return (setq answer (sp-mul* answer
1986 (sign-of-permutation tem)))))))
1987 answer)
1988
1989 (defun sp-product-of-pivots (sp-mat &aux entry (answer 1) col )
1990 (loop for i below (sp-number-of-rows sp-mat)
1991 when (setq col (aref (sp-column-used-in-row sp-mat) i))
1992 do
1993 (setq entry ( sp-entry sp-mat i col))
1994 (cond (entry
1995 (setq answer (sp-mul* answer entry)))
1996 (t (format t "~%Zero entry in (~A ,~A)" i col)))
1997 else
1998 do (format t "~%There is no pivot in row ~A" i)
1999 finally
2000 (return answer)))
2001
2002 (defun gcd-and-cofs (b a )
2003 (let ((bp 0 ) (b 1)(ap 1)( a 0)
2004 (rp a ) (rpp b)r q)
2005 (loop
2006 do
2007 (setq q (quotient rpp rp))
2008 ; (print (list r q a b))
2009 (setq r (- rpp (* q rp)))
2010 (cond ((eql r 0) (return (list a b rp))))
2011 (setq b (prog1 (- bp (* b q))
2012 (setq bp b)))
2013 (setq a (prog1 (- ap (* a q))
2014 (setq ap a)))
2015 (setq rpp rp rp r))))
2016
2017 ;(defun test (m n) (setq ans (gcd-and-cofs m n))
2018 ; (setq gc (+ (* m (nth 0 ans)) (* n (nth 1 ans))))
2019 ; (- gc (gcd m n)))
2020
2021 (defun convert-to-sparse-matrix (row-list-matrix &key re-use-sparse-matrix
2022 (old-index-base 1)
2023 (rows-to-use nil ruse) (columns-to-use nil cuse)
2024 (rows-not-to-use nil rnuse)
2025 (renumber-columns 0)
2026 (columns-not-to-use nil cnuse)
2027 &aux sp ar ar-rows jj)
2028 "converts a macsyma matrix or a list of lists to the corresponding sparse matrix
2029 and sets up the type of entries. It can be used to pick out a submatrix. Note
2030 If you don't renumber the columns the determinant of submatrices may have wrong
2031 sign"
2032
2033 (setq row-list-matrix (if ($matrixp row-list-matrix) (cdr row-list-matrix)
2034 row-list-matrix))
2035 (loop for v in row-list-matrix
2036 for i from old-index-base
2037 do (setq v (if ($listp v) (cdr v) v))
2038 when (and (or (null ruse) (member i rows-to-use))
2039 ( or (null rnuse) (not (member i rows-not-to-use))))
2040 collecting
2041 (setq ar (make-array (* 2 (loop for u in v when (not ($zerop u)) counting u))
2042 :fill-pointer 0 :adjustable t))
2043 into rows
2044 and
2045 do
2046 (setq jj (1- renumber-columns))
2047 (loop for w
2048 in v
2049 for j from old-index-base
2050 when
2051 (and (or (null cuse) (member j columns-to-use))
2052 ( or (null cnuse) (not (member j columns-not-to-use))))
2053 do
2054 (cond (renumber-columns (incf jj))
2055 (t (setq jj j)))
2056
2057 and when
2058 (not ($zerop w) )
2059 do
2060 (vector-push-extend jj ar)
2061 (vector-push-extend w ar ))
2062 finally
2063 (setq ar-rows (make-array (length rows) :fill-pointer (length rows)))
2064 (fillarray ar-rows rows)
2065
2066 (cond (re-use-sparse-matrix (setq sp re-use-sparse-matrix))
2067 (t (setq sp (make-sparse-matrix))))
2068
2069 (setf (sp-constants-column sp) nil)
2070 (sp-set-rows sp ar-rows)
2071 (sp-choose-type-of-entries sp)
2072 (return sp)))
2073
2074
2075
2076 ;;((5 3 1 1 4) (1 1 2 1 5) (4 5 0 1 4) (1 1 1 1 0) (3 1 5 2 4))
2077 ;;the sp-determinant function does not work in float mode.
2078 ;;I think the problem is that it allows pivots which should be zero entries
2079 ;;since the $zerop function only recognize 3.1111e-50 as non zero yet its really
2080 ;;zero.
2081
2082 ;;the following matrix has determinant -36 according to math:determinant.
2083 ;((1 1 3 5 3) (3 3 3 1 3) (2 1 4 3 4) (0 0 0 4 4) (4 3 0 0 4))
2084 ;;Macsyma makes it 152 and my program makes it +152.
2085 ;;((5 5 ) (3 4)) gave a determinant of 20 with their program but 5 with mine.
2086 ;;Note that the :integer type for the determinant program is bad since the
2087 ;;addition of new rows makes trouble.
2088
2089 ;(defun test-determinant (n &aux mat)
2090 ; (setq mat (loop for i below n
2091 ; collecting
2092 ; (loop for j below n collecting (random 6))))
2093 ; (dett mat))
2094
2095 ;;;I had to wrap a floating point $zerop around it :
2096 ;;; (function-let ($zerop float-zerop)
2097 ;;; (sp-determinant sp))
2098 ;;;this should be somehow automatic but I hates to change $zerop to check
2099 ;;;for floating point since I use it so rarely..
2100 ;;;these agreed (since I outlawed the use of sp-determinant coerces away from rational
2101 ;;;type
2102 ;
2103 ;(compare-functions
2104 ;(defun dett (mat &aux a1)
2105 ; (show (setq a1 ($determinant ($coerce_matrix mat))))
2106 ; (iassert (equal (listarray (te1 mat)) (apply 'append mat)))
2107 ; (show mat)
2108 ; (convert-to-sparse-matrix mat :re-use-sparse-matrix *sparse-matrix*)
2109 ; (iassert (equal a1 (sp-determinant *sparse-matrix*))))
2110 ;(defun dett (mat &aux )
2111 ; (setq mat (copy-list mat))
2112 ; (setf (car mat) (mapcar 'float (car mat)))
2113 ; (convert-to-sparse-matrix mat :re-use-sparse-matrix *sparse-matrix*)
2114 ; (round (sp-determinant *sparse-matrix*)))
2115 ;)
2116
2117 (defun make-test (n)
2118 (cons '($matrix)
2119 (loop for i from 1 to (expt n 2)
2120 collecting
2121 (cons '(mlist)
2122 (loop for
2123 j from 1 to (expt n 2)
2124 when (eql i j) collecting 4.
2125 else
2126 when (eql (abs (- i j)) 1) collecting -1
2127 else
2128 when (eql (abs (- i j)) n) collecting -1
2129 else collecting 0)))))
2130
2131
2132 ;;for n=10 the 100 by 100 matrix has determinant
2133 ;;6265706907310779461620940495418760520632027297067936
2134 ;;the only factors I could find were: (613 1 23 1 13 1 2 5)
2135
2136 (defvar *sparse-matrix* (make-sparse-matrix))
2137
2138 (defun $coerce_matrix(mat &aux rows)
2139 (cond (($matrixp mat) mat)
2140 ((arrayp mat)
2141 (let ((dims (array-dimensions mat)))
2142 (case (length dims)
2143 (1 (setq rows (cons '(mlist) (listarray mat))))
2144 (2 (setq rows (loop for i below (car dims)
2145 collecting
2146 (cons '(mlist) (loop for j below (second dims)
2147 collecting (aref mat i j)))))))
2148 (cons '($matrix) rows)))
2149 (t (loop for v in mat
2150 collecting (cons '(mlist) v) into tem
2151 finally (return (cons '($matrix ) tem))))))
2152
2153
2154 (defun sp-float (sp-mat &aux tem)
2155 (loop for i below (sp-number-of-rows sp-mat)
2156 do (let ((row (aref (sp-rows sp-mat) i)))
2157 (loop for j from 1 below (fill-pointer row) by 2
2158 when (numberp (setq tem (aref row j)))
2159 do (setf (aref row j) (float tem))))
2160 finally (sp-set-type-of-entries sp-mat :float)))
2161
2162 (defun float-zerop (n)
2163 (cond ((numberp n) (or (zerop n)(and (floatp n) (< (abs n) .0001))))
2164 ((atom n) nil)
2165 ((affine-polynomialp n) nil)
2166 ((rational-functionp n)(float-zerop (car n)))
2167 (t
2168 (let (tem type-of-n)
2169 (setf type-of-n (caar n))
2170 (cond ((member type-of-n '(mrat rat) :test #'eq)
2171 (equal (cdr n) (rzero)))
2172 (t (and (numberp (setq tem ($ratsimp n)))(zerop tem))))))))
2173
2174 (defun sp-make-sparse-matrix (list-of-rows &key (sp (make-sparse-matrix)) (type-of-entries :any-macsyma) &aux rows)
2175 "Takes a List-of-rows which is a list or array of arrays with fill-pointer and alternating column number and entry
2176 and converts to a sparse matrix."
2177 (cond ((arrayp list-of-rows)(setq rows list-of-rows))
2178 (t
2179 (setq rows (MAKE-ARRAY (length list-of-rows) :fill-pointer (length list-of-rows)))
2180 (fillarray rows list-of-rows)))
2181 (sp-set-rows sp rows)
2182 (sp-set-type-of-entries sp type-of-entries)
2183 sp)
2184
2185 (defun sp-quotient-space-basis ( list-subspace-rows all-rows &aux sp-sub sp)
2186 "takes list-subspace-rows and finds a basis for the space spanned by all-rows
2187 modulo list-subspace-rows. Returns two values: the reduce sparse matrix for the
2188 subspace and then the basis for the quotient space. Note all-rows need not
2189 contain list-subspace-rows."
2190 ; (declare (values sp-quotient sp-sub))
2191 (setq sp-sub (sp-make-sparse-matrix list-subspace-rows))
2192 (sp-reduce sp-sub)
2193 (cond ((arrayp all-rows)(setq all-rows (listarray all-rows))))
2194 (loop for v in all-rows do (sp-reduce-row-with-respect-to-rows sp-sub v)
2195 when (not (zerop (fill-pointer v)))
2196 collecting v into real-rows
2197 finally
2198 (setq all-rows real-rows))
2199 (setq sp (sp-make-sparse-matrix all-rows))
2200 (sp-reduce sp)
2201 (values sp sp-sub))
2202
2203 (defun $fast_determinant(matrix &aux sp)
2204 (cond ((null *sparse-matrix*) (setq *sparse-matrix* (make-sparse-matrix))))
2205 (setq sp (convert-to-sparse-matrix matrix :re-use-sparse-matrix *sparse-matrix*))
2206 (new-disrep (sp-determinant sp)))
Maxima, a Computer Algebra System