684 - Integral Determinant

[lowai]

What's wrong with this code?
It stores fractions to avoid float computation.
[cpp]
#include <stdio.h>
#include <math.h>

#define max_matrix_size 31

typedef struct {
long long a, b;
} tfraction;

typedef struct {
int size;
tfraction matrix[max_matrix_size][max_matrix_size];
} tmatrix;

tmatrix res;

long gcd(long long a, long long b)
{
return (a%b==0) ? b : gcd(b, a%b);
}

tfraction fracdiv(tfraction frac1, tfraction frac2)
{
tfraction t;
t.a = frac1.a * frac2.b;
t.b = frac1.b * frac2.a;
long long k = gcd(t.a, t.b);
t.a /= k;
t.b /= k;
return t;
}

tfraction fracadd(tfraction frac1, tfraction frac2)
{
tfraction t;
t.b = frac1.b * frac2.b;
t.a = frac1.a * frac2.b + frac2.a * frac1.b;
long long k = gcd(t.a, t.b);
t.a /= k;
t.b /= k;
return t;
}

tfraction fracmulti(tfraction frac1, tfraction frac2)
{
tfraction t;
t.a = frac1.a * frac2.a;
t.b = frac1.b * frac2.b;
long long k = gcd(t.a, t.b);
t.a /= k;
t.b /= k;
return t;
}

void swapline(int ln1, int ln2)
{
int i;
tfraction t;
for (i = 1; i <= res.size; i++) {
t = res.matrix[ln1];
res.matrix[ln1] = res.matrix[ln2];
res.matrix[ln2] = t;
}
}

tfraction det()
{
tfraction d, mik;
int i, j, k, ik;
float max;
tfraction tmp;

d.a = 1;
d.b = 1;
for (k = 1; k <= res.size - 1; k++) {
max = -1e10;
for (i = k; i <= res.size; i++)
if (res.matrix[k].a != 0)
if (res.matrix[k].a / res.matrix[k].b > max) {
max = res.matrix[k].a / res.matrix[k].b;
ik = i;
}
if (res.matrix[ik][k].a == 0) {
d.a = 0;
return d; /* the matrix is singular */
}
if (ik != k) {
swapline(ik, k);
d.a = -d.a;
}
d = fracmulti(d, res.matrix[k][k]);
for (i = k + 1; i <= res.size; i++) {
mik = fracdiv(res.matrix[k], res.matrix[k][k]);
res.matrix[i][k] = mik;
for (j = k; j <= res.size; j++) {
tmp = fracmulti(mik, res.matrix[k][j]);
tmp.a = -tmp.a;
res.matrix[i][j] = fracadd(tmp, res.matrix[i][j]);
}
}
}
d = fracmulti(d, res.matrix[res.size][res.size]);
return d;
}

int main()
{
int i, j;
long k;
tfraction detans;

while (scanf("%d", &res.size) == 1) {
if (res.size == 0) break;
for (i = 1; i <= res.size; i++)
for (j = 1; j <= res.size; j++) {
scanf("%ld", &k);
res.matrix[i][j].a = k;
res.matrix[i][j].b = 1;
}
detans = det();
printf("%lld\n", detans.a / detans.b);
}
printf("*\n");
return 0;
}
[/cpp][/cpp]

[Dominik Michniewski]

Could anyone help me and post some special cases to this problem ?
I get WA or RTE and I don't know , what I'm doing wrong ....

I test my data with a few matrixes which contains zeros at arbitrary posiotions and so on ....

Thanks for help
Dominik Michniewski

[Caesum]

I'd be interested to hear how people have solved this one - is floating point the way to go? I am currently using a fairly complicated integer algorithm that includes computation of egcd's to speed up the column additions as I can reduce two columns headed by x y to d 0 from this where d=gcd(x,y) and it runs very fast on 30*30 matrices but it always comes back as TLE on the judge. A quick play with floating point gave a wa which is no surprise when dealing with determinants and floating point errors. If I can get my program running faster then will integer overflow be a problem as well ? Or have people implemented fractional algorithms ?

[Ivan Golubev]

I've got accepted with floating point and actually I think that it's very hard to get precision error with floating point calculations in this problem!

[Dominik Michniewski]

I use fractional algorithm, but I have some mistake within it I dont know where (I know in which loop, but I don't find where it's ...)

I don't know - my algorithm run with time less than 0.5 sec before crash ... maybe I got TLE too ... but I don't think so ... when I send to judge my program without one loop (this in error) program end solving problem with 0.25 sec ....

Best regards
Dominik

PS. Some special cases - someone could help me ? I try to make floating version of this problem ) Thanks Ivan for suggestion

[Ivan Golubev]

It's a bit hard to create correct IO for this problem (without final 2^31 overflow), so only one small test:
Input:

Code: Select all

30
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 9 0 0 0 0 0 0 0 0 0 0 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 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 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 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
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 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 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 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 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 3 4 2 1 2 0 1 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 9 5 4 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 2 0 3 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 7 5 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 3 0 1 0 0 0 7 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 2 0 0 0
0 0 0 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 0 0 0 0 0 0 0 1 0 0 1 0 0 0 7 0
0 0 0 0 0 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 0 0 0 0 0 2 0 0 0 0 1 0 3 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 1 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 5 1 2 9 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 7 4 1

Answer is -1725324300.

My algorithm is straightforward -- floating point and recursion. Also a small adjustment before outputting: [c] double d = getdet(0, n);
if (d < 0)
d -= 0.5;
else
d += 0.5;
printf("%d\n", (int)(d));[/c]

[Adrian Kuegel]

I solved this problem with fractions and long long (but with a very bad runtime and a lot of WA and TLE before success). I heared from my teammate that it is solvable with double by using the strategie to select as pivot element always the maximum value.

[Caesum]

My current version of this uses doubles and rounds as Ivan suggests:
[c]
if(det>=0)
det=floor(det+0.5);
else
det=-floor(-det+0.5);
printf("%.0lf\n",det);
[/c]

It is however WA.
I have tried several strategies for the reduction method the current one uses:
1. pick a row, normally the first
2. Ensure each column has positive value in that row (ie multiply column by -1 and det by -1 if it is -ve).
3. Choose the smallest value in the row (I have tried the largest too).
4. For each column subtract q times the column from 3, so that each value in the row becomes zero.
5. Then exchange rows and columns to get the remaining value to the bottom right and multiply det by it, and continue with smaller matrix.

It passes the tests that I have come up with, which include 30*30 random matrices tested in Excel, and various sizes in between, the matrix Ivan gives above, etc.

My first version of the program just used ints and an extended gcd algorithm to reduce the values but it was too slow...... (isnt the point of this question to do that ???). I'm not sure about using a fractional algorithm because you may get some large intermediate values in numerators/denominators, and Ive tried this with double values in several different ways without success.

I think my idea is now to scrap it and do it again in a couple of months, as this has worked with other problems

[Ivan Golubev]

Caesum, (a stupid question but sometimes it happens) do you print asterisk after last line of output?

[Caesum]

yes

[Ivan Golubev]

OK, I'll try to describe my algorithm. Main idea -- to leave coefficients only in main diagonal and to zero upper triangle (sorry, don't know how this triangle must be named correctly). So, for given matrix above bottom-right part of the matrix after calculations will be:

Code: Select all

0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
0.00  9.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
0.00  1.00  1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
0.00  7.00  0.00  1.11  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
1.00  3.00  0.00 -1.67 -8.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
0.00  0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
0.00  0.00  0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00  0.00  0.00 
0.00  0.00  0.00  0.00  1.00 -0.69 -0.25  0.88  0.00  0.00  0.00  0.00  0.00 
0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00 
0.00  0.00  0.00  0.00  2.00 -1.38 -0.50 -0.25  1.43  4.43  0.00  0.00  0.00 
0.00  0.00  0.00  0.00  0.00  0.00  0.00  2.00 -1.43 -3.43  1.06  0.00  0.00 
0.00  0.00  0.00  0.00  1.00 -0.69 -0.25 -0.12  5.71  2.71  2.03 28.09  0.00 
0.00  0.00  0.00  0.00  0.00  4.00  0.00  0.00  0.00 -8.00  6.48 79.91 -4.69 

Algorithm:
1. For each row:
1.1 Find first column (starting from 'current row number') which doesn't has zero (i.e. fabs(matrix[row][column] > 1e-8).
1.2 If there no such column then determinant == 0.
1.3 If column number not equal row number swap this column with column at 'current row number' and remember that we're need to switch sign of result.
1.4 For each column...
1.4.1 ...greater than 'row number' find coefficient that will produce zero at matrix[row][current column] after multiply this coefficient with matrix[row][row] and substract result from matrix[row][current column].
1.4.2 Adjust all remaining values in current column.

[Den]

My task is "Integral Determinant". My number is 684.

I send my program but I have " complier error" . Why??? I have done all as said at "How to submit Pascal programs"… and have a compile error. I know that my program is work … may be not effectively… but it's work … and I confirm it's correctly counting determinants of all matrixes…

So if my program will not accepted I will have been banished from University .. if it will have been happened I will have only one way … way to Russian Army… way to Hell where I have finished my days..

May be I wrong and it's not correctly counting determinants … so show me the way how to make it correctly …. If you know how to dispose from this stupid error .. so show me how I can do it… if you won't people see your comments or program's code… send me to mail pushoc@freemail.ru
I put down my program below..

You know that in the web we exist without skin color, without nationality, without religious bias.. why you should not help me if you can ??? I'm just one alone person living in 9 billion's world. Is there somebody who help me and don't turn his back .. Believe me my life in your hands…

Sincerely Den.

program PXXX(input, output); { ISO states this line in this way }
{$IFNDEF ONLINE_JUDGE}
{$M 65500, 0 , 400000}
{$ENDIF}

const size=30; {The size of the matrix}
type Matrix=array[1..size,1..size] of integer;

var n:integer;
A:Matrix;

{//////////////////////////////////////////////////////
Input a matrix
//////////////////////////////////////////////////////}

function cut_number(search:string;var k:integer):integer;
var j,abba,Code:integer;
str:string;
begin

j:=1;

while search[k]=' ' do k:=k+1;
while (search[k]<>' ')and(k<>length(search)+1) do
begin
str[j]:=search[k];
j:=j+1;
k:=k+1;
end;

str[0]:=chr(j-1);
Val(str,abba,Code);
cut_number:=abba;

end;

procedure read_string(var f:text;i:integer;n:integer;var A:Matrix);
var s:string;
j,k:integer;

begin
k:=1;
readln(f,s);
for j:=1 to n
do A[i,j]:=cut_number(s,k);
end;

procedure read_matrix(var data:text;var n:integer;var A:Matrix);
var s:string;
k:integer;
i:integer;
begin

k:=1;
readln(data,s);
n:=cut_number(s,k);

if n=0 then exit;

for i:=1 to n do
begin
read_string(data,i,n,A);
end;

end;

{//////////////////////////////////////////////////////
Calculating the determinant of the matrix
///////////////////////////////////////////////////////}

{Function of stopping a processing}

Function finish(A:Matrix;n:integer):boolean;
var i,sum:integer;

begin

sum:=0;
for i:=1 to n do
if A[1,i]<>0 then sum:=sum+1;

if (sum=1)or(sum=0) then finish:=true else finish:=false;

end;

{Exchange the rows}

procedure Exchange(x,y:integer; var A:Matrix;n:integer);
var i,buffer:integer;
begin

for i:=1 to n do
begin

buffer:=A[i,x];
A[i,x]:=A[i,y];
A[i,y]:=buffer;

end;

end;

{Sort the first string of the matrix using buble sort}

Procedure Sort(var A:Matrix;n:integer;var s:integer);
var i,j:integer;
begin
for i:=2 to n do
for j:=n downto i do
if abs(A[1,j-1])>abs(A[1,j]) then
begin
Exchange(j-1,j,A,n);
s:=-s;
end;

end;

{Get sign of a number}

function sign(a:integer):integer;
begin

if a>0 then sign:=1 else
if a=0 then sign:=0 else
if a<0 then sign:=-1;

end;

{Subtract the rows using mathematics rules}

Procedure Subtract(x,y:integer;var A:Matrix;n:integer);
var i:integer;
c:integer;

begin

c:=sign(A[1,y])*sign(A[1,x])*abs(A[1,y] div A[1,x]);

for i:=1 to n do
begin
A[i,y]:=A[i,y]-c*A[i,x] ;
end;

end;

{Process the matrix}

Procedure Processing(var A:Matrix;n:integer);
var i,j:integer;

begin

j:=1;
while(A[1,j]=0)
do j:=j+1;

for i:=j+1 to n do Subtract(j,i,A,n);

end;

{Make current minor from the matrix}

procedure GetMatrix(a:Matrix; var b:Matrix; m,i,j:integer);
var ki,kj,di,dj:integer;
begin
di:=0;
for ki:=1 to m-1 do
begin

if (ki=i) then di:=1;
dj:=0;

for kj:=1 to m-1 do
begin

if (kj=j) then dj:=1;
b[ki,kj]:=a[ki+di,kj+dj];

end;

end;
end;

{Main function}

Function Determinant(a:Matrix;n:integer):integer;
var i,j,d,k,d_sign:integer;
b:Matrix;
begin

{set raw data}
d:=0;
k:=1;
d_sign:=1;

{calculating}

if (n=1) then d:=a[1,1]
else
if (n=2) then d:=a[1,1]*a[2,2]-a[2,1]*a[1,2]
else { n>2 }
begin

while (finish(a,n)=false) do
begin
Sort(A,n,d_sign);
Processing(A,n);
end;

for i:=1 to n do
begin
if a[1,i]<>0 then
begin
GetMatrix(a,b,n,1,i);
d:=d+k*a[1,i]*Determinant(b,n-1);

end;
k:=-k;
end;
end;

{last breath}
Determinant:=d_sign*d;
end;


begin

while not eof(input) do begin
if eoln(input) then begin{ this should follow the while sentence with nothing between, just like this }
readln(input); { skip rest of line }
{writeln(output);} { reproduce the eoln character }
continue; { re-test for EOF }
end;
read_matrix(input,n,A);
if n=0 then exit;
writeln(output,Determinant(A,n));

end;

end.

[little joey]

Well, I haven't solved this problem yet, but I can give you some hints:

The judge always gives you a reason for a compiler error, so if you read that and follow the hints it gives, you can resolve most of these.
In this case it objects against setting the stringlength by accessing the zeroth element and advises you to use the setlength procedure instead:
change [pascal]str[0]:=chr(j-1);[/pascal]to[pascal]setlength(str,j-1);[/pascal].

Now you've got only Wrong Answer to deal with
One obvious mistake is that you don't print a star after the last line of output. But even correcting this, won't get you Accepted.

As I said, I haven't solved the problem, so this is all the help I can give.

[Den]

little joey thanks alot... may be you are litlle
...but kind and helpful ...

Thanks.Den.

[friggstad]

I tried solving this using a modular approach. Lets say that that the answer is 'a'. If we have a prime 'p' where we know that -p/2 < a < p/2 and we compute the answer mod p, we can recover the integer answer by simply printing the residue of the answer mod p that is in the range -p/2...p/2. I computed the determinant mod p by the standard method using row operations, except the arithmetic was done mod p.

This method works because p is large enough so that no two different answers can map to the same residue. Also, the determinant is simply an expression using *, + and - so working mod p and recovering the answer is allowed by the natural ring homomorphism from the integers to the residues mod p.

We are guaranteed no overflow in the answer, so the answer is in the range of [2^31..2^31) and an appropriate prime should be at least 2^32. However, I took a bit of a gamble and used the prime just previous to 2^32 (to avoid overflow in 64 bits when multiplying) and it still worked.

This involves a bit more coding than using doubles, but the approach feels more "pure" to me than relying on floating point arithmetic.