From f5650d8e23e7681ba452a7c10baf1d2e0ebc5e05 Mon Sep 17 00:00:00 2001 From: Ben Edgington Date: Thu, 4 Feb 2021 14:15:15 +0000 Subject: [PATCH] Implement polynomial division --- src/Makefile | 4 +-- src/poly.c | 57 ++++++++++++++++++++++++++++++ src/poly.h | 20 +++++++++++ src/poly_test.c | 93 +++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 172 insertions(+), 2 deletions(-) create mode 100644 src/poly.c create mode 100644 src/poly.h create mode 100644 src/poly_test.c diff --git a/src/Makefile b/src/Makefile index bde9b84..cbd9a2d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,5 +1,5 @@ -TESTS = blst_util_test fft_common_test fft_fr_test fft_g1_test -LIB_SRC = blst_util.c fft_common.c fft_fr.c fft_g1.c +TESTS = blst_util_test fft_common_test fft_fr_test fft_g1_test poly_test +LIB_SRC = poly.c blst_util.c fft_common.c fft_fr.c fft_g1.c LIB_OBJ = $(LIB_SRC:.c=.o) CFLAGS = diff --git a/src/poly.c b/src/poly.c new file mode 100644 index 0000000..2724821 --- /dev/null +++ b/src/poly.c @@ -0,0 +1,57 @@ +/* + * Copyright 2021 Benjamin Edgington + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "poly.h" + +static void poly_factor_div(blst_fr *out, const blst_fr *a, const blst_fr *b) { + blst_fr_eucl_inverse(out, b); + blst_fr_mul(out, out, a); +} + +// Call this to find out how much space to allocate for the result +uint64_t poly_long_div_length(const uint64_t len_dividend, const uint64_t len_divisor) { + return len_dividend - len_divisor + 1; +} + +// `out` must have been pre-allocated to the correct size, and the length is provided +// as a check +C_KZG_RET poly_long_div(blst_fr *out, const uint64_t len_out, const blst_fr *dividend, const uint64_t len_dividend, const blst_fr *divisor, const uint64_t len_divisor) { + uint64_t a_pos = len_dividend - 1; + uint64_t b_pos = len_divisor - 1; + uint64_t diff = a_pos - b_pos; + blst_fr a[len_dividend]; + + ASSERT(len_out == diff + 1, C_KZG_BADARGS); + + for (uint64_t i = 0; i < len_dividend; i++) { + a[i] = dividend[i]; + } + + while (true) { + poly_factor_div(&out[diff], &a[a_pos], &divisor[b_pos]); + for (uint64_t i = 0; i <= b_pos; i++) { + blst_fr tmp; + // a[diff + i] -= b[i] * quot + blst_fr_mul(&tmp, &out[diff], &divisor[i]); + blst_fr_sub(&a[diff + i], &a[diff + i], &tmp); + } + if (diff == 0) break; + --diff; + --a_pos; + } + + return C_KZG_SUCCESS; +} diff --git a/src/poly.h b/src/poly.h new file mode 100644 index 0000000..07838f2 --- /dev/null +++ b/src/poly.h @@ -0,0 +1,20 @@ +/* + * Copyright 2021 Benjamin Edgington + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "c_kzg.h" + +uint64_t poly_long_div_length(const uint64_t len_dividend, const uint64_t len_divisor); +C_KZG_RET poly_long_div(blst_fr *out, const uint64_t len_out, const blst_fr *dividend, const uint64_t len_dividend, const blst_fr *divisor, const uint64_t len_divisor); diff --git a/src/poly_test.c b/src/poly_test.c new file mode 100644 index 0000000..d8da6b6 --- /dev/null +++ b/src/poly_test.c @@ -0,0 +1,93 @@ +/* + * Copyright 2021 Benjamin Edgington + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "../inc/acutest.h" +#include "debug_util.h" +#include "blst_util.h" +#include "poly.h" + +void poly_div_length() { + TEST_CHECK(3 == poly_long_div_length(4, 2)); +} + +void poly_div_0() { + // Calculate (x^2 - 1) / (x + 1) = x - 1 + blst_fr dividend[3]; + blst_fr divisor[2]; + blst_fr expected[2], actual[2]; + + // Set up dividend + fr_from_uint64(÷nd[0], 1); + fr_negate(÷nd[0], ÷nd[0]); + fr_from_uint64(÷nd[1], 0); + fr_from_uint64(÷nd[2], 1); + + // Set up divisor + fr_from_uint64(&divisor[0], 1); + fr_from_uint64(&divisor[1], 1); + + // Set up result + fr_from_uint64(&expected[0], 1); + fr_negate(&expected[0], &expected[0]); + fr_from_uint64(&expected[1], 1); + + TEST_CHECK(poly_long_div(actual, 2, dividend, 3, divisor, 2) == C_KZG_SUCCESS); + TEST_CHECK(fr_equal(expected + 0, actual + 0)); + TEST_CHECK(fr_equal(expected + 1, actual + 1)); +} + +void poly_div_1() { + // Calculate (12x^3 - 11x^2 + 9x + 18) / (4x + 3) = 3x^2 - 5x + 6 + blst_fr dividend[4]; + blst_fr divisor[2]; + blst_fr expected[3], actual[3]; + + // Set up dividend + fr_from_uint64(÷nd[0], 18); + fr_from_uint64(÷nd[1], 9); + fr_from_uint64(÷nd[2], 11); + fr_negate(÷nd[2], ÷nd[2]); + fr_from_uint64(÷nd[3], 12); + + // Set up divisor + fr_from_uint64(&divisor[0], 3); + fr_from_uint64(&divisor[1], 4); + + // Set up result + fr_from_uint64(&expected[0], 6); + fr_from_uint64(&expected[1], 5); + fr_negate(&expected[1], &expected[1]); + fr_from_uint64(&expected[2], 3); + + TEST_CHECK(poly_long_div(actual, 3, dividend, 4, divisor, 2) == C_KZG_SUCCESS); + TEST_CHECK(fr_equal(expected + 0, actual + 0)); + TEST_CHECK(fr_equal(expected + 1, actual + 1)); + TEST_CHECK(fr_equal(expected + 2, actual + 2)); +} + +void poly_wrong_size(void) { + blst_fr dividend[1], divisor[1], result[1]; + TEST_CHECK(poly_long_div(result, 5, dividend, 20, divisor, 7) == C_KZG_BADARGS); +} + +TEST_LIST = + { + {"poly_div_length", poly_div_length}, + {"poly_div_0", poly_div_0}, + {"poly_div_1", poly_div_1}, + {"poly_wrong_size", poly_wrong_size}, + { NULL, NULL } /* zero record marks the end of the list */ + };