crypto: much faster pippenger/bucket-style multi-scalar multiplication for range proofs + tests

src/crypto/crypto-sugar.h

@@ -764,6 +764,24 @@ namespace crypto
       return *this;
     }
 
+    point_t modify_mul_pow_2(size_t power)
+    {
+      if (power > 0)
+      {
+        ge_p1p1 p1;
+        ge_p2 p2;
+        ge_p3_to_p2(&p2, &m_p3);
+        for (size_t i = 1; i < power; ++i)
+        {
+          ge_p2_dbl(&p1, &p2);
+          ge_p1p1_to_p2(&p2, &p1);
+        }
+        ge_p2_dbl(&p1, &p2);
+        ge_p1p1_to_p3(&m_p3, &p1);
+      }
+      return *this;
+    }
+
     // returns a * this + G
     point_t mul_plus_G(const scalar_t& a) const
     {

src/crypto/msm.cpp

@@ -0,0 +1,25 @@
+// Copyright (c) 2023-2023 Zano Project
+// Copyright (c) 2023-2023 sowle (val@zano.org, crypto.sowle@gmail.com)
+// Distributed under the MIT/X11 software license, see the accompanying
+// file COPYING or http://www.opensource.org/licenses/mit-license.php.
+//
+#include "epee/include/misc_log_ex.h"
+//#include "zarcanum.h"
+#include "msm.h"
+//#include "../currency_core/crypto_config.h"    // TODO: move it to the crypto
+//#include "../common/crypto_stream_operators.h" // TODO: move it to the crypto
+
+#if 0
+#  define DBG_VAL_PRINT(x) std::cout << std::setw(30) << std::left << #x ": " << x << std::endl
+#  define DBG_PRINT(x)     std::cout << x << std::endl
+#else
+#  define DBG_VAL_PRINT(x) (void(0))
+#  define DBG_PRINT(x)     (void(0))
+#endif
+
+namespace crypto
+{
+
+
+
+} // namespace crypto

src/crypto/msm.h

@@ -0,0 +1,173 @@
+// Copyright (c) 2023-2023 Zano Project (https://zano.org/)
+// Copyright (c) 2023-2023 sowle (val@zano.org, crypto.sowle@gmail.com)
+// Distributed under the MIT/X11 software license, see the accompanying
+// file COPYING or http://www.opensource.org/licenses/mit-license.php.
+#pragma once
+
+// This file contains Multi-Scalar Multiplication routines
+
+#include "epee/include/misc_log_ex.h"
+#include "crypto-sugar.h"
+
+namespace crypto
+{
+
+  template<typename CT>
+  bool msm_and_check_zero_naive(const scalar_vec_t& g_scalars, const scalar_vec_t& h_scalars, const point_t& summand)
+  {
+    CHECK_AND_ASSERT_MES(g_scalars.size() <= CT::c_bpp_mn_max, false, "g_scalars oversized");
+    CHECK_AND_ASSERT_MES(h_scalars.size() <= CT::c_bpp_mn_max, false, "h_scalars oversized");
+
+    point_t result = summand;
+
+    for (size_t i = 0; i < g_scalars.size(); ++i)
+      result += g_scalars[i] * CT::get_generator(false, i);
+
+    for (size_t i = 0; i < h_scalars.size(); ++i)
+      result += h_scalars[i] * CT::get_generator(true, i);
+
+    if (!result.is_zero())
+    {
+      LOG_PRINT_L0("msm result is non zero: " << result);
+      return false;
+    }
+    return true;
+  }
+
+
+  // https://eprint.iacr.org/2022/999.pdf
+  // "Pippenger algorithm [1], and its variant that is widely used in the ZK space is called the bucket method"
+  template<typename CT>
+  bool msm_and_check_zero_pippenger_v3(const scalar_vec_t& g_scalars, const scalar_vec_t& h_scalars, const point_t& summand, uint8_t c)
+  {
+    // TODO: with c = 8 and with direct access got much worse result than with c = 7 and get_bits() for N = 128..256, consider checking again for bigger datasets (N>256) 
+    // TODO: consider preparing a cached generators' points
+
+    CHECK_AND_ASSERT_MES(g_scalars.size() <= CT::c_bpp_mn_max, false, "g_scalars oversized");
+    CHECK_AND_ASSERT_MES(h_scalars.size() <= CT::c_bpp_mn_max, false, "h_scalars oversized");
+    CHECK_AND_ASSERT_MES(c < 10, false, "c is too big");
+
+    size_t C = 1ull << c;
+
+    // k_max * c + (c-1) >= max_bit_idx
+    // 
+    //                 max_bit_idx - (c - 1)            max_bit_idx - (c - 1) + (c - 1)              max_bit_idx
+    // k_max = ceil ( --------------------- ) = floor ( ------------------------------ )  =  floor ( ----------- )
+    //                           c                                    c                                   c      
+    const size_t b = 253; // the maximum number of bits in x  https://eprint.iacr.org/2022/999.pdf  TODO: we may also scan for maximum bit used in all the scalars if all the scalars are small
+    const size_t max_bit_idx = b - 1;
+    const size_t k_max = max_bit_idx / c;
+    const size_t K = k_max + 1;
+
+    std::vector<point_t> buckets(C * K);
+    std::vector<bool> buckets_inited(C * K);
+    std::vector<point_t> Sk(K);
+    std::vector<bool> Sk_inited(K);
+    std::vector<point_t> Gk(K);
+    std::vector<bool> Gk_inited(K);
+
+    // first loop, calculate partial bucket sums
+    for (size_t n = 0; n < g_scalars.size(); ++n)
+    {
+      for (size_t k = 0; k < K; ++k)
+      {
+        uint64_t l = g_scalars[n].get_bits((uint8_t)(k * c), c); // l in [0; 2^c-1]
+        if (l != 0)
+        {
+          size_t bucket_id = l * K + k;
+          if (buckets_inited[bucket_id])
+            buckets[bucket_id] += CT::get_generator(false, n);
+          else
+          {
+            buckets[bucket_id] =  CT::get_generator(false, n);
+            buckets_inited[bucket_id] = true;
+          }
+        }
+      }
+    }
+    // still the first loop (continued)
+    for (size_t n = 0; n < h_scalars.size(); ++n)
+    {
+      for (size_t k = 0; k < K; ++k)
+      {
+        uint64_t l = h_scalars[n].get_bits((uint8_t)(k * c), c); // l in [0; 2^c-1]
+        if (l != 0)
+        {
+          size_t bucket_id = l * K + k;
+          if (buckets_inited[bucket_id])
+            buckets[bucket_id] += CT::get_generator(true, n);
+          else
+          {
+            buckets[bucket_id] =  CT::get_generator(true, n);
+            buckets_inited[bucket_id] = true;
+          }
+        }
+      }
+    }
+
+    // the second loop
+    for (size_t l = C - 1; l > 0; --l)
+    {
+      for (size_t k = 0; k < K; ++k)
+      {
+        size_t bucket_id = l * K + k;
+        if (buckets_inited[bucket_id])
+        {
+          if (Sk_inited[k])
+            Sk[k] += buckets[bucket_id];
+          else
+          {
+            Sk[k] = buckets[bucket_id];
+            Sk_inited[k] = true;
+          }
+        }
+
+        if (Sk_inited[k])
+        {
+          if (Gk_inited[k])
+            Gk[k] += Sk[k];
+          else
+          {
+            Gk[k] = Sk[k];
+            Gk_inited[k] = true;
+          }
+        }
+      }
+    }
+
+    // the third loop: Horner’s rule
+    point_t result = Gk_inited[K - 1] ? Gk[K - 1] : c_point_0;
+    for (size_t k = K - 2; k != SIZE_MAX; --k)
+    {
+      result.modify_mul_pow_2(c);
+      if (Gk_inited[k])
+        result += Gk[k];
+    }
+
+    result += summand;
+
+    if (!result.is_zero())
+    {
+      LOG_PRINT_L0("multiexp result is non zero: " << result);
+      return false;
+    }
+
+    return true;
+  }
+
+
+
+  
+
+  // Just switcher
+
+  template<typename CT>
+  bool msm_and_check_zero(const scalar_vec_t& g_scalars, const scalar_vec_t& h_scalars, const point_t& summand)
+  {
+    //return msm_and_check_zero_naive<CT>(g_scalars, h_scalars, summand);
+    return msm_and_check_zero_pippenger_v3<CT>(g_scalars, h_scalars, summand, 7);
+  }
+
+
+
+} // namespace crypto

src/crypto/range_proof_bpp.h

@@ -716,7 +716,7 @@ namespace crypto
 
     point_t GH_exponents = c_point_0;
     CT::calc_pedersen_commitment(G_scalar, H_scalar, GH_exponents);
-    bool result = multiexp_and_check_being_zero<CT>(g_scalars, h_scalars, summand + GH_exponents);
+    bool result = msm_and_check_zero<CT>(g_scalars, h_scalars, summand + GH_exponents);
     if (result)
       DBG_PRINT(ENDL << " . . . . bpp_verify() -- SUCCEEDED!!!" << ENDL);
     return result;

src/crypto/range_proof_bppe.h

@@ -734,7 +734,7 @@ namespace crypto
 
     point_t GH_exponents = c_point_0;
     CT::calc_pedersen_commitment_2(G_scalar, H_scalar, H2_scalar, GH_exponents);
-    bool result = multiexp_and_check_being_zero<CT>(g_scalars, h_scalars, summand + GH_exponents);
+    bool result = msm_and_check_zero<CT>(g_scalars, h_scalars, summand + GH_exponents);
     if (result)
       DBG_PRINT(ENDL << " . . . . bppe_verify() -- SUCCEEDED!!!" << ENDL);
     return result;

src/crypto/range_proofs.h

@@ -41,7 +41,7 @@ namespace crypto
   ////////////////////////////////////////
   struct bpp_ct_generators_HGX
   {
-    // NOTE! This notation follows the original BP+ whitepaper, see mapping to Zano's generators below
+    // NOTE! This notation follows the original BP+ whitepaper, see mapping to Zano's generators in range_proofs.cpp
     static const point_t& bpp_G;
     static const point_t& bpp_H;
     static const point_t& bpp_H2;
@@ -49,7 +49,7 @@ namespace crypto
 
   struct bpp_ct_generators_UGX
   {
-    // NOTE! This notation follows the original BP+ whitepaper, see mapping to Zano's generators below
+    // NOTE! This notation follows the original BP+ whitepaper, see mapping to Zano's generators in range_proofs.cpp
     static const point_t& bpp_G;
     static const point_t& bpp_H;
     static const point_t& bpp_H2;
@@ -138,32 +138,11 @@ namespace crypto
 
   typedef bpp_crypto_trait_zano<bpp_ct_generators_HGX, 128, 16> bpp_crypto_trait_Zarcanum;
 
-  
-  // efficient multiexponentiation (naive stub implementation atm, TODO)
-  template<typename CT>
-  bool multiexp_and_check_being_zero(const scalar_vec_t& g_scalars, const scalar_vec_t& h_scalars, const point_t& summand)
-  {
-    CHECK_AND_ASSERT_MES(g_scalars.size() <= CT::c_bpp_mn_max, false, "g_scalars oversized");
-    CHECK_AND_ASSERT_MES(h_scalars.size() <= CT::c_bpp_mn_max, false, "h_scalars oversized");
-
-    point_t result = summand;
-
-    for (size_t i = 0; i < g_scalars.size(); ++i)
-      result += g_scalars[i] * CT::get_generator(false, i);
-
-    for (size_t i = 0; i < h_scalars.size(); ++i)
-      result += h_scalars[i] * CT::get_generator(true, i);
-
-    if (!result.is_zero())
-    {
-      LOG_PRINT_L0("multiexp result is non zero: " << result);
-      return false;
-    }
-    return true;
-  }
-
 
 } // namespace crypto
 
+#include "epee/include/profile_tools.h" // <- remove this, sowle
+
+#include "msm.h"
 #include "range_proof_bpp.h"
 #include "range_proof_bppe.h"

tests/functional_tests/crypto_tests.cpp

@@ -1619,6 +1619,7 @@ TEST(crypto, schnorr_sig)
   return true;
 }
 
+
 TEST(crypto, point_negation)
 {
   ASSERT_EQ(c_point_0, -c_point_0);
@@ -1694,6 +1695,73 @@ TEST(crypto, scalar_get_bits)
   return true;
 }
 
+template<typename CT>
+bool crypto_msm_runner(size_t N, size_t low_bits_to_clear, size_t high_bits_to_clear)
+{
+  scalar_vec_t g_scalars, h_scalars;
+  g_scalars.resize_and_make_random(N);
+  h_scalars.resize_and_make_random(N);
+  if (N > 4)
+  {
+    g_scalars[0] = c_scalar_Lm1; // always include the max and the min
+    h_scalars[0] = c_scalar_Lm1;
+    g_scalars[1] = 0;
+    h_scalars[1] = 0;
+  }
+
+  point_t sum = c_point_0;
+  for(size_t i = 0; i < N; ++i)
+  {
+    for(size_t bit_index = 0; bit_index < low_bits_to_clear; ++bit_index)
+    {
+      g_scalars[i].clear_bit(bit_index);
+      h_scalars[i].clear_bit(bit_index);
+    }
+    for(size_t bit_index = 256 - high_bits_to_clear; bit_index < 256; ++bit_index)
+    {
+      g_scalars[i].clear_bit(bit_index);
+      h_scalars[i].clear_bit(bit_index);
+    }
+    sum += g_scalars[i] * CT::get_generator(false, i) + h_scalars[i] * CT::get_generator(true, i);
+  }
+
+  //TIME_MEASURE_START(t);
+  bool r = msm_and_check_zero<CT>(g_scalars, h_scalars, -sum);
+  //TIME_MEASURE_FINISH(t);
+  return r;
+}
+
+TEST(crypto, msm)
+{
+  // test the default msm_and_check_zero correctness
+  bool r = false;
+
+  for(size_t N = 1; N <= 128; ++N)
+  {
+    std::cout << "N = " << N << ENDL;
+    r = crypto_msm_runner<bpp_crypto_trait_Zarcanum>(N, 0, 0);
+    ASSERT_TRUE(r);
+    r = crypto_msm_runner<bpp_crypto_trait_ZC_out>(N, 0, 0);
+    ASSERT_TRUE(r);
+  }
+
+  for(size_t i = 0; i <= 128; ++i)
+  {
+    std::cout << "i = " << i << ENDL;
+    r = crypto_msm_runner<bpp_crypto_trait_Zarcanum>(128, i, 0);
+    ASSERT_TRUE(r);
+    r = crypto_msm_runner<bpp_crypto_trait_Zarcanum>(128, 0, i);
+    ASSERT_TRUE(r);
+    r = crypto_msm_runner<bpp_crypto_trait_ZC_out>(256, i, 0);
+    ASSERT_TRUE(r);
+    r = crypto_msm_runner<bpp_crypto_trait_ZC_out>(256, 0, i);
+    ASSERT_TRUE(r);
+  }
+
+  return true;
+}
+
+
 //
 // test's runner
 //