initial improt: some basic C FFI seems to work

.gitignore

@@ -0,0 +1,3 @@
+*
+!*/
+!*.*

LICENSE.md

@@ -0,0 +1,5 @@
+Licensed and distributed under either of
+[MIT license](http://opensource.org/licenses/MIT) or
+[Apache License, Version 2.0](http://www.apache.org/licenses/LICENSE-2.0)
+at your option. These files may not be copied, modified, or distributed except
+according to those terms.

README.md

@@ -0,0 +1,42 @@
+Nim/C implementation of Poseidon2 over the Goldilocks field
+===========================================================
+
+Experimental implementation of the [Poseidon2][1] cryptographic hash function,
+specialized to the Goldilocks field `p=2^64-2^32+1` and `t=12`. 
+Uses a C implementation internally.
+
+Installation
+------------
+
+Use the [Nimble][2] package manager to add `poseidon2-goldilocks` to an existing
+project. Add the following to its `.nimble` file:
+
+```nim
+requires "poseidon2-goldilocks >= 0.0.1 & < 0.0.1"
+```
+
+Usage
+-----
+
+Hashing bytes into a field element with the sponge construction:
+```nim
+import poseidon2_goldilocks
+
+let input = [1'u8, 2'u8, 3'u8] # some bytes that you want to hash
+let digest: F = Sponge.digest(input) # a field element
+```
+
+Converting a hash digest (4 field elements) into bytes:
+```nim
+let output: array[32, byte] = digest.toBytes
+```
+
+Combining field elements, useful for constructing a binary Merkle tree:
+```nim
+let left  = Sponge.digest([1'u8, 2'u8, 3'u8])
+let right = Sponge.digest([4'u8, 5'u8, 6'u8])
+let combination = compress(left, right)
+```
+
+[1]: https://eprint.iacr.org/2023/323.pdf
+[2]: https://github.com/nim-lang/nimble

cbits/goldilocks.c

@@ -0,0 +1,417 @@
+
+#include <stdint.h>
+#include <stdio.h>      // for testing only
+
+#include "goldilocks.h"
+
+//==============================================================================
+// *** Goldilocks field ***
+
+uint64_t goldilocks_neg(uint64_t x) {
+  return (x==0) ? 0 : (GOLDILOCKS_PRIME - x);
+}
+
+uint64_t goldilocks_add(uint64_t x, uint64_t y) {
+  uint64_t z = x + y;
+  return ( (z >= GOLDILOCKS_PRIME) || (z<x) ) ? (z - GOLDILOCKS_PRIME) : z;
+}
+
+uint64_t goldilocks_sub(uint64_t x, uint64_t y) {
+  uint64_t z = x - y;
+  return (z > x) ? (z + GOLDILOCKS_PRIME) : z;
+}
+
+uint64_t goldilocks_sub_safe(uint64_t x, uint64_t y) {
+  return goldilocks_add( x , goldilocks_neg(y) );
+}
+
+/*
+
+// add together 3 field elements
+uint64_t goldilocks_add3( uint64_t x0, uint64_t x1, uint64_t x2 ) {
+  uint64_t x01 = goldilocks_add( x0 , x1  );
+  return goldilocks_add( x01, x2 );
+}
+
+//--------------------------------------
+
+uint64_t goldilocks_div_by_2(uint64_t x) {
+  return (x & 1) ? (x/2 + 0x7fffffff80000001) : (x/2);
+}
+
+uint64_t goldilocks_div_by_3(uint64_t x) {
+  uint64_t m = x % 3;
+  uint64_t r;
+  switch(m) {
+    case 0:
+      r = (x/3);
+      break;
+    case 1:
+      r = (x/3 + 0xaaaaaaaa00000001);      // (x+2*p) / 3 = x/3 + (2*p+1)/3
+      break;
+    case 2:
+      r = (x/3 + 0x5555555500000001);      // (x+p) / 3 = x/3 + (p+1)/3 
+      break;
+  }
+  return r;
+}
+
+uint64_t goldilocks_div_by_4(uint64_t x) {
+  return goldilocks_div_by_2(goldilocks_div_by_2(x));
+}
+
+*/
+
+//--------------------------------------
+
+uint64_t goldilocks_rdc(__uint128_t x) {
+  // x = n0 + 2^64 * n1 + 2^96 * n2
+  uint64_t n0 = (uint64_t)x;
+  uint64_t n1 = (x >> 64) & 0xffffffff;
+  uint64_t n2 = (x >> 96);
+  
+  uint64_t mid = (n1 << 32) - n1;     // (2^32 - 1) * n1
+  uint64_t tmp = n0 + mid;
+  if (tmp < n0) { tmp -= GOLDILOCKS_PRIME; }
+
+  uint64_t res = tmp - n2;
+  if (res > tmp) { res += GOLDILOCKS_PRIME; }
+  return (res >= GOLDILOCKS_PRIME) ? (res - GOLDILOCKS_PRIME) : res;
+}
+
+// reduce to 64-bit, but it can be still bigger than `p`
+uint64_t goldilocks_rdc_to_uint64(__uint128_t x) {
+  // x = n0 + 2^64 * n1 + 2^96 * n2
+  uint64_t n0 = (uint64_t)x;
+  uint64_t n1 = (x >> 64) & 0xffffffff;
+  uint64_t n2 = (x >> 96);
+  
+  uint64_t mid = (n1 << 32) - n1;     // (2^32 - 1) * n1
+  uint64_t tmp = n0 + mid;
+  if (tmp < n0) { tmp -= GOLDILOCKS_PRIME; }
+
+  uint64_t res = tmp - n2;
+  if (res > tmp) { res += GOLDILOCKS_PRIME; }
+  return res;
+}
+
+// we assume x < 2^96
+uint64_t goldilocks_rdc_small(__uint128_t x) {
+  // x = n0 + 2^64 * n1
+  uint64_t n0 = (uint64_t)x;
+  uint64_t n1 = (x >> 64);
+
+  uint64_t mid = (n1 << 32) - n1;     // (2^32 - 1) * n1
+  uint64_t tmp = n0 + mid;
+  if (tmp < n0) { tmp -= GOLDILOCKS_PRIME; }
+
+  uint64_t res = tmp;
+  return (res >= GOLDILOCKS_PRIME) ? (res - GOLDILOCKS_PRIME) : res;
+}
+
+//--------------------------------------
+
+uint64_t goldilocks_mul(uint64_t x, uint64_t y) {
+  __uint128_t z = (__uint128_t)x * (__uint128_t)y;
+  return goldilocks_rdc(z); 
+}
+
+uint64_t goldilocks_mul_add128(uint64_t x, uint64_t y, __uint128_t z) {
+  __uint128_t w = (__uint128_t)x * (__uint128_t)y + z;
+  return goldilocks_rdc(w); 
+}
+
+uint64_t goldilocks_sqr(uint64_t x) {
+  __uint128_t z = (__uint128_t)x * (__uint128_t)x;
+  return goldilocks_rdc(z); 
+}
+
+uint64_t goldilocks_sqr_add(uint64_t x, uint64_t y) {
+  __uint128_t z = (__uint128_t)x * x + y;
+  return goldilocks_rdc(z); 
+}
+
+// only reduce to uint64, not to [0..p-1]
+uint64_t goldilocks_sqr_add_to_uint64(uint64_t x, uint64_t y) {
+  __uint128_t z = (__uint128_t)x * x + y;
+  return goldilocks_rdc_to_uint64(z); 
+}
+
+uint64_t goldilocks_mul_small(uint64_t x, uint32_t y) {
+  __uint128_t z = (__uint128_t)x * (__uint128_t)y;
+  return goldilocks_rdc_small(z); 
+}
+
+//==============================================================================
+// *** debugging ***
+
+void debug_print_state(const char *msg, int n, uint64_t *state) {
+  printf("-----------------\n");
+  printf("%s\n",msg);
+  for(int i=0;i<n;i++) {
+    printf(" - 0x%016llx = %llu\n",state[i],state[i]);
+  }
+}
+
+//==============================================================================
+// *** Poseidon2 ***
+//
+// compatible with <https://github.com/HorizenLabs/poseidon2>
+// NOT compatible with <https://extgit.iaik.tugraz.at/krypto/zkfriendlyhashzoo>
+// (presumably they use different constants or whatever)
+//
+
+#include "poseidon2_constants.inc"
+
+/* 
+poseidon2 test vector (permutation of [0..11]) 
+----------------------------------------------
+from <https://github.com/HorizenLabs/poseidon2/blob/main/plain_implementations/src/poseidon2/poseidon2.rs#L284>
+
+0x01eaef96bdf1c0c1
+0x1f0d2cc525b2540c
+0x6282c1dfe1e0358d
+0xe780d721f698e1e6
+0x280c0b6f753d833b
+0x1b942dd5023156ab
+0x43f0df3fcccb8398
+0xe8e8190585489025
+0x56bdbf72f77ada22
+0x7911c32bf9dcd705
+0xec467926508fbe67
+0x6a50450ddf85a6ed
+*/
+
+uint64_t goldilocks_poseidon2_sbox(uint64_t x0, uint64_t rc) {
+  uint64_t x  = goldilocks_add( x0 , rc );
+  uint64_t x2 = goldilocks_sqr( x       );
+  uint64_t x4 = goldilocks_sqr( x2      );
+  uint64_t x6 = goldilocks_mul( x4 , x2 );
+  uint64_t x7 = goldilocks_mul( x6 , x  );
+  return x7;
+}
+
+// remark: (p-1)^2 + 12*(p-1) does not overflow in 2^128
+void goldilocks_poseidon2_internal_diffusion(uint64_t *inp, uint64_t *out) {
+  __uint128_t s0 = inp[0];
+  __uint128_t s1 = inp[6];
+  s0 += inp[1]; s1 += inp[7];
+  s0 += inp[2]; s1 += inp[8];
+  s0 += inp[3]; s1 += inp[9];
+  s0 += inp[4]; s1 += inp[10];
+  s0 += inp[5]; s1 += inp[11];
+//  uint64_t s = goldilocks_rdc_small( s0 + s1 );
+  __uint128_t s = s0 + s1;
+
+  for(int i=0; i<12; i++) { 
+    out[i] = goldilocks_mul_add128( inp[i] , internal_diag[i] , s );
+  }
+}
+
+//--------------------------------------
+
+/*
+
+// multiplies a vector of size 4 by the 4x4 MDS matrix on the left:
+//
+//       [ 5 7 1 3 ]
+//  M4 = [ 4 6 1 1 ]
+//       [ 1 3 5 7 ]
+//       [ 1 1 4 6 ]
+//
+void goldilocks_mul_by_M4(uint64_t *inp, uint64_t *out) {
+  uint64_t a = inp[0];
+  uint64_t b = inp[1];
+  uint64_t c = inp[2];
+  uint64_t d = inp[3];
+
+  uint64_t a2 = goldilocks_add( a  , a  );
+  uint64_t a4 = goldilocks_add( a2 , a2 );
+  uint64_t a5 = goldilocks_add( a4 , a  );
+
+  uint64_t b2 = goldilocks_add( b  , b  );
+  uint64_t b3 = goldilocks_add( b2 , b  );
+  uint64_t b6 = goldilocks_add( b3 , b3 );
+  uint64_t b7 = goldilocks_add( b6 , b  );
+
+  uint64_t c2 = goldilocks_add( c  , c  );
+  uint64_t c4 = goldilocks_add( c2 , c2 );
+  uint64_t c5 = goldilocks_add( c4 , c  );
+
+  uint64_t d2 = goldilocks_add( d  , d  );
+  uint64_t d3 = goldilocks_add( d2 , d  );
+  uint64_t d6 = goldilocks_add( d3 , d3 );
+  uint64_t d7 = goldilocks_add( d6 , d  );
+
+  out[0] = goldilocks_add( goldilocks_add( a5 , b7 ) , goldilocks_add( c  , d3 ) );
+  out[1] = goldilocks_add( goldilocks_add( a4 , b6 ) , goldilocks_add( c  , d  ) );
+  out[2] = goldilocks_add( goldilocks_add( a  , b3 ) , goldilocks_add( c5 , d7 ) );
+  out[3] = goldilocks_add( goldilocks_add( a  , b  ) , goldilocks_add( c4 , d6 ) );
+}
+
+// returns 2*a + b + c
+uint64_t goldilocks_weighted_add_211(uint64_t a, uint64_t b, uint64_t c) {
+  uint64_t a2 = goldilocks_add( a , a );
+  uint64_t bc = goldilocks_add( b , c );
+  return goldilocks_add( a2 , bc );
+}
+
+// multiplies by 12x12 block-circulant matrix [2*M4, M4, M4]
+void goldilocks_poseidon2_external_diffusion(uint64_t *inp, uint64_t *out) {
+  uint64_t us[4];
+  uint64_t vs[4];
+  uint64_t ws[4];
+
+  goldilocks_mul_by_M4( inp + 0 , us );
+  goldilocks_mul_by_M4( inp + 4 , vs );
+  goldilocks_mul_by_M4( inp + 8 , ws );
+
+  out[0]  = goldilocks_weighted_add_211( us[0] , vs[0] , ws[0] );
+  out[1]  = goldilocks_weighted_add_211( us[1] , vs[1] , ws[1] );
+  out[2]  = goldilocks_weighted_add_211( us[2] , vs[2] , ws[2] );
+  out[3]  = goldilocks_weighted_add_211( us[3] , vs[3] , ws[3] );
+
+  out[4]  = goldilocks_weighted_add_211( vs[0] , ws[0] , us[0] );
+  out[5]  = goldilocks_weighted_add_211( vs[1] , ws[1] , us[1] );
+  out[6]  = goldilocks_weighted_add_211( vs[2] , ws[2] , us[2] );
+  out[7]  = goldilocks_weighted_add_211( vs[3] , ws[3] , us[3] );
+
+  out[ 8] = goldilocks_weighted_add_211( ws[0] , us[0] , vs[0] );
+  out[ 9] = goldilocks_weighted_add_211( ws[1] , us[1] , vs[1] );
+  out[10] = goldilocks_weighted_add_211( ws[2] , us[2] , vs[2] );
+  out[11] = goldilocks_weighted_add_211( ws[3] , us[3] , vs[3] );
+}
+
+*/
+
+//--------------------------------------
+
+// multiplies a vector of size 4 by the 4x4 MDS matrix on the left
+void uint64_mul_by_M4(uint64_t *inp, uint64_t *out) {
+  uint64_t a = inp[0];
+  uint64_t b = inp[1];
+  uint64_t c = inp[2];
+  uint64_t d = inp[3];
+
+  uint64_t a4 = a << 2;
+  uint64_t a5 = a4 + a;
+
+  uint64_t b2 = b  + b  ;
+  uint64_t b3 = b2 + b  ;
+  uint64_t b6 = b3 + b3 ;
+  uint64_t b7 = b6 + b  ;
+
+  uint64_t c4 = c << 2 ;
+  uint64_t c5 = c4 + c ;
+
+  uint64_t d2 = d  + d  ;
+  uint64_t d3 = d2 + d  ;
+  uint64_t d6 = d3 + d3 ;
+  uint64_t d7 = d6 + d  ;
+
+  out[0] = a5 + b7 + c  + d3 ;
+  out[1] = a4 + b6 + c  + d  ;
+  out[2] = a  + b3 + c5 + d7 ;
+  out[3] = a  + b  + c4 + d6 ;
+}
+
+// multiplies by 12x12 block-circulant matrix [2*M4, M4, M4]
+void uint64_mul_by_poseidon2_circulant12(uint64_t *inp, uint64_t *out) {
+  uint64_t us[4];
+  uint64_t vs[4];
+  uint64_t ws[4];
+
+  uint64_mul_by_M4( inp + 0 , us );
+  uint64_mul_by_M4( inp + 4 , vs );
+  uint64_mul_by_M4( inp + 8 , ws );
+
+  out[0]  = 2*us[0] + vs[0] + ws[0];
+  out[1]  = 2*us[1] + vs[1] + ws[1];
+  out[2]  = 2*us[2] + vs[2] + ws[2];
+  out[3]  = 2*us[3] + vs[3] + ws[3];
+
+  out[4]  = us[0] + 2*vs[0] + ws[0];
+  out[5]  = us[1] + 2*vs[1] + ws[1];
+  out[6]  = us[2] + 2*vs[2] + ws[2];
+  out[7]  = us[3] + 2*vs[3] + ws[3];
+
+  out[ 8] = us[0] + vs[0] + 2*ws[0];
+  out[ 9] = us[1] + vs[1] + 2*ws[1];
+  out[10] = us[2] + vs[2] + 2*ws[2];
+  out[11] = us[3] + vs[3] + 2*ws[3];
+}
+
+void goldilocks_poseidon2_external_diffusion_split(uint64_t *inp, uint64_t *out) {
+  uint64_t inp_lo[12];
+  uint64_t inp_hi[12];
+  uint64_t out_lo[12];
+  uint64_t out_hi[12];
+
+  for(int i=0; i<12; i++) { 
+    uint64_t x = inp[i];
+    inp_lo[i] = x & 0xffffffff;
+    inp_hi[i] = x >> 32;
+  }
+
+  uint64_mul_by_poseidon2_circulant12(inp_lo, out_lo);
+  uint64_mul_by_poseidon2_circulant12(inp_hi, out_hi);
+
+  for(int i=0; i<12; i++) {
+    __uint128_t x = (((__uint128_t)out_hi[i]) << 32) + out_lo[i];
+    out[i] = goldilocks_rdc_small(x);
+  }
+}
+
+//--------------------------------------
+
+// 0 <= round_idx < 22
+void goldilocks_poseidon2_internal_round(int round_idx, uint64_t *state) {
+  state[0] = goldilocks_poseidon2_sbox( state[0] , internal_round_consts[round_idx] );
+  goldilocks_poseidon2_internal_diffusion( state, state );
+}
+
+void goldilocks_poseidon2_external_round(const uint64_t *rcs, uint64_t *state) {
+  for (int i=0; i<12; i++) {
+    state[i] = goldilocks_poseidon2_sbox( state[i] , rcs[i] );
+  }
+  goldilocks_poseidon2_external_diffusion_split( state, state );
+}
+
+void goldilocks_poseidon2_permutation(uint64_t *state) {
+  goldilocks_poseidon2_external_diffusion_split( state, state );
+  goldilocks_poseidon2_external_round( intial_round_consts + 0  , state );
+  goldilocks_poseidon2_external_round( intial_round_consts + 12 , state );
+  goldilocks_poseidon2_external_round( intial_round_consts + 24 , state );
+  goldilocks_poseidon2_external_round( intial_round_consts + 36 , state );
+  for(int idx=0; idx<22; idx++) {
+    goldilocks_poseidon2_internal_round( idx, state );
+  }
+  goldilocks_poseidon2_external_round( final_round_consts + 0  , state );
+  goldilocks_poseidon2_external_round( final_round_consts + 12 , state );
+  goldilocks_poseidon2_external_round( final_round_consts + 24 , state );
+  goldilocks_poseidon2_external_round( final_round_consts + 36 , state );
+}
+
+//------------------------------------------------------------------------------
+
+// compression function: input is two 4-element vector of field elements, 
+// and the output is a vector of 4 field elements
+void goldilocks_poseidon2_keyed_compress(uint64_t key, const uint64_t *x, const uint64_t *y, uint64_t *out) {
+  uint64_t state[12];
+  for(int i=0; i<4; i++) {
+    state[i  ] = x[i];
+    state[i+4] = y[i];
+    state[i+8] = 0;
+  }
+  state[8] = key;
+  goldilocks_poseidon2_permutation(state);
+  for(int i=0; i<4; i++) {
+    out[i] = state[i];
+  }
+}
+
+void goldilocks_poseidon2_compress(const uint64_t *x, const uint64_t *y, uint64_t *out) {
+  goldilocks_poseidon2_keyed_compress(0, x, y, out);
+}
+
+//------------------------------------------------------------------------------

cbits/goldilocks.h

@@ -0,0 +1,32 @@
+
+#include <stdint.h>
+
+//------------------------------------------------------------------------------
+
+#define GOLDILOCKS_PRIME 0xffffffff00000001
+
+//------------------------------------------------------------------------------
+
+uint64_t goldilocks_neg(uint64_t x);
+uint64_t goldilocks_add(uint64_t x, uint64_t y);
+uint64_t goldilocks_sub(uint64_t x, uint64_t y);
+uint64_t goldilocks_mul(uint64_t x, uint64_t y);
+uint64_t goldilocks_mul_small(uint64_t x, uint32_t y);
+
+uint64_t goldilocks_div_by_2(uint64_t x);
+uint64_t goldilocks_div_by_3(uint64_t x);
+uint64_t goldilocks_div_by_4(uint64_t x);
+
+uint64_t goldilocks_add3(uint64_t x, uint64_t y, uint64_t z);
+
+//uint64_t goldilocks_rdc(__uint128_t x);
+
+//------------------------------------------------------------------------------
+
+void goldilocks_poseidon2_permutation(uint64_t *state);
+void goldilocks_monolith_permutation (uint64_t *state);
+
+void monolith_print_sbox_table();
+void monolith_print_sbox_table_c_format();
+
+//------------------------------------------------------------------------------

cbits/poseidon2_constants.inc

@@ -0,0 +1,165 @@
+
+// based on the reference implementation:
+// <https://github.com/HorizenLabs/poseidon2/blob/main/plain_implementations/src/poseidon2/poseidon2_instance_goldilocks.rs>
+
+#include <stdint.h>
+
+//------------------------------------------------------------------------------
+
+// internal matrices are the constant 1 matrix + diagonal of these
+// (in particular the diagonals of the resulting matrix will be these plus one!)
+const uint64_t internal_diag[12] = 
+  { 0xc3b6c08e23ba9300
+  , 0xd84b5de94a324fb6
+  , 0x0d0c371c5b35b84f
+  , 0x7964f570e7188037
+  , 0x5daf18bbd996604b
+  , 0x6743bc47b9595257
+  , 0x5528b9362c59bb70
+  , 0xac45e25b7127b68b
+  , 0xa2077d7dfbb606b5
+  , 0xf3faac6faee378ae
+  , 0x0c6388b51545e883
+  , 0xd27dbb6944917b60
+  };
+
+//------------------------------------------------------------------------------
+
+const uint64_t intial_round_consts[4*12] = 
+  { 0x13dcf33aba214f46
+  , 0x30b3b654a1da6d83
+  , 0x1fc634ada6159b56
+  , 0x937459964dc03466
+  , 0xedd2ef2ca7949924
+  , 0xede9affde0e22f68
+  , 0x8515b9d6bac9282d
+  , 0x6b5c07b4e9e900d8
+  , 0x1ec66368838c8a08
+  , 0x9042367d80d1fbab
+  , 0x400283564a3c3799
+  , 0x4a00be0466bca75e
+  // -----------------
+  , 0x7913beee58e3817f
+  , 0xf545e88532237d90
+  , 0x22f8cb8736042005
+  , 0x6f04990e247a2623
+  , 0xfe22e87ba37c38cd
+  , 0xd20e32c85ffe2815
+  , 0x117227674048fe73
+  , 0x4e9fb7ea98a6b145
+  , 0xe0866c232b8af08b
+  , 0x00bbc77916884964
+  , 0x7031c0fb990d7116
+  , 0x240a9e87cf35108f
+  // -----------------
+  , 0x2e6363a5a12244b3
+  , 0x5e1c3787d1b5011c
+  , 0x4132660e2a196e8b
+  , 0x3a013b648d3d4327
+  , 0xf79839f49888ea43
+  , 0xfe85658ebafe1439
+  , 0xb6889825a14240bd
+  , 0x578453605541382b
+  , 0x4508cda8f6b63ce9
+  , 0x9c3ef35848684c91
+  , 0x0812bde23c87178c
+  , 0xfe49638f7f722c14
+  // -----------------
+  , 0x8e3f688ce885cbf5
+  , 0xb8e110acf746a87d
+  , 0xb4b2e8973a6dabef
+  , 0x9e714c5da3d462ec
+  , 0x6438f9033d3d0c15
+  , 0x24312f7cf1a27199
+  , 0x23f843bb47acbf71
+  , 0x9183f11a34be9f01
+  , 0x839062fbb9d45dbf
+  , 0x24b56e7e6c2e43fa
+  , 0xe1683da61c962a72
+  , 0xa95c63971a19bfa7
+  };
+
+//------------------------------------------------------------------------------
+
+const uint64_t final_round_consts[4*12] = 
+  { 0xc68be7c94882a24d
+  , 0xaf996d5d5cdaedd9
+  , 0x9717f025e7daf6a5
+  , 0x6436679e6e7216f4
+  , 0x8a223d99047af267
+  , 0xbb512e35a133ba9a
+  , 0xfbbf44097671aa03
+  , 0xf04058ebf6811e61
+  , 0x5cca84703fac7ffb
+  , 0x9b55c7945de6469f
+  , 0x8e05bf09808e934f
+  , 0x2ea900de876307d7
+  //------------------
+  , 0x7748fff2b38dfb89
+  , 0x6b99a676dd3b5d81
+  , 0xac4bb7c627cf7c13
+  , 0xadb6ebe5e9e2f5ba
+  , 0x2d33378cafa24ae3
+  , 0x1e5b73807543f8c2
+  , 0x09208814bfebb10f
+  , 0x782e64b6bb5b93dd
+  , 0xadd5a48eac90b50f
+  , 0xadd4c54c736ea4b1
+  , 0xd58dbb86ed817fd8
+  , 0x6d5ed1a533f34ddd
+  //------------------
+  , 0x28686aa3e36b7cb9
+  , 0x591abd3476689f36
+  , 0x047d766678f13875
+  , 0xa2a11112625f5b49
+  , 0x21fd10a3f8304958
+  , 0xf9b40711443b0280
+  , 0xd2697eb8b2bde88e
+  , 0x3493790b51731b3f
+  , 0x11caf9dd73764023
+  , 0x7acfb8f72878164e
+  , 0x744ec4db23cefc26
+  , 0x1e00e58f422c6340
+  //------------------
+  , 0x21dd28d906a62dda
+  , 0xf32a46ab5f465b5f
+  , 0xbfce13201f3f7e6b
+  , 0xf30d2e7adb5304e2
+  , 0xecdf4ee4abad48e9
+  , 0xf94e82182d395019
+  , 0x4ee52e3744d887c5
+  , 0xa1341c7cac0083b2
+  , 0x2302fb26c30c834a
+  , 0xaea3c587273bf7d3
+  , 0xf798e24961823ec7
+  , 0x962deba3e9a2cd94
+  };
+
+//------------------------------------------------------------------------------
+
+const uint64_t internal_round_consts[22] = 
+  { 0x4adf842aa75d4316
+  , 0xf8fbb871aa4ab4eb
+  , 0x68e85b6eb2dd6aeb
+  , 0x07a0b06b2d270380
+  , 0xd94e0228bd282de4
+  , 0x8bdd91d3250c5278
+  , 0x209c68b88bba778f
+  , 0xb5e18cdab77f3877
+  , 0xb296a3e808da93fa
+  , 0x8370ecbda11a327e
+  , 0x3f9075283775dad8
+  , 0xb78095bb23c6aa84
+  , 0x3f36b9fe72ad4e5f
+  , 0x69bc96780b10b553
+  , 0x3f1d341f2eb7b881
+  , 0x4e939e9815838818
+  , 0xda366b3ae2a31604
+  , 0xbc89db1e7287d509
+  , 0x6102f411f9ef5659
+  , 0x58725c5e7ac1f0ab
+  , 0x0df5856c798883e7
+  , 0xf7bb62a8da4c961b
+  };
+
+//------------------------------------------------------------------------------

poseidon2/goldilocks.nim

@@ -0,0 +1,8 @@
+{. compile: "../cbits/goldilocks.c" .}
+
+import ./types
+
+func neg* (x: F   ): F {. header: "../cbits/goldilocks.h", importc: "goldilocks_neg", cdecl .}
+func `+`* (x, y: F): F {. header: "../cbits/goldilocks.h", importc: "goldilocks_add", cdecl .}
+func `-`* (x, y: F): F {. header: "../cbits/goldilocks.h", importc: "goldilocks_sub", cdecl .}
+func `*`* (x, y: F): F {. header: "../cbits/goldilocks.h", importc: "goldilocks_mul", cdecl .}

poseidon2/types.nim

@@ -0,0 +1,29 @@
+
+#-------------------------------------------------------------------------------
+
+type F* = distinct uint64
+
+func `==`* (x, y: F): bool =
+  return (uint64(x) == uint64(y))
+
+func fromF* (x: F): uint64 = 
+  return uint64(x)
+
+func toF* (x: uint64): F = 
+  return F(x)
+
+#-------------------------------------------------------------------------------
+
+const zero* : F = toF(0)
+const one*  : F = toF(1)
+const two*  : F = toF(2)
+
+#-------------------------------------------------------------------------------
+
+type F4*  = array[4 , F]
+type F12* = array[12, F]
+
+type Digest* = distinct F4
+type State*  = distinct F12
+
+#-------------------------------------------------------------------------------

poseidon2_goldilocks.nim

@@ -0,0 +1,17 @@
+# import poseidon2/types
+# import poseidon2/io
+# import poseidon2/sponge
+# import poseidon2/compress
+# import poseidon2/merkle
+# import poseidon2/spongemerkle
+
+# export sponge
+# export compress
+# export merkle
+# export spongemerkle
+# export fromBytes
+# export toBytes
+# export toF
+# export elements
+# export types
+

poseidon2_goldilocks.nimble

@@ -0,0 +1,5 @@
+version = "0.0.1"
+author = "Balazs Komuves"
+description = "Poseidon2 hash function over the goldilocks field"
+license = "MIT"
+

reference/Goldilocks.hs

@@ -0,0 +1,93 @@
+
+-- | Reference (slow) implementation of the Goldilocks prime field
+
+{-# LANGUAGE BangPatterns, NumericUnderscores #-}
+module Goldilocks where
+
+--------------------------------------------------------------------------------
+
+import Prelude hiding ( div )
+import qualified Prelude
+
+import Data.Bits
+import Data.Ratio
+
+import Text.Printf
+
+--------------------------------------------------------------------------------
+
+type F = Goldilocks
+
+newtype Goldilocks 
+  = Goldilocks Integer 
+  deriving Eq
+
+instance Show Goldilocks where
+  show (Goldilocks k) = printf "0x%016x" k
+
+--------------------------------------------------------------------------------
+
+instance Num Goldilocks where
+  fromInteger = mkGoldilocks
+  negate = neg
+  (+)    = add
+  (-)    = sub
+  (*)    = mul
+  abs    = id
+  signum _ = Goldilocks 1
+
+instance Fractional Goldilocks where
+  fromRational y = fromInteger (numerator y) `div` fromInteger (denominator y)
+  recip  = inv
+  (/)    = div
+
+--------------------------------------------------------------------------------
+
+-- | @p = 2^64 - 2^32 + 1@
+goldilocksPrime :: Integer
+goldilocksPrime = 0x_ffff_ffff_0000_0001
+
+modp :: Integer -> Integer
+modp a = mod a goldilocksPrime
+
+mkGoldilocks :: Integer -> Goldilocks
+mkGoldilocks = Goldilocks . modp
+
+--------------------------------------------------------------------------------
+
+neg :: Goldilocks -> Goldilocks
+neg (Goldilocks k) = mkGoldilocks (negate k) 
+
+add :: Goldilocks -> Goldilocks -> Goldilocks
+add (Goldilocks a) (Goldilocks b) = mkGoldilocks (a+b) 
+
+sub :: Goldilocks -> Goldilocks -> Goldilocks
+sub (Goldilocks a) (Goldilocks b) = mkGoldilocks (a-b) 
+
+sqr :: Goldilocks -> Goldilocks
+sqr x = mul x x
+
+mul :: Goldilocks -> Goldilocks -> Goldilocks
+mul (Goldilocks a) (Goldilocks b) = mkGoldilocks (a*b) 
+
+inv :: Goldilocks -> Goldilocks
+inv x = pow x (goldilocksPrime - 2)
+
+div :: Goldilocks -> Goldilocks -> Goldilocks
+div a b = mul a (inv b)
+
+--------------------------------------------------------------------------------
+
+pow :: Goldilocks -> Integer -> Goldilocks
+pow x e 
+  | e == 0    = 1
+  | e <  0    = pow (inv x) (negate e)
+  | otherwise = go 1 x e
+  where
+    go !acc _  0     = acc
+    go !acc !s !expo = case expo .&. 1 of
+      0 -> go acc     (sqr s) (shiftR expo 1)
+      _ -> go (acc*s) (sqr s) (shiftR expo 1)
+
+--------------------------------------------------------------------------------
+

reference/README.md

@@ -0,0 +1,6 @@
+
+Haskell reference implementation
+--------------------------------
+
+Used for generating test cases (as Nim doesn't support integers larger than 64 bit)
+

reference/TestGen.hs

@@ -0,0 +1,87 @@
+
+-- | Generate test cases for Nim
+
+module TestGen where
+
+--------------------------------------------------------------------------------
+
+import System.IO
+
+import Goldilocks
+
+--------------------------------------------------------------------------------
+
+centered :: Integer -> Integer -> [Integer]
+centered center width = [center-width .. center+width]
+
+mkTestFieldElems :: Integer -> [F]
+mkTestFieldElems width = map fromInteger $ concat
+  [ centered (0   ) width
+  , centered (2^16) width
+  , centered (2^31) width
+  , centered (2^32) width
+  , centered (2^33) width
+  , centered (2^48) width
+  , centered (2^63) width
+  ]
+
+testFieldElems :: [F]
+testFieldElems = mkTestFieldElems 7
+
+testFieldPairs :: [(F,F)]
+testFieldPairs = [ (x,y) | x<-list, y<-list ] where
+  list = mkTestFieldElems 3
+
+--------------------------------------------------------------------------------
+
+nimShow :: F -> String
+nimShow x = show x ++ "'u64"
+
+nimShowPair :: (F,F) -> String
+nimShowPair (x,y) = "( " ++ nimShow x ++ " , " ++ nimShow y ++ " )"
+
+nimShowTriple :: (F,F,F) -> String
+nimShowTriple (x,y,z) = "( " ++ nimShow x ++ " , " ++ nimShow y ++ " , " ++ nimShow z ++ " )"
+
+showPairs :: [(F,F)] -> [String]
+showPairs xys = zipWith (++) prefix (map nimShowPair xys) where
+  prefix = "  [ " : repeat "  , "
+
+showTriples :: [(F,F,F)] -> [String]
+showTriples xyzs = zipWith (++) prefix (map nimShowTriple xyzs) where
+  prefix = "  [ " : repeat "  , "
+
+----------------------------------------
+
+unary :: String -> (F -> F) -> [F] -> String
+unary varname f xs = unlines (header : stuff ++ footer) where
+  header = "const " ++ varname ++ "* : array[" ++ show (length xs) ++ ", tuple[x:uint64, y:uint64]] = " 
+  footer = ["  ]",""]
+  stuff  = showPairs [ (x, f x) | x<-xs ]
+
+binary :: String -> (F -> F -> F) -> [(F,F)] -> String
+binary varname f xys = unlines (header : stuff ++ footer) where
+  header = "const " ++ varname ++ "* : array[" ++ show (length xys) ++ ", tuple[x:uint64, y:uint64, z:uint64]] = " 
+  footer = ["  ]",""]
+  stuff  = showTriples [ (x, y, f x y) | (x,y)<-xys ]
+
+--------------------------------------------------------------------------------
+
+printTests :: IO ()
+printTests = hPrintTests stdout
+
+hPrintTests :: Handle -> IO ()
+hPrintTests h = hPutStrLn h $ unlines 
+  [ unary  "testcases_neg" negate testFieldElems
+  , binary "testcases_add" (+)    testFieldPairs
+  , binary "testcases_sub" (-)    testFieldPairs
+  , binary "testcases_mul" (*)    testFieldPairs
+  ]
+
+writeTests :: IO ()
+writeTests = withFile "fieldTestCases.nim" WriteMode $ \h -> do
+  hPutStrLn   h "# generated by TestGen.hs\n"
+  -- hPutStrLn   h "import poseidon2/types\n"
+  hPrintTests h
+
+--------------------------------------------------------------------------------

tests/nim.cfg

@@ -0,0 +1 @@
+--path:".."

tests/poseidon2/reference.nim

@@ -0,0 +1,5 @@
+
+# import poseidon2/types
+# 
+# func add(x, y: F): F = 
+#   (uint)

tests/poseidon2/testField.nim

@@ -0,0 +1,58 @@
+
+import std/unittest
+# import std/sequtils
+
+import poseidon2/types
+import poseidon2/goldilocks
+
+import ./fieldTestCases
+
+#-------------------------------------------------------------------------------
+
+suite "field":
+
+  test "negation":
+    var ok = true
+    for (x0,y0) in testcases_neg:
+      let x = toF(x0)
+      let y = toF(y0)
+      if neg(x) != y: 
+        ok = false
+        break
+    check ok
+    
+  test "addition":
+    var ok = true
+    for (x0,y0,z0) in testcases_add:
+      let x = toF(x0)
+      let y = toF(y0)
+      let z = toF(z0)
+      if x + y != z: 
+        ok = false
+        break
+    check ok
+
+  test "subtraction":
+    var ok = true
+    for (x0,y0,z0) in testcases_sub:
+      let x = toF(x0)
+      let y = toF(y0)
+      let z = toF(z0)
+      if x - y != z: 
+        ok = false
+        break
+    check ok
+
+  test "multiplication":
+    var ok = true
+    for (x0,y0,z0) in testcases_mul:
+      let x = toF(x0)
+      let y = toF(y0)
+      let z = toF(z0)
+      if x * y != z: 
+        ok = false
+        break
+    check ok
+
+#-------------------------------------------------------------------------------
+

tests/test.nim

@@ -0,0 +1,11 @@
+#import ./poseidon2/testPermutation
+#import ./poseidon2/testSponge
+#import ./poseidon2/testCompress
+#import ./poseidon2/testMerkle
+#import ./poseidon2/testSpongeMerkle
+#import ./poseidon2/testIo
+#import ./poseidon2/testReadme
+
+import ./poseidon2/testField
+
+{.warning[UnusedImport]: off.}