/* 
 * integer_math.c
 * implementations for integer math routines needed for the 
 * ECDSA and RSA  algorithm. implemented using the array of  words big 
 * integer notation
 * 
 */


#include "bb_nn.h"
#include "binary_field.h"
#include "poly_math.h"

/*
 * multiply digits a = b*c
 */
static void bigint_digit_mult(bigint_digit a[2], bigint_digit b, 
			      bigint_digit c){
  bigint_digit t, u;
  bigint_half_digit bHigh, bLow, cHigh, cLow;
  
  bHigh = (bigint_half_digit) HIGHER_HALF (b);
  bLow = (bigint_half_digit) LOWER_HALF (b);
  cHigh = (bigint_half_digit) HIGHER_HALF (c);
  cLow = (bigint_half_digit) LOWER_HALF (c);

  a[0] = (bigint_digit) bLow * (bigint_digit)cLow;
  t = (bigint_digit)bLow * (bigint_digit) cHigh;
  u = (bigint_digit)bHigh * (bigint_digit) cLow;
  a[1] = (bigint_digit) bHigh * (bigint_digit) cHigh;
  
  if((t += u) < u)
    a[1] += TO_HIGHER_HALF(1);
  u = TO_HIGHER_HALF(t);

  if((a[0] += u) < u)
    a[1]++;
  a[1] += HIGHER_HALF(t);

}
  
/*
 * computes a = b + c*d 
 * return carry
 * a, b, d, arrays. c digit.
 */
static bigint_digit bigint_add_digit_mult(bigint_digit *a, bigint_digit *b, bigint_digit c, bigint_digit *d, int digits){

  bigint_digit carry, t[2];
  int i;
  
  /* check if multiplier is zero */
  if (c==0)
    return (0);

  carry = 0;
  for(i=0; i < digits; i++){
    bigint_digit_mult(t, c, d[i]);
    if((a[i] = b[i]+ carry) < carry)
      carry = 1;
    else
      carry = 0;
    if((a[i] += t[0]) < t[0])
      carry++;
    carry += t[1];
  }
  return carry;
}

/*
 * computes a = b - c*d for 
 * a[] b[] d[] and c is a digit
 * return borrow
 */
static bigint_digit bigint_sub_digit_mult(bigint_digit *a, bigint_digit *b, 
					  bigint_digit c, bigint_digit *d, 
					  int digits){

  bigint_digit borrow, t[2];
  int i;
  
  if(c ==0)
    return 0;

  borrow = 0;
  for(i = 0; i < digits; i++){
    bigint_digit_mult(t, c, d[i]);
    if((a[i] = b[i] - borrow) > (MAX_BIGINT_DIGIT - borrow))
      borrow = 1;
    else
      borrow = 0;
    if((a[i] -= t[0]) > (MAX_BIGINT_DIGIT - t[0]))
      borrow++;
    borrow += t[1];
  }

  return borrow;
}

/* 
 * returns how many bits are significant 
 */
static int bigint_digit_bits(bigint_digit a){
  int i;
  for(i=0; i< BIGINT_DIGIT_BITS; i++, a >>=1)
    if(a == 0)
      break;
  
  return i;
}

/*
 * a = b/ c where all are digits 
 * length of b is b[2]
 * b[1] < c and HIGHER_HALF(c) > 0
 */

void bigint_digit_div(bigint_digit *a, bigint_digit b[2], bigint_digit c){
  
  bigint_digit t[2], u, v;
  bigint_half_digit aHigh, aLow, cHigh, cLow;

  cHigh = (bigint_half_digit) HIGHER_HALF(c);
  cLow = (bigint_half_digit) LOWER_HALF(c);
  
  t[0] = b[0];
  t[1] = b[1];

  /* under estimate higher half of quotient, subtract */
  if(cHigh == LOW_MASK)
    aHigh = (bigint_half_digit) HIGHER_HALF (t[1]);

  else
    aHigh = (bigint_half_digit) (t[1]/(cHigh + 1));

  u = (bigint_digit) aHigh * (bigint_digit) cLow;
  v = (bigint_digit) aHigh * (bigint_digit) cHigh;

  if((t[0] -= TO_HIGHER_HALF (u)) > (MAX_BIGINT_DIGIT - TO_HIGHER_HALF (u)))
    t[1]--;
  t[1] -= HIGHER_HALF (u);
  t[1] -= v;

  /* Correct estimate.
   */
  while ((t[1] > cHigh) ||
         ((t[1] == cHigh) && (t[0] >= TO_HIGHER_HALF (cLow)))) {
    if ((t[0] -= TO_HIGHER_HALF (cLow)) > MAX_BIGINT_DIGIT - TO_HIGHER_HALF (cLow))
      t[1]--;
    t[1] -= cHigh;
    aHigh++;
  }

  /* Underestimate low half of quotient and subtract.
   */
  if (cHigh == LOW_MASK)
    aLow = (bigint_half_digit)LOWER_HALF (t[1]);
  else
    aLow = 
      (bigint_half_digit)((TO_HIGHER_HALF (t[1]) + HIGHER_HALF (t[0])) / (cHigh + 1));
  u = (bigint_digit)aLow * (bigint_digit)cLow;
  v = (bigint_digit)aLow * (bigint_digit)cHigh;
  if ((t[0] -= u) > (MAX_BIGINT_DIGIT - u))
    t[1]--;
  if ((t[0] -= TO_HIGHER_HALF (v)) > (MAX_BIGINT_DIGIT - TO_HIGHER_HALF (v)))
    t[1]--;
  t[1] -= HIGHER_HALF (v);

  /* Correct estimate.
   */
  while ((t[1] > 0) || ((t[1] == 0) && t[0] >= c)) {
    if ((t[0] -= c) > (MAX_BIGINT_DIGIT - c))
      t[1]--;
    aLow++;
  }
  
  *a = TO_HIGHER_HALF (aHigh) + aLow;
}

#ifdef OLD_IMP
/*
 * decodes character string b[len] to number a[digits]
 * assumes b has a string that will fit in a.
 */

void bigint_decode(bigint_digit *a, int digits, u8 *b, int len){
  bigint_digit t;
  int i, j, u;
  
  for(i =0, j=len-1; i < digits && j >=0; i++){
    t = 0;
    for(u=0; j >= 0 && u < BIGINT_DIGIT_BITS; j--, u+=8)
      t |= ((bigint_digit)b[j] ) << u;
    a[i] = t;
  }
  for(; i < digits; i++)
    a[i] = 0;
}
#endif

/*
 * encodes character string a[len] to number b[digits]
 */

 
void bigint_encode (u8 *a, int len, bigint_digit *b, int digits)
{
  bigint_digit t;
  int j;
  unsigned int i, u;

  for (i = 0, j = len - 1; i < digits && j >= 0; i++) {
    t = b[i];
    for (u = 0; j >= 0 && u < BIGINT_DIGIT_BITS; j--, u += 8)
      a[j] = (u8)(t >> u);
  }

  for (; j >= 0; j--)
    a[j] = 0;
}

/* 
 * copies b to a 
 */
void bigint_copy(bigint_digit  *a, bigint_digit *b, int digits){
  int i;
  for(i=0; i< digits; i++){
    a[i] = b[i];
  }
}

/* 
 * zeros 
 */
void bigint_zero(bigint_digit *a, int digits){
  int i;
  for(i=0; i< digits; i++){
    a[i] = 0;
  }
}

#ifdef OLD_IMP
/*
 * raises to two power
 * a = 2^b
 * b < digits * BIGINT_DIGIT_BITS
 */

void bigint_two_exp(bigint_digit *a, int b, int digits){
  bigint_zero(a, digits);
  if(b >= digits * BIGINT_DIGIT_BITS)
    return;
  a[b/BIGINT_DIGIT_BITS] = (bigint_digit) 1<<(b % BIGINT_DIGIT_BITS);
}
#endif
/* 
 * a = b+c
 * a,b,c all bigint_digit  arrays
 */
bigint_digit bigint_add(bigint_digit *a, bigint_digit *b, bigint_digit *c, int digits){
  bigint_digit ai, carry;
  int i;
  carry = 0;
  
  for(i = 0; i < digits; i++){
    if ((ai = b[i] + carry) < carry)
      ai = c[i];
    else if ((ai += c[i]) < c[i])
      carry = 1;
    else
      carry = 0;
    a[i] = ai;
  }

  return (carry);
}

/*
 * a = b - c all arrays, returns borrow
 */
bigint_digit bigint_sub(bigint_digit *a, bigint_digit *b, bigint_digit *c, int digits){
  bigint_digit ai, borrow;
  int i;
  borrow = 0;
  
  for (i = 0; i < digits; i++) {
    if ((ai = b[i] - borrow) > (MAX_BIGINT_DIGIT - borrow))
      ai = MAX_BIGINT_DIGIT - c[i];
    else if ((ai -= c[i]) > (MAX_BIGINT_DIGIT - c[i]))
      borrow = 1;
    else
      borrow = 0;
    a[i] = ai;
  }

  return (borrow);
}


/*
 * return length in digits
 */
int bigint_digits(bigint_digit *a, int digits){
  int i;
  for(i= digits - 1; i >= 0; i--)
    if (a[i])
      break;

  return (i + 1);
}

/* Computes a = b * c.
   Lengths: a[2*digits], b[digits], c[digits].
 */
void bigint_mult (bigint_digit *a, bigint_digit *b, bigint_digit *c, 
		  int digits)
{
  bigint_digit  t[2*MAX_BIGINT_DIGITS];
  int bDigits, cDigits, i;

  bigint_zero (t, 2 * digits);
  
  bDigits = bigint_digits (b, digits);
  cDigits = bigint_digits (c, digits);

  for (i = 0; i < bDigits; i++){
    t[i+cDigits] += bigint_add_digit_mult(&t[i], &t[i], b[i], c, cDigits);
  }
  bigint_copy(a, t, 2 * digits);
}

/*
 * a = b << c (shifts b left c bits)
 *return carry 
 */

bigint_digit  bigint_left_shift (bigint_digit *a, bigint_digit  *b, 
				 int c, int digits){
  bigint_digit bi, carry;
  int i, t;
  
  if (c >= BIGINT_DIGIT_BITS)
    return (0);
  
  t = BIGINT_DIGIT_BITS - c;
  carry = 0;
  for (i = 0; i < digits; i++) {
    bi = b[i];
    a[i] = (bi << c) | carry;
    carry = c ? (bi >> t) : 0;
  }
  return (carry);
}

/*
 * shifts b right c times and returns in a 
 */
bigint_digit bigint_right_shift(bigint_digit *a, bigint_digit *b, 
				int c, int digits){
  bigint_digit  bi, carry;
  int i;
  unsigned int t;
  
  if (c >= BIGINT_DIGIT_BITS)
    return (0);
  
  t = BIGINT_DIGIT_BITS - c;
  carry = 0;

  for (i = digits - 1; i >= 0; i--) {
    bi = b[i];
    a[i] = (bi >> c) | carry;
    carry = c ? (bi << t) : 0;
  }
  
  return (carry);
}


/*
 * compare and return sign of a-b
 */
int bigint_cmp(bigint_digit *a, bigint_digit *b, int digits){
  int i;
  
  for (i = digits - 1; i >= 0; i--) {
    if (a[i] > b[i])
      return (1);
    if (a[i] < b[i])
      return (-1);
  }

  return (0);
}


/*
 * a = c/d and b = c % d
 * a[cDigits], b[dDigits], c[cDigits], d[dDigits], 
 * cDigits < 2*MAX_BIGINT_DIGITS
 */
void bigint_div (bigint_digit *a, bigint_digit *b, bigint_digit *c, 
		 int cDigits, bigint_digit  *d, int dDigits){
  bigint_digit  ai, cc[2*MAX_BIGINT_DIGITS+1], dd[MAX_BIGINT_DIGITS], t;
  int i, ddDigits, shift;
  
  ddDigits = bigint_digits (d, dDigits);
  if (ddDigits == 0)
    return;
  
  /* Normalize operands.
   */
  shift = BIGINT_DIGIT_BITS - bigint_digit_bits (d[ddDigits-1]);
  bigint_zero (cc, ddDigits);
  cc[cDigits] = bigint_left_shift(cc, c, shift, cDigits);
  bigint_left_shift (dd, d, shift, ddDigits);
  t = dd[ddDigits-1];
  
  bigint_zero(a, cDigits);

  for (i = cDigits-ddDigits; i >= 0; i--) {
    /* Underestimate quotient digit and subtract.
     */
    if (t == MAX_BIGINT_DIGIT)
      ai = cc[i+ddDigits];
    else
      bigint_digit_div(&ai, &cc[i+ddDigits-1], t + 1);
    cc[i+ddDigits] -= bigint_sub_digit_mult(&cc[i], &cc[i], ai, dd, ddDigits);

    /* Correct estimate.
     */
    while (cc[i+ddDigits] || (bigint_cmp(&cc[i], dd, ddDigits) >= 0)) {
      ai++;
      cc[i+ddDigits] -= bigint_sub(&cc[i], &cc[i], dd, ddDigits);
    }
    
    a[i] = ai;
  }
  
  bigint_zero(b, dDigits);
  bigint_right_shift(b, cc, shift, ddDigits);
}

/*
 * compute a = b mod c 
 * a[cDigits], b[bDigits], c[cDigits]
 * c > 0, bDigits < 2*MAX_BIGINT_DIGITS, cDigits < MAX_BIGINT_DIGITS
 */
void bigint_mod(bigint_digit *a, bigint_digit *b, int bDigits, bigint_digit *c, int cDigits){
  bigint_digit  t[2 * MAX_BIGINT_DIGITS];

  bigint_div (t, a, b, bDigits, c, cDigits);
}

/*
 * a = b*c mod d 
 */
void bigint_mod_mult(bigint_digit *a,  bigint_digit *b, bigint_digit *c, bigint_digit *d, int digits)
{
  bigint_digit  t[2*MAX_BIGINT_DIGITS];

  bigint_mult(t, b, c, digits);
  bigint_mod(a, t, 2 * digits, d, digits);
}


/* a = b^c mod d.
 * Lengths: a[dDigits], b[dDigits], c[cDigits], d[dDigits].
 * Assumes d > 0, cDigits > 0, dDigits < MAX_BIGINT_DIGITS.
 */
void bigint_mod_exp (bigint_digit *a, bigint_digit *b, bigint_digit *c, int cDigits, bigint_digit *d, int dDigits)
{
  bigint_digit  bPower[3][MAX_BIGINT_DIGITS], ci, t[MAX_BIGINT_DIGITS];
  int i;
  unsigned int ciBits, j, s;
  bigint_digit exp;
  unsigned int need[4];
  bigint_digit setbits = 0;
  
  /* zero and power one come for free */

  need[0] = 1;
  need[1] = 1;
  need[2] = 0;
  need[3] = 0;

  /* check input to see if power 3 is needed, then turn on 2.*/
  exp = c[0];
  
  for(i = 0; i< 16; i++){
      /* take last two bits */
      setbits = exp & 0x00000003;
      /* you need to compute that power */
      need[setbits]++;
      /* shift right two bits */
      exp = exp >> 2;
  }
  if(need[3] >0){
      need[2] = 1; /* need to compute anyway */
  }
    
  /* Store b, b^2 mod d, and b^3 mod d.
   */
  
  bigint_copy(bPower[0], b, dDigits);
  if(need[2] > 0){
      bigint_mod_mult(bPower[1], bPower[0], b, d, dDigits);
  }
  if(need[3] > 0){
      bigint_mod_mult(bPower[2], bPower[1], b, d, dDigits);
  }

  
  BIGINT_ASSIGN_DIGIT (t, 1, dDigits);

  cDigits = bigint_digits (c, cDigits);
  for (i = cDigits - 1; i >= 0; i--) {
    ci = c[i];
    ciBits = BIGINT_DIGIT_BITS;
    
    /* Scan past leading zero bits of most significant digit.
     */
    if (i == (int)(cDigits - 1)) {
      while (! BIGINT_DIGIT_2MSBS (ci)) {
        ci <<= 2;
        ciBits -= 2;
      }
    }

    for (j = 0; j < ciBits; j += 2, ci <<= 2) {
      /* Compute t = t^4 * b^s mod d, where s = two MSB's of ci.
       */
      bigint_mod_mult(t, t, t, d, dDigits);
      bigint_mod_mult(t, t, t, d, dDigits);
      if ((s = BIGINT_DIGIT_2MSBS (ci)) != 0)
        bigint_mod_mult(t, t, bPower[s-1], d, dDigits);
    }
  }
  
  bigint_copy(a, t, dDigits);
}

/* 
 * return nonzero iff a is zero
 */
int bigint_iszero(bigint_digit *a, int digits){
  int i;
    
  for (i = 0; i < digits; i++)
    if (a[i])
      return (0);
    
  return (1);
}

/*
 * return length of a in bits 
 */
int bigint_bits(bigint_digit *a, int digits){
  if ((digits = bigint_digits (a, digits)) == 0)
    return (0);

  return ((digits - 1) * BIGINT_DIGIT_BITS + bigint_digit_bits (a[digits-1]));
}



void
bsl_rsa_verify(char *result, unsigned long *certsign, unsigned long *certpublickey, unsigned long *certexponent,  int num_bits){
 /* grab from cert */
    unsigned long eDigits, nDigits;
    bigint_digit bign[MAX_BIGINT_DIGITS];
    bigint_digit bigm[MAX_BIGINT_DIGITS];
    bigint_digit bige[MAX_BIGINT_DIGITS];
    bigint_digit  bigc[MAX_BIGINT_DIGITS];
    int outlen;
    int i;
    int num_words = num_bits/BIGINT_DIGIT_BITS;

    for(i=0; i< MAX_BIGINT_DIGITS; i++){
        bign[i] = 0;
        bigm[i] = 0;
        bige[i] = 0;
        bigc[i] = 0;
    }
    
    for(i =0; i < num_words; i++){
        bign[num_words -1 -i] = certpublickey[i];
    }
    for(i =0; i< num_words; i++){
        bigm[num_words -1 -i] = certsign[i];
    }
    bige[0] = *certexponent;
    nDigits = bigint_digits(bign, MAX_BIGINT_DIGITS);
    eDigits = 1;
    
                
    bigint_mod_exp(bigc, bigm, bige, eDigits, bign, nDigits);
    
    outlen = (num_bits + 7)/8;
    bigint_encode(result, outlen, bigc, nDigits);
    
}



void 
field_to_bigint(field_2n *a, bigint_digit *b, int digits){
  short int i;

  for(i= 0; i < digits; i++){
    b[digits-1-i] = a->e[i];    
  }
}


void bigint_to_field(bigint_digit *a, field_2n *b, int digits){
  short int i;
  null(b);
  for(i= 0; i < digits; i++){
    b->e[i] = a[digits-1-i];
  }
}


/* compute a = 1/b mod c
 * this is only for ecc, so use the smaller MAX_ECC_DIGITS
 * gcd(b,c)=1
 */
void bigint_mod_inv (bigint_digit *a, bigint_digit *b, 
		     bigint_digit *c, int digits){
  bigint_digit q[MAX_ECC_DIGITS], t1[MAX_ECC_DIGITS], t3[MAX_ECC_DIGITS],
    u1[MAX_ECC_DIGITS], u3[MAX_ECC_DIGITS], v1[MAX_ECC_DIGITS],
    v3[MAX_ECC_DIGITS], w[2*MAX_ECC_DIGITS];
  int u1Sign;

  /* Apply extended Euclidean algorithm, modified to avoid negative
     numbers.
   */
  BIGINT_ASSIGN_DIGIT (u1, 1, digits);
  bigint_zero(v1, digits);
  bigint_copy(u3, b, digits);
  bigint_copy(v3, c, digits);
  u1Sign = 1;

  while (! bigint_iszero(v3, digits)) {
    bigint_div(q, t3, u3, digits, v3, digits);
    bigint_mult(w, q, v1, digits);
    bigint_add(t1, u1, w, digits);
    bigint_copy(u1, v1, digits);
    bigint_copy(v1, t1, digits);
    bigint_copy(u3, v3, digits);
    bigint_copy (v3, t3, digits);
    u1Sign = -u1Sign;
  }
  
  /* Negate result if sign is negative.
    */
  if (u1Sign < 0)
    bigint_sub(a, c, u1, digits);
  else
    bigint_copy (a, u1, digits);
}