Implement linear combinations of G1 points

src/blst_util.c

@@ -15,6 +15,7 @@
  */
 
 #include "blst_util.h"
+#include "debug_util.h"
 
 bool fr_is_one(const blst_fr *fr_p) {
     uint64_t a[4];
@@ -22,18 +23,19 @@ bool fr_is_one(const blst_fr *fr_p) {
     return a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 0;
 }
 
-void fr_from_uint64(blst_fr *a, uint64_t n) {
+void fr_from_uint64(blst_fr *a, const uint64_t n) {
     uint64_t vals[] = {n, 0, 0, 0};
     blst_fr_from_uint64(a, vals);
 }
 
-bool fr_equal(blst_fr *aa, blst_fr *bb) {
+bool fr_equal(const blst_fr *aa, const blst_fr *bb) {
     uint64_t a[4], b[4];
     blst_uint64_from_fr(a, aa);
     blst_uint64_from_fr(b, bb);
     return a[0] == b[0] && a[1] == b[1] && a[2] == b[2] && a[3] == b[3];
 }
 
+// TODO: Is there really no better way to do this?
 void p1_mul(blst_p1 *out, const blst_p1 *a, const blst_fr *b) {
     blst_scalar s;
     blst_scalar_from_fr(&s, b);
@@ -45,3 +47,13 @@ void p1_sub(blst_p1 *out, const blst_p1 *a, const blst_p1 *b) {
     blst_p1_cneg(&bneg, true);
     blst_p1_add_or_double(out, a, &bneg);
 }
+
+// TODO: would be good to have an optimised multiexp for this
+void linear_combination_g1(blst_p1 *out, const blst_p1 *p, const blst_fr *coeffs, const uint64_t len) {
+    blst_p1 tmp;
+    blst_p1_from_affine(out, &identity_g1_affine);
+    for (uint64_t i = 0; i < len; i++) {
+        p1_mul(&tmp, &p[i], &coeffs[i]);
+        blst_p1_add_or_double(out, out, &tmp);
+    }
+}

src/blst_util.h

@@ -20,8 +20,12 @@
 static const blst_fr one =
     {0x00000001fffffffeL, 0x5884b7fa00034802L, 0x998c4fefecbc4ff5L, 0x1824b159acc5056fL};
 
+// The G1 identity/infinity in affine representation
+static const blst_p1_affine identity_g1_affine = {{0,0,0,0,0,0},{0,0,0,0,0,0}};
+
 bool fr_is_one(const blst_fr *fr_p);
-void fr_from_uint64(blst_fr *a, uint64_t n);
-bool fr_equal(blst_fr *aa, blst_fr *bb);
+void fr_from_uint64(blst_fr *a, const uint64_t n);
+bool fr_equal(const blst_fr *aa, const blst_fr *bb);
 void p1_mul(blst_p1 *out, const blst_p1 *a, const blst_fr *b);
 void p1_sub(blst_p1 *out, const blst_p1 *a, const blst_p1 *b);
+void linear_combination_g1(blst_p1 *out, const blst_p1 *p, const blst_fr *coeffs, const uint64_t len);

src/blst_util_test.c

@@ -70,7 +70,30 @@ void p1_sub_works(void) {
     TEST_CHECK(blst_p1_is_equal(&tmp, &res));
 }
 
+void identity_g1_is_infinity(void) {
+    blst_p1 identity_g1;
+    blst_p1_from_affine(&identity_g1, &identity_g1_affine);
+    TEST_CHECK(blst_p1_is_inf(&identity_g1));
+}
 
+void g1_linear_combination(void) {
+    int len = 255;
+    blst_fr coeffs[len], tmp;
+    blst_p1 p[len], res, exp, g1_gen;
+    for (int i = 0; i < len; i++) {
+        fr_from_uint64(coeffs + i, i + 1);
+        blst_p1_from_affine(p + i, &BLS12_381_G1);
+    }
+
+    // Expected result
+    fr_from_uint64(&tmp, len * (len + 1) / 2);
+    blst_p1_from_affine(&g1_gen, &BLS12_381_G1);
+    p1_mul(&exp, &g1_gen, &tmp);
+
+    // Test result
+    linear_combination_g1(&res, p, coeffs, len);
+    TEST_CHECK(blst_p1_is_equal(&exp, &res));
+}
 
 TEST_LIST =
     {
@@ -79,5 +102,7 @@ TEST_LIST =
      {"fr_equal_works", fr_equal_works},
      {"p1_mul_works", p1_mul_works},
      {"p1_sub_works", p1_sub_works},
+     {"identity_g1_is_infinity", identity_g1_is_infinity},
+     {"g1_linear_combination", g1_linear_combination},
      { NULL, NULL }     /* zero record marks the end of the list */
     };

src/debug_util.c

@@ -45,6 +45,15 @@ void print_fr(const blst_fr *a) {
     print_bytes_as_hex_le(b.b, 0, 32);
 }
 
+//
+// Fp Utilities
+//
+
+void print_limbs(const blst_fp *fp) {
+    printf("(%08lx, %08lx, %08lx, %08lx, %08lx, %08lx)",
+           fp->l[0], fp->l[1], fp->l[2], fp->l[3], fp->l[4], fp->l[5]);
+}
+
 //
 // G1 and G2 utilities
 //
@@ -57,7 +66,7 @@ void print_p1_bytes(byte p1[96]) {
     printf("]\n");
 }
 
-/* "Pretty" print a point in G1 */
+/* "Pretty" print serialisation of a point in G1 */
 void print_p1(const blst_p1 *p1) {
     byte *p1_bytes = (byte *)malloc(96);
     blst_p1_serialize(p1_bytes, p1);
@@ -65,7 +74,7 @@ void print_p1(const blst_p1 *p1) {
     free(p1_bytes);
 }
 
-/* "Pretty" print an affine point in G1 */
+/* "Pretty" print serialisation of an affine point in G1 */
 void print_p1_affine(const blst_p1_affine *p1) {
     byte *p1_bytes = (byte *)malloc(96);
     blst_p1_affine_serialize(p1_bytes, p1);
@@ -73,7 +82,27 @@ void print_p1_affine(const blst_p1_affine *p1) {
     free(p1_bytes);
 }
 
-/* "Pretty" print an affine point in G2 */
+/* "Pretty" print internals of a point in G1 */
+void print_p1_limbs(const blst_p1 *p1) {
+    printf("x = ");
+    print_limbs(&p1->x);
+    printf(", y = ");
+    print_limbs(&p1->y);
+    printf(", z = ");
+    print_limbs(&p1->z);
+    printf("\n");
+}
+
+/* "Pretty" print internals of an affine point in G1 */
+void print_p1_affine_limbs(const blst_p1_affine *p1) {
+    printf("x = ");
+    print_limbs(&p1->x);
+    printf(", y = ");
+    print_limbs(&p1->y);
+    printf("\n");
+}
+
+/* "Pretty" print serialisation of an affine point in G2 */
 void print_p2_affine(const blst_p2_affine *p2) {
     byte *p2_hex = (byte *)malloc(192);
     blst_p2_affine_serialize(p2_hex, p2);

src/debug_util.h

@@ -16,6 +16,7 @@
 
 #include <stdlib.h>
 #include <stdio.h>
+#include <assert.h>
 #include "c_kzg.h"
 
 // General Utilities
@@ -25,8 +26,13 @@ void print_bytes_as_hex_le(byte *bytes, int start, int len);
 // Fr utilities
 void print_fr(const blst_fr *a);
 
+// Fp Utilities
+void print_limbs(const blst_fp *fp);
+
 // G1 and G2 utilities
 void print_p1_bytes(byte p1[96]);
 void print_p1(const blst_p1 *p1);
 void print_p1_affine(const blst_p1_affine *p1);
+void print_p1_limbs(const blst_p1 *p1);
+void print_p1_affine_limbs(const blst_p1_affine *p1);
 void print_p2_affine(const blst_p2_affine *p2);