Skip to content

Instantly share code, notes, and snippets.

@xdsopl

xdsopl/cauchy.cc

Last active April 19, 2023 13:59
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save xdsopl/8371feede32a1e8ceadf64beff7b5b06 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save xdsopl/8371feede32a1e8ceadf64beff7b5b06 to your computer and use it in GitHub Desktop.
Download ZIP
Generating and inverting Cauchy matrices using Galois field values
Raw
cauchy.cc
#include <iostream>
#include "galois_field.hh"
#if 0
typedef CODE::GaloisField<8, 0b100011101, uint8_t> GF;
#else
typedef CODE::GaloisField<16, 0b10001000000001011, uint16_t> GF;
#endif
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
// The Art of Computer Programming
// Volume 1: Fundamental Algorithms
// Donald E. Knuth - 1997
// Paragraph 1.2.3: Sums and Products
// Page 37, Cauchy's matrix:
// $a_{ij} = \frac{1}{x_i + y_j}$
IndexType cauchy_matrix(ValueType row, int j, int n) {
ValueType col(ValueType::N - j);
return rcp(index(row + col));
}
// Page 38, Exercise 41:
// $b_{ij} = \frac{\prod_{k=1}^{n}{(x_j + y_k)(x_k + y_i)}}{(x_j + y_i)\prod_{k \ne j}^{n}{(x_j - x_k)}\prod_{k \ne i}^{n}{(y_i - y_k)}}$
IndexType inverse_cauchy_matrix(const ValueType *rows, int i, int j, int n) {
ValueType col_i(ValueType::N - i);
IndexType prod_xy(0);
for (int k = 0; k < n; k++) {
ValueType col_k(ValueType::N - k);
prod_xy *= index(rows[j] + col_k) * index(rows[k] + col_i);
}
IndexType prod_x(0);
for (int k = 0; k < n; k++) {
if (k != j) {
prod_x *= index(rows[j] + rows[k]);
}
}
IndexType prod_y(0);
for (int k = 0; k < n; k++) {
if (k != i) {
ValueType col_k(ValueType::N - k);
prod_y *= index(col_i + col_k);
}
}
return prod_xy / (index(rows[j] + col_i) * prod_x * prod_y);
}
void matrix_multiply(IndexType *A, IndexType *B, ValueType *result, int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
ValueType sum(0);
for (int k = 0; k < n; k++) {
sum += value(A[n*i+k] * B[n*k+j]);
}
*result++ = sum;
}
}
}
bool matrices_are_equal(ValueType *A, ValueType *B, int n) {
for (int i = 0; i < n * n; i++) {
if (*A++ != *B++) {
return false;
}
}
return true;
}
void print_matrix(ValueType *matrix, int n) {
if (1) {
std::cout << std::endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
std::cout << " " << (int)matrix[n*i+j].v;
}
std::cout << std::endl;
}
std::cout << std::endl;
}
}
void print_matrix(IndexType *matrix, int n) {
if (1) {
std::cout << std::endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
std::cout << " " << (int)matrix[n*i+j].i;
}
std::cout << std::endl;
}
std::cout << std::endl;
}
}
void fisher_yates_shuffle(ValueType *vector, int n) {
for (int i = 0; i < n - 1; i++) {
int j = i + rand() % (n - i);
std::swap(vector[i], vector[j]);
}
}
int main(int argc, char **argv) {
if (argc != 2)
return 1;
int n = std::atoi(argv[1]);
if (n < 0 || n > ValueType::Q / 2)
return 1;
int m = ValueType::Q - n;
srand(time(0));
GF *instance = new GF();
ValueType *rows = new ValueType[m];
for (int i = 0; i < m; i++) {
rows[i] = ValueType(i);
}
fisher_yates_shuffle(rows, m);
IndexType *matrix = new IndexType[n*n];
IndexType *inverse = new IndexType[n*n];
ValueType *product = new ValueType[n*n];
ValueType *identity = new ValueType[n*n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
identity[n*i+j] = ValueType(i == j);
}
}
print_matrix(identity, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
matrix[n*i+j] = cauchy_matrix(rows[i], j, n);
}
}
print_matrix(matrix, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
inverse[n*i+j] = inverse_cauchy_matrix(rows, i, j, n);
}
}
print_matrix(inverse, n);
matrix_multiply(matrix, inverse, product, n);
print_matrix(product, n);
if (matrices_are_equal(product, identity, n)) {
std::cerr << "The inverse is correct." << std::endl;
} else {
std::cerr << "The inverse is incorrect." << std::endl;
}
delete[] matrix;
delete[] inverse;
delete[] product;
delete[] identity;
delete instance;
return 0;
}
Raw
large.cc
#include <iostream>
#include "galois_field.hh"
typedef CODE::GaloisFieldReference<32, 0x100400007, uint32_t> GFR;
// The Art of Computer Programming
// Volume 1: Fundamental Algorithms
// Donald E. Knuth - 1997
// Paragraph 1.2.3: Sums and Products
// Page 37, Cauchy's matrix:
// $a_{ij} = \frac{1}{x_i + y_j}$
GFR cauchy_matrix(GFR row, int j, int n) {
GFR col(GFR::N - j);
return rcp(row + col);
}
// Page 38, Exercise 41:
// $b_{ij} = \frac{\prod_{k=1}^{n}{(x_j + y_k)(x_k + y_i)}}{(x_j + y_i)\prod_{k \ne j}^{n}{(x_j - x_k)}\prod_{k \ne i}^{n}{(y_i - y_k)}}$
GFR inverse_cauchy_matrix(const GFR *rows, int i, int j, int n) {
GFR col_i(GFR::N - i);
GFR prod_xy(1);
for (int k = 0; k < n; k++) {
GFR col_k(GFR::N - k);
prod_xy *= (rows[j] + col_k) * (rows[k] + col_i);
}
GFR prod_x(1);
for (int k = 0; k < n; k++) {
if (k != j) {
prod_x *= rows[j] + rows[k];
}
}
GFR prod_y(1);
for (int k = 0; k < n; k++) {
if (k != i) {
GFR col_k(GFR::N - k);
prod_y *= col_i + col_k;
}
}
return prod_xy / ((rows[j] + col_i) * prod_x * prod_y);
}
void matrix_multiply(GFR *A, GFR *B, GFR *result, int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
GFR sum(0);
for (int k = 0; k < n; k++) {
sum += A[n*i+k] * B[n*k+j];
}
*result++ = sum;
}
}
}
bool matrices_are_equal(GFR *A, GFR *B, int n) {
for (int i = 0; i < n * n; i++) {
if (*A++ != *B++) {
return false;
}
}
return true;
}
void print_matrix(GFR *matrix, int n) {
if (1) {
std::cout << std::endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
std::cout << " " << (int)matrix[n*i+j].v;
}
std::cout << std::endl;
}
std::cout << std::endl;
}
}
void fisher_yates_shuffle(GFR *vector, int n) {
for (int i = 0; i < n - 1; i++) {
int j = i + rand() % (n - i);
std::swap(vector[i], vector[j]);
}
}
int main(int argc, char **argv) {
if (argc != 3)
return 1;
int n = std::atoi(argv[1]);
if (n < 0)
return 1;
int m = std::atoi(argv[2]);
if (m < 0)
return 1;
srand(time(0));
GFR *rows = new GFR[m];
for (int i = 0; i < m; i++) {
rows[i] = GFR(i);
}
fisher_yates_shuffle(rows, m);
GFR *matrix = new GFR[n*n];
GFR *inverse = new GFR[n*n];
GFR *product = new GFR[n*n];
GFR *identity = new GFR[n*n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
identity[n*i+j] = GFR(i == j);
}
}
print_matrix(identity, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
matrix[n*i+j] = cauchy_matrix(rows[i], j, n);
}
}
print_matrix(matrix, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
inverse[n*i+j] = inverse_cauchy_matrix(rows, i, j, n);
}
}
print_matrix(inverse, n);
matrix_multiply(matrix, inverse, product, n);
print_matrix(product, n);
if (matrices_are_equal(product, identity, n)) {
std::cerr << "The inverse is correct." << std::endl;
} else {
std::cerr << "The inverse is incorrect." << std::endl;
}
delete[] matrix;
delete[] inverse;
delete[] product;
delete[] identity;
return 0;
}
Raw
reference.cc
#include <iostream>
#include "galois_field.hh"
#if 0
typedef CODE::GaloisFieldReference<8, 0b100011101, uint8_t> GFR;
#else
typedef CODE::GaloisFieldReference<16, 0b10001000000001011, uint16_t> GFR;
#endif
// The Art of Computer Programming
// Volume 1: Fundamental Algorithms
// Donald E. Knuth - 1997
// Paragraph 1.2.3: Sums and Products
// Page 37, Cauchy's matrix:
// $a_{ij} = \frac{1}{x_i + y_j}$
GFR cauchy_matrix(GFR row, int j, int n) {
GFR col(GFR::N - j);
return rcp(row + col);
}
// Page 38, Exercise 41:
// $b_{ij} = \frac{\prod_{k=1}^{n}{(x_j + y_k)(x_k + y_i)}}{(x_j + y_i)\prod_{k \ne j}^{n}{(x_j - x_k)}\prod_{k \ne i}^{n}{(y_i - y_k)}}$
GFR inverse_cauchy_matrix(const GFR *rows, int i, int j, int n) {
GFR col_i(GFR::N - i);
GFR prod_xy(1);
for (int k = 0; k < n; k++) {
GFR col_k(GFR::N - k);
prod_xy *= (rows[j] + col_k) * (rows[k] + col_i);
}
GFR prod_x(1);
for (int k = 0; k < n; k++) {
if (k != j) {
prod_x *= rows[j] + rows[k];
}
}
GFR prod_y(1);
for (int k = 0; k < n; k++) {
if (k != i) {
GFR col_k(GFR::N - k);
prod_y *= col_i + col_k;
}
}
return prod_xy / ((rows[j] + col_i) * prod_x * prod_y);
}
void matrix_multiply(GFR *A, GFR *B, GFR *result, int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
GFR sum(0);
for (int k = 0; k < n; k++) {
sum += A[n*i+k] * B[n*k+j];
}
*result++ = sum;
}
}
}
bool matrices_are_equal(GFR *A, GFR *B, int n) {
for (int i = 0; i < n * n; i++) {
if (*A++ != *B++) {
return false;
}
}
return true;
}
void print_matrix(GFR *matrix, int n) {
if (1) {
std::cout << std::endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
std::cout << " " << (int)matrix[n*i+j].v;
}
std::cout << std::endl;
}
std::cout << std::endl;
}
}
void fisher_yates_shuffle(GFR *vector, int n) {
for (int i = 0; i < n - 1; i++) {
int j = i + rand() % (n - i);
std::swap(vector[i], vector[j]);
}
}
int main(int argc, char **argv) {
if (argc != 2)
return 1;
int n = std::atoi(argv[1]);
if (n < 0 || n > GFR::Q / 2)
return 1;
int m = GFR::Q - n;
srand(time(0));
GFR *rows = new GFR[m];
for (int i = 0; i < m; i++) {
rows[i] = GFR(i);
}
fisher_yates_shuffle(rows, m);
GFR *matrix = new GFR[n*n];
GFR *inverse = new GFR[n*n];
GFR *product = new GFR[n*n];
GFR *identity = new GFR[n*n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
identity[n*i+j] = GFR(i == j);
}
}
print_matrix(identity, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
matrix[n*i+j] = cauchy_matrix(rows[i], j, n);
}
}
print_matrix(matrix, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
inverse[n*i+j] = inverse_cauchy_matrix(rows, i, j, n);
}
}
print_matrix(inverse, n);
matrix_multiply(matrix, inverse, product, n);
print_matrix(product, n);
if (matrices_are_equal(product, identity, n)) {
std::cerr << "The inverse is correct." << std::endl;
} else {
std::cerr << "The inverse is incorrect." << std::endl;
}
delete[] matrix;
delete[] inverse;
delete[] product;
delete[] identity;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment