Topic: How are you guys representing 256-bit ints in C/C++? (Read 336 times)

#define MultiplyWordsLoHi(low, high, a, b) asm  ( "mulx  %2, %0, %1;" : "=r"(low), "=r"(high) :  "gr" (a), "d" (b) : "cc");
// compute a * b = (low, high)

#define AccAdd4WordsBy4_wc(a0, a1, a2, a3, b0, b1, b2)  asm  ("addq %4, %0; adcx %5, %1; adcx %6, %2; adcq $0, %3;" : "+r"(a0), "+r"(a1), "+r"(a2), "+r"(a3) : "r"(b0), "r"(b1), "r"(b2) : "cc");
// (a0, a1, a2, a3) = (a0, a1, a2, a3) + (b0, b1, b2, 0) without carry

#define MulAcc(c, a0, a1, a, b) asm  ("mulx %3, %3, %4; addq %3, %1; adcq %4, %2; adcq $0, %0;" : "+r"(c), "+r"(a0), "+r"(a1), "=a"(a), "=d"(b) : "a"(a), "d"(b) : "cc");

#define MulAcc_11(a0, a1, c0, a, b)  asm ("mulx %3, %0, %1; addq %2, %0; adcq $0, %1;" : "+r"(a0), "+r"(a1): "r"(c0), "r"(a), "d"(b) : "cc");


//compute u * v  = r mod p
void mul(uint64_t *r, const uint64_t *u, const uint64_t *v) {
  
  uint64_t u0 = u[0];
  uint64_t u1 = u[1];
  uint64_t u2 = u[2];
  uint64_t u3 = u[3];


  uint64_t v0 = v[0];
  uint64_t v1 = v[1];
  uint64_t v2 = v[2];
  uint64_t v3 = v[3];

  uint64_t r0, r1, r2, r3, r4, r5, r6, r7;
  uint64_t z1, z2, z3, z4, z5, z6, z7, z8, z44, z66;


  z2 = z3 = z4 = z5 = z6 = z7 = z8 = r1 = r2 = r3 = r4 = r5 = r6 = r7 = 0;

  MultiplyWordsLoHi(r0, z1, u0, v0) //x1 --> r0 ok
  MultiplyWordsLoHi(z2, z3, u1, v0)
  MultiplyWordsLoHi(z4, z5, u2, v0)
  MultiplyWordsLoHi(z6, z7, u3, v0)
  MultiplyWordsLoHi(z66, z8, u3, v1)//
  AccAdd4WordsBy4_wc(z2, z4, z6, z7, z1, z3, z5)


  MulAcc_11(r1, z1, z2, u0, v1) //x1 --> r1 ok
  MultiplyWordsLoHi(z2, z3, u1, v1)
  MultiplyWordsLoHi(z44, z5, u2, v1)
  AccAdd4WordsBy4_wc(z1, z3, z5, z8, z4, z6, z7)
  AccAdd4WordsBy4_wc(z2, z44, z66, z8, z1, z3, z5)

  
  MulAcc_11(r2, z1, z2, u0, v2) //x1 --> r2 ok
  MultiplyWordsLoHi(z2, z3, u1, v2)
  MultiplyWordsLoHi(z4, z5, u2, v2)
  MultiplyWordsLoHi(z6, z7, u3, v2)
  AccAdd4WordsBy4_wc(z1, z3, z5, z7, z44, z66, z8)
  AccAdd4WordsBy4_wc(z2, z4, z6, z7, z1, z3, z5)

  MulAcc_11(r3, z1, z2, u0, v3) //x1 --> r3 ok
  MultiplyWordsLoHi(r4, z3, u1, v3)
  MultiplyWordsLoHi(r5, z5, u2, v3)
  MultiplyWordsLoHi(r6, r7, u3, v3)
  AccAdd4WordsBy4_wc(z1, z3, z5, r7, z4, z6, z7)
  AccAdd4WordsBy4_wc(r4, r5, r6, r7, z1, z3, z5) //r4, r5, r6, r7 ok
  

  //Reduction
  
  uint64_t p = 0x1000003d1;
  MultiplyWordsLoHi(z3, z4, r5, p)
  MultiplyWordsLoHi(z5, z6, r6, p)
  MultiplyWordsLoHi(z7, z8, r7, p)
  

  MulAcc_11(z1, z2, r0, r4, p)
  AccAdd4WordsBy4_wc(z2, z4, z6, z8, r1, r2, r3)

  
  uint64_t c = 0;
  AccAdd4WordsBy4_wc(z3, z5, z7, z8, z2, z4, z6)
  MulAcc(c, z1, z3, p, z8)
  
  r[0] = z1;
  r[1] = z3;
 
  
  if(c == 1){

  asm (
     "addq $1, %0; adcq $0, %1; \n"
       : "=r" (z5), "=r" (z7)
       : : "cc");

  }
  
  r[2] = z5;
  r[3] = z7;
  
  
}

##############################################################
#compute r = (a - b) mod p

void sub (uint64_t *r, const uint64_t *a, const uint64_t *b) {
  
  uint64_t a0, a1, a2, a3, b0, b1, b2, b3, p;

  a0 = a[0];
  a1 = a[1];
  a2 = a[2];
  a3 = a[3];

  b0 = b[0];
  b1 = b[1];
  b2 = b[2];
  b3 = b[3];
  
  size_t borrow = 0;

  //compute:   (a0, a1, a2, a3) = (a0, a1, a2, a3) - (b0, b1, b2, b3)  with borrow //

  asm("subq %5, %1\n\t"
   "sbbq %6, %2\n\t"
  "sbbq %7, %3\n\t"
"sbbq %8, %4\n\t"
"sbbq $0, %0"
: "+r" (borrow), "+r" (a0),  "+r" (a1), "+r" (a2), "+r" (a3)
: "r" (b0), "r" (b1), "r" (b2), "r" (b3)
: "cc");
  
  if(borrow == 0){   //if a >= b
    r[0] = a0;
    r[1] = a1;
    r[2] = a2;
    r[3] = a3;
    return;
  }
  
  //if a < b

  p = 0x1000003d1;

  //compute:  (a0, a1, a2, a3) = (a0, a1, a2, a3) - p
  asm("subq %4, %0\n\t"
   "sbbq $0, %1\n\t"
  "sbbq $0, %2\n\t"
"sbbq $0, %3"
: "+r" (a0),  "+r" (a1), "+r" (a2), "+r" (a3)
: "r" (p)
: "cc");

    r[0] = a0;
    r[1] = a1;
    r[2] = a2;
    r[3] = a3;
    
    return;
 
}
##################################################################

in the main function:

//a = a0 + a1*(2^64) + a2*(2^128) + a3*(2^196)

uint64_t  a =  {0x59f2815b16f81798, 0x029bfcdb2dce28d9, 0x55a06295ce870b07, 0x79be667ef9dcbbac};
uint64_t  b =  {0x9c47d08ffb10d4b8, 0xfd17b448a6855419, 0x5da4fbfc0e1108a8, 0x483ada7726a3c465};

uint64_t c [4];

uint64_t* ptra = &a[0];
uint64_t* ptrb = &b[0];
uint64_t* ptrc = &c[0];

sub(ptrc, ptra, ptrb);  //compute c = a - b mod p

import time


def is_odd_rec(n):
    if n == 0:
        return False
    else:
        return not is_odd_rec(n - 1)


def is_odd_mod(n):
    return n % 2 == 1


def is_odd_bit(n):
    return n & 1


start_time = time.time()
for i in range(1000000):
    is_odd_rec(40)
print(time.time() - start_time)
# output 8.055624008178711

start_time = time.time()
for i in range(1000000):
    is_odd_mod(40)
print(time.time() - start_time)
# output 0.24014711380004883

start_time = time.time()
for i in range(1000000):
    is_odd_bit(40)
print(time.time() - start_time)
# output 0.23915982246398926