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
