読者です 読者をやめる 読者になる 読者になる

桁数が爆発的に増える数値計算実験(GaucheとMIT Schemeを比較)

前回までのScheme数値計算

nibosiiwasi.hatenablog.com

nibosiiwasi.hatenablog.com


この後、ちょっとパソコンを新しく買ったりしてて更新遅くなりました。新しいパソコンのWindows 10にMIT Schemeをインストールしました。
そこで、MIT SchemeGaucheを比較してみました。

ヒルベルト行列の逆行列を正確数を使って計算してみました。


<OS:Mac 処理系:Gauche 0.9.2>

ヒルベルト行列のサイズを大きくして逆行列を計算すると、
CPU使用率100%のまま止まらないです。手動で強制終了するしかない。
heapのサイズの決め方がよくわからないです。



<OS:Windows 10 処理系:MIT Scheme 9.2>

オプションで --heap 256 とする。

(あと、MITは引数の小文字と大文字を判別してくれない所があったようなのでソースを少しだけ変更した)

ヒルベルト行列のサイズを大きくして逆行列を計算すると、

サイズ52x52では逆行列計算完了。

サイズ53x53では異常終了。
こんなメッセージが出てくる。

f:id:nibosiiwasi:20161002224407p:plain

これだと、ちゃんとout of memory の警告が自動的に出てくるので安心して使える。


すみません。私にわかるのは、ここまでです。Gaucheの開発者の方に、この状況を報告したいのですが、どうすればいいのでしょうか?

追記2016年10月3日19:58

コメントありがとうございます。コードを追記致します。
この比較実験には、下記のコードを使用しました。

(define (my-map-item2 op a b)
  (define (my-map-item2-pre op a b tmp)
    (if (and (null? a) (null? b))
	tmp
	(my-map-item2-pre op
			  (cdr a)
			  (cdr b)
			  (append tmp (list (op (car a) (car b)))))))
  (my-map-item2-pre op a b '() ))

  
(define (add-mat A1 A2)
  (define (add-each a1 a2) (my-map-item2 + a1 a2))
  (my-map-item2 add-each A1 A2))

(define (scalar-mat a AA)
  (define (scalar-each a a1) (map (lambda (x) (* a x)) a1))
  (map (lambda (x) (scalar-each a x)) AA))

(define (sub-mat A1 A2) (add-mat A1 (scalar-mat -1 A2)))

(define (l-accum2 op1 op2 ini a b)
  (if (and (null? a) (null? b))
      ini
      (l-accum2 op1 op2 (op1 ini
			     (op2 (car a) (car b)))
		(cdr a) (cdr b))))

(define (sum-of-mul-each a1 a2) (l-accum2 + * 0 a1 a2))

(define (col1sep-mat A)
  (define (col1sep-mat-pre  A tmp1 tmp2)
    (if (null? A)
	(cons tmp1 tmp2)
	(col1sep-mat-pre (cdr A)
				(append tmp1 (list (list (caar A))))
				(append tmp2 (list (cdr (car A)))))))
  (col1sep-mat-pre A '() '() ))

(define (col2line col)
  (define (col2line-p col tmp)
    (if (null? col)
	tmp
	(col2line-p (cdr col)
		    (append tmp (list (caar col))))))
  (col2line-p col '() ))

(define (transposed-mat A)
  (define (transposed-mat-p A tmp)
    (if (null? (car A))
	tmp
	(transposed-mat-p (cdr (col1sep-mat A))
			  (append tmp
				  (list (col2line (car (col1sep-mat A))))))))
  (transposed-mat-p A '() ))

(define (mul-oneline-mat oneline A2)
  (map (lambda (x) (sum-of-mul-each oneline x)) (transposed-mat A2)))

(define (mul-mat A1 A2)
  (define (mul-mat-pre A1 A2 tmp)
    (if (null? A1)
	tmp
	(mul-mat-pre (cdr A1) A2
		     (append tmp
			     (list (mul-oneline-mat (car A1) A2))))))
  (mul-mat-pre A1 A2 '() ))

(define (eye-line len loc1)
  (define (eye-line-p len loc1 tmp)
    (if (<= len 0)
	tmp
	(eye-line-p (- len 1)
		    (- loc1 1)
		    (append tmp (list (if (= loc1 1) 1 0))))))
  (eye-line-p len loc1 '() ))

(define (eye-mat N)
  (define (eye-mat-pre N i tmp)
    (if (> i N)
	tmp
	(eye-mat-pre N (+ i 1) (append tmp (list (eye-line N i))))))
  (eye-mat-pre N 1 '()))


(define (x-th-item x list-a)
  (if (null? list-a)
      (error "x-th-item")
      (if (= x 1)
          (car list-a)
          (x-th-item (- x 1) (cdr list-a)))))

(define (extract-from-list x1 x2 list-a)
  (define (extract-from-list-p x1 x2 list-a tmp count)
    (cond ((> count x2) tmp)
          ((and (>= count x1)
                (<= count x2))
           (extract-from-list-p x1 x2 (cdr list-a)
                                (append tmp (list (car list-a)))
                                (+ count 1)))
          ((< count x1) (extract-from-list-p x1 x2 (cdr list-a) tmp (+ count 1)))))
  (if (and (> x1 0)
           (>= x2 x1)
           (>= (length list-a) x2))
      (extract-from-list-p x1 x2 list-a '() 1)
      (error "EXTRACT-FROM-LIST: (extract-from-list x1 x2 list-a) WRONG X1 X2!")))

(define (size-mat A)
  (define (N-mat A j)
    (if (null? (car A))
        j
        (N-mat (list (cdr (car A))) (+ j 1))))
  (define (M-mat A i)
    (if (null? A)
        i
        (M-mat (cdr A) (+ i 1))))
  (cons (M-mat A 0) (N-mat A 0)))


(define (x-th-scaled-list-to-each x list-a A)
  (define (x-th-scaled-list-to-each-p x list-a A tmp)
    (if (null? A)
        tmp
        (x-th-scaled-list-to-each-p x list-a (cdr A)
                                    (append tmp
                                            (scalar-mat (/ (x-th-item x (car A))
                                                           (x-th-item x list-a))
                                                        (list list-a))))))
  (x-th-scaled-list-to-each-p x list-a A '()))


(define (gyakugyouretu-hakidashi AE)
  (define (hakidashi-ahead i tmp-mat-up tmp-mat-low b-line size-mat-AE)
    (if (> i size-mat-AE)
        tmp-mat-up
        (begin (set! b-line (map (lambda (a) (/ a (x-th-item i (car tmp-mat-low))))
                                 (car tmp-mat-low))) ;; normalize                                                       
               (set! tmp-mat-up (append tmp-mat-up (list b-line)))
               (if (not (null? (cdr tmp-mat-low))) ;;or (< i size-mat-AE)                                               
                   (set! tmp-mat-low (sub-mat (cdr tmp-mat-low)
                                              (x-th-scaled-list-to-each i b-line (cdr tmp-mat-low)))))                                                                                            
               (hakidashi-ahead (+ i 1) tmp-mat-up tmp-mat-low b-line size-mat-AE))))                
  (define (hakidashi-back i tmp-mat-up tmp-mat-low size-mat-AE)                                                                                           
    (if (> i size-mat-AE)
        (reverse tmp-mat-up)
        (begin (set! tmp-mat-up (append tmp-mat-up (list (car tmp-mat-low))))
               (if (not (null? (cdr tmp-mat-low))) ;;or (< i size-mat-AE)                                               
                   (set! tmp-mat-low (sub-mat (cdr tmp-mat-low)
                                              (x-th-scaled-list-to-each (- (+ size-mat-AE 1) i)
                                                                        (car tmp-mat-low)
                                                                        (cdr tmp-mat-low)))))
               (hakidashi-back (+ i 1) tmp-mat-up tmp-mat-low size-mat-AE))))
  (let ((N (car (size-mat AE))))
    (hakidashi-back 1 '() (reverse (hakidashi-ahead 1 '() AE '() N))
                    N)))

(define (connect-mat-side A B)
  (define (connect-mat-side-p A B tmp)
    (if (null? A)
        tmp
        (connect-mat-side-p (cdr A) (cdr B)
                            (append tmp (list (append (car A) (car B)))))))
  (if (= (car (size-mat A))
         (car (size-mat B)))
      (connect-mat-side-p A B '())
      (error "CONNECT-MAT-SIDE: (connect-mat-side A B) WRONG SIZE MATRIXES !")))

(define (extract-mat A i j k l)
  (define (extract-mat-p D k l tmp)
    (if (null? D)
        tmp
        (extract-mat-p (cdr D) k l
                       (append tmp (list (extract-from-list k l (car D)))))))

  (extract-mat-p (extract-from-list i j A) k l '()))


(define (gyakugyouretu A)
  (let ((N (car (size-mat A))))
    (extract-mat (gyakugyouretu-hakidashi (connect-mat-side A (eye-mat N)))
                 1 N (+ N 1) (* 2 N))))


(define (hilbert-matrix-r n)
  (define (hilbert-matrix-r-p op ini count n)
    (if (= count n)
        '()
        (cons (op ini n)
              (hilbert-matrix-r-p op (+ ini 1) (+ count 1) n))))
  (define (hilbert-matrix-r-each-line ini n)
    (if (= n 0)
        '()
        (cons (/ 1 ini) (hilbert-matrix-r-each-line (+ ini 1) (- n 1)))))
  (hilbert-matrix-r-p hilbert-matrix-r-each-line 1 0 n))


以下がテスト用のコードです。

52x52のヒルベルト行列をhilbとし、その逆行列を計算するコード。

(define hilb (hilbert-matrix-r 52))

(gyakugyouretu hilb)

52の部分をいろいろな自然数に変えて実験しました。



行列と逆行列の積が単位行列になっているか検算するコード

(mul-mat hilb (gyakugyouretu hilb))