diff options
Diffstat (limited to 'src/maths')
| -rw-r--r-- | src/maths/complex.c | 122 | ||||
| -rw-r--r-- | src/maths/format.c | 265 | ||||
| -rw-r--r-- | src/maths/matrix.c | 192 | ||||
| -rw-r--r-- | src/maths/matrix.rs | 310 | ||||
| -rw-r--r-- | src/maths/simd.c | 424 | ||||
| -rw-r--r-- | src/maths/vector.c | 136 | ||||
| -rw-r--r-- | src/maths/vector.rs | 258 |
7 files changed, 596 insertions, 1111 deletions
diff --git a/src/maths/complex.c b/src/maths/complex.c index 44e255f..419d0b0 100644 --- a/src/maths/complex.c +++ b/src/maths/complex.c @@ -4,126 +4,126 @@ struct PsiComplex psi_new_complex(double real, double imaginary) { - return (struct PsiComplex){ - real, - imaginary, - }; + return (struct PsiComplex){ + real, + imaginary, + }; } struct PsiComplex psi_new_complex_from_real(double real) { - return (struct PsiComplex){ - real, - 0.0, - }; + return (struct PsiComplex){ + real, + 0.0, + }; } struct PsiComplex psi_conjugate_complex(struct PsiComplex z) { - return (struct PsiComplex){ - z.real, - -z.imaginary, - }; + return (struct PsiComplex){ + z.real, + -z.imaginary, + }; } struct PsiComplex psi_neg_complex(struct PsiComplex z) { - return (struct PsiComplex){ - -z.real, - -z.imaginary, - }; + return (struct PsiComplex){ + -z.real, + -z.imaginary, + }; } double psi_phase_complex(struct PsiComplex z) { - return atan2(z.imaginary, z.real); + return atan2(z.imaginary, z.real); } double psi_norm2_complex(struct PsiComplex z) { - return z.real * z.real + z.imaginary * z.imaginary; + return z.real * z.real + z.imaginary * z.imaginary; } double psi_abs_complex(struct PsiComplex z) { - return sqrt(psi_norm2_complex(z)); + return sqrt(psi_norm2_complex(z)); } struct PsiComplex psi_sqrt_complex(struct PsiComplex z) { - double r = psi_abs_complex(z); - double half_theta = psi_phase_complex(z) / 2.0; - double sqrt_r = sqrt(r); + double r = psi_abs_complex(z); + double half_theta = psi_phase_complex(z) / 2.0; + double sqrt_r = sqrt(r); - return (struct PsiComplex){ - sqrt_r * cos(half_theta), - sqrt_r * sin(half_theta), - }; + return (struct PsiComplex){ + sqrt_r * cos(half_theta), + sqrt_r * sin(half_theta), + }; } struct PsiComplex psi_add_complex(struct PsiComplex a, struct PsiComplex b) { - return (struct PsiComplex){ - a.real + b.real, - a.imaginary + b.imaginary, - }; + return (struct PsiComplex){ + a.real + b.real, + a.imaginary + b.imaginary, + }; } struct PsiComplex psi_sub_complex(struct PsiComplex a, struct PsiComplex b) { - return (struct PsiComplex){ - a.real - b.real, - a.imaginary - b.imaginary, - }; + return (struct PsiComplex){ + a.real - b.real, + a.imaginary - b.imaginary, + }; } struct PsiComplex psi_mul_complex(struct PsiComplex a, struct PsiComplex b) { - // (a + bi)(c + di) = (ac - bd) + (ad + bc)i - return (struct PsiComplex){ - a.real * b.real - a.imaginary * b.imaginary, - a.real * b.imaginary + a.imaginary * b.real, - }; + // (a + bi)(c + di) = (ac - bd) + (ad + bc)i + return (struct PsiComplex){ + a.real * b.real - a.imaginary * b.imaginary, + a.real * b.imaginary + a.imaginary * b.real, + }; } struct PsiComplex psi_div_complex(struct PsiComplex a, struct PsiComplex b) { - // (a + bi) / (c + di) = ((ac + bd) + (bc - ad)i) / (c² + d²) - double denom = b.real * b.real + b.imaginary * b.imaginary; - return (struct PsiComplex){ - (a.real * b.real + a.imaginary * b.imaginary) / denom, - (a.imaginary * b.real - a.real * b.imaginary) / denom, - }; + // (a + bi) / (c + di) = ((ac + bd) + (bc - ad)i) / (c² + d²) + double denom = b.real * b.real + b.imaginary * b.imaginary; + return (struct PsiComplex){ + (a.real * b.real + a.imaginary * b.imaginary) / denom, + (a.imaginary * b.real - a.real * b.imaginary) / denom, + }; } struct PsiComplex psi_add_complex_real(struct PsiComplex a, double b) { - return (struct PsiComplex){ - a.real + b, - a.imaginary, - }; + return (struct PsiComplex){ + a.real + b, + a.imaginary, + }; } struct PsiComplex psi_sub_complex_real(struct PsiComplex a, double b) { - return (struct PsiComplex){ - a.real - b, - a.imaginary, - }; + return (struct PsiComplex){ + a.real - b, + a.imaginary, + }; } struct PsiComplex psi_mul_complex_real(struct PsiComplex a, double b) { - return (struct PsiComplex){ - a.real * b, - a.imaginary * b, - }; + return (struct PsiComplex){ + a.real * b, + a.imaginary * b, + }; } struct PsiComplex psi_div_complex_real(struct PsiComplex a, double b) { - return (struct PsiComplex){ - a.real / b, - a.imaginary / b, - }; + return (struct PsiComplex){ + a.real / b, + a.imaginary / b, + }; } diff --git a/src/maths/format.c b/src/maths/format.c index f3db83d..52b2b5a 100644 --- a/src/maths/format.c +++ b/src/maths/format.c @@ -12,136 +12,189 @@ static const double INV_SQRT_32 = 0.1767766952966369; static bool approx_eq(double a, double b) { - return fabs(a - b) < EPSILON; + return fabs(a - b) < EPSILON; } static bool real_symbolic(double v, char* out, size_t cap) { - double abs_v = fabs(v); - const char* sign = v < 0.0 ? "-" : ""; + double abs_v = fabs(v); + const char* sign = v < 0.0 ? "-" : ""; - if (approx_eq(abs_v, 0.0)) - { - snprintf(out, cap, "%s", "0"); - return true; - } + if (approx_eq(abs_v, 0.0)) + { + snprintf(out, cap, "%s", "0"); + return true; + } - const char* sym = NULL; - if (approx_eq(abs_v, 1.0)) - sym = "1"; - else if (approx_eq(abs_v, 0.5)) - sym = "½"; - else if (approx_eq(abs_v, 0.25)) - sym = "¼"; - else if (approx_eq(abs_v, 0.75)) - sym = "¾"; - else if (approx_eq(abs_v, 0.125)) - sym = "⅛"; - else if (approx_eq(abs_v, SQRT_2)) - sym = "√2"; - else if (approx_eq(abs_v, INV_SQRT_2)) - sym = "¹⁄√2"; - else if (approx_eq(abs_v, INV_SQRT_8)) - sym = "¹⁄√8"; - else if (approx_eq(abs_v, INV_SQRT_32)) - sym = "¹⁄√32"; - else if (approx_eq(abs_v, 2.0)) - sym = "2"; - else if (approx_eq(abs_v, 1.0 / 3.0)) - sym = "⅓"; - else if (approx_eq(abs_v, 2.0 / 3.0)) - sym = "⅔"; - else - return false; + const char* sym = NULL; + if (approx_eq(abs_v, 1.0)) + sym = "1"; + else if (approx_eq(abs_v, 0.5)) + sym = "½"; + else if (approx_eq(abs_v, 0.25)) + sym = "¼"; + else if (approx_eq(abs_v, 0.75)) + sym = "¾"; + else if (approx_eq(abs_v, 0.125)) + sym = "⅛"; + else if (approx_eq(abs_v, SQRT_2)) + sym = "√2"; + else if (approx_eq(abs_v, INV_SQRT_2)) + sym = "¹⁄√2"; + else if (approx_eq(abs_v, INV_SQRT_8)) + sym = "¹⁄√8"; + else if (approx_eq(abs_v, INV_SQRT_32)) + sym = "¹⁄√32"; + else if (approx_eq(abs_v, 2.0)) + sym = "2"; + else if (approx_eq(abs_v, 1.0 / 3.0)) + sym = "⅓"; + else if (approx_eq(abs_v, 2.0 / 3.0)) + sym = "⅔"; + else + return false; - snprintf(out, cap, "%s%s", sign, sym); - return true; + snprintf(out, cap, "%s%s", sign, sym); + return true; } char* psi_format_amplitude(struct PsiComplex c, char* out, size_t cap) { - double re = c.real; - double im = c.imaginary; + double re = c.real; + double im = c.imaginary; - bool re_zero = approx_eq(fabs(re), 0.0); - bool im_zero = approx_eq(fabs(im), 0.0); + bool re_zero = approx_eq(fabs(re), 0.0); + bool im_zero = approx_eq(fabs(im), 0.0); - if (re_zero && im_zero) - { - snprintf(out, cap, "%s", "0"); - return out; - } + if (re_zero && im_zero) + { + snprintf(out, cap, "%s", "0"); + return out; + } - if (im_zero) - { - if (!real_symbolic(re, out, cap)) - snprintf(out, cap, "%.4f", re); + if (im_zero) + { + if (!real_symbolic(re, out, cap)) + snprintf(out, cap, "%.4f", re); - return out; - } + return out; + } - if (re_zero) - { - if (approx_eq(fabs(im), 1.0)) - { - snprintf(out, cap, "%s", im > 0.0 ? "i" : "-i"); - return out; - } + if (re_zero) + { + if (approx_eq(fabs(im), 1.0)) + { + snprintf(out, cap, "%s", im > 0.0 ? "i" : "-i"); + return out; + } - char sym[32]; - if (real_symbolic(im, sym, sizeof sym)) - snprintf(out, cap, "%si", sym); - else - snprintf(out, cap, "%.4fi", im); + char sym[32]; + if (real_symbolic(im, sym, sizeof sym)) + snprintf(out, cap, "%si", sym); + else + snprintf(out, cap, "%.4fi", im); - return out; - } + return out; + } - char re_str[32]; - if (!real_symbolic(re, re_str, sizeof re_str)) - snprintf(re_str, sizeof re_str, "%.4f", re); + char re_str[32]; + if (!real_symbolic(re, re_str, sizeof re_str)) + snprintf(re_str, sizeof re_str, "%.4f", re); - char im_str[32]; - if (approx_eq(fabs(im), 1.0)) - { - snprintf(im_str, sizeof im_str, "%s", im > 0.0 ? "+i" : "-i"); - } - else - { - const char* sign = im > 0.0 ? "+" : "-"; - char sym[32]; - if (real_symbolic(fabs(im), sym, sizeof sym)) - snprintf(im_str, sizeof im_str, "%s%si", sign, sym); - else - snprintf(im_str, sizeof im_str, "%s%.4fi", sign, fabs(im)); - } + char im_str[32]; + if (approx_eq(fabs(im), 1.0)) + { + snprintf(im_str, sizeof im_str, "%s", im > 0.0 ? "+i" : "-i"); + } + else + { + const char* sign = im > 0.0 ? "+" : "-"; + char sym[32]; + if (real_symbolic(fabs(im), sym, sizeof sym)) + snprintf(im_str, sizeof im_str, "%s%si", sign, sym); + else + snprintf(im_str, sizeof im_str, "%s%.4fi", sign, fabs(im)); + } - snprintf(out, cap, "%s%s", re_str, im_str); - return out; + snprintf(out, cap, "%s%s", re_str, im_str); + return out; } char* psi_format_probability(double p, char* out, size_t cap) { - if (approx_eq(p, 0.0)) - snprintf(out, cap, "%s", "0"); - else if (approx_eq(p, 1.0)) - snprintf(out, cap, "%s", "1"); - else if (approx_eq(p, 0.5)) - snprintf(out, cap, "%s", "½"); - else if (approx_eq(p, 0.25)) - snprintf(out, cap, "%s", "¼"); - else if (approx_eq(p, 0.75)) - snprintf(out, cap, "%s", "¾"); - else if (approx_eq(p, 0.125)) - snprintf(out, cap, "%s", "⅛"); - else if (approx_eq(p, 0.0625)) - snprintf(out, cap, "%s", "¹⁄₁₆"); - else if (approx_eq(p, 1.0 / 3.0)) - snprintf(out, cap, "%s", "⅓"); - else if (approx_eq(p, 2.0 / 3.0)) - snprintf(out, cap, "%s", "⅔"); - else - snprintf(out, cap, "%.4f", p); + if (approx_eq(p, 0.0)) + snprintf(out, cap, "%s", "0"); + else if (approx_eq(p, 1.0)) + snprintf(out, cap, "%s", "1"); + else if (approx_eq(p, 0.5)) + snprintf(out, cap, "%s", "½"); + else if (approx_eq(p, 0.25)) + snprintf(out, cap, "%s", "¼"); + else if (approx_eq(p, 0.75)) + snprintf(out, cap, "%s", "¾"); + else if (approx_eq(p, 0.125)) + snprintf(out, cap, "%s", "⅛"); + else if (approx_eq(p, 0.0625)) + snprintf(out, cap, "%s", "¹⁄₁₆"); + else if (approx_eq(p, 1.0 / 3.0)) + snprintf(out, cap, "%s", "⅓"); + else if (approx_eq(p, 2.0 / 3.0)) + snprintf(out, cap, "%s", "⅔"); + else + snprintf(out, cap, "%.4f", p); - return out; + return out; +} + +void psi_print_matrix(struct PsiMatrix m, FILE* out) +{ + int re_width = 0; + int im_width = 0; + for (size_t k = 0; k < m.rows * m.cols; k++) + { + int rl = snprintf(NULL, 0, "%.2f", m.data[k].real); + int il = snprintf(NULL, 0, "%.2f", fabs(m.data[k].imaginary)); + if (rl > re_width) + re_width = rl; + if (il > im_width) + im_width = il; + } + + for (size_t i = 0; i < m.rows; i++) + { + fputs(i == 0 ? "┌" : (i == m.rows - 1 ? "└" : "│"), out); + + for (size_t j = 0; j < m.cols; j++) + { + struct PsiComplex e = m.data[i * m.cols + j]; + const char* sign = e.imaginary > 0.0 ? "+" : "-"; + fprintf(out, "%*.2f %s %*.2fi", re_width, e.real, sign, im_width, fabs(e.imaginary)); + if (j != m.cols - 1) + fputs(", ", out); + } + + fputs(i == 0 ? "┐" : (i == m.rows - 1 ? "┘" : "│"), out); + if (i != m.rows - 1) + fputc('\n', out); + } +} + +void psi_print_vector(struct PsiVector v, FILE* out) +{ + if (v.kind == PSI_COLUMN_VECTOR) + { + struct PsiMatrix m = psi_matrix_from_vector(v); + psi_print_matrix(m, out); + psi_free_matrix(&m); + return; + } + + fputc('[', out); + for (size_t i = 0; i < v.size; i++) + { + fprintf(out, "%g + %gi", v.data[i].real, v.data[i].imaginary); + if (i != v.size - 1) + fputs(", ", out); + } + fputc(']', out); } diff --git a/src/maths/matrix.c b/src/maths/matrix.c index 2415395..8b41065 100644 --- a/src/maths/matrix.c +++ b/src/maths/matrix.c @@ -6,181 +6,181 @@ struct PsiMatrix psi_new_matrix(size_t rows, size_t cols) { - struct PsiComplex* data = calloc(rows * cols, sizeof(struct PsiComplex)); - assert(data != NULL || rows * cols == 0); + struct PsiComplex* data = calloc(rows * cols, sizeof(struct PsiComplex)); + assert(data != NULL || rows * cols == 0); - return (struct PsiMatrix){ - data, - rows, - cols, - }; + return (struct PsiMatrix){ + data, + rows, + cols, + }; } struct PsiMatrix psi_new_matrix_from(const struct PsiComplex* data, size_t rows, size_t cols) { - struct PsiMatrix m = psi_new_matrix(rows, cols); - memcpy(m.data, data, rows * cols * sizeof(struct PsiComplex)); + struct PsiMatrix m = psi_new_matrix(rows, cols); + memcpy(m.data, data, rows * cols * sizeof(struct PsiComplex)); - return m; + return m; } struct PsiMatrix psi_clone_matrix(struct PsiMatrix m) { - return psi_new_matrix_from(m.data, m.rows, m.cols); + return psi_new_matrix_from(m.data, m.rows, m.cols); } struct PsiMatrix psi_identity_matrix(size_t size) { - struct PsiMatrix m = psi_new_matrix(size, size); - for (size_t i = 0; i < size; i++) - m.data[i * size + i] = psi_new_complex(1.0, 0.0); + struct PsiMatrix m = psi_new_matrix(size, size); + for (size_t i = 0; i < size; i++) + m.data[i * size + i] = psi_new_complex(1.0, 0.0); - return m; + return m; } void psi_free_matrix(struct PsiMatrix* m) { - free(m->data); - m->data = NULL; - m->rows = 0; - m->cols = 0; + free(m->data); + m->data = NULL; + m->rows = 0; + m->cols = 0; } struct PsiComplex psi_get_matrix(struct PsiMatrix m, size_t row, size_t col) { - assert(row < m.rows && col < m.cols); - return m.data[row * m.cols + col]; + assert(row < m.rows && col < m.cols); + return m.data[row * m.cols + col]; } void psi_set_matrix(struct PsiMatrix* m, size_t row, size_t col, struct PsiComplex value) { - assert(row < m->rows && col < m->cols); - m->data[row * m->cols + col] = value; + assert(row < m->rows && col < m->cols); + m->data[row * m->cols + col] = value; } struct PsiMatrix psi_dot_matrix(struct PsiMatrix a, struct PsiMatrix b) { - assert(a.cols == b.rows); + assert(a.cols == b.rows); - struct PsiMatrix result = psi_new_matrix(a.rows, b.cols); - for (size_t i = 0; i < a.rows; i++) - for (size_t j = 0; j < b.cols; j++) - { - struct PsiComplex sum = psi_new_complex(0.0, 0.0); - for (size_t k = 0; k < a.cols; k++) - sum = psi_add_complex( - sum, psi_mul_complex(a.data[i * a.cols + k], b.data[k * b.cols + j])); + struct PsiMatrix result = psi_new_matrix(a.rows, b.cols); + for (size_t i = 0; i < a.rows; i++) + for (size_t j = 0; j < b.cols; j++) + { + struct PsiComplex sum = psi_new_complex(0.0, 0.0); + for (size_t k = 0; k < a.cols; k++) + sum = psi_add_complex( + sum, psi_mul_complex(a.data[i * a.cols + k], b.data[k * b.cols + j])); - result.data[i * result.cols + j] = sum; - } + result.data[i * result.cols + j] = sum; + } - return result; + return result; } struct PsiMatrix psi_kronecker_matrix(struct PsiMatrix a, struct PsiMatrix b) { - struct PsiMatrix result = psi_new_matrix(a.rows * b.rows, a.cols * b.cols); - for (size_t i = 0; i < a.rows; i++) - for (size_t j = 0; j < a.cols; j++) - { - struct PsiComplex a_val = a.data[i * a.cols + j]; - for (size_t k = 0; k < b.rows; k++) - for (size_t l = 0; l < b.cols; l++) - { - size_t row = i * b.rows + k; - size_t col = j * b.cols + l; - result.data[row * result.cols + col] = - psi_mul_complex(a_val, b.data[k * b.cols + l]); - } - } + struct PsiMatrix result = psi_new_matrix(a.rows * b.rows, a.cols * b.cols); + for (size_t i = 0; i < a.rows; i++) + for (size_t j = 0; j < a.cols; j++) + { + struct PsiComplex a_val = a.data[i * a.cols + j]; + for (size_t k = 0; k < b.rows; k++) + for (size_t l = 0; l < b.cols; l++) + { + size_t row = i * b.rows + k; + size_t col = j * b.cols + l; + result.data[row * result.cols + col] = + psi_mul_complex(a_val, b.data[k * b.cols + l]); + } + } - return result; + return result; } struct PsiMatrix psi_transpose_matrix(struct PsiMatrix m) { - struct PsiMatrix result = psi_new_matrix(m.cols, m.rows); - for (size_t i = 0; i < m.rows; i++) - for (size_t j = 0; j < m.cols; j++) - result.data[j * result.cols + i] = m.data[i * m.cols + j]; + struct PsiMatrix result = psi_new_matrix(m.cols, m.rows); + for (size_t i = 0; i < m.rows; i++) + for (size_t j = 0; j < m.cols; j++) + result.data[j * result.cols + i] = m.data[i * m.cols + j]; - return result; + return result; } struct PsiMatrix psi_add_matrix(struct PsiMatrix a, struct PsiMatrix b) { - assert(a.rows == b.rows && a.cols == b.cols); + assert(a.rows == b.rows && a.cols == b.cols); - struct PsiMatrix result = psi_new_matrix(a.rows, a.cols); - for (size_t i = 0; i < a.rows * a.cols; i++) - result.data[i] = psi_add_complex(a.data[i], b.data[i]); + struct PsiMatrix result = psi_new_matrix(a.rows, a.cols); + for (size_t i = 0; i < a.rows * a.cols; i++) + result.data[i] = psi_add_complex(a.data[i], b.data[i]); - return result; + return result; } struct PsiMatrix psi_sub_matrix(struct PsiMatrix a, struct PsiMatrix b) { - assert(a.rows == b.rows && a.cols == b.cols); + assert(a.rows == b.rows && a.cols == b.cols); - struct PsiMatrix result = psi_new_matrix(a.rows, a.cols); - for (size_t i = 0; i < a.rows * a.cols; i++) - result.data[i] = psi_sub_complex(a.data[i], b.data[i]); + struct PsiMatrix result = psi_new_matrix(a.rows, a.cols); + for (size_t i = 0; i < a.rows * a.cols; i++) + result.data[i] = psi_sub_complex(a.data[i], b.data[i]); - return result; + return result; } struct PsiMatrix psi_scale_matrix(struct PsiMatrix m, struct PsiComplex scalar) { - struct PsiMatrix result = psi_new_matrix(m.rows, m.cols); - for (size_t i = 0; i < m.rows * m.cols; i++) - result.data[i] = psi_mul_complex(m.data[i], scalar); + struct PsiMatrix result = psi_new_matrix(m.rows, m.cols); + for (size_t i = 0; i < m.rows * m.cols; i++) + result.data[i] = psi_mul_complex(m.data[i], scalar); - return result; + return result; } struct PsiMatrix psi_matrix_from_vector(struct PsiVector v) { - if (v.kind == PSI_ROW_VECTOR) - return psi_new_matrix_from(v.data, 1, v.size); + if (v.kind == PSI_ROW_VECTOR) + return psi_new_matrix_from(v.data, 1, v.size); - return psi_new_matrix_from(v.data, v.size, 1); + return psi_new_matrix_from(v.data, v.size, 1); } struct PsiVector psi_vector_from_matrix(struct PsiMatrix m, enum PsiVectorKind kind) { - return psi_new_vector_from(m.data, m.rows * m.cols, kind); + return psi_new_vector_from(m.data, m.rows * m.cols, kind); } struct PsiVector psi_mul_vector_matrix(struct PsiVector v, struct PsiMatrix m) { - if (v.kind == PSI_COLUMN_VECTOR) - { - assert(m.cols == v.size); + if (v.kind == PSI_COLUMN_VECTOR) + { + assert(m.cols == v.size); - struct PsiVector result = psi_new_vector(m.rows, PSI_COLUMN_VECTOR); - for (size_t i = 0; i < m.rows; i++) - { - struct PsiComplex sum = psi_new_complex(0.0, 0.0); - for (size_t j = 0; j < m.cols; j++) - sum = psi_add_complex(sum, psi_mul_complex(m.data[i * m.cols + j], v.data[j])); + struct PsiVector result = psi_new_vector(m.rows, PSI_COLUMN_VECTOR); + for (size_t i = 0; i < m.rows; i++) + { + struct PsiComplex sum = psi_new_complex(0.0, 0.0); + for (size_t j = 0; j < m.cols; j++) + sum = psi_add_complex(sum, psi_mul_complex(m.data[i * m.cols + j], v.data[j])); - result.data[i] = sum; - } + result.data[i] = sum; + } - return result; - } + return result; + } - assert(v.size == m.rows); + assert(v.size == m.rows); - struct PsiVector result = psi_new_vector(m.cols, PSI_ROW_VECTOR); - for (size_t j = 0; j < m.cols; j++) - { - struct PsiComplex sum = psi_new_complex(0.0, 0.0); - for (size_t i = 0; i < m.rows; i++) - sum = psi_add_complex(sum, psi_mul_complex(v.data[i], m.data[i * m.cols + j])); + struct PsiVector result = psi_new_vector(m.cols, PSI_ROW_VECTOR); + for (size_t j = 0; j < m.cols; j++) + { + struct PsiComplex sum = psi_new_complex(0.0, 0.0); + for (size_t i = 0; i < m.rows; i++) + sum = psi_add_complex(sum, psi_mul_complex(v.data[i], m.data[i * m.cols + j])); - result.data[j] = sum; - } + result.data[j] = sum; + } - return result; + return result; } diff --git a/src/maths/matrix.rs b/src/maths/matrix.rs deleted file mode 100644 index c9d96a4..0000000 --- a/src/maths/matrix.rs +++ /dev/null @@ -1,310 +0,0 @@ -use super::Float; -use core::{fmt, ops}; - -#[macro_export] -macro_rules! matrix { - ( $( $( $x:expr ),* );* ) => {{ - let mut data = Vec::new(); - let mut rows = 0; - let mut cols = 0; - - $( - let row_data = $( $x )*; - if cols == 0 { - cols = row_data.len(); - } - assert_eq!(cols, row_data.len(), "All rows must have the same number of columns."); - data.extend(row_data); - rows += 1; - )* - - $crate::Matrix::new(rows, cols, data) - }}; -} - -macro_rules! impl_matrix_ops { - ($($trait:ident, $method:ident, $other:ty, $output:ty, $scale_fn:ident),* $(,)?) => { - $( - impl<T: Float> core::ops::$trait<$other> for Matrix<T> { - type Output = $output; - - fn $method(self, other: $other) -> Self::Output { - self.$scale_fn(other) - } - } - )* - }; - ($($trait:ident, $method:ident, $other:ty, $scale_fn:ident),* $(,)?) => { - $( - impl<T: Float> core::ops::$trait<$other> for Matrix<T> { - fn $method(&mut self, other: $other) { - *self = self.$scale_fn(other); - } - } - )* - }; -} - -#[derive(Clone)] -pub struct Matrix<T: Float> { - pub data: Vec<T>, - pub rows: usize, - pub cols: usize, -} - -impl<T: Float> Matrix<T> { - pub fn new(rows: usize, cols: usize, data: Vec<T>) -> Self { - Matrix { data, rows, cols } - } - - pub fn get(&self, row: usize, col: usize) -> T { - self.data[row * self.cols + col] - } - - pub fn set(&mut self, row: usize, col: usize, value: T) { - self.data[row * self.cols + col] = value; - } - - pub fn dot(&self, other: &Self) -> Option<Matrix<T>> { - if self.cols != other.rows { - return None; - } - - let mut result = Matrix::new( - self.rows, - other.cols, - vec![T::zero(); self.rows * other.cols], - ); - for i in 0..self.rows { - for j in 0..other.cols { - let mut sum = T::zero(); - for k in 0..self.cols { - sum += self.get(i, k) * other.get(k, j) ; - } - result.set(i, j, sum); - } - } - Some(result) - } - - pub fn kronecker(&self, other: &Self) -> Matrix<T> { - let new_rows = self.rows * other.rows; - let new_cols = self.cols * other.cols; - - let mut result = Matrix::new(new_rows, new_cols, vec![T::zero(); new_rows * new_cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let self_val = self.get(i, j); - for k in 0..other.rows { - for l in 0..other.cols { - let result_row = i * other.rows + k; - let result_col = j * other.cols + l; - result.set(result_row, result_col, self_val * other.get(k, l)); - } - } - } - } - - result - } - - pub fn transpose(&self) -> Matrix<T> { - let mut result = Matrix::new(self.cols, self.rows, vec![T::zero(); self.cols * self.rows]); - - for i in 0..self.rows { - for j in 0..self.cols { - let value = self.get(i, j); - result.set(j, i, value); - } - } - - result - } - - pub fn add_to(&self, other: &Self) -> Option<Matrix<T>> { - if self.rows != other.rows || self.cols != other.cols { - return None; - } - - let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let sum = self.get(i, j) + other.get(i, j); - result.set(i, j, sum); - } - } - Some(result) - } - - pub fn subtract(&self, other: &Self) -> Option<Matrix<T>> { - if self.rows != other.rows || self.cols != other.cols { - return None; - } - - let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let diff = self.get(i, j) - other.get(i, j); - result.set(i, j, diff); - } - } - Some(result) - } - - pub fn scale(&self, scalar: T) -> Matrix<T> { - let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let scaled_value = self.get(i, j) * scalar; - result.set(i, j, scaled_value); - } - } - result - } -} - -impl<T: Float> ops::Index<(usize, usize)> for Matrix<T> { - type Output = T; - - fn index(&self, index: (usize, usize)) -> &Self::Output { - &self.data[index.0 * self.cols + index.1] - } -} - -impl<T: Float> ops::IndexMut<(usize, usize)> for Matrix<T> { - fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output { - &mut self.data[index.0 * self.cols + index.1] - } -} - -impl<T: Float> ops::AddAssign<&Matrix<T>> for Matrix<T> { - fn add_assign(&mut self, other: &Matrix<T>) { - if let Some(result) = self.add_to(other) { - *self = result; - } - } -} - -impl<T: Float> ops::SubAssign<&Matrix<T>> for Matrix<T> { - fn sub_assign(&mut self, other: &Matrix<T>) { - if let Some(result) = self.subtract(other) { - *self = result; - } - } -} - -impl_matrix_ops! { - Add, add, &Matrix<T>, Option<Matrix<T>>, add_to, - Sub, sub, &Matrix<T>, Option<Matrix<T>>, subtract, - Mul, mul, T, Matrix<T>, scale, - Div, div, T, Matrix<T>, scale, -} - -impl_matrix_ops! { - MulAssign, mul_assign, T, scale, - DivAssign, div_assign, T, scale, -} - -impl<T: Float + fmt::Debug> fmt::Debug for Matrix<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - for i in 0..self.rows { - for j in 0..self.cols { - write!(f, "{:?} ", self.get(i, j))?; - } - writeln!(f)?; - } - Ok(()) - } -} - -impl<T: Float + fmt::Display> fmt::Display for Matrix<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let elements: Vec<String> = self.data.iter().map(ToString::to_string).collect(); - let is_complex = elements.iter().any(|element| element.contains("i")); - - let normalized: Vec<(f64, f64)> = self - .data - .iter() - .map(|element| { - let element_string = element.to_string(); - - if is_complex { - let element_string = element_string.trim_end_matches('i').trim(); - let element_split: Vec<&str> = element_string.split_whitespace().collect(); - let real = element_split[0].parse::<f64>().unwrap(); - let imaginary = element_split - .get(2) - .map_or(0.0, |&s| s.parse::<f64>().unwrap()); - (real, imaginary) - } else { - (element_string.parse::<f64>().unwrap(), 0.0) - } - }) - .collect(); - - let max_widths = normalized - .iter() - .fold((0, 0), |(max_0, max_1), &(real, imag)| { - let new_max_0 = max_0.max(format!("{:.2}", real).len()); - let new_max_1 = if is_complex { - max_1.max(format!("{:.2}", imag.abs()).len()) - } else { - max_1 - }; - (new_max_0, new_max_1) - }); - - let aligned: Vec<String> = normalized - .iter() - .map(|&(real, imag)| { - if is_complex { - format!( - "{:>rewidth$.2} {} {:>imwidth$.2}i", - real, - if imag > 0.0 { "+" } else { "-" }, - imag.abs(), - rewidth = max_widths.0, - imwidth = max_widths.1, - ) - } else { - format!("{:>width$.2}", real, width = max_widths.0) - } - }) - .collect(); - - for i in 0..self.rows { - if i == 0 { - write!(f, "┌")?; - } else if i == self.rows - 1 { - write!(f, "└")?; - } else { - write!(f, "│")?; - } - - for j in 0..self.cols { - write!(f, "{}", aligned[i + j * self.rows])?; - if j != self.cols - 1 { - write!(f, ", ")?; - } - } - - if i == 0 { - write!(f, "┐")?; - } else if i == self.rows - 1 { - write!(f, "┘")?; - } else { - write!(f, "│")?; - } - - if i != self.rows - 1 { - writeln!(f)?; - } - } - - Ok(()) - } -} diff --git a/src/maths/simd.c b/src/maths/simd.c index 9e83dcc..fc25de8 100644 --- a/src/maths/simd.c +++ b/src/maths/simd.c @@ -12,144 +12,144 @@ enum PsiSimdCapability psi_detect_simd(void) { #if defined(__x86_64__) || defined(__i386__) - if (__builtin_cpu_supports("avx512f") && __builtin_cpu_supports("avx512dq")) - return PSI_SIMD_AVX512; - if (__builtin_cpu_supports("avx2") && __builtin_cpu_supports("fma")) - return PSI_SIMD_AVX2; + if (__builtin_cpu_supports("avx512f") && __builtin_cpu_supports("avx512dq")) + return PSI_SIMD_AVX512; + if (__builtin_cpu_supports("avx2") && __builtin_cpu_supports("fma")) + return PSI_SIMD_AVX2; - return PSI_SIMD_NONE; + return PSI_SIMD_NONE; #elif defined(__aarch64__) - return PSI_SIMD_NEON; + return PSI_SIMD_NEON; #else - return PSI_SIMD_NONE; + return PSI_SIMD_NONE; #endif } const char* psi_simd_name(enum PsiSimdCapability cap) { - switch (cap) - { - case PSI_SIMD_NONE: return "Scalar"; - case PSI_SIMD_AVX2: return "AVX2+FMA"; - case PSI_SIMD_AVX512: return "AVX-512"; - case PSI_SIMD_NEON: return "NEON"; - } + switch (cap) + { + case PSI_SIMD_NONE: return "Scalar"; + case PSI_SIMD_AVX2: return "AVX2+FMA"; + case PSI_SIMD_AVX512: return "AVX-512"; + case PSI_SIMD_NEON: return "NEON"; + } - return "Scalar"; + return "Scalar"; } static void apply_pair(struct PsiComplex* state, size_t i, size_t j, struct PsiComplex g00, struct PsiComplex g01, struct PsiComplex g10, struct PsiComplex g11) { - struct PsiComplex s0 = state[i]; - struct PsiComplex s1 = state[j]; + struct PsiComplex s0 = state[i]; + struct PsiComplex s1 = state[j]; - state[i] = psi_new_complex(s0.real * g00.real - s0.imaginary * g00.imaginary + - s1.real * g01.real - s1.imaginary * g01.imaginary, - s0.real * g00.imaginary + s0.imaginary * g00.real + - s1.real * g01.imaginary + s1.imaginary * g01.real); - state[j] = psi_new_complex(s0.real * g10.real - s0.imaginary * g10.imaginary + - s1.real * g11.real - s1.imaginary * g11.imaginary, - s0.real * g10.imaginary + s0.imaginary * g10.real + - s1.real * g11.imaginary + s1.imaginary * g11.real); + state[i] = psi_new_complex(s0.real * g00.real - s0.imaginary * g00.imaginary + + s1.real * g01.real - s1.imaginary * g01.imaginary, + s0.real * g00.imaginary + s0.imaginary * g00.real + + s1.real * g01.imaginary + s1.imaginary * g01.real); + state[j] = psi_new_complex(s0.real * g10.real - s0.imaginary * g10.imaginary + + s1.real * g11.real - s1.imaginary * g11.imaginary, + s0.real * g10.imaginary + s0.imaginary * g10.real + + s1.real * g11.imaginary + s1.imaginary * g11.real); } static void apply_scalar(struct PsiComplex* state, const struct PsiComplex gate[2][2], size_t target, size_t num_qubits) { - size_t target_bit = num_qubits - 1 - target; - size_t step = (size_t)1 << target_bit; - size_t dim = (size_t)1 << num_qubits; + size_t target_bit = num_qubits - 1 - target; + size_t step = (size_t)1 << target_bit; + size_t dim = (size_t)1 << num_qubits; - for (size_t i = 0; i < dim; i++) - { - if ((i >> target_bit) & 1) - continue; + for (size_t i = 0; i < dim; i++) + { + if ((i >> target_bit) & 1) + continue; - apply_pair(state, i, i | step, gate[0][0], gate[0][1], gate[1][0], gate[1][1]); - } + apply_pair(state, i, i | step, gate[0][0], gate[0][1], gate[1][0], gate[1][1]); + } } static size_t (*build_pairs(size_t dim, size_t target_bit, size_t step, size_t* out_count))[2] { - size_t (*pairs)[2] = malloc((dim / 2) * sizeof(*pairs)); - assert(pairs != NULL || dim == 0); + size_t (*pairs)[2] = malloc((dim / 2) * sizeof(*pairs)); + assert(pairs != NULL || dim == 0); - size_t n = 0; - for (size_t i = 0; i < dim; i++) - if (((i >> target_bit) & 1) == 0) - { - pairs[n][0] = i; - pairs[n][1] = i | step; - n++; - } + size_t n = 0; + for (size_t i = 0; i < dim; i++) + if (((i >> target_bit) & 1) == 0) + { + pairs[n][0] = i; + pairs[n][1] = i | step; + n++; + } - *out_count = n; - return pairs; + *out_count = n; + return pairs; } #if defined(__aarch64__) static void apply_neon(struct PsiComplex* state, const struct PsiComplex gate[2][2], size_t target, size_t num_qubits) { - size_t target_bit = num_qubits - 1 - target; - size_t step = (size_t)1 << target_bit; - size_t dim = (size_t)1 << num_qubits; + size_t target_bit = num_qubits - 1 - target; + size_t step = (size_t)1 << target_bit; + size_t dim = (size_t)1 << num_qubits; - struct PsiComplex g00 = gate[0][0]; - struct PsiComplex g01 = gate[0][1]; - struct PsiComplex g10 = gate[1][0]; - struct PsiComplex g11 = gate[1][1]; + struct PsiComplex g00 = gate[0][0]; + struct PsiComplex g01 = gate[0][1]; + struct PsiComplex g10 = gate[1][0]; + struct PsiComplex g11 = gate[1][1]; - size_t np; - size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np); - size_t chunks = np / 2; + size_t np; + size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np); + size_t chunks = np / 2; - float64x2_t g00_re = vdupq_n_f64(g00.real); - float64x2_t g00_im = vdupq_n_f64(g00.imaginary); - float64x2_t g01_re = vdupq_n_f64(g01.real); - float64x2_t g01_im = vdupq_n_f64(g01.imaginary); - float64x2_t g10_re = vdupq_n_f64(g10.real); - float64x2_t g10_im = vdupq_n_f64(g10.imaginary); - float64x2_t g11_re = vdupq_n_f64(g11.real); - float64x2_t g11_im = vdupq_n_f64(g11.imaginary); + float64x2_t g00_re = vdupq_n_f64(g00.real); + float64x2_t g00_im = vdupq_n_f64(g00.imaginary); + float64x2_t g01_re = vdupq_n_f64(g01.real); + float64x2_t g01_im = vdupq_n_f64(g01.imaginary); + float64x2_t g10_re = vdupq_n_f64(g10.real); + float64x2_t g10_im = vdupq_n_f64(g10.imaginary); + float64x2_t g11_re = vdupq_n_f64(g11.real); + float64x2_t g11_im = vdupq_n_f64(g11.imaginary); - for (size_t chunk = 0; chunk < chunks; chunk++) - { - size_t i0 = pairs[chunk * 2][0]; - size_t j0 = pairs[chunk * 2][1]; - size_t i1 = pairs[chunk * 2 + 1][0]; - size_t j1 = pairs[chunk * 2 + 1][1]; + for (size_t chunk = 0; chunk < chunks; chunk++) + { + size_t i0 = pairs[chunk * 2][0]; + size_t j0 = pairs[chunk * 2][1]; + size_t i1 = pairs[chunk * 2 + 1][0]; + size_t j1 = pairs[chunk * 2 + 1][1]; - double s0r[2] = { state[i0].real, state[i1].real }; - double s0i[2] = { state[i0].imaginary, state[i1].imaginary }; - double s1r[2] = { state[j0].real, state[j1].real }; - double s1i[2] = { state[j0].imaginary, state[j1].imaginary }; + double s0r[2] = { state[i0].real, state[i1].real }; + double s0i[2] = { state[i0].imaginary, state[i1].imaginary }; + double s1r[2] = { state[j0].real, state[j1].real }; + double s1i[2] = { state[j0].imaginary, state[j1].imaginary }; - float64x2_t s0_re = vld1q_f64(s0r); - float64x2_t s0_im = vld1q_f64(s0i); - float64x2_t s1_re = vld1q_f64(s1r); - float64x2_t s1_im = vld1q_f64(s1i); + float64x2_t s0_re = vld1q_f64(s0r); + float64x2_t s0_im = vld1q_f64(s0i); + float64x2_t s1_re = vld1q_f64(s1r); + float64x2_t s1_im = vld1q_f64(s1i); - float64x2_t new0_re = vaddq_f64(vfmsq_f64(vmulq_f64(s0_re, g00_re), s0_im, g00_im), - vfmsq_f64(vmulq_f64(s1_re, g01_re), s1_im, g01_im)); - float64x2_t new0_im = vaddq_f64(vfmaq_f64(vmulq_f64(s0_re, g00_im), s0_im, g00_re), - vfmaq_f64(vmulq_f64(s1_re, g01_im), s1_im, g01_re)); - float64x2_t new1_re = vaddq_f64(vfmsq_f64(vmulq_f64(s0_re, g10_re), s0_im, g10_im), - vfmsq_f64(vmulq_f64(s1_re, g11_re), s1_im, g11_im)); - float64x2_t new1_im = vaddq_f64(vfmaq_f64(vmulq_f64(s0_re, g10_im), s0_im, g10_re), - vfmaq_f64(vmulq_f64(s1_re, g11_im), s1_im, g11_re)); + float64x2_t new0_re = vaddq_f64(vfmsq_f64(vmulq_f64(s0_re, g00_re), s0_im, g00_im), + vfmsq_f64(vmulq_f64(s1_re, g01_re), s1_im, g01_im)); + float64x2_t new0_im = vaddq_f64(vfmaq_f64(vmulq_f64(s0_re, g00_im), s0_im, g00_re), + vfmaq_f64(vmulq_f64(s1_re, g01_im), s1_im, g01_re)); + float64x2_t new1_re = vaddq_f64(vfmsq_f64(vmulq_f64(s0_re, g10_re), s0_im, g10_im), + vfmsq_f64(vmulq_f64(s1_re, g11_re), s1_im, g11_im)); + float64x2_t new1_im = vaddq_f64(vfmaq_f64(vmulq_f64(s0_re, g10_im), s0_im, g10_re), + vfmaq_f64(vmulq_f64(s1_re, g11_im), s1_im, g11_re)); - state[i0] = psi_new_complex(vgetq_lane_f64(new0_re, 0), vgetq_lane_f64(new0_im, 0)); - state[j0] = psi_new_complex(vgetq_lane_f64(new1_re, 0), vgetq_lane_f64(new1_im, 0)); - state[i1] = psi_new_complex(vgetq_lane_f64(new0_re, 1), vgetq_lane_f64(new0_im, 1)); - state[j1] = psi_new_complex(vgetq_lane_f64(new1_re, 1), vgetq_lane_f64(new1_im, 1)); - } + state[i0] = psi_new_complex(vgetq_lane_f64(new0_re, 0), vgetq_lane_f64(new0_im, 0)); + state[j0] = psi_new_complex(vgetq_lane_f64(new1_re, 0), vgetq_lane_f64(new1_im, 0)); + state[i1] = psi_new_complex(vgetq_lane_f64(new0_re, 1), vgetq_lane_f64(new0_im, 1)); + state[j1] = psi_new_complex(vgetq_lane_f64(new1_re, 1), vgetq_lane_f64(new1_im, 1)); + } - for (size_t p = chunks * 2; p < np; p++) - apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11); + for (size_t p = chunks * 2; p < np; p++) + apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11); - free(pairs); + free(pairs); } #endif @@ -158,167 +158,167 @@ __attribute__((target("avx2,fma"))) static void apply_avx2(struct PsiComplex* st const struct PsiComplex gate[2][2], size_t target, size_t num_qubits) { - size_t target_bit = num_qubits - 1 - target; - size_t step = (size_t)1 << target_bit; - size_t dim = (size_t)1 << num_qubits; + size_t target_bit = num_qubits - 1 - target; + size_t step = (size_t)1 << target_bit; + size_t dim = (size_t)1 << num_qubits; - struct PsiComplex g00 = gate[0][0]; - struct PsiComplex g01 = gate[0][1]; - struct PsiComplex g10 = gate[1][0]; - struct PsiComplex g11 = gate[1][1]; + struct PsiComplex g00 = gate[0][0]; + struct PsiComplex g01 = gate[0][1]; + struct PsiComplex g10 = gate[1][0]; + struct PsiComplex g11 = gate[1][1]; - size_t np; - size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np); - size_t chunks = np / 2; + size_t np; + size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np); + size_t chunks = np / 2; - for (size_t chunk = 0; chunk < chunks; chunk++) - { - size_t i0 = pairs[chunk * 2][0]; - size_t j0 = pairs[chunk * 2][1]; - size_t i1 = pairs[chunk * 2 + 1][0]; - size_t j1 = pairs[chunk * 2 + 1][1]; + for (size_t chunk = 0; chunk < chunks; chunk++) + { + size_t i0 = pairs[chunk * 2][0]; + size_t j0 = pairs[chunk * 2][1]; + size_t i1 = pairs[chunk * 2 + 1][0]; + size_t j1 = pairs[chunk * 2 + 1][1]; - __m256d s0_re = - _mm256_set_pd(state[j1].real, state[i1].real, state[j0].real, state[i0].real); - __m256d s0_im = _mm256_set_pd(state[j1].imaginary, state[i1].imaginary, state[j0].imaginary, - state[i0].imaginary); + __m256d s0_re = + _mm256_set_pd(state[j1].real, state[i1].real, state[j0].real, state[i0].real); + __m256d s0_im = _mm256_set_pd(state[j1].imaginary, state[i1].imaginary, state[j0].imaginary, + state[i0].imaginary); - __m256d g_re_0 = _mm256_set_pd(g01.real, g00.real, g01.real, g00.real); - __m256d g_im_0 = _mm256_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary); - __m256d g_re_1 = _mm256_set_pd(g11.real, g10.real, g11.real, g10.real); - __m256d g_im_1 = _mm256_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary); + __m256d g_re_0 = _mm256_set_pd(g01.real, g00.real, g01.real, g00.real); + __m256d g_im_0 = _mm256_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary); + __m256d g_re_1 = _mm256_set_pd(g11.real, g10.real, g11.real, g10.real); + __m256d g_im_1 = _mm256_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary); - __m256d prod0_re = _mm256_fmsub_pd(s0_re, g_re_0, _mm256_mul_pd(s0_im, g_im_0)); - __m256d prod0_im = _mm256_fmadd_pd(s0_re, g_im_0, _mm256_mul_pd(s0_im, g_re_0)); - __m256d prod1_re = _mm256_fmsub_pd(s0_re, g_re_1, _mm256_mul_pd(s0_im, g_im_1)); - __m256d prod1_im = _mm256_fmadd_pd(s0_re, g_im_1, _mm256_mul_pd(s0_im, g_re_1)); + __m256d prod0_re = _mm256_fmsub_pd(s0_re, g_re_0, _mm256_mul_pd(s0_im, g_im_0)); + __m256d prod0_im = _mm256_fmadd_pd(s0_re, g_im_0, _mm256_mul_pd(s0_im, g_re_0)); + __m256d prod1_re = _mm256_fmsub_pd(s0_re, g_re_1, _mm256_mul_pd(s0_im, g_im_1)); + __m256d prod1_im = _mm256_fmadd_pd(s0_re, g_im_1, _mm256_mul_pd(s0_im, g_re_1)); - double res0_re[4]; - double res0_im[4]; - double res1_re[4]; - double res1_im[4]; + double res0_re[4]; + double res0_im[4]; + double res1_re[4]; + double res1_im[4]; - _mm256_storeu_pd(res0_re, prod0_re); - _mm256_storeu_pd(res0_im, prod0_im); - _mm256_storeu_pd(res1_re, prod1_re); - _mm256_storeu_pd(res1_im, prod1_im); + _mm256_storeu_pd(res0_re, prod0_re); + _mm256_storeu_pd(res0_im, prod0_im); + _mm256_storeu_pd(res1_re, prod1_re); + _mm256_storeu_pd(res1_im, prod1_im); - state[i0] = psi_new_complex(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); - state[j0] = psi_new_complex(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); - state[i1] = psi_new_complex(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); - state[j1] = psi_new_complex(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); - } + state[i0] = psi_new_complex(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); + state[j0] = psi_new_complex(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); + state[i1] = psi_new_complex(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); + state[j1] = psi_new_complex(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); + } - for (size_t p = chunks * 2; p < np; p++) - apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11); + for (size_t p = chunks * 2; p < np; p++) + apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11); - free(pairs); + free(pairs); } __attribute__((target("avx512f,avx512dq"))) static void apply_avx512(struct PsiComplex* state, const struct PsiComplex gate[2][2], size_t target, size_t num_qubits) { - size_t target_bit = num_qubits - 1 - target; - size_t step = (size_t)1 << target_bit; - size_t dim = (size_t)1 << num_qubits; + size_t target_bit = num_qubits - 1 - target; + size_t step = (size_t)1 << target_bit; + size_t dim = (size_t)1 << num_qubits; - struct PsiComplex g00 = gate[0][0]; - struct PsiComplex g01 = gate[0][1]; - struct PsiComplex g10 = gate[1][0]; - struct PsiComplex g11 = gate[1][1]; + struct PsiComplex g00 = gate[0][0]; + struct PsiComplex g01 = gate[0][1]; + struct PsiComplex g10 = gate[1][0]; + struct PsiComplex g11 = gate[1][1]; - size_t np; - size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np); - size_t chunks = np / 4; + size_t np; + size_t (*pairs)[2] = build_pairs(dim, target_bit, step, &np); + size_t chunks = np / 4; - for (size_t chunk = 0; chunk < chunks; chunk++) - { - size_t base = chunk * 4; - size_t i0 = pairs[base][0], j0 = pairs[base][1]; - size_t i1 = pairs[base + 1][0], j1 = pairs[base + 1][1]; - size_t i2 = pairs[base + 2][0], j2 = pairs[base + 2][1]; - size_t i3 = pairs[base + 3][0], j3 = pairs[base + 3][1]; + for (size_t chunk = 0; chunk < chunks; chunk++) + { + size_t base = chunk * 4; + size_t i0 = pairs[base][0], j0 = pairs[base][1]; + size_t i1 = pairs[base + 1][0], j1 = pairs[base + 1][1]; + size_t i2 = pairs[base + 2][0], j2 = pairs[base + 2][1]; + size_t i3 = pairs[base + 3][0], j3 = pairs[base + 3][1]; - __m512d s0_re = - _mm512_set_pd(state[j3].real, state[i3].real, state[j2].real, state[i2].real, - state[j1].real, state[i1].real, state[j0].real, state[i0].real); - __m512d s0_im = _mm512_set_pd(state[j3].imaginary, state[i3].imaginary, state[j2].imaginary, - state[i2].imaginary, state[j1].imaginary, state[i1].imaginary, - state[j0].imaginary, state[i0].imaginary); + __m512d s0_re = + _mm512_set_pd(state[j3].real, state[i3].real, state[j2].real, state[i2].real, + state[j1].real, state[i1].real, state[j0].real, state[i0].real); + __m512d s0_im = _mm512_set_pd(state[j3].imaginary, state[i3].imaginary, state[j2].imaginary, + state[i2].imaginary, state[j1].imaginary, state[i1].imaginary, + state[j0].imaginary, state[i0].imaginary); - __m512d g_re_0 = _mm512_set_pd(g01.real, g00.real, g01.real, g00.real, g01.real, g00.real, - g01.real, g00.real); - __m512d g_im_0 = _mm512_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary, - g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary); - __m512d g_re_1 = _mm512_set_pd(g11.real, g10.real, g11.real, g10.real, g11.real, g10.real, - g11.real, g10.real); - __m512d g_im_1 = _mm512_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary, - g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary); + __m512d g_re_0 = _mm512_set_pd(g01.real, g00.real, g01.real, g00.real, g01.real, g00.real, + g01.real, g00.real); + __m512d g_im_0 = _mm512_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary, + g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary); + __m512d g_re_1 = _mm512_set_pd(g11.real, g10.real, g11.real, g10.real, g11.real, g10.real, + g11.real, g10.real); + __m512d g_im_1 = _mm512_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary, + g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary); - __m512d prod0_re = _mm512_fmsub_pd(s0_re, g_re_0, _mm512_mul_pd(s0_im, g_im_0)); - __m512d prod0_im = _mm512_fmadd_pd(s0_re, g_im_0, _mm512_mul_pd(s0_im, g_re_0)); - __m512d prod1_re = _mm512_fmsub_pd(s0_re, g_re_1, _mm512_mul_pd(s0_im, g_im_1)); - __m512d prod1_im = _mm512_fmadd_pd(s0_re, g_im_1, _mm512_mul_pd(s0_im, g_re_1)); + __m512d prod0_re = _mm512_fmsub_pd(s0_re, g_re_0, _mm512_mul_pd(s0_im, g_im_0)); + __m512d prod0_im = _mm512_fmadd_pd(s0_re, g_im_0, _mm512_mul_pd(s0_im, g_re_0)); + __m512d prod1_re = _mm512_fmsub_pd(s0_re, g_re_1, _mm512_mul_pd(s0_im, g_im_1)); + __m512d prod1_im = _mm512_fmadd_pd(s0_re, g_im_1, _mm512_mul_pd(s0_im, g_re_1)); - double res0_re[8]; - double res0_im[8]; - double res1_re[8]; - double res1_im[8]; + double res0_re[8]; + double res0_im[8]; + double res1_re[8]; + double res1_im[8]; - _mm512_storeu_pd(res0_re, prod0_re); - _mm512_storeu_pd(res0_im, prod0_im); - _mm512_storeu_pd(res1_re, prod1_re); - _mm512_storeu_pd(res1_im, prod1_im); + _mm512_storeu_pd(res0_re, prod0_re); + _mm512_storeu_pd(res0_im, prod0_im); + _mm512_storeu_pd(res1_re, prod1_re); + _mm512_storeu_pd(res1_im, prod1_im); - state[i0] = psi_new_complex(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); - state[j0] = psi_new_complex(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); - state[i1] = psi_new_complex(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); - state[j1] = psi_new_complex(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); - state[i2] = psi_new_complex(res0_re[4] + res0_re[5], res0_im[4] + res0_im[5]); - state[j2] = psi_new_complex(res1_re[4] + res1_re[5], res1_im[4] + res1_im[5]); - state[i3] = psi_new_complex(res0_re[6] + res0_re[7], res0_im[6] + res0_im[7]); - state[j3] = psi_new_complex(res1_re[6] + res1_re[7], res1_im[6] + res1_im[7]); - } + state[i0] = psi_new_complex(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); + state[j0] = psi_new_complex(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); + state[i1] = psi_new_complex(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); + state[j1] = psi_new_complex(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); + state[i2] = psi_new_complex(res0_re[4] + res0_re[5], res0_im[4] + res0_im[5]); + state[j2] = psi_new_complex(res1_re[4] + res1_re[5], res1_im[4] + res1_im[5]); + state[i3] = psi_new_complex(res0_re[6] + res0_re[7], res0_im[6] + res0_im[7]); + state[j3] = psi_new_complex(res1_re[6] + res1_re[7], res1_im[6] + res1_im[7]); + } - for (size_t p = chunks * 4; p < np; p++) - apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11); + for (size_t p = chunks * 4; p < np; p++) + apply_pair(state, pairs[p][0], pairs[p][1], g00, g01, g10, g11); - free(pairs); + free(pairs); } #endif void psi_apply_single_qubit_gate_simd(struct PsiComplex* state, const struct PsiComplex gate[2][2], size_t target, size_t num_qubits) { - enum PsiSimdCapability cap = psi_detect_simd(); + enum PsiSimdCapability cap = psi_detect_simd(); #if defined(__x86_64__) || defined(__i386__) - if (cap == PSI_SIMD_AVX512) - { - apply_avx512(state, gate, target, num_qubits); - return; - } - if (cap == PSI_SIMD_AVX2) - { - apply_avx2(state, gate, target, num_qubits); - return; - } + if (cap == PSI_SIMD_AVX512) + { + apply_avx512(state, gate, target, num_qubits); + return; + } + if (cap == PSI_SIMD_AVX2) + { + apply_avx2(state, gate, target, num_qubits); + return; + } #elif defined(__aarch64__) - if (cap == PSI_SIMD_NEON) - { - apply_neon(state, gate, target, num_qubits); - return; - } + if (cap == PSI_SIMD_NEON) + { + apply_neon(state, gate, target, num_qubits); + return; + } #endif - (void)cap; - apply_scalar(state, gate, target, num_qubits); + (void)cap; + apply_scalar(state, gate, target, num_qubits); } void psi_apply_single_qubit_gate_simd_parallel(struct PsiComplex* state, const struct PsiComplex gate[2][2], size_t target, size_t num_qubits) { - apply_scalar(state, gate, target, num_qubits); + apply_scalar(state, gate, target, num_qubits); } diff --git a/src/maths/vector.c b/src/maths/vector.c index 1c5c4f9..5cd49e0 100644 --- a/src/maths/vector.c +++ b/src/maths/vector.c @@ -6,155 +6,155 @@ static enum PsiVectorKind flip_kind(enum PsiVectorKind kind) { - if (kind == PSI_ROW_VECTOR) - return PSI_COLUMN_VECTOR; + if (kind == PSI_ROW_VECTOR) + return PSI_COLUMN_VECTOR; - return PSI_ROW_VECTOR; + return PSI_ROW_VECTOR; } struct PsiVector psi_new_vector(size_t size, enum PsiVectorKind kind) { - struct PsiComplex* data = calloc(size, sizeof(struct PsiComplex)); - assert(data != NULL || size == 0); + struct PsiComplex* data = calloc(size, sizeof(struct PsiComplex)); + assert(data != NULL || size == 0); - return (struct PsiVector){ - data, - size, - kind, - }; + return (struct PsiVector){ + data, + size, + kind, + }; } struct PsiVector psi_new_vector_from(const struct PsiComplex* data, size_t size, enum PsiVectorKind kind) { - struct PsiVector v = psi_new_vector(size, kind); - memcpy(v.data, data, size * sizeof(struct PsiComplex)); + struct PsiVector v = psi_new_vector(size, kind); + memcpy(v.data, data, size * sizeof(struct PsiComplex)); - return v; + return v; } struct PsiVector psi_clone_vector(struct PsiVector v) { - return psi_new_vector_from(v.data, v.size, v.kind); + return psi_new_vector_from(v.data, v.size, v.kind); } void psi_free_vector(struct PsiVector* v) { - free(v->data); - v->data = NULL; - v->size = 0; + free(v->data); + v->data = NULL; + v->size = 0; } struct PsiComplex psi_get_vector(struct PsiVector v, size_t index) { - assert(index < v.size); - return v.data[index]; + assert(index < v.size); + return v.data[index]; } void psi_set_vector(struct PsiVector* v, size_t index, struct PsiComplex value) { - assert(index < v->size); - v->data[index] = value; + assert(index < v->size); + v->data[index] = value; } struct PsiComplex psi_dot_vector(struct PsiVector a, struct PsiVector b) { - assert(a.size == b.size); + assert(a.size == b.size); - struct PsiComplex sum = psi_new_complex(0.0, 0.0); - for (size_t i = 0; i < a.size; i++) - sum = psi_add_complex(sum, psi_mul_complex(a.data[i], b.data[i])); + struct PsiComplex sum = psi_new_complex(0.0, 0.0); + for (size_t i = 0; i < a.size; i++) + sum = psi_add_complex(sum, psi_mul_complex(a.data[i], b.data[i])); - return sum; + return sum; } struct PsiComplex psi_norm_vector(struct PsiVector v) { - struct PsiComplex sum = psi_new_complex(0.0, 0.0); - for (size_t i = 0; i < v.size; i++) - sum = psi_add_complex(sum, psi_mul_complex(v.data[i], v.data[i])); + struct PsiComplex sum = psi_new_complex(0.0, 0.0); + for (size_t i = 0; i < v.size; i++) + sum = psi_add_complex(sum, psi_mul_complex(v.data[i], v.data[i])); - return psi_sqrt_complex(sum); + return psi_sqrt_complex(sum); } struct PsiComplex psi_sum_vector(struct PsiVector v) { - struct PsiComplex sum = psi_new_complex(0.0, 0.0); - for (size_t i = 0; i < v.size; i++) - sum = psi_add_complex(sum, v.data[i]); + struct PsiComplex sum = psi_new_complex(0.0, 0.0); + for (size_t i = 0; i < v.size; i++) + sum = psi_add_complex(sum, v.data[i]); - return sum; + return sum; } static int less_than(struct PsiComplex a, struct PsiComplex b) { - if (a.real != b.real) - return a.real < b.real; + if (a.real != b.real) + return a.real < b.real; - return a.imaginary < b.imaginary; + return a.imaginary < b.imaginary; } struct PsiComplex psi_max_vector(struct PsiVector v) { - if (v.size == 0) - return psi_new_complex(0.0, 0.0); + if (v.size == 0) + return psi_new_complex(0.0, 0.0); - struct PsiComplex best = v.data[0]; - for (size_t i = 1; i < v.size; i++) - if (less_than(best, v.data[i])) - best = v.data[i]; + struct PsiComplex best = v.data[0]; + for (size_t i = 1; i < v.size; i++) + if (less_than(best, v.data[i])) + best = v.data[i]; - return best; + return best; } struct PsiComplex psi_min_vector(struct PsiVector v) { - if (v.size == 0) - return psi_new_complex(0.0, 0.0); + if (v.size == 0) + return psi_new_complex(0.0, 0.0); - struct PsiComplex best = v.data[0]; - for (size_t i = 1; i < v.size; i++) - if (less_than(v.data[i], best)) - best = v.data[i]; + struct PsiComplex best = v.data[0]; + for (size_t i = 1; i < v.size; i++) + if (less_than(v.data[i], best)) + best = v.data[i]; - return best; + return best; } struct PsiVector psi_add_vector(struct PsiVector a, struct PsiVector b) { - assert(a.size == b.size); + assert(a.size == b.size); - struct PsiVector result = psi_new_vector(a.size, a.kind); - for (size_t i = 0; i < a.size; i++) - result.data[i] = psi_add_complex(a.data[i], b.data[i]); + struct PsiVector result = psi_new_vector(a.size, a.kind); + for (size_t i = 0; i < a.size; i++) + result.data[i] = psi_add_complex(a.data[i], b.data[i]); - return result; + return result; } struct PsiVector psi_sub_vector(struct PsiVector a, struct PsiVector b) { - assert(a.size == b.size); + assert(a.size == b.size); - struct PsiVector result = psi_new_vector(a.size, a.kind); - for (size_t i = 0; i < a.size; i++) - result.data[i] = psi_sub_complex(a.data[i], b.data[i]); + struct PsiVector result = psi_new_vector(a.size, a.kind); + for (size_t i = 0; i < a.size; i++) + result.data[i] = psi_sub_complex(a.data[i], b.data[i]); - return result; + return result; } struct PsiVector psi_scale_vector(struct PsiVector v, struct PsiComplex scalar) { - struct PsiVector result = psi_new_vector(v.size, v.kind); - for (size_t i = 0; i < v.size; i++) - result.data[i] = psi_mul_complex(v.data[i], scalar); + struct PsiVector result = psi_new_vector(v.size, v.kind); + for (size_t i = 0; i < v.size; i++) + result.data[i] = psi_mul_complex(v.data[i], scalar); - return result; + return result; } struct PsiVector psi_transpose_vector(struct PsiVector v) { - struct PsiVector result = psi_clone_vector(v); - result.kind = flip_kind(v.kind); + struct PsiVector result = psi_clone_vector(v); + result.kind = flip_kind(v.kind); - return result; + return result; } diff --git a/src/maths/vector.rs b/src/maths/vector.rs deleted file mode 100644 index 11f3d29..0000000 --- a/src/maths/vector.rs +++ /dev/null @@ -1,258 +0,0 @@ -use super::{Float, Matrix}; -use core::{fmt, ops}; - -#[macro_export] -macro_rules! row_vector { - ($($x:expr),*) => { - RowVector::new(vec![$($x),*]) - }; - ($($x:expr,)*) => { - RowVector::new(vec![$($x),*]) - }; -} - -#[macro_export] -macro_rules! column_vector { - ($($x:expr),*) => { - ColumnVector::new(vec![$($x),*]) - }; - ($($x:expr,)*) => { - ColumnVector::new(vec![$($x),*]) - }; -} - -pub trait Vector<T: Float> { - fn new(data: Vec<T>) -> Self; - fn get(&self, index: usize) -> T; - fn set(&mut self, index: usize, value: T); - fn size(&self) -> usize; - - fn dot(&self, other: &Self) -> T; - fn norm(&self) -> T; - - fn max(&self) -> T; - fn min(&self) -> T; - fn sum(&self) -> T; - - fn from_matrix(matrix: &Matrix<T>) -> Self; -} - -pub trait VectorMatrix<T: Float> { - fn to_matrix(&self) -> Matrix<T>; -} - -#[derive(Clone)] -pub struct VectorImpl<T: Float, const ROWS: usize, const COLS: usize>(Vec<T>); -pub type RowVector<T> = VectorImpl<T, 1, 0>; -pub type ColumnVector<T> = VectorImpl<T, 0, 1>; - -impl<T: Float> ColumnVector<T> { - pub fn mul_matrix(&self, matrix: &Matrix<T>) -> Option<ColumnVector<T>> { - if matrix.cols != self.size() { - return None; - } - - let mut result = ColumnVector::new(vec![T::zero(); matrix.rows]); - - for i in 0..matrix.rows { - let mut sum = T::zero(); - for j in 0..matrix.cols { - sum += matrix.get(i, j) * self.get(j) ; - } - result.set(i, sum); - } - - Some(result) - } - - pub fn transpose(&self) -> RowVector<T> { - RowVector::new(self.0.clone()) - } -} - -impl<T: Float> RowVector<T> { - pub fn mul_matrix(&self, matrix: &Matrix<T>) -> Option<RowVector<T>> { - if self.size() != matrix.rows { - return None; - } - - let mut result = RowVector::new(vec![T::zero(); matrix.cols]); - - for j in 0..matrix.cols { - let mut sum = T::zero(); - for i in 0..matrix.rows { - sum += self.get(i) * matrix.get(i, j) ; - } - result.set(j, sum); - } - - Some(result) - } - - pub fn transpose(&self) -> ColumnVector<T> { - ColumnVector::new(self.0.clone()) - } -} - -impl<T: Float> VectorMatrix<T> for RowVector<T> { - fn to_matrix(&self) -> Matrix<T> { - Matrix::new(1, self.size(), self.0.clone()) - } -} - -impl<T: Float> VectorMatrix<T> for ColumnVector<T> { - fn to_matrix(&self) -> Matrix<T> { - Matrix::new(self.size(), 1, self.0.clone()) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> Vector<T> for VectorImpl<T, ROWS, COLS> { - fn from_matrix(matrix: &Matrix<T>) -> Self { - Self::new(matrix.data.clone()) - } - - fn new(data: Vec<T>) -> Self { - Self(data) - } - - fn get(&self, index: usize) -> T { - self.0[index] - } - - fn set(&mut self, index: usize, value: T) { - self.0[index] = value; - } - - fn size(&self) -> usize { - self.0.len() - } - - fn dot(&self, other: &Self) -> T { - self.0 - .iter() - .zip(other.0.iter()) - .map(|(a, b)| *a * *b) - .fold(T::zero(), |acc, x| acc + x) - } - - fn norm(&self) -> T { - self.0 - .iter() - .map(|x| *x * *x) - .fold(T::zero(), |acc, x| acc + x) - .sqrt() - } - - fn max(&self) -> T { - *self - .0 - .iter() - .max_by(|a, b| a.partial_cmp(b).unwrap()) - .unwrap_or(&T::zero()) - } - - fn min(&self) -> T { - *self - .0 - .iter() - .min_by(|a, b| a.partial_cmp(b).unwrap()) - .unwrap_or(&T::zero()) - } - - fn sum(&self) -> T { - self.0.iter().fold(T::zero(), |acc, x| acc + *x) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> VectorImpl<T, ROWS, COLS> { - pub fn add_to(&self, other: &Self) -> Option<VectorImpl<T, ROWS, COLS>> { - if self.size() != other.size() { - return None; - } - - let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]); - - for i in 0..self.size() { - let sum = self.get(i) + other.get(i); - result.set(i, sum); - } - - Some(result) - } - - pub fn subtract(&self, other: &Self) -> Option<VectorImpl<T, ROWS, COLS>> { - if self.size() != other.size() { - return None; - } - - let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]); - - for i in 0..self.size() { - let sum = self.get(i) - other.get(i); - result.set(i, sum); - } - - Some(result) - } - - pub fn scale(&self, scalar: T) -> VectorImpl<T, ROWS, COLS> { - let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]); - - for i in 0..self.size() { - let product = self.get(i) * scalar; - result.set(i, product); - } - - result - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::Index<usize> - for VectorImpl<T, ROWS, COLS> -{ - type Output = T; - - fn index(&self, index: usize) -> &Self::Output { - &self.0[index] - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::IndexMut<usize> - for VectorImpl<T, ROWS, COLS> -{ - fn index_mut(&mut self, index: usize) -> &mut Self::Output { - &mut self.0[index] - } -} - -impl<T: Float + fmt::Debug> fmt::Debug for RowVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "RowVector({:?})", self.0) - } -} - -impl<T: Float + fmt::Debug> fmt::Debug for ColumnVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "ColumnVector({:?})", self.0) - } -} - -impl<T: Float + fmt::Display> fmt::Display for RowVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "[{}]", - self.0 - .iter() - .map(|x| x.to_string()) - .collect::<Vec<String>>() - .join(", ") - ) - } -} - -impl<T: Float + fmt::Display> fmt::Display for ColumnVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "{}", self.to_matrix()) - } -} |