// Function: To calculate the number of distinct sequences defined by the elements of powers >= 0 of n X n (0,1) matrices. // // Developed by Chris Gribble in 2014 as a console project under Microsoft Visual Studio Express 2013. // // The version of the program uploaded to the OEIS compiles with g++ (gcc 4.5.2). using namespace std; #include #include #include #include #include #include #include #include #include static const int MaxDS = 400000; // Maximum number of distinct sequences. static const int MaxN = 10; // Maximum value of N. static const int MaxP = 15; // Maximum power to which M is raised. static const int MinN = 1; // Minimum value of N. static int A[MaxN * MaxN]; // Vector of binary numbers. static int DS[MaxDS][MaxP + 1]; // Set of distinct sequences. static int Freq[MaxDS]; // Frequency with which each value of FreqSeq occurs. static int FreqSeq[MaxDS]; // Frequency with which each sequence occurs. static int M[MaxN][MaxN]; // N X N (0,1) matrix. static int MaxFreqSeq; // Maximum value of FreqSeq. static int N; // User-specified dimension of N X N (0,1) matrix. static int N2; // N X N. static int NumDistinctSeqs; // Number of distinct sequences. static int P[MaxN][MaxN]; // Power of M. static int Q[MaxN][MaxN]; // Temporary power of M. static int S[MaxP + 1][MaxN][MaxN]; // Sequences of all element values of successive powers of M. fstream OpFile; // Output file. int main() { __int64 i; // Loop counter. int j, k, l, p; // Loop counters. bool duplicate; // Create output file. OpFile.open("Matrices.txt", ios_base::out | ios_base::trunc); // Get matrix dimension. cout << "Enter N "; cin >> N; if ((N < MinN) || (N > MaxN)) { OpFile << "Value of N " << N << " out of range." << endl; return 1; } OpFile << "N = " << N << endl; N2 = N * N; for (i = 0; i < pow(2, N2); i++) { if (i % 10000 == 0) { cout << i << " " << NumDistinctSeqs << endl; } // Form binary matrix. for (j = 0; j < N; j++) { for (k = 0; k < N; k++) { M[j][k] = A[j * N + k]; P[j][k] = M[j][k]; if (j == k) // Record 1st term of each sequence. { S[0][j][k] = 1; } else { S[0][j][k] = 0; } S[1][j][k] = P[j][k]; // Record 2nd term of each sequence. } } // Raise M to successive powers. for (p = 2; p <= MaxP; p++) { for (j = 0; j < N; j++) { for (k = 0; k < N; k++) { Q[j][k] = 0; } } for (j = 0; j < N; j++) { for (k = 0; k < N; k++) { for (l = 0; l < N; l++) { Q[j][k] += P[j][l] * M[l][k]; } } } for (j = 0; j < N; j++) { for (k = 0; k < N; k++) { P[j][k] = Q[j][k]; S[p][j][k] = P[j][k]; // Record pth term of each sequence. } } } /* for (j = 0; j < N; j++) { for (k = 0; k < N; k++) { OpFile << j << " " << k << " "; for (p = 1; p < MaxP; p++) { OpFile << S[p][j][k] << ", "; } OpFile << endl; } } OpFile << endl; */ // Record distinct sequences and frequencies of occurrence. for (j = 0; j < N; j++) { for (k = 0; k < N; k++) { if (NumDistinctSeqs == 0) { for (p = 0; p <= MaxP; p++) { DS[NumDistinctSeqs][p] = S[p][j][k]; } FreqSeq[NumDistinctSeqs] = 1; NumDistinctSeqs++; } else { duplicate = false; for (l = 0; l < NumDistinctSeqs; l++) { for (p = 0; p <= MaxP; p++) { if (S[p][j][k] != DS[l][p]) { break; } } if (p == MaxP + 1) { duplicate = true; FreqSeq[l]++; break; } } if (duplicate == false) { for (p = 0; p <= MaxP; p++) { DS[NumDistinctSeqs][p] = S[p][j][k]; } FreqSeq[NumDistinctSeqs] = 1; NumDistinctSeqs++; } } } } // Produce next binary number. A[N2 - 1] = A[N2 - 1] + 1; for (j = 1; j <= N2; j++) { if (A[N2 - j] == 2) { A[N2 - j - 1] = A[N2 - j - 1] + 1; A[N2 - j] = 0; } else { break; } } } for (i = 0; i < NumDistinctSeqs; i++) { OpFile << FreqSeq[i] << " "; for (p = 0; p < MaxP; p++) { OpFile << DS[i][p] << ", "; } OpFile << DS[i][MaxP] << endl; Freq[i] = FreqSeq[i]; } OpFile << endl << "NumDistinctSeqs = " << NumDistinctSeqs << endl << endl; int c, d, position, swap; for (c = 0; c < (NumDistinctSeqs - 1); c++) { position = c; for (d = c + 1; d < NumDistinctSeqs; d++) { if (Freq[position] > Freq[d]) position = d; } if (position != c) { swap = Freq[c]; Freq[c] = Freq[position]; Freq[position] = swap; } } j = 1; for (i = 1; i <= NumDistinctSeqs; i++) { if (Freq[i] == Freq[i - 1]) { j++; } else { OpFile << Freq[i - 1] << " " << j << endl; j = 1; } } return 0; }