Thanks for the tips. I'm not very familiar with such algorithms in linear algebra. Do you recommend any sources, related to J or otherwise, on solving such a equations, and maybe where to start on the xAx = y problem?
On 14 December 2010 03:10, Marshall Lochbaum <[email protected]> wrote: > If it helps, an algorithm for reducing a symmetric matrix is given below. It > factors matrix y into s.r.s(T), so that x(T).y.x = x(T).s.r.s(T).x = > (s(T)x)(T).r.(s(T)x) . > A further factoring is simple: L = s.(sqrt(r)) so y=L.L(T) and we have (for > scalar z) x(T).y.x = z -> (x(T).L).(x(T).L)(T) = z . > However, it does not do integer calculations. For the system Ax=y I applied > Euclid's method for GCD to the rows multiple times; I hope something similar > is possible for this system. > > NB. y is a symmetric matrix > NB. returns r;s with r diagonal, s lower triangular > NB. such that y=s.r.s(T) > symreduce=:3 :0 > if. 1 1-:$y do. y;(1 1$1) return. end. > n=. #y > I=. =/~i.n NB. identity matrix > i=. <:+:I NB. 1 on diagonal, _1 elsewhere--i*p = (%.p) for Gaussian p > > m=.n {. % (<0 0){y > NB. p*i = %.p is a gaussian matrix so that pi.y.pi(T) has only > NB. a diagonal element in the first row and column > p=.(I * -. n{.1) + m *"1 y > y=.pi +/ .* y +/ .* |:pi=.p*i > 'r0 s0'=.symreduce 1 1}.y > > r=.(n{. (<0 0){y) ,. 0, r0 NB. diagonal matrix r > s=.({."1 p) ,. 0, s0 NB. lower triangular s > r;s > ) > > Marshall > > -----Original Message----- > From: [email protected] > [mailto:[email protected]] On Behalf Of Justin Paston-Cooper > Sent: Sunday, December 12, 2010 1:14 PM > To: Programming forum > Subject: Re: [Jprogramming] Solving diophantine equations > > Thanks for the code. I'm getting a syntax error on backsub, though. > What should the original definition of r be? Do you have any tips for > solving xAx = y? > > On 12 December 2010 15:45, Marshall Lochbaum > <[email protected]> wrote: >> A working back-substitution routine is: >> >> backsub=:3 :0 >> n=.# (-. 0...@{.) y >> firstind=. i.&1@:~:&0 >> for_i. i.&.<: n do. >> 'ind pivot'=.(firstind ([ , {) ]) i{r >> mult=.n {. - pivot <....@%~ ind ([ {"1 {.) r >> r=. (mult */ i{r) + r >> end. >> ) >> >> backsub@:rowreduceint puts an integer matrix A into an upper >> triangular form so that: >> Zero rows are at the bottom. >> The first element of each row is positive, and these first values do >> not decrease for later rows. >> If a row has first element e, all other numbers in that column are >> zero if they are below e and between zero and (e-1) if they are above it. >> >> This is about the closest you get to reduced row echelon for integer >> matrices. >> >> I can't really post representative timings because for random matrices >> using this process the entries become giant (ie. 10^173 for 100 1...@$10). >> >> Marshall >> >> -----Original Message----- >> From: Marshall Lochbaum [mailto:[email protected]] >> Sent: Sunday, December 12, 2010 10:17 AM >> To: 'Programming forum' >> Subject: RE: [Jprogramming] Solving diophantine equations >> >> A while ago I made a integer row-reduction routine that does the >> equivalent of forward elimination for integer matrices: >> >> rowreduceint=:3 :0 >> if. 0=*/$y do. y return. end. >> y=.(* <:@+:@(>&0)@{."1) y >> >> n=.+/*{."1 y >> if. n=0 do. 0,.rowreduceint }."1 y return. end. >> while. n>1 do. >> i=.(i.<./) (+ _*=&0) {."1 y >> y=.y - (i{y)*/~ 0 i} (<....@% i&{) {."1 y >> n=.+/*{."1 y >> end. >> y=. \:~ y >> ({.y) , 0,. rowreduceint 1 1}. y >> ) >> >> You could run this on the extended matrix A|y as a start to solving >> the integer equation Ax=y . Then you would have to find a good way to >> do back-substitution, because it will not work perfectly with integer > matrices. >> >> Marshall >> >> >> -----Original Message----- >> From: [email protected] >> [mailto:[email protected]] On Behalf Of Justin >> Paston-Cooper >> Sent: Sunday, December 12, 2010 3:18 AM >> To: Programming forum >> Subject: [Jprogramming] Solving diophantine equations >> >> Hello, >> >> Could anyone please tell me what facilities there are for J for >> solving diophantine equations? Specifically, I have square integer >> matrices A of rank n, and I want to find integer vectors x of length n >> such that xA = y (another integer vector) and x.x = z (an integer). >> I'll also have cases where it is useful to find scalars z with xAx = z. >> Have such packages been written for J? Should I be looking somewhere else? >> I've been working with Mathematica, and it hasn't been particularly fast. >> >> Thanks in advance, >> >> Justin >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forums.htm >> >> ---------------------------------------------------------------------- >> For information about J forums see http://www.jsoftware.com/forums.htm >> > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
