-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpow.ijs
341 lines (315 loc) · 10.4 KB
/
pow.ijs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
NB. Raise matrix to an integer power(s)
NB.
NB. gepow Raise a general matrix to integer power(s)
NB. dipow Raise a diagonalizable matrix to integer
NB. power(s)
NB. hepow Raise a Hermitian (symmetric) matrix to
NB. integer power(s)
NB.
NB. testgepow Test gepow by square matrix
NB. testdipow Test dipow by diagonalizable matrix
NB. testhepow Test hepow by Hermitian (symmetric) matrix
NB. testpow Adv. to make verb to test xxpow by matrix of
NB. generator and shape given
NB.
NB. Copyright 2010,2011,2013,2017,2018,2020,2021,2023,2024,
NB. 2025 Igor Zhuravlov
NB.
NB. This file is part of mt
NB.
NB. mt is free software: you can redistribute it and/or
NB. modify it under the terms of the GNU Lesser General
NB. Public License as published by the Free Software
NB. Foundation, either version 3 of the License, or (at your
NB. option) any later version.
NB.
NB. mt is distributed in the hope that it will be useful, but
NB. WITHOUT ANY WARRANTY; without even the implied warranty
NB. of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
NB. See the GNU Lesser General Public License for more
NB. details.
NB.
NB. You should have received a copy of the GNU Lesser General
NB. Public License along with mt. If not, see
NB. <http://www.gnu.org/licenses/>.
NB. =========================================================
NB. Configuration
coclass 'mt'
NB. =========================================================
NB. Local definitions
NB. =========================================================
NB. Interface
NB. ---------------------------------------------------------
NB. gepow
NB.
NB. Description:
NB. Raise a general matrix A to integer power(s)
NB.
NB. Syntax:
NB. P=. p gepow A
NB. where
NB. A - n×n-matrix, a general matrix
NB. p - sh-array of non-negative integers, power(s)
NB. P - sh×n×n-array if r>0, the matrix A in power(s) p
NB. n×n-array if r=0, the matrix A in power p
NB. sh - r-vector of non-negative integers, the shape of p
NB. r ≥ 0, the rank of p
NB.
NB. References:
NB. [1] Roger Hui. Linear Recurrences and Matrix Powers.
NB. 2006-08-09 09:20:34
NB. http://code.jsoftware.com/wiki/Essays/Linear_Recurrences
gepow=: 4 : 0
mpi3=. mp/^:(3 = #@$) NB. apply mp/ only to 3-rank arrays (stiff rank)
pl=. i. >: <. 2 ^. >./ , x NB. powers list: 2^i
pc=. mp~^:pl y NB. powers cache: A^2^i
pb=. (<@I.@|."1@#:) x NB. pl bits boxed array of shape sh
pc mpi3@({~ >)"3 0 pb NB. extract and mp A's powers for each pl atom
)
NB. ---------------------------------------------------------
NB. dipow
NB.
NB. Description:
NB. Raise a diagonalizable matrix to integer power(s)
NB.
NB. Syntax:
NB. P=. p dipow iLl ; vl ; Ll
NB. P=. p dipow (ct Rl) ; vl ; ct iRl
NB. P=. p dipow (ct iLu) ; vu ; ct Lu
NB. P=. p dipow Ru ; vu ; iRu
NB. where
NB. Ll - n×n-matrix, rows are left eigenvectors of A,
NB. output of geevlvx
NB. Lu - n×n-matrix, columns are left eigenvectors of A,
NB. output of geevuvx
NB. Rl - n×n-matrix, rows are right eigenvectors of A,
NB. output of geevlxv
NB. Ru - n×n-matrix, columns are right eigenvectors of A,
NB. output of geevuxv
NB. vl - n-vector, eigenvalues of A, output of geevlxx
NB. vu - n-vector, eigenvalues of A, output of geevuxx
NB. iLl -:%. Ll
NB. iLu -:%. Lu
NB. iRl -:%. Rl
NB. iRu -:%. Ru
NB. p - sh-array of positive integers, power(s)
NB. P - sh×n×n-array if r>0,
NB. n×n-array if r=0, a matrix A in power(s) p
NB. sh - r-vector of non-negative integers, the shape of p
NB. r ≥ 0, the rank of p
NB.
NB. Assertions (with appropriate comparison tolerance):
NB. - lower case:
NB. ((-: ~.) v) +. ((-: ct) A) NB. A must be normal (diagonalizable)
NB. (L mp A) -: V mp L NB. verify...
NB. (A mp ct R) -: (ct R) mp V NB. ...eigendecomposition
NB. W -: L mp (ct R) NB. LAPACK doesn't guarantee (W -: idmat # A)
NB. iL -: (ct R) %"1 w NB. see [1]
NB. (ct iR) -: L % w NB. see [1]
NB. A -: iL mp V mp L
NB. A -: iL mp v * L
NB. A -: (ct R) mp V mp ct iR
NB. A -: (ct R) mp v * ct iR
NB. F -: diagmat"1 f
NB. P -: iL mp"2 F mp"2 L
NB. P -: iL mp"2 f *"1 2 L
NB. P -: (ct R) mp"2 F mp"2 ct iR
NB. P -: (ct R) mp f *"1 2 ct iR
NB. P -: p dipow iL ; v ; L
NB. P -: p dipow (ct R) ; v ; (ct iR)
NB. where
NB. 'v LR'=. geevuvv A NB. use definition from ggevuxx application notes
NB. 'L R'=. LR
NB. iL=. %. L
NB. iR=. %. R
NB. V=. diagmat v
NB. w=. (ct L) mp"1 |: R
NB. W=. diagmat w
NB. P=. p gepow A
NB. f=. v ^1 0 p
NB. F=. p gepow V
NB. - upper case:
NB. ((-: ~.) v) +. ((-: ct) A) NB. A must be normal (diagonalizable)
NB. ((ct L) mp A) -: V mp ct L NB. verify...
NB. (A mp R) -: R mp V NB. ...eigendecomposition
NB. W -: (ct L) mp R NB. LAPACK doesn't guarantee (W -: idmat # A)
NB. (ct iL) -: R %"1 w NB. see [1]
NB. iR -: (ct L) % w NB. see [1]
NB. A -: (ct iL) mp V mp ct L
NB. A -: (ct iL) mp v * ct L
NB. A -: R mp V mp iR
NB. A -: R mp v * iR
NB. F -: diagmat"1 f
NB. P -: (ct iL) mp"2 F mp"2 ct L
NB. P -: (ct iL) mp"2 f *"1 2 ct L
NB. P -: R mp"2 F mp"2 iR
NB. P -: R mp f *"1 2 iR
NB. P -: p dipow (ct iL) ; v ; ct L
NB. P -: p dipow R ; v ; iR
NB. where
NB. 'v LR'=. geevuvv A NB. use definition from ggevuxx application notes
NB. 'L R'=. LR
NB. iL=. %. L
NB. iR=. %. R
NB. V=. diagmat v
NB. w=. (ct L) mp"1 |: R
NB. W=. diagmat w
NB. P=. p gepow A
NB. f=. v ^1 0 p
NB. F=. p gepow V
NB.
NB. References:
NB. [1] Julien Langou. LAPACK/ScaLAPACK Development - DGEEVX
NB. and left eigenvectors. Dec 22, 2006 5:15 pm.
NB. http://icl.cs.utk.edu/lapack-forum/viewtopic.php?p=985#p985
dipow=: (0 {:: ]) mp"2 ([ ^"1 0~ 1 {:: ]) (*"1 2) 2 {:: ]
NB. ---------------------------------------------------------
NB. hepow
NB.
NB. Description:
NB. Raise a Hermitian (symmetric) matrix to integer
NB. power(s)
NB.
NB. Syntax:
NB. P=. p hepow vl ; iRl
NB. P=. p hepow vu ; Ru
NB. where
NB. Rl - n×n-matrix, rows are eigenvectors of A, output of
NB. heevlv
NB. Ru - n×n-matrix, columns are eigenvectors of A, output
NB. of heevuv
NB. vl - n-vector, eigenvalues of A, output of heevlx
NB. vu - n-vector, eigenvalues of A, output of heevux
NB. iRl -:%. Rl
NB. p - sh-array of positive integers, power(s)
NB. P - sh×n×n-array if r>0,
NB. n×n-array if r=0, a matrix A in power(s) p
NB. sh - r-vector of non-negative integers, the shape of p
NB. r ≥ 0, the rank of p
NB.
NB. Assertions (with appropriate comparison tolerance):
NB. - lower case:
NB. (-: ct) A NB. A must be Hermitian (symmetric)
NB. (R mp A) -: V mp R NB. verify eigendecomposition
NB. iR -: ct R
NB. A -: iR mp V mp R
NB. A -: iR mp v * R
NB. F -: diagmat"1 f
NB. P -: iR mp"2 F mp"2 R
NB. P -: iR mp"2 f *"1 2 R
NB. P -: p hepow v ; R
NB. where
NB. 'v R'=. heevlv A NB. use definition from ggevlxx application notes
NB. iR=. %. R
NB. V=. diagmat v
NB. P=. p gepow A
NB. f=. v ^1 0 p
NB. F=. p gepow V
NB. - upper case:
NB. (-: ct) A NB. A must be Hermitian (symmetric)
NB. (A mp R) -: R mp V NB. verify eigendecomposition
NB. iR -: ct R
NB. A -: R mp V mp iR
NB. A -: R mp v * iR
NB. F -: diagmat"1 f
NB. P -: R mp"2 F mp"2 iR
NB. P -: R mp"2 f *"1 2 iR
NB. P -: p hepow v ; R
NB. where
NB. 'v R'=. heevuv A NB. use definition from ggevuxx application notes
NB. iR=. %. R
NB. V=. diagmat v
NB. P=. p gepow A
NB. f=. v ^1 0 p
NB. F=. p gepow V
hepow=: (1 {:: ]) mp"2 ([ ^"1 0~ 0 {:: ]) *"1 2 ct@(1 {:: ])
NB. =========================================================
NB. Test suite
NB. ---------------------------------------------------------
NB. testgepow
NB.
NB. Description:
NB. Test gepow by square matrix
NB.
NB. Syntax:
NB. log=. testgepow A
NB. where
NB. A - n×n-matrix
NB. log - 6-vector of boxes, test log
NB.
NB. Notes:
NB. - fixed powers vector (p -: 5 7) is used
testgepow=: 'gepow' tdyad ((5 7"_)`]`]`geconi`nan`nan)
NB. ---------------------------------------------------------
NB. testdipow
NB.
NB. Description:
NB. Test dipow by diagonalizable matrix
NB.
NB. Syntax:
NB. log=. testdipow A
NB. where
NB. A - n×n-matrix, the diagonalizable
NB. log - 6-vector of boxes, test log
NB.
NB. Notes:
NB. - fixed powers vector (p -: 5 7) is used
testdipow=: 3 : 0
NB. use for a while the definition from ggevlxx application notes
geevlvv=. {.&.>`(((* *@+@((i. >./)"1@sorim{"0 1 ])) % normsr)"2&.>)"0@ggevlvv@(,: idmat@c)
try.
'v LR'=. geevlvv y NB. eigendecomposition
'L R'=. LR
assert ((-: ~.) v) +. (-: ct) y NB. A must be normal (diagonalizable)
iRh=. L ([ % (mp"1 +)) R NB. reconstruct R^_1^H , see [1] in dipow
catch.
R=. v=. iRh=. _.
end.
log=. ('dipow' tdyad ((5 7"_)`]`]`geconi`nan`nan)) (ct R) ; v ; iRh
)
NB. ---------------------------------------------------------
NB. testhepow
NB.
NB. Description:
NB. Test hepow by Hermitian (symmetric) matrix
NB.
NB. Syntax:
NB. log=. testhepow A
NB. where
NB. A - n×n-matrix, the Hermitian (symmetric)
NB. log - 6-vector of boxes, test log
testhepow=: 3 : 0
NB. use for a while the definition from ggevlxx application notes
heevlv=. (9 o. {.)&.>`((%%:@diag@(mp ct))&.>)"0@ggevlvn@(,: idmat@c)
try.
'v R'=. heevlv y NB. eigendecomposition
catch.
v=. R=. _.
end.
log=. ('hepow' tdyad ((5 7"_)`]`]`heconi`nan`nan)) v ; ct R
)
NB. ---------------------------------------------------------
NB. testpow
NB.
NB. Description:
NB. Adv. to make verb to test xxpow by matrix of
NB. generator and shape given
NB.
NB. Syntax:
NB. log=. (mkmat testpow) (m,n)
NB. where
NB. mkmat - monad to generate a matrix; is called as:
NB. mat=. mkmat (m,n)
NB. (m,n) - 2-vector of integers, the shape of matrix mat
NB. log - 6-vector of boxes, test log
NB.
NB. Application:
NB. - test by random square real matrix with elements
NB. distributed uniformly with support (0,1):
NB. log=. ?@$&0 testpow_mt_ 150 150
NB. - test by random square real matrix with elements with
NB. limited value's amplitude:
NB. log=. _1 1 0 4 _6 4&gemat_mt_ testpow_mt_ 150 150
NB. - test by random square complex matrix:
NB. log=. (gemat_mt_ j. gemat_mt_) testpow_mt_ 150 150
testpow=: 1 : 'nolog_mt_`(lcat_mt_@(testgepow_mt_@u`(testdipow_mt_@(u dimat_mt_ u))`(testhepow_mt_@(u hemat_mt_))`:0))@.(=/)'