Thanks, Cliff
Yes, I also wrote a function some years ago, "m_gauss_jordan", merely
converted from a
gauss-jordan I found somewhere - sorry not to attribute it... I
couldn't locate it under
addons/math, but perhaps it's there.
I think the only changes were to replace the primitives * % + - and (+/
. *) by their modular
equivalents, and not to worry about checking for tolerance as
everything's finite integer.
m_gauss_jordan=: 3 : 0
0 m_gauss_jordan y
:
M =. x
NB. modular equivalents of %, *, +, -, and (+/ . *)
div =. M mdivide
by =. M mtimes
add =. M madd
sub =. add -
mby =. add/ . by
m =. y
'r c'=. $m
rws =. i.r
i =. j=. 0
max =. i.>./
while. (i<r) *. j<c do.
k=. max col=. | i}. j{"1 m
col=. | i}. j{"1 m
NB. if all col < tol, set to 0: original comment in "standard" gauss-jordan
if. 0 = k{col do. NB. no need for tolerance here!
m=. 0 (<(i}.rws);j) } m
NB. otherwise sort and pivot:
else.
if. k do.
m=. (<i,i+k) C. m
end.
col =. j{"1 m
m =. m sub (col sub i=i.#m) by"0/ (i{m) div i{col
i =. >:i
end.
j=. >:j
end.
m
)
It can do some of your Md calculation:
a
4 1 3
1 0 3
1 1 1
[aa =: (,.=@i.@#) a
4 1 3 1 0 0
1 0 3 0 1 0
1 1 1 0 0 1
[ai =: 3}."1 ] 17 m_gauss_jordan aa
15 7 2
7 12 11
12 15 5
It doesn't manage to work with your A =: a, 1 , though:
17 m_gauss_jordan A,.=i.4
1 0 0 15 7 0 2
0 1 0 7 12 0 11
0 0 1 12 15 0 5
0 0 0 0 0 1 16
Interesting,
Mike
On 05/04/2023 20:09, Clifford Reiter wrote:
load modular_matrix_divide.ijs
open script for function/adverb definitions
gcd's, gauss-jordan, modulo divide and modular matrix divide
I think this is a draft of the modular matrix system solver that I imagined
a couple decades ago. It may still have some rough edges/mistakes.
gcd2x 51 119
17 _2 1
_2 1 +/ . * 51 119
17
NB. more general than integer
gcd2x 53r3 20
1r3 17 _15
NB. md is meant to be like % mod m
7 md 6
6
m md i.&.<:m=:7
1 4 5 2 3 6
m md i.&.<:m=:8
1 _ 3 _ 5 _ 7
5 m md i.&.<:m
5 _ 7 _ 1 _ 3
NB. Md is meant to be like %. mod m
-/ . * a=:?.3 3$5
_7
a
4 1 3
1 0 3
1 1 1
]ai=: 17 Md a
15 7 2
7 12 11
12 15 5
17|ai mp a
1 0 0
0 1 0
0 0 1
1 2 3 (17) Md a
1 13 6
NB. Non square case includes:
NB. Can take an independent direction b and
NB. make one, B, orthogonal to the columns of A --
NB. as in Gram-Schmidt even if not "least squares fit"
$A=:a,1
4 3
17 Md A NB. Left inverse
15 7 1 1
7 12 14 14
12 15 11 11
17|(17 Md A)mp A
1 0 0
0 1 0
0 0 1
b=:1 2 3 4
s=:b 17 Md A
17| (B=:17|b-A mp s) mp A
0 0 0
B
0 0 8 9
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
