icicle/examples/c++/pedersen-commitment/example.cu

#include <iostream>
#include <iomanip>
#include <chrono>
#include <cassert>
#include <nvml.h>

#include "api/bn254.h"
#include "msm/msm.cuh"
using namespace bn254;

typedef point_field_t T;

// modular power
T modPow(T base, T exp) {
  T r = T::one();
  T b = base;
  T e = exp;
  while (e != T::zero()) {
      // If exp is odd, multiply the base with result
      if (T::is_odd(e)) {
          r = r * b;
      }
      // Now exp must be even, divide it by 2
      e =T::div2(e);
      b = b * b;
  }
  return r;
}

// Check if y2 is a quadratic residue using Euler's Criterion
bool quadratic_residue(T y2) {
  return modPow(y2, T::div2(T::zero() - T::one())) == T::one();
}

// modular square root adapted from:
// https://github.com/ShahjalalShohag/code-library/blob/main/Number%20Theory/Tonelli%20Shanks%20Algorithm.cpp
bool mySQRT(T a, T *result) {
  if (a == T::zero()) {
    *result = T::zero();
    return true;
  }
  if (modPow(a, T::div2(T::zero() - T::one())) != T::one() ) {
    return false; // solution does not exist
  }
  // TODO: consider special cases
  // if (p % 4 == 3) return power(a, (p + 1) / 4, p); 
  T s = T::zero() - T::one(); // p - 1, 
  T n = T::one() + T::one(); //2;
  T r = T::zero(); 
  T m;
  while (T::is_even(s)) {
    r = r + T::one();
    s = T::div2(s); //s /= 2;
  }
  // find a non-square mod p
  while (modPow(n, T::div2((T::zero() - T::one())) ) != T::zero() - T::one()) {
    n = n + T::one();
  }
  T x = modPow(a, T::div2(s + T::one()));
  T b = modPow(a, s);
  T g = modPow(n, s);
  for (;; r = m) {
    T t = b;
    for (m = T::zero(); T::lt(m,r) /* m < r*/ && t != T::one(); m = m + T::one()) t =  t * t;
    if (m == T::zero() ) {
      *result = x;
      return true;
    }
    T gs = modPow(g, modPow(T::one() + T::one(), r - m - T::one()) );
    g = gs * gs ;
    x = x * gs ;
    b =  b * g ;
  }
}

void point_near_x(T x, affine_t *point) {
  const T wb = T { weierstrass_b };
  T y2;
  while (y2 = x*x*x + wb, quadratic_residue(y2) == false)
  {
    x = x + T::one();
  };
  T y;
  bool found = mySQRT(y2, &y);
  assert(y*y == y2);
  point->x = x;
  point->y = y;
}

static int seed = 0;
static HOST_INLINE T rand_host_seed()
  {
    std::mt19937_64 generator(seed++);
    std::uniform_int_distribution<unsigned> distribution;
    
    T value;
    for (unsigned i = 0; i <  T::TLC-1 ; i++)
    // TODO: use the full range of limbs: for (unsigned i = 0; i <  T::TLC ; i++)
      value.limbs_storage.limbs[i] = distribution(generator);
    // while (lt(Field{get_modulus()}, value))
    //   value = value - Field{get_modulus()};
    return value;
  }

using FpMilliseconds = std::chrono::duration<float, std::chrono::milliseconds::period>;
#define START_TIMER(timer) auto timer##_start = std::chrono::high_resolution_clock::now();
#define END_TIMER(timer, msg) printf("%s: %.0f ms\n", msg, FpMilliseconds(std::chrono::high_resolution_clock::now() - timer##_start).count());

int main(int argc, char** argv)
{
  const unsigned N = pow(2, 10);
  std::cout << "Commitment vector size: " << N << "+1 for salt (a.k.a blinding factor)" << std::endl;
  T* xs = new T[N+1];
  
  std::cout << "Generating random points transparently using publicly chosen seed" << std::endl;
  std::cout << "Public seed prevents committer from knowing the discrete logs of points used in the commitment" << std::endl;
  seed = 1234;
  std::cout << "Using seed: " << seed << std::endl;
  std::cout << "Generating random field values" << std::endl;
  START_TIMER(gen);
  
  for (unsigned i = 0; i < N; i++) {
    xs[i] = rand_host_seed();
  }
  END_TIMER(gen, "Time to generate field values");
  std::cout << "xs[0]: " << xs[0]  << std::endl;
  std::cout << "xs[1]: " << xs[1]  << std::endl;
  
  // affine_t points[N];
  affine_t* points = new affine_t[N+1];
  std::cout << "Generating point about random field values" << std::endl;
  START_TIMER(points);
  for (unsigned i = 0; i < N+1; i++) {
    point_near_x(xs[i], &points[i]);
  }
  END_TIMER(points, "Time to generate points");
  
  std::cout << "Generating commitment vector" << std::endl;
  projective_t result;
  scalar_t* scalars = new scalar_t[N+1];
  scalar_t::rand_host_many(scalars, N);

  std::cout << "Generating salt" << std::endl;
  scalars[N] = scalar_t::rand_host();

  std::cout << "Executing MSM" << std::endl;
  auto config = msm::default_msm_config();
  START_TIMER(msm);
  bn254_msm_cuda(scalars, points, N+1, config, &result);
  END_TIMER(msm, "Time to execute MSM");

  std::cout << "Computed commitment: " << result << std::endl;

  std::cout << "Cleaning up..." << std::endl;
  delete[] xs;
  delete[] scalars;
  delete[] points;
  return 0;
}