ringct: add bos coster multiexp

4 files changed, 455 insertions, 51 deletions

diff --git a/src/ringct/CMakeLists.txt b/src/ringct/CMakeLists.txt
index c8dcdca26..29f166a3b 100644
--- a/src/ringct/CMakeLists.txt
+++ b/src/ringct/CMakeLists.txt
@@ -30,11 +30,13 @@ set(ringct_basic_sources
   rctOps.cpp
   rctTypes.cpp
   rctCryptoOps.c
+  multiexp.cc
   bulletproofs.cc)
 
 set(ringct_basic_private_headers
   rctOps.h
   rctTypes.h
+  multiexp.h
   bulletproofs.h)
 
 monero_private_headers(ringct_basic
diff --git a/src/ringct/bulletproofs.cc b/src/ringct/bulletproofs.cc
index aa56b0589..057f19029 100644
--- a/src/ringct/bulletproofs.cc
+++ b/src/ringct/bulletproofs.cc
@@ -38,6 +38,7 @@ extern "C"
 #include "crypto/crypto-ops.h"
 }
 #include "rctOps.h"
+#include "multiexp.h"
 #include "bulletproofs.h"
 
 #undef MONERO_DEFAULT_LOG_CATEGORY
@@ -58,6 +59,7 @@ static rct::key inner_product(const rct::keyV &a, const rct::keyV &b);
 static constexpr size_t maxN = 64;
 static constexpr size_t maxM = 16;
 static rct::key Hi[maxN*maxM], Gi[maxN*maxM];
+static ge_p3 Hi_p3[maxN*maxM], Gi_p3[maxN*maxM];
 static ge_dsmp Gprecomp[maxN*maxM], Hprecomp[maxN*maxM];
 static const rct::key TWO = { {0x02, 0x00, 0x00,0x00 , 0x00, 0x00, 0x00,0x00 , 0x00, 0x00, 0x00,0x00 , 0x00, 0x00, 0x00,0x00 , 0x00, 0x00, 0x00,0x00 , 0x00, 0x00, 0x00,0x00 , 0x00, 0x00, 0x00,0x00 , 0x00, 0x00, 0x00,0x00  } };
 static const rct::keyV oneN = vector_dup(rct::identity(), maxN);
@@ -110,9 +112,12 @@ static void init_exponents()
   {
     Hi[i] = get_exponent(rct::H, i * 2);
     rct::precomp(Hprecomp[i], Hi[i]);
+    CHECK_AND_ASSERT_THROW_MES(ge_frombytes_vartime(&Hi_p3[i], Hi[i].bytes) == 0, "ge_frombytes_vartime failed");
     Gi[i] = get_exponent(rct::H, i * 2 + 1);
     rct::precomp(Gprecomp[i], Gi[i]);
+    CHECK_AND_ASSERT_THROW_MES(ge_frombytes_vartime(&Gi_p3[i], Gi[i].bytes) == 0, "ge_frombytes_vartime failed");
   }
+  MINFO("cache size: " << (sizeof(Hi)+sizeof(Hprecomp)+sizeof(Hi_p3))*2/1024 << " kB");
   init_done = true;
 }
 
@@ -121,6 +126,26 @@ static rct::key vector_exponent(const rct::keyV &a, const rct::keyV &b)
 {
   CHECK_AND_ASSERT_THROW_MES(a.size() == b.size(), "Incompatible sizes of a and b");
   CHECK_AND_ASSERT_THROW_MES(a.size() <= maxN*maxM, "Incompatible sizes of a and maxN");
+#if 1
+  std::vector<MultiexpData> multiexp_data;
+  multiexp_data.reserve(a.size()*2);
+  for (size_t i = 0; i < a.size(); ++i)
+  {
+    if (!(a[i] == rct::zero()))
+    {
+      multiexp_data.resize(multiexp_data.size() + 1);
+      multiexp_data.back().scalar = a[i];
+      multiexp_data.back().point = Gi_p3[i];
+    }
+    if (!(b[i] == rct::zero()))
+    {
+      multiexp_data.resize(multiexp_data.size() + 1);
+      multiexp_data.back().scalar = b[i];
+      multiexp_data.back().point = Hi_p3[i];
+    }
+  }
+  return bos_coster_heap_conv_robust(multiexp_data);
+#else
   ge_p3 res_p3 = ge_p3_identity;
   for (size_t i = 0; i < a.size(); ++i)
   {
@@ -129,6 +154,7 @@ static rct::key vector_exponent(const rct::keyV &a, const rct::keyV &b)
   rct::key res;
   ge_p3_tobytes(res.bytes, &res_p3);
   return res;
+#endif
 }
 
 /* Compute a custom vector-scalar commitment */
@@ -138,6 +164,26 @@ static rct::key vector_exponent_custom(const rct::keyV &A, const rct::keyV &B, c
   CHECK_AND_ASSERT_THROW_MES(a.size() == b.size(), "Incompatible sizes of a and b");
   CHECK_AND_ASSERT_THROW_MES(a.size() == A.size(), "Incompatible sizes of a and A");
   CHECK_AND_ASSERT_THROW_MES(a.size() <= maxN*maxM, "Incompatible sizes of a and maxN");
+#if 1
+  std::vector<MultiexpData> multiexp_data;
+  multiexp_data.reserve(a.size()*2);
+  for (size_t i = 0; i < a.size(); ++i)
+  {
+    if (!(a[i] == rct::zero()))
+    {
+      multiexp_data.resize(multiexp_data.size() + 1);
+      multiexp_data.back().scalar = a[i];
+      CHECK_AND_ASSERT_THROW_MES(ge_frombytes_vartime(&multiexp_data.back().point, A[i].bytes) == 0, "ge_frombytes_vartime failed");
+    }
+    if (!(b[i] == rct::zero()))
+    {
+      multiexp_data.resize(multiexp_data.size() + 1);
+      multiexp_data.back().scalar = b[i];
+      CHECK_AND_ASSERT_THROW_MES(ge_frombytes_vartime(&multiexp_data.back().point, B[i].bytes) == 0, "ge_frombytes_vartime failed");
+    }
+  }
+  return bos_coster_heap_conv_robust(multiexp_data);
+#else
   ge_p3 res_p3 = ge_p3_identity;
   for (size_t i = 0; i < a.size(); ++i)
   {
@@ -174,6 +220,7 @@ static rct::key vector_exponent_custom(const rct::keyV &A, const rct::keyV &B, c
   rct::key res;
   ge_p3_tobytes(res.bytes, &res_p3);
   return res;
+#endif
 }
 
 /* Given a scalar, construct a vector of powers */
@@ -924,7 +971,7 @@ bool bulletproof_VERIFY(const Bulletproof &proof)
 
   PERF_TIMER_START_BP(VERIFY_line_61);
   // PAPER LINE 61
-  rct::key L61Left;
+  rct::key L61Left, L61Right;
   rct::addKeys2(L61Left, proof.taux, proof.t, rct::H);
 
   const rct::keyV zpow = vector_powers(z, M+3);
@@ -939,36 +986,61 @@ bool bulletproof_VERIFY(const Bulletproof &proof)
   }
   PERF_TIMER_STOP(VERIFY_line_61);
 
-  PERF_TIMER_START_BP(VERIFY_line_61rl);
-  sc_muladd(tmp.bytes, z.bytes, ip1y.bytes, k.bytes);
-  rct::key L61Right = rct::scalarmultKey(rct::H, tmp);
-  ge_p3 L61Right_p3;
-  CHECK_AND_ASSERT_THROW_MES(ge_frombytes_vartime(&L61Right_p3, L61Right.bytes) == 0, "ge_frombytes_vartime failed");
-  for (size_t j = 0; j+1 < proof.V.size(); j += 2)
+  // multiexp is slower for small numbers of calcs
+  if (M >= 16)
   {
-    CHECK_AND_ASSERT_MES(j+2+1 < zpow.size(), false, "invalid zpow index");
-    ge_dsmp precomp0, precomp1;
-    rct::precomp(precomp0, j < proof.V.size() ? proof.V[j] : rct::identity());
-    rct::precomp(precomp1, j+1 < proof.V.size() ? proof.V[j+1] : rct::identity());
-    rct::addKeys3acc_p3(&L61Right_p3, zpow[j+2], precomp0, zpow[j+2+1], precomp1);
+    PERF_TIMER_START_BP(VERIFY_line_61rl_new);
+    sc_muladd(tmp.bytes, z.bytes, ip1y.bytes, k.bytes);
+    std::vector<MultiexpData> multiexp_data;
+    multiexp_data.reserve(3+M);
+    multiexp_data.push_back({tmp, rct::H});
+    for (size_t j = 0; j < M; j++)
+    {
+      if (!(zpow[j+2] == rct::zero()))
+        multiexp_data.push_back({zpow[j+2], j < proof.V.size() ? proof.V[j] : rct::identity()});
+    }
+    if (!(x == rct::zero()))
+      multiexp_data.push_back({x, proof.T1});
+    rct::key xsq;
+    sc_mul(xsq.bytes, x.bytes, x.bytes);
+    if (!(xsq == rct::zero()))
+      multiexp_data.push_back({xsq, proof.T2});
+    L61Right = bos_coster_heap_conv_robust(multiexp_data);
+    PERF_TIMER_STOP(VERIFY_line_61rl_new);
   }
-  for (size_t j = proof.V.size() & 0xfffffffe; j < M; j++)
+  else
   {
-    CHECK_AND_ASSERT_MES(j+2 < zpow.size(), false, "invalid zpow index");
-    // faster equivalent to:
-    // tmp = rct::scalarmultKey(j < proof.V.size() ? proof.V[j] : rct::identity(), zpow[j+2]);
-    // rct::addKeys(L61Right, L61Right, tmp);
-    if (j < proof.V.size())
-      addKeys_acc_p3(&L61Right_p3, zpow[j+2], proof.V[j]);
-  }
+    PERF_TIMER_START_BP(VERIFY_line_61rl_old);
+    sc_muladd(tmp.bytes, z.bytes, ip1y.bytes, k.bytes);
+    L61Right = rct::scalarmultKey(rct::H, tmp);
+    ge_p3 L61Right_p3;
+    CHECK_AND_ASSERT_THROW_MES(ge_frombytes_vartime(&L61Right_p3, L61Right.bytes) == 0, "ge_frombytes_vartime failed");
+    for (size_t j = 0; j+1 < proof.V.size(); j += 2)
+    {
+      CHECK_AND_ASSERT_MES(j+2+1 < zpow.size(), false, "invalid zpow index");
+      ge_dsmp precomp0, precomp1;
+      rct::precomp(precomp0, j < proof.V.size() ? proof.V[j] : rct::identity());
+      rct::precomp(precomp1, j+1 < proof.V.size() ? proof.V[j+1] : rct::identity());
+      rct::addKeys3acc_p3(&L61Right_p3, zpow[j+2], precomp0, zpow[j+2+1], precomp1);
+    }
+    for (size_t j = proof.V.size() & 0xfffffffe; j < M; j++)
+    {
+      CHECK_AND_ASSERT_MES(j+2 < zpow.size(), false, "invalid zpow index");
+      // faster equivalent to:
+      // tmp = rct::scalarmultKey(j < proof.V.size() ? proof.V[j] : rct::identity(), zpow[j+2]);
+      // rct::addKeys(L61Right, L61Right, tmp);
+      if (j < proof.V.size())
+        addKeys_acc_p3(&L61Right_p3, zpow[j+2], proof.V[j]);
+    }
 
-  addKeys_acc_p3(&L61Right_p3, x, proof.T1);
+    addKeys_acc_p3(&L61Right_p3, x, proof.T1);
 
-  rct::key xsq;
-  sc_mul(xsq.bytes, x.bytes, x.bytes);
-  addKeys_acc_p3(&L61Right_p3, xsq, proof.T2);
-  ge_p3_tobytes(L61Right.bytes, &L61Right_p3);
-  PERF_TIMER_STOP(VERIFY_line_61rl);
+    rct::key xsq;
+    sc_mul(xsq.bytes, x.bytes, x.bytes);
+    addKeys_acc_p3(&L61Right_p3, xsq, proof.T2);
+    ge_p3_tobytes(L61Right.bytes, &L61Right_p3);
+    PERF_TIMER_STOP(VERIFY_line_61rl_old);
+  }
 
   if (!(L61Right == L61Left))
   {
@@ -998,7 +1070,6 @@ bool bulletproof_VERIFY(const Bulletproof &proof)
   PERF_TIMER_START_BP(VERIFY_line_24_25);
   // Basically PAPER LINES 24-25
   // Compute the curvepoints from G[i] and H[i]
-  ge_p3 inner_prod_p3 = ge_p3_identity;
   rct::key yinvpow = rct::identity();
   rct::key ypow = rct::identity();
 
@@ -1009,6 +1080,9 @@ bool bulletproof_VERIFY(const Bulletproof &proof)
     winv[i] = invert(w[i]);
   PERF_TIMER_STOP(VERIFY_line_24_25_invert);
 
+  std::vector<MultiexpData> multiexp_data;
+  multiexp_data.clear();
+  multiexp_data.reserve(MN*2);
   for (size_t i = 0; i < MN; ++i)
   {
     // Convert the index to binary IN REVERSE and construct the scalar exponent
@@ -1040,9 +1114,10 @@ bool bulletproof_VERIFY(const Bulletproof &proof)
     sc_muladd(tmp.bytes, z.bytes, ypow.bytes, tmp.bytes);
     sc_mulsub(h_scalar.bytes, tmp.bytes, yinvpow.bytes, h_scalar.bytes);
 
-    // Now compute the basepoint's scalar multiplication
-    // Each of these could be written as a multiexp operation instead
-    addKeys3acc_p3(&inner_prod_p3, g_scalar, Gprecomp[i], h_scalar, Hprecomp[i]);
+    if (!(g_scalar == rct::zero()))
+      multiexp_data.push_back({g_scalar, Gi_p3[i]});
+    if (!(h_scalar == rct::zero()))
+      multiexp_data.push_back({h_scalar, Hi_p3[i]});
 
     if (i != MN-1)
     {
@@ -1050,36 +1125,64 @@ bool bulletproof_VERIFY(const Bulletproof &proof)
       sc_mul(ypow.bytes, ypow.bytes, y.bytes);
     }
   }
-  rct::key inner_prod;
-  ge_p3_tobytes(inner_prod.bytes, &inner_prod_p3);
+
+  rct::key inner_prod = bos_coster_heap_conv_robust(multiexp_data);
   PERF_TIMER_STOP(VERIFY_line_24_25);
 
-  PERF_TIMER_START_BP(VERIFY_line_26);
-  // PAPER LINE 26
   rct::key pprime;
-  sc_sub(tmp.bytes, rct::zero().bytes, proof.mu.bytes);
-  rct::addKeys(pprime, P, rct::scalarmultBase(tmp));
-  ge_p3 pprime_p3;
-  CHECK_AND_ASSERT_MES(ge_frombytes_vartime(&pprime_p3, pprime.bytes) == 0, false, "ge_frombytes_vartime failed");
+  // multiexp does not seem to give any speedup here
+  if(0)
+  {
+    PERF_TIMER_START_BP(VERIFY_line_26_new);
+    // PAPER LINE 26
+    std::vector<MultiexpData> multiexp_data;
+    multiexp_data.reserve(1+2*rounds);
+
+    sc_sub(tmp.bytes, rct::zero().bytes, proof.mu.bytes);
+    rct::addKeys(pprime, P, rct::scalarmultBase(tmp));
+    for (size_t i = 0; i < rounds; ++i)
+    {
+      sc_mul(tmp.bytes, w[i].bytes, w[i].bytes);
+      sc_mul(tmp2.bytes, winv[i].bytes, winv[i].bytes);
+      if (!(tmp == rct::zero()))
+        multiexp_data.push_back({tmp, proof.L[i]});
+      if (!(tmp2 == rct::zero()))
+        multiexp_data.push_back({tmp2, proof.R[i]});
+    }
+    sc_mul(tmp.bytes, proof.t.bytes, x_ip.bytes);
+    if (!(tmp == rct::zero()))
+      multiexp_data.push_back({tmp, rct::H});
+    addKeys(pprime, pprime, bos_coster_heap_conv_robust(multiexp_data));
+    PERF_TIMER_STOP(VERIFY_line_26_new);
+  }
 
-  for (size_t i = 0; i < rounds; ++i)
   {
-    sc_mul(tmp.bytes, w[i].bytes, w[i].bytes);
-    sc_mul(tmp2.bytes, winv[i].bytes, winv[i].bytes);
+    PERF_TIMER_START_BP(VERIFY_line_26_old);
+    // PAPER LINE 26
+    sc_sub(tmp.bytes, rct::zero().bytes, proof.mu.bytes);
+    rct::addKeys(pprime, P, rct::scalarmultBase(tmp));
+    ge_p3 pprime_p3;
+    CHECK_AND_ASSERT_MES(ge_frombytes_vartime(&pprime_p3, pprime.bytes) == 0, false, "ge_frombytes_vartime failed");
+
+    for (size_t i = 0; i < rounds; ++i)
+    {
+      sc_mul(tmp.bytes, w[i].bytes, w[i].bytes);
+      sc_mul(tmp2.bytes, winv[i].bytes, winv[i].bytes);
 #if 1
-    ge_dsmp cacheL, cacheR;
-    rct::precomp(cacheL, proof.L[i]);
-    rct::precomp(cacheR, proof.R[i]);
-    addKeys3acc_p3(&pprime_p3, tmp, cacheL, tmp2, cacheR);
+      ge_dsmp cacheL, cacheR;
+      rct::precomp(cacheL, proof.L[i]);
+      rct::precomp(cacheR, proof.R[i]);
+      addKeys3acc_p3(&pprime_p3, tmp, cacheL, tmp2, cacheR);
 #else
-    rct::addKeys(pprime, pprime, rct::scalarmultKey(proof.L[i], tmp));
-    rct::addKeys(pprime, pprime, rct::scalarmultKey(proof.R[i], tmp2));
+      rct::addKeys(pprime, pprime, rct::scalarmultKey(proof.L[i], tmp));
+      rct::addKeys(pprime, pprime, rct::scalarmultKey(proof.R[i], tmp2));
 #endif
+    }
+    sc_mul(tmp.bytes, proof.t.bytes, x_ip.bytes);
+    addKeys_acc_p3(&pprime_p3, tmp, rct::H);
+    ge_p3_tobytes(pprime.bytes, &pprime_p3);
+    PERF_TIMER_STOP(VERIFY_line_26_old);
   }
-  sc_mul(tmp.bytes, proof.t.bytes, x_ip.bytes);
-  addKeys_acc_p3(&pprime_p3, tmp, rct::H);
-  ge_p3_tobytes(pprime.bytes, &pprime_p3);
-  PERF_TIMER_STOP(VERIFY_line_26);
 
   PERF_TIMER_START_BP(VERIFY_step2_check);
   sc_mul(tmp.bytes, proof.a.bytes, proof.b.bytes);
diff --git a/src/ringct/multiexp.cc b/src/ringct/multiexp.cc
new file mode 100644
index 000000000..2545325ae
--- /dev/null
+++ b/src/ringct/multiexp.cc
@@ -0,0 +1,239 @@
+// Copyright (c) 2017, The Monero Project
+// 
+// All rights reserved.
+// 
+// Redistribution and use in source and binary forms, with or without modification, are
+// permitted provided that the following conditions are met:
+// 
+// 1. Redistributions of source code must retain the above copyright notice, this list of
+//    conditions and the following disclaimer.
+// 
+// 2. Redistributions in binary form must reproduce the above copyright notice, this list
+//    of conditions and the following disclaimer in the documentation and/or other
+//    materials provided with the distribution.
+// 
+// 3. Neither the name of the copyright holder nor the names of its contributors may be
+//    used to endorse or promote products derived from this software without specific
+//    prior written permission.
+// 
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
+// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
+// THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
+// THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Adapted from Python code by Sarang Noether
+
+#include "misc_log_ex.h"
+#include "common/perf_timer.h"
+extern "C"
+{
+#include "crypto/crypto-ops.h"
+}
+#include "rctOps.h"
+#include "multiexp.h"
+
+#undef MONERO_DEFAULT_LOG_CATEGORY
+#define MONERO_DEFAULT_LOG_CATEGORY "multiexp.boscoster"
+
+//#define MULTIEXP_PERF(x) x
+#define MULTIEXP_PERF(x)
+
+namespace rct
+{
+
+static inline bool operator<(const rct::key &k0, const rct::key&k1)
+{
+  for (int n = 31; n >= 0; --n)
+  {
+    if (k0.bytes[n] < k1.bytes[n])
+      return true;
+    if (k0.bytes[n] > k1.bytes[n])
+      return false;
+  }
+  return false;
+}
+
+static inline rct::key div2(const rct::key &k)
+{
+  rct::key res;
+  int carry = 0;
+  for (int n = 31; n >= 0; --n)
+  {
+    int new_carry = (k.bytes[n] & 1) << 7;
+    res.bytes[n] = k.bytes[n] / 2 + carry;
+    carry = new_carry;
+  }
+  return res;
+}
+
+rct::key bos_coster_heap_conv(std::vector<MultiexpData> &data)
+{
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(bos_coster, 1000000));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(setup, 1000000));
+  size_t points = data.size();
+  CHECK_AND_ASSERT_THROW_MES(points > 1, "Not enough points");
+  std::vector<size_t> heap(points);
+  for (size_t n = 0; n < points; ++n)
+    heap[n] = n;
+
+  auto Comp = [&](size_t e0, size_t e1) { return data[e0].scalar < data[e1].scalar; };
+  std::make_heap(heap.begin(), heap.end(), Comp);
+  MULTIEXP_PERF(PERF_TIMER_STOP(setup));
+
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(loop, 1000000));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(pop, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(pop));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(add, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(add));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(sub, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(sub));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(push, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(push));
+  while (heap.size() > 1)
+  {
+    MULTIEXP_PERF(PERF_TIMER_RESUME(pop));
+    std::pop_heap(heap.begin(), heap.end(), Comp);
+    size_t index1 = heap.back();
+    heap.pop_back();
+    std::pop_heap(heap.begin(), heap.end(), Comp);
+    size_t index2 = heap.back();
+    heap.pop_back();
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(pop));
+
+    MULTIEXP_PERF(PERF_TIMER_RESUME(add));
+    ge_cached cached;
+    ge_p3_to_cached(&cached, &data[index1].point);
+    ge_p1p1 p1;
+    ge_add(&p1, &data[index2].point, &cached);
+    ge_p1p1_to_p3(&data[index2].point, &p1);
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(add));
+
+    MULTIEXP_PERF(PERF_TIMER_RESUME(sub));
+    sc_sub(data[index1].scalar.bytes, data[index1].scalar.bytes, data[index2].scalar.bytes);
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(sub));
+
+    MULTIEXP_PERF(PERF_TIMER_RESUME(push));
+    if (!(data[index1].scalar == rct::zero()))
+    {
+      heap.push_back(index1);
+      std::push_heap(heap.begin(), heap.end(), Comp);
+    }
+
+    heap.push_back(index2);
+    std::push_heap(heap.begin(), heap.end(), Comp);
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(push));
+  }
+  MULTIEXP_PERF(PERF_TIMER_STOP(push));
+  MULTIEXP_PERF(PERF_TIMER_STOP(sub));
+  MULTIEXP_PERF(PERF_TIMER_STOP(add));
+  MULTIEXP_PERF(PERF_TIMER_STOP(pop));
+  MULTIEXP_PERF(PERF_TIMER_STOP(loop));
+
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(end, 1000000));
+  //return rct::scalarmultKey(data[index1].point, data[index1].scalar);
+  std::pop_heap(heap.begin(), heap.end(), Comp);
+  size_t index1 = heap.back();
+  heap.pop_back();
+  ge_p2 p2;
+  ge_scalarmult(&p2, data[index1].scalar.bytes, &data[index1].point);
+  rct::key res;
+  ge_tobytes(res.bytes, &p2);
+  return res;
+}
+
+rct::key bos_coster_heap_conv_robust(std::vector<MultiexpData> &data)
+{
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(bos_coster, 1000000));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(setup, 1000000));
+  size_t points = data.size();
+  CHECK_AND_ASSERT_THROW_MES(points > 1, "Not enough points");
+  std::vector<size_t> heap(points);
+  for (size_t n = 0; n < points; ++n)
+    heap[n] = n;
+
+  auto Comp = [&](size_t e0, size_t e1) { return data[e0].scalar < data[e1].scalar; };
+  std::make_heap(heap.begin(), heap.end(), Comp);
+  MULTIEXP_PERF(PERF_TIMER_STOP(setup));
+
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(loop, 1000000));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(pop, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(pop));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(div, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(div));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(add, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(add));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(sub, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(sub));
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(push, 1000000)); MULTIEXP_PERF(PERF_TIMER_PAUSE(push));
+  while (heap.size() > 1)
+  {
+    MULTIEXP_PERF(PERF_TIMER_RESUME(pop));
+    std::pop_heap(heap.begin(), heap.end(), Comp);
+    size_t index1 = heap.back();
+    heap.pop_back();
+    std::pop_heap(heap.begin(), heap.end(), Comp);
+    size_t index2 = heap.back();
+    heap.pop_back();
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(pop));
+
+    ge_cached cached;
+    ge_p1p1 p1;
+
+    MULTIEXP_PERF(PERF_TIMER_RESUME(div));
+    while (1)
+    {
+      rct::key s1_2 = div2(data[index1].scalar);
+      if (!(data[index2].scalar < s1_2))
+       break;
+      if (data[index1].scalar.bytes[0] & 1)
+      {
+        data.resize(data.size()+1);
+        data.back().scalar = rct::identity();
+        data.back().point = data[index1].point;
+        heap.push_back(data.size() - 1);
+        std::push_heap(heap.begin(), heap.end(), Comp);
+      }
+      data[index1].scalar = div2(data[index1].scalar);
+      ge_p3_to_cached(&cached, &data[index1].point);
+      ge_add(&p1, &data[index1].point, &cached);
+      ge_p1p1_to_p3(&data[index1].point, &p1);
+    }
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(div));
+
+    MULTIEXP_PERF(PERF_TIMER_RESUME(add));
+    ge_p3_to_cached(&cached, &data[index1].point);
+    ge_add(&p1, &data[index2].point, &cached);
+    ge_p1p1_to_p3(&data[index2].point, &p1);
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(add));
+
+    MULTIEXP_PERF(PERF_TIMER_RESUME(sub));
+    sc_sub(data[index1].scalar.bytes, data[index1].scalar.bytes, data[index2].scalar.bytes);
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(sub));
+
+    MULTIEXP_PERF(PERF_TIMER_RESUME(push));
+    if (!(data[index1].scalar == rct::zero()))
+    {
+      heap.push_back(index1);
+      std::push_heap(heap.begin(), heap.end(), Comp);
+    }
+
+    heap.push_back(index2);
+    std::push_heap(heap.begin(), heap.end(), Comp);
+    MULTIEXP_PERF(PERF_TIMER_PAUSE(push));
+  }
+  MULTIEXP_PERF(PERF_TIMER_STOP(push));
+  MULTIEXP_PERF(PERF_TIMER_STOP(sub));
+  MULTIEXP_PERF(PERF_TIMER_STOP(add));
+  MULTIEXP_PERF(PERF_TIMER_STOP(pop));
+  MULTIEXP_PERF(PERF_TIMER_STOP(loop));
+
+  MULTIEXP_PERF(PERF_TIMER_START_UNIT(end, 1000000));
+  //return rct::scalarmultKey(data[index1].point, data[index1].scalar);
+  std::pop_heap(heap.begin(), heap.end(), Comp);
+  size_t index1 = heap.back();
+  heap.pop_back();
+  ge_p2 p2;
+  ge_scalarmult(&p2, data[index1].scalar.bytes, &data[index1].point);
+  rct::key res;
+  ge_tobytes(res.bytes, &p2);
+  return res;
+}
+
+}
diff --git a/src/ringct/multiexp.h b/src/ringct/multiexp.h
new file mode 100644
index 000000000..108db7c39
--- /dev/null
+++ b/src/ringct/multiexp.h
@@ -0,0 +1,60 @@
+// Copyright (c) 2017, The Monero Project
+// 
+// All rights reserved.
+// 
+// Redistribution and use in source and binary forms, with or without modification, are
+// permitted provided that the following conditions are met:
+// 
+// 1. Redistributions of source code must retain the above copyright notice, this list of
+//    conditions and the following disclaimer.
+// 
+// 2. Redistributions in binary form must reproduce the above copyright notice, this list
+//    of conditions and the following disclaimer in the documentation and/or other
+//    materials provided with the distribution.
+// 
+// 3. Neither the name of the copyright holder nor the names of its contributors may be
+//    used to endorse or promote products derived from this software without specific
+//    prior written permission.
+// 
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
+// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
+// THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
+// THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Adapted from Python code by Sarang Noether
+
+#pragma once
+
+#ifndef MULTIEXP_H
+#define MULTIEXP_H
+
+#include <vector>
+#include "crypto/crypto.h"
+#include "rctTypes.h"
+
+namespace rct
+{
+
+struct MultiexpData {
+  rct::key scalar;
+  ge_p3 point;
+
+  MultiexpData() {}
+  MultiexpData(const rct::key &s, const ge_p3 &p): scalar(s), point(p) {}
+  MultiexpData(const rct::key &s, const rct::key &p): scalar(s)
+  {
+    CHECK_AND_ASSERT_THROW_MES(ge_frombytes_vartime(&point, p.bytes) == 0, "ge_frombytes_vartime failed");
+  }
+};
+
+rct::key bos_coster_heap_conv(std::vector<MultiexpData> &data);
+rct::key bos_coster_heap_conv_robust(std::vector<MultiexpData> &data);
+
+}
+
+#endif