/*********************************************************************** cramer4.c (June,2000) We want to study the problem of solving the system below, later referred to as system (1): [ m0 + u00 u01 u02 u03] [a0] [-1] [ u10 m1 + u11 u12 u13] [a1] = [-1] [ u20 u21 m2 + u22 u23] [a2] [-1] [ u30 u31 u32 m3 + u33] [a3] [-1] Most of the uij are equal to -1, but in each row we choose one uij = 0. Further restrictions: m1 = 2, 3, ... 10 m2 = 2, 3, 4 m3 = 1, 2 m0 is unknown. The coefficient matrix of system (1) will be denoted as C in the comments below. All the numbers a0, a1, a2 and a3 are to be positive integers and any three of them should be relatively prime. We solve by Cramer's Rule...Let d_i denote the determinant obtained by replacing column i of C by -1's. The determinant of C depends on m0 and we get things like a0 = d0/(a*m0 + b). The fact that we want this to be a positive integer puts restrictions on m0. We then let m0 run through all possible values and look for 4-tuples [a0,a1,a2,a3] that satisfy the system (1), are positive integers, and any three are relatively prime. ***************************************************************************/ #include #include main() { int i,j; int M[4][4]; /* diagonal matrix with m_i = M[i][i] */ int m[4]; int U[4][4]; /* matrix whose entries are mostly -1 */ int z0, z1, z2, z3; /* where to put 0 in ith row of U */ /* if j(i) = k then zi = k */ int C[4][4]; /* the coefficient matrix of system (1) */ int R[4][4]; /* a copy of C for row reduction */ int badmatrix; /* don't let C[3][3] = 0 */ int nocondition; /* if det C does not depend on m0 */ int D, a, b; /* det C = D = a*m0 + b */ int d0, d1, d2, d3; /* det's from Cramer's rule */ int S[3][3]; /* a 3x3 matrix, usually a subdet of di */ int signchange; /* did we switch two rows in reduction? */ int GiveUp; /* Is there a problem with this 4tuple? */ int alpha[4]; /* Series case (a = 0, d0 nonzero) */ int shift[4]; /* */ int multiple; /* is shift a multiple of (d1,d2,d3)? */ int A[4]; /* candidate 4-tuple */ int B[4]; /* same as A, but throw away signs */ int t1, t2, t3, t4; /* gcd of three of the A[i]'s */ int absb, L; /* how far do we let m0 run? */ int temp; /* used for sorting --- write the ai's */ /* in increasing order */ int gcd(int a, int b){ int r; if (a > b){ r = a; a = b; b = r; } r = b%a; while (r > 0){ b = a; a = r; r = b%a; } return a; } int det( int S[3][3] ){ int s1,s2,s3; s1 = S[0][0]*(S[1][1]*S[2][2] - S[1][2]*S[2][1]); s2 = S[1][0]*(S[0][2]*S[2][1] - S[0][1]*S[2][2]); s3 = S[2][0]*(S[0][1]*S[1][2] - S[0][2]*S[1][1]); return s1+s2+s3; } FILE *fp; FILE *mfp; FILE *nfp; FILE *kfp; /****************************************** Initialize the matrix M: *******************************************/ for(i = 0; i<4; i++){ for(j=0; j<4; j++){ M[i][j] = 0; } } /************************************************************** All solutions to the system (1) where det C is nonzero and where d_0 (replace col 0 of C by -1's and take det) is also nonzero. If both d_0 and det C are zero then we have to check it out. These examples are written in the file "zerozero" Normally, det C depends on the unknown m0 --- det C will be a*m0 + b. If a = 0, then these examples will be written to the file "series" and we expect these to lead to infinite series of solutions to (1). They may also produce individual solutions that do not appear elsewhere. Solutions in this case which do not fall into an infinite series will be printed to the file "multiples". **************************************************************/ fp = fopen("ordered4tuples", "w"); mfp = fopen("zerozero", "w"); nfp = fopen("series","w"); fprintf(nfp, "If det C does not depend on m0 then we get: \n\n"); kfp = fopen("multiples","w"); /************************************** Let m[3] run from 1 to 2 Let m[2] run from 2 to 4. Let m[1] run from 2 to 10. In each case m[i] should equal M[i,i] **************************************/ for (m[3]=1; m[3]<3; m[3]++){ M[3][3] = m[3]; for(m[2]=2; m[2]<5; m[2]++){ M[2][2]=m[2]; for(m[1]=2; m[1]<11; m[1]++){ M[1][1] = m[1]; /* at this point we have chosen a matrix M */ /* now we need to define the entries of U */ /* Start by defining U to have all entries */ /* equal to -1. Then choose one entry in */ /* each row to be 0. */ for(z0=0;z0<4; z0++){ for(z1=0;z1<4; z1++){ for(z2=0; z2<4; z2++){ for(z3=0; z3<4; z3++){ for(i=0; i<4; i++){ for(j=0; j<4; j++) { U[i][j] = -1; } } /* Reinitialize U to be all -1's */ U[0][z0] = 0; /* put a 0 in the 0th row */ U[1][z1] = 0; /* 1st */ U[2][z2] = 0; /* 2nd */ U[3][z3] = 0; /* 3rd */ /***************************************** Now the coefficient matrix C of the system (1) can be formed by adding M and U. Save a copy as R *********************************************/ for(i=0; i<4; i++){ for(j=0; j < 4; j++){ C[i][j] = M[i][j] + U[i][j]; R[i][j] = C[i][j]; } } /*There is one bad case. It cannot happen */ /* in our problem that the 3,3 entry of C */ /* will be zero. This would correspond to */ /* choosing m3 = 1 and NOT having j = 3 for */ /* i = 3. But this would mean that two of */ /* the a_i's sum to 1, an impossibility */ badmatrix = 0; if(C[3][3] == 0) badmatrix = 1; if(badmatrix == 0){ /********************************************************* Now we want to solve system (1) using Cramer's rule. Replace column 0 of C by -1's and do three row operations to get (-1,0,0,0) in column 0. Let S be the 3x3 submatrix obtained by deleting row 0 and column 0 of the resulting matrix. Then d_0 will be (-1)*det S. ***********************************************************/ for(i=1; i<4; i++){ for(j=1; j<4; j++){ S[i-1][j-1] = C[i][j] - C[0][j]; } } d0 = (-1)* det(S); /**************************************************** Next we want to compute the determinant of C, where the 0,0 entry is the unknown m0 + U[0][0]. This will be of the form a*m0 + b where a is the det of the 3x3 submatrix obtained by deleting col 0 and row 0 in the coefficient matrix C. The other term will be more complicated and its computation is explained below. First we compute a. The U[i][j] are 0 or -1. So the first column has 0's and/or -1's in rows 1,2,3. We want to do row operations to get (*,-1,0,0) or (*,0 0 0) in column 0. If we have a -1 in row 1, then use it to get 0's below. If not, we we look at row 2. If there is a -1 then we switch row 1 and row 2 (changing the sign of the determinant) and row reduce. If there is a 0 in rows 1 and 2 but a -1 in row 3 then we switch rows 1 and 3. etc. *****************************************************/ signchange = 1; /* Assume that we don't need to switch rows */ if(R[1][0]==-1){ if(R[2][0]==-1){ for(j=0; j<4; j++){ R[2][j]=C[2][j] - C[1][j]; } } if(R[3][0]==-1){ for(j=0; j<4; j++){ R[3][j] = C[3][j] - C[1][j]; } } } /* end of the case where column 0 started as (*,-1,*,*) */ /* At this point column 0 looks like (*,-1,0,0) or we started with (*,0,*,*) in column 0. In that case, if we have (*,0,-1,*) we need to switch row 1 and row 2 and make sure we have a 0 in row 3 */ if(R[2][0]==-1){ for(j=0; j<4; j++){ R[1][j] = C[2][j]; R[2][j] = C[1][j]; signchange = -1; if(C[3][0] ==-1) R[3][j] = C[3][j] - C[2][j]; } } /* end of the case where we start with (*,0,-1,0) */ /* Again at this point column 0 looks like (*,-1,0,0) or we started with (*,0,0,*) in column 0. If the last entry is -1 then we need to switch rows 1 and 3. If the last entry is 0, then we are done. */ if(R[3][0]==-1){ for(j=0; j<4; j++){ R[1][j] = C[3][j]; R[3][j] = C[1][j]; signchange = -1; } } /* Now pull out the 3x3 submatrix of this reduced matrix consisting of cols 1,2,3 and rows 1,2,3 Take the determinant and adjust the sign if necessary. det C will be of the form a*m_0 + b where m_0 is unknown. */ for(i = 1; i<4; i++){ for(j=1; j<4; j++){ S[i-1][j-1]= R[i][j]; } } a = signchange*det(S); /************************************************************ Now det C is of the form a*m0 + b. We need to compute b. If we are in the case where the reduced matrix had (*,0,0,0) as its column 0, then b is simple. We just get b = U[0][0]*a. Otherwise, first column 0 looks like (*,-1,0,0) and we have to compute the subdeterminant we get by deleting row 1 and column 0 and add that to U[0][0]*a. *************************************************************/ if(R[1][0] == 0){ b = U[0][0]*a; } else { for(j=1; j<4; j++){ S[0][j-1] = R[0][j]; S[1][j-1] = R[2][j]; S[2][j-1] = R[3][j]; } b = U[0][0]*a + signchange*det(S); } /********************************************************** Now we have det C = a*m0 + b. Since a0 = d0/det C = d0/(a*m0 + b) and a0 is to be a positive integer we get restrictions on how m0 can be chosen since both a0 and m0 must be integers larger than or equal to 1. There are three cases to consider. Case 1: If d0 is zero, then the system will give no 4-tuples unless we choose m0 so that det C is also 0, that is we must choose m0 to be -b/a. So we check if -b/a is a positive integer. If not, then we will get no solutions to this system and we should move on. Otherwise, if -b/a is a integer, then we should set m0 equal to this value and save the matrix C (where now the 00 entry is m0 + u00) for later consideration. We may be in a degenerate case and we'll need to think about it more carefully. ************************************************************/ if(d0 == 0 && a != 0){ if (-b%a == 0){ /* otherwise a0 = 0 and that is no good for us */ if( -b/a > 0){ m[0] = -b/a; C[0][0] = C[0][0] + m[0]; /* mfp points to the file "zerozero" */ fprintf(mfp, "d0 = 0 and det C = 0 if\n"); fprintf(mfp, "[m0, m1, m2, m3] = [%5d, %5d, %5d, %5d]\n", m[0], m[1], m[2], m[3]); for(i=0; i<4; i++){ for(j=0; j<4; j++) fprintf(mfp, "%5d", C[i][j]); fprintf(mfp, "\n"); } /* end of for */ } /* end of if -b/a is postive */ } /* end of if -b/a is an integer */ } /* end of the case where d0 is zero but a is not. */ /************************************************************** Case 2: If a = 0 then there is no condition on m[0]. We expect this case to lead to the infinite series. Now det C = b. Compute d0. If d0/b is not a positive integer, then give up. Otherwise, compute d1,d2 and d3. These depend on m0. Let d_i = m0*alpha_i + beta_i. The beta[i] don't matter for us but the alpha_i do. For m0 between 1 and |b| we compute d1,d2,d3. If di/b is an integer for every i, then a[i] = m0*alpha[i] + beta[i] ---------------------- b is an integer and when we increase m0 by |b| we will increase a[i] by (|b|*alpha[i]/b). Let shift[i] denote this quantity. For m0 between 1 and |b|, find all 4-tuples [a0,a1,a2,a3] = [d0/b,d1/b,d2/b,d3/b] such that d0/b will be a positive integer and so that any three are relatively prime. Then [a0,a1,a2,a3] = [ a0, a1 + k*shift[1], a2 + k*shift[2], a3 + k*shift[3]] will be another solution in the same series. If the shift vector is a multiple of the vector [a1,a2,a3] then we will not get a series. In that case, for k > 0, we will get 4-tuples that fail the relatively prime condition. However the k=0 case may be a new solution. These will be written to the file "multiples". ****************************************************************/ if(a == 0) nocondition = 1; else nocondition = 0; if(nocondition == 1){ if(b != 0){ if( b < 0 ) absb = -1*b; else absb = b; for( m[0] = 1; m[0] <= absb; m[0]++){ /* initialize A */ for(i=0; i<4; i++) A[i] = 1; GiveUp = 0; D = b; if(d0 % D != 0) GiveUp = 1; else A[0] = d0/D; /***************************************************************** Compute alpha[1], shift[1] and d1 ******************************************************************/ for(i=0; i<4; i++) for(j=0; j<4; j++) R[i][j] = C[i][j]; /* Reinitialize R */ R[0][0] = m[0] + U[0][0]; for(i=0; i<4; i++) R[i][1] = -1; /* Set col 1 to be all -1 */ /* Now alpha[1] will be the determinant of the 3x3 matrix obtained by deleting column 1 and row 0 of the current R */ for(i=1; i<4; i++){ for(j=1; j<4; j++){ S[i-1][j-1] = R[i][j]; } } alpha[1] = det(S); shift[1] = absb*alpha[1]/b; for(i = 1; i < 4; i++){ /* subtract row 0 from rows 1,2,3 */ for(j=0; j<4; j++){ R[i][j] = R[i][j] - R[0][j]; } } for(i=1; i<4; i++) S[i-1][0] = R[i][0]; for(i=1; i<4; i++) for(j=2; j<4; j++) S[i-1][j-1] = R[i][j]; d1 = det(S); if(d1 == 0) GiveUp = 1; if(d1 % D != 0) GiveUp = 1; else A[1] = d1/D; /********************************************************* Compute alpha[2], d2 and shift[2] **********************************************************/ for(i=0; i<4; i++) for(j=0; j<4; j++) R[i][j] = C[i][j]; /* Reinitialize R */ for(i=0; i<4; i++) R[i][2] = -1; /* Set col 2 to be all -1 */ R[0][0] = m[0] + U[0][0]; for(i=1; i<4; i++){ for(j=1; j<4; j++){ S[i-1][j-1] = R[i][j]; } } alpha[2] = det(S); shift[2] = absb*alpha[2]/b; for(i = 1; i < 4; i++){ /* Subtract row 0 from rows 1,2,3 */ for(j=0; j<4; j++){ R[i][j] = R[i][j] - R[0][j]; } } for(i=1; i<4; i++){ for(j=0; j<2; j++) S[i-1][j] = R[i][j]; S[i-1][2] = R[i][3]; } d2 = -1*det(S); if(d2 == 0) GiveUp = 1; if(d2 % D != 0) GiveUp = 1; else A[2] = d2/D; /********************************************************* Compute alpha[3], d3 and shift[3] **********************************************************/ for(i=0; i<4; i++) for(j=0; j<4; j++) R[i][j] = C[i][j]; /* Reinitialize R */ R[0][0] = m[0] + U[0][0]; for(i=0; i<4; i++) R[i][3] = -1; /* Set col 3 to be all -1 */ for(i=1; i<4; i++){ for(j=1; j<4; j++){ S[i-1][j-1] = R[i][j]; } } alpha[3] = det(S); shift[3] = absb*alpha[3]/b; for(i = 1; i < 4; i++){ /* subtract row 0 from rows 1,2,3 */ for(j=0; j<4; j++){ R[i][j] = R[i][j] - R[0][j]; } } for(i=1; i<4; i++){ for(j=0; j<3; j++) S[i-1][j] = R[i][j]; } d3 = det(S); if(d3 == 0) GiveUp = 1; if(d3 % D != 0) GiveUp = 1; else A[3] = d3/D; /************************************************************* Now the candidate A is computed. All we have checked so far is that we can find integer solutions using Cramer's rule. Now we have to check the relatively prime condition. It could happen that some of the A[i]'s are negative. The function gcd needs positive integers as input. So save a copy as B, throwing away any signs *************************************************************/ if(GiveUp==0){ for(i=0; i<4; i++){ if(A[i]>0) B[i] = A[i]; else B[i] = -1*A[i]; } t1 = gcd(gcd(B[0], B[1]),B[2]); t2 = gcd(gcd(B[0], B[1]),B[3]); t3 = gcd(gcd(B[1], B[2]),B[3]); t4 = gcd(gcd(B[0], B[2]),B[3]); if(t1*t2*t3*t4 != 1) GiveUp = 1; } /************************************************************ We also need to check if [alpha_1, alpha_2, alpha_3] is a multiple of [d1/b, d2/b, d3/b]. If it is then we don't get a series, but it could be a new solution. These should be printed to the file "multiples". ************************************************************/ multiple = 0; if (alpha[1]*d2 - d1*alpha[2] == 0){ if (alpha[2]*d3 - alpha[3]*d2 == 0){ if (alpha[1]*d3 - alpha[3]*d1 == 0){ multiple = 1; } } } if(GiveUp == 0 && multiple == 0){ /* nfp points to the file "series" */ fprintf(nfp, "\n\n"); fprintf(nfp, "[%5d, %5d + %3dk, %5d + %3dk, %5d + %3dk]\n", A[0], A[1], shift[1],A[2],shift[2], A[3],shift[3]); if(shift[1] <= 0){ fprintf(nfp, "watch out\n"); } if(shift[2] <= 0){ fprintf(nfp, "watch out\n"); } if(shift[3] <= 0){ fprintf(nfp, "watch out\n"); } } if(GiveUp == 0 && multiple == 1){ /* kfp points to the file "multiples" */ fprintf(kfp, "\n\n"); fprintf(kfp,"[%5d, %5d, %5d, %5d]", A[0], A[1], A[2], A[3]); } }/* end of m[0] loop */ } /* end of if b is not 0 */ } /* end of no condition case. */ /****************************************************************** That takes care of the bad cases for the problem. Now we consider what happens when det C depends on m0 (that is, a is not 0) and d0 is not zero so we will get a nonzero value for a0. There will be a limit on the possible size of m0. If d0 is positive, then look at d0 d0 d0 ----------, ---------, ----------,........ a + b 2a + b 3a + b (m0 = 1) (m0 = 2) (m0 = 3) and so on until the denominator gets larger than d0. If any of these is a positive integer, then we have a candidate for a0. Set m0 to this value and then compute d_i/(m0*a + b). If these are also positive integers, then we have a candidate 4-tuple [a0,a1,a2,a3]. They may not be listed in increasing order and we still must check the rel. prime condition. If that is satisfied, these 4-tuples will be written to the file "ordered4tuples" in increasing order. To get started: If d0 > 0 then we need d0 > (m0*a + b), which leads to the condition d0 - b m0 < -------- and so -b < m0*a < d0 - b. a If d0 < 0 then we need d0 - b < m0*a < -b. Depending on whether a is positive or negative we get d0 - b -b -b ---------- < m0 < --- or --- < m0 < d0 - b a a a *******************************************************************/ if(d0 != 0){ if(a != 0){ L = (d0 - b)/a; if ( -b/a > L) L = -b/a; for( m[0] = 1; m[0] <= L; m[0]++){ GiveUp = 0; D = a*m[0] + b; if(D==0){ GiveUp = 1; } if(D!=0){ if(d0 % D != 0){ GiveUp = 1; } else{ A[0] = d0/D; if (A[0]<1) { GiveUp = 1; } } /*************************************************************** if this value for m0 is working so far, if it gives a legal value for a0, then we should compute a1,a2,a3. We have to adjust the matrix C to reflect the current value of m0 to compute the other determinants. *****************************************************************/ C[0][0] = m[0] + U[0][0]; /************************************** Compute d1 **************************************/ for(i=0; i<4; i++) for(j=0; j<4; j++) R[i][j] = C[i][j]; /* Reinitialize R */ for(i=0; i<4; i++) R[i][1] = -1; /* replace col 1 by -1's */ for(i = 1; i < 4; i++){ for(j=0; j<4; j++){ R[i][j] = R[i][j] - R[0][j]; /* row reduce */ } } for(i=1; i<4; i++) /* S = the 0,1 minor of R */ S[i-1][0] = R[i][0]; for(i=1; i<4; i++) for(j=2; j<4; j++) S[i-1][j-1] = R[i][j]; d1 = det(S); if(d1 % D != 0){ GiveUp = 1; } else{ A[1] = d1/D; if(A[1]<1){ GiveUp=1; } } /************************************** Compute d2 **************************************/ for(i=0; i<4; i++) for(j=0; j<4; j++) R[i][j] = C[i][j]; /* Reinitialize R */ for(i=0; i<4; i++) R[i][2] = -1; for(i = 1; i < 4; i++){ for(j=0; j<4; j++){ R[i][j] = R[i][j] - R[0][j]; } } for(i=1; i<4; i++){ for(j=0; j<2; j++) S[i-1][j] = R[i][j]; S[i-1][2] = R[i][3]; } d2 = -1*det(S); if(d2 % D != 0){ GiveUp = 1; } else{ A[2] = d2/D; if(A[2]<1){ GiveUp=1; } } /************************************** Compute d3 **************************************/ for(i=0; i<4; i++) for(j=0; j<4; j++) R[i][j] = C[i][j]; /* Reinitialize R */ for(i=0; i<4; i++) R[i][3] = -1; for(i = 1; i < 4; i++){ for(j=0; j<4; j++){ R[i][j] = R[i][j] - R[0][j]; } } for(i=1; i<4; i++){ for(j=0; j<3; j++) S[i-1][j] = R[i][j]; } d3 = det(S); if(d3 % D != 0){ GiveUp = 1; } else{ A[3] = d3/D; if(A[3]<1){ GiveUp=1; } } if(GiveUp == 0){ /* Test for relative primeness */ t1 = gcd(gcd(A[0], A[1]),A[2]); t2 = gcd(gcd(A[0], A[1]),A[3]); t3 = gcd(gcd(A[1], A[2]),A[3]); t4 = gcd(gcd(A[0], A[2]),A[3]); if(t1*t2*t3*t4 != 1){ GiveUp = 1; } } /***************************************************************** The 4-tuples produced so far will not be in increasing order. Next, we reorder the A[i]'s to get A[0] <= A[1] <= A[2] <= A[3] when we record them in the file "ordered4tuples. *******************************************************************/ if(GiveUp==0){ for(i=1; i<4; i++){ if(A[i] < A[0]){ temp = A[0]; A[0] = A[i]; A[i] = temp; } } for(i=2; i<4; i++){ if(A[i] < A[1]){ temp = A[1]; A[1] = A[i]; A[i] = temp; } } if(A[3]