今日は、逆行列を計算するプログラムを書いてみた。このページの後ろの方で、動作確認してみる。

まず逆行列を単純にガウスの消去法（＝掃き出し法）で計算するプログラム。

ページ最後の追記をご参照下さい。

( 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 ))) ( set! tmp-mat-up ( append tmp-mat-up ( list b-line ))) ( if ( not ( null? ( cdr tmp-mat-low ))) ( set! tmp-mat-low ( sub-mat ( cdr tmp-mat-low ) ( x-th-scaled-list-to-each i b-line ( cdr tmp-mat-low ))))) tml is null. ( 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 ) rsed . and output-matrix too. ( 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 ))) ( 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 ))))

Schemeでは、分数の計算が簡単にできる。Schemeは途中の演算で小数点を使わなければ自動的に分数計算する。以前のブログ

http://nibosiiwasi.hatenablog.com/entry/2016/08/04/231338

のソースを一部変更して、分数計算するようにする。

( define ( sub-mat A1 A2 ) ( add-mat A1 ( scalar-mat -1 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 '() ))

それでは、ためしに以前の日記

http://nibosiiwasi.hatenablog.com/entry/2016/08/20/221037

で書いたヒルベルト行列の逆行列を計算してみる。

hilbert-matrix-r はヒルベルト行列を生成する。

( 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 ))

実験として、20x20のヒルベルト行列の逆行列を分数を用いて計算してみる。結構大きいが、結果はすぐに返ってきた。（スペースを取るので実行結果は下に書く。）はやかった。自分用で使う分には、このレベルで十分強力な気がする。

逆行列の計算にかかった時間を印字してみた。

;(time (gyakugyouretu hil20)) ; real 0.029 ; user 0.030 ; sys 0.000

逆行列の右下の方を見ると、かなり桁数が大きいので、メモリーはだいぶ食ってると思う。メモリーがどういう風に使われたのかどうすれば確認できるんだろう？

(define hil20 (hilbert-matrix-r 20)) hil20 => ((1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20) (1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21) (1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22) (1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23) (1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24) (1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25) (1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26) (1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27) (1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28) (1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29) (1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30) (1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31) (1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32) (1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33) (1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34) (1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35) (1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36) (1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36 1/37) (1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36 1/37 1/38) (1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36 1/37 1/38 1/39)) (define gh20 (gyakugyouretu hil20)) gh20 => ((400 -79800 5266800 -171609900 3294910080 -41186376000 356948592000 -2237302782000 10440746316000 -37006645275600 100927214388000 -213323430411000 350069219136000 -444318624288000 431623806451200 -314725692204000 166619484108000 -60440401098000 13431200244000 -1378465288200) (-79800 21226800 -1576089900 54777880080 -1095557601600 14085740592000 -124619677182000 793496720016000 -3749272002075600 13423319513604000 -36914128662411000 78568660369836000 -129700645689888000 165464255684851200 -161454280100652000 118188754060608000 -62787775594698000 22846471615044000 -5091096452488200 523816809516000) (5266800 -1576089900 124826320080 -4519175106600 92965887907200 -1220177278782000 10966531592016000 -70700557753425600 337434480186804000 -1218166245859563000 3373383450072636000 -7222704706855638000 11984339661745651200 -15357151230750252000 15043739981143104000 -11050648504666848000 5888832426829044000 -2148710655394888200 480017665520316000 -49500688499262000) (-171609900 54777880080 -4519175106600 168285473017200 -3533994933361200 47119932444816000 -428791385247825600 2792314957736304000 -13438015734105963000 48851589962162988000 -136086572037454038000 292867300483909351200 -488112167473182252000 627936850324010304000 -617257652189248224000 454821427928919744000 -243045200549516488200 88904324471894316000 -19906187129228862000 2057028610969332000) (3294910080 -1095557601600 92965887907200 -3533994933361200 75391891911705600 -1017790540808025600 9355448405407104000 -61430929070198688000 297703733186347488000 -1088692576299632304000 3048339213638970451200 -6589514260887960402000 11025592488805999104000 -14233235274010900224000 14034489776092380672000 -10369928556779370163200 5555318869703234016000 -2036717251537942512000 456976817575340832000 -47311658052294636000) (-41186376000 14085740592000 -1220177278782000 47119932444816000 -1017790540808025600 13878961920109440000 -128637415574347680000 850582094818135680000 -4146587712238411440000 15241696068194852256000 -42867270191798021970000 93028436624300617440000 -156195893591418320640000 202261764420154897920000 -199991479309316424576000 148141836525419573760000 -79542065634387214320000 29222464913370479520000 -6569041752645524460000 681287875953042758400) (356948592000 -124619677182000 10966531592016000 -428791385247825600 9355448405407104000 -128637415574347680000 1200615878693911680000 -7986020779125829440000 39131501817716564256000 -144478577313097037010000 407939512413450457440000 -888364144801191698640000 1496192243875691281920000 -1942836614902487880576000 1925843874830454458880000 -1429793179798367704320000 769290895652382623520000 -283160273442983396460000 63763498612345890758400 -6623632127321249040000) (-2237302782000 793496720016000 -70700557753425600 2792314957736304000 -61430929070198688000 850582094818135680000 -7986020779125829440000 53392253209012688256000 -262789996263109325010000 974059652089259255520000 -2759835680919567890640000 6028662370412383621920000 -10181740892252025672576000 13254337256186068922880000 -13168270131145899644160000 9796684058920137899520000 -5281025000511636836460000 1947218076313218066758400 -439187362891157921040000 45689544061930248480000) (10440746316000 -3749272002075600 337434480186804000 -13438015734105963000 297703733186347488000 -4146587712238411440000 39131501817716564256000 -262789996263109325010000 1298491746241246076520000 -4829712441609243808620000 13726551149836798192920000 -30067953572431763314326000 50908704461260128362880000 -66422303749750640852160000 66127617397710930821760000 -49289566671441943806960000 26616366002578649655758400 -9829706635234995048540000 2220336112394187267480000 -231303316813521882930000) (-37006645275600 13423319513604000 -1218166245859563000 48851589962162988000 -1088692576299632304000 15241696068194852256000 -144478577313097037010000 974059652089259255520000 -4829712441609243808620000 18019628875711681578360000 -51355942295778292498326000 112776921688485978803880000 -191379018622885297364160000 250215634844739665861760000 -249577329653809207530480000 186351072808177541622758400 -100790844667884487656540000 37278274729144512543480000 -8431990712544592122930000 879523723192157283240000) (100927214388000 -36914128662411000 3373383450072636000 -136086572037454038000 3048339213638970451200 -42867270191798021970000 407939512413450457440000 -2759835680919567890640000 13726551149836798192920000 -51355942295778292498326000 146731263702223692852360000 -322952093926118939302020000 549174575178714331566720000 -719369950178626539352560000 718782709402970517687782400 -537551171562050600834880000 291173551262777408785560000 -107840723323596625572210000 24423697236336059942280000 -2550618797257256121396000) (-213323430411000 78568660369836000 -7222704706855638000 292867300483909351200 -6589514260887960402000 93028436624300617440000 -888364144801191698640000 6028662370412383621920000 -30067953572431763314326000 112776921688485978803880000 -322952093926118939302020000 712281693323218846365720000 -1213516958995113590104560000 1592364028560688696695782400 -1593613921832714355013440000 1193570920878823097170560000 -647405655744540206724210000 240083137818163818458280000 -54438604501155664185396000 5691463431896356634520000) (350069219136000 -129700645689888000 11984339661745651200 -488112167473182252000 11025592488805999104000 -156195893591418320640000 1496192243875691281920000 -10181740892252025672576000 50908704461260128362880000 -191379018622885297364160000 549174575178714331566720000 -1213516958995113590104560000 2071068943351660527111782400 -2721989792411433669565440000 2728162104865881364961280000 -2046121578649411023720960000 1111255684956145642193280000 -412587318324548191720896000 93657814195536626555520000 -9801964799377058648340000) (-444318624288000 165464255684851200 -15357151230750252000 627936850324010304000 -14233235274010900224000 202261764420154897920000 -1942836614902487880576000 13254337256186068922880000 -66422303749750640852160000 250215634844739665861760000 -719369950178626539352560000 1592364028560688696695782400 -2721989792411433669565440000 3582789983174023804385280000 -3595846797981509313876480000 2700329396185347990497280000 -1468304109175782969832896000 545758754722107075515520000 -124016939248194856280340000 12991953343553024480640000) (431623806451200 -161454280100652000 15043739981143104000 -617257652189248224000 14034489776092380672000 -199991479309316424576000 1925843874830454458880000 -13168270131145899644160000 66127617397710930821760000 -249577329653809207530480000 718782709402970517687782400 -1593613921832714355013440000 2728162104865881364961280000 -3595846797981509313876480000 3613560329006048768624640000 -2716862025141584814928896000 1478936989492394959739520000 -550283540316104136728340000 125167374677213361440640000 -13124524296038259424320000) (-314725692204000 118188754060608000 -11050648504666848000 454821427928919744000 -10369928556779370163200 148141836525419573760000 -1429793179798367704320000 9796684058920137899520000 -49289566671441943806960000 186351072808177541622758400 -537551171562050600834880000 1193570920878823097170560000 -2046121578649411023720960000 2700329396185347990497280000 -2716862025141584814928896000 2044949911396891796183040000 -1114337939999478146748180000 415028663403391672145280000 -94489096570053223832640000 9916307245895573787264000) (166619484108000 -62787775594698000 5888832426829044000 -243045200549516488200 5555318869703234016000 -79542065634387214320000 769290895652382623520000 -5281025000511636836460000 26616366002578649655758400 -100790844667884487656540000 291173551262777408785560000 -647405655744540206724210000 1111255684956145642193280000 -1468304109175782969832896000 1478936989492394959739520000 -1114337939999478146748180000 607820694545169898226280000 -226587340130160526888140000 51631542054350511594264000 -5422980525099141914910000) (-60440401098000 22846471615044000 -2148710655394888200 88904324471894316000 -2036717251537942512000 29222464913370479520000 -283160273442983396460000 1947218076313218066758400 -9829706635234995048540000 37278274729144512543480000 -107840723323596625572210000 240083137818163818458280000 -412587318324548191720896000 545758754722107075515520000 -550283540316104136728340000 415028663403391672145280000 -226587340130160526888140000 84541831107387625158264000 -19279944336904242362910000 2026580957476495940520000) (13431200244000 -5091096452488200 480017665520316000 -19906187129228862000 456976817575340832000 -6569041752645524460000 63763498612345890758400 -439187362891157921040000 2220336112394187267480000 -8431990712544592122930000 24423697236336059942280000 -54438604501155664185396000 93657814195536626555520000 -124016939248194856280340000 125167374677213361440640000 -94489096570053223832640000 51631542054350511594264000 -19279944336904242362910000 4400227536350517776520000 -462861082880434258020000) (-1378465288200 523816809516000 -49500688499262000 2057028610969332000 -47311658052294636000 681287875953042758400 -6623632127321249040000 45689544061930248480000 -231303316813521882930000 879523723192157283240000 -2550618797257256121396000 5691463431896356634520000 -9801964799377058648340000 12991953343553024480640000 -13124524296038259424320000 9916307245895573787264000 -5422980525099141914910000 2026580957476495940520000 -462861082880434258020000 48722219250572027160000)) (mul-mat hil20 gh20) => ((1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1))