#include <stdio.h>
#include <ctype.h>
#include <assert.h>
#include <string.h>

#define isEven(x) (0==(x & 0x01))
#define isOdd(x)  (0!=(x & 0x01))
#define swap(x,y) ((x^=y),(y^=x),(x^=y))

#ifndef TRUE
#define TRUE (0==0)
#define FALSE (!TRUE)
#endif

int FBN_shift(unsigned char *, int, char *);
int die(char *);  /* print message and quit */


/*
 *  Fixed Big Number library
 *    http://www.notatla.demon.co.uk/CRYPTO/FBN_funcs.c
 *
 *  This is integer arithmetic code for up to 1024 bits
 *  (with intermediate results up to 2048 bits). 
 *  This is what I mean by 'Fixed Big Number' - not arbitrary sizes.
 *  It has been tested against 'bc' for some calculations of each
 *  type, but not right up to the size limits.   This code is simple
 *  and not especially fast.
 *
 *  FBN_powmod() is effectively RSA.  This code contains no
 *  primality testing which you would need to generate RSA keys.
 *
 *  Parts of this code are derived from the book "Applied Cryptography"
 *  2nd ed by Bruce Schneier.  
 *
 *  Distribution and use is free; no GPL, Berkeley or other licences
 *  apply.   The RSA patent in the U.S is due to expire September 2000.
 *  More details at http://cwis.kub.nl/~frw/people/koops/lawsurvy.htm.
 *
 *  Enail any bug reports to me please.
 *       Antonomasia <ant notatla.demon.co.uk>   13Feb1999
 */

int FBN_byte_sbc(unsigned char *alter, unsigned char deduct, int scale)
{
    int i;

    i = 255 - scale;
    if ((i > 256) || (i < 0)) {
	die("bad scale for FBN_byte_sbc");
    }
    for (; (i >= 0) && (0 != deduct); i--) {
	if (alter[i] >= deduct) {
	    alter[i] -= deduct;
	    deduct = 0;
	} else {
	    alter[i] = alter[i] + 256 - deduct;
	    deduct = 1;
	}
    }

    if (0 != deduct) {
	die("subtraction failed (outstanding borrow at end)");
    }
    return 0;
}
/*
 * ############################################ 
 */
int FBN_short_mult(unsigned char *left, int right)
{
    int i, tmp;
    unsigned char carry = 0;

    for (i = 255; i > 0; i--) {
	tmp = carry + (left[i] * right);
	left[i] = tmp % 256;
	tmp /= 256;
	if (tmp > 255) {
	    die("improper carry detected");
	}
	carry = tmp;
    }

    if (0 != tmp) {
	die("carry attempted past end of number");
    }
    return 0;
}

/*
 * ############################################ 
 */
int FBN_add(unsigned char *left, unsigned char *right)
{
    int i, tmp;
    unsigned char carry = 0;

    for (i = 255; i > 0; i--) {
	tmp = carry + left[i] + right[i];
	left[i] = tmp % 256;
	tmp /= 256;
	if (tmp > 255) {
	    die("improper carry detected");
	}
	carry = tmp;
    }

    if (0 != tmp) {
	die("carry attempted past end of number");
    }
    return 0;
}
/*
 * ############################################ 
 */
int FBN_cmp(unsigned char *left, unsigned char *right)
{
    int i;

    for (i = 0; i < 256; i++) {
	if (left[i] < right[i])
	    return -1;
	if (left[i] > right[i])
	    return +1;
    }
    return 0;
}
/*
 * ############################################ 
 */
int FBN_iszero(unsigned char *polo)
{
/*
 * is this number round, with a hole ? 
 */
    int i;
    for (i = 0; i < 256; i++)
	if (0 != *(polo + i))
	    return FALSE;
    return TRUE;
}

/*
 * ############################################ 
 */
int FBN_isone(unsigned char *stick)
{
/*
 * is this number long and vertical ?
 */
    int i;
    for (i = 0; i < 255; i++) {
	if (0 != *(stick + i))
	    return FALSE;
    }
    if (1 != stick[255])
	return FALSE;
    return TRUE;
}

/*
 * ############################################ 
 */

void FBN_swap(unsigned char *x, unsigned char *y)
{
    int i;
    unsigned char tmp;

    for (i = 0; i < 256; i++) {
	tmp = x[i];
	x[i] = y[i];
	y[i] = tmp;
    }
}
/*
 * ############################################ 
 */

int FBN_readfile(unsigned char *v, char *fname)
{
    FILE *fp;
    char decnum[400];
    unsigned char scale[256];
    int i, j;

    fp = fopen(fname, "r");
    if (NULL == fp) {
	perror("fopen in FBN_readfile:");
	die("could not open input file");
    }
    memset(decnum, 0, sizeof(decnum));
    memset(scale, 0, sizeof(scale));
    scale[255] = 1;		/*
				 * one in LSB 
				 */

    fgets(decnum, sizeof(decnum), fp);
    fclose(fp);
    i = strlen(decnum);

/*
 * string can end with NULL and newline 
 */
    for (; (decnum[i] < '0') || (decnum[i] > '9'); i--);

    for (; i >= 0; i--) {
	for (j = 0; j < (decnum[i] - '0'); j++) {
	    FBN_add(v, scale);
	}
	FBN_short_mult(scale, 10);
    }

    return 0;
}
/*
 * ############################################ 
 */

int FBN_print(unsigned char *v)
{
    int i, killzero = 1;

    printf("0x");
    for (i = 0; i < 256; i++) {
	if (0 != v[i])
	    killzero = 0;
	if ((0 == killzero) || (255 == i))
	    printf("%02x", v[i]);
    }
    printf("\n");

    return 0;
}

/*
 * ############################################ 
 */

int FBN_sbc(unsigned char *decreasing, unsigned char *reduction)
{
    int j;
    for (j = 0; j < 256; j++) {
	FBN_byte_sbc(decreasing, reduction[255 - j], j);
    }
    return 0;
}

/*
 * ############################################ 
 */
void FBNExtBinEuclid(unsigned char *u,
		     unsigned char *v,
		     unsigned char *u1,
		     unsigned char *u2,
		     unsigned char *u3)
{
/*
 * based on Schneier's code in AC2, page 247 
 */
    /*
     * warning; u and v will be swapped if u < v 
     */
    int k;
    unsigned char t1[256], t2[256], t3[256];
    memset(t1, 0, 256);
    memset(t2, 0, 256);
    memset(t3, 0, 256);

    if (-1 == FBN_cmp(u, v))
	FBN_swap(u, v);

    for (k = 0; isEven(u[255]) && isEven(v[255]); ++k) {
	FBN_shift(u, 1, "right");
	FBN_shift(v, 1, "right");
    }

    memset(u1, 0, 256);
    u1[255] = 1;
    memset(u2, 0, 256);
    memcpy(u3, u, 256);
    memcpy(t1, v, 256);
    memcpy(t2, u, 256);
    FBN_byte_sbc(t2, 1, 0);
    memcpy(t3, v, 256);

    do {
	do {
	    if (isEven(u3[255])) {
		if (isOdd(u1[255]) || isOdd(u2[255])) {
		    FBN_add(u1, v);
		    FBN_add(u2, u);
		}
		FBN_shift(u1, 1, "right");
		FBN_shift(u2, 1, "right");
		FBN_shift(u3, 1, "right");
	    }
	    if (isEven(t3[255]) ||
		(-1 == FBN_cmp(u3, t3))) {
		FBN_swap(u1, t1);
		FBN_swap(u2, t2);
		FBN_swap(u3, t3);
	    }
	} while (isEven(u3[255]));
	while ((-1 == FBN_cmp(u1, t1)) ||
	       (-1 == FBN_cmp(u2, t2))) {
	    FBN_add(u1, v);
	    FBN_add(u2, u);
	}
	FBN_sbc(u1, t1);
	FBN_sbc(u2, t2);
	FBN_sbc(u3, t3);
    } while (!FBN_iszero(t3));

    while ((FBN_cmp(u1, v) >= 0) &&
	   (FBN_cmp(u2, u) >= 0)) {
	FBN_sbc(u1, v);
	FBN_sbc(u2, u);
    }

    FBN_shift(u, k, "left");
    FBN_shift(v, k, "left");
    FBN_shift(u3, k, "left");
}

/*
 * ############################################ 
 */
int FBN_bitshift_left(unsigned char *x, int by)
{
/*
 * increases value of the number 
 */
    int i, tmp;
    unsigned char carry;

    if ((by < 1) | (by > 7)) {
	fprintf(stderr, "In FBN_bitshift_left by=%d - quitting\n", by);
	die("Value of by is outside valid range (1-7)");
    }
    carry = 0;
    for (i = 255; i >= 0; i--) {

	tmp = (x[i] << by);
	x[i] = tmp % 256;
	tmp /= 256;
	x[i] += carry;
	carry = tmp;
    }

    return 0;
}
/*
 * ############################################ 
 */
int FBN_bitshift_right(unsigned char *x, int by)
{
/*
 * decreases value of the number 
 */
    int i, tmp;
    unsigned char carry;

    if ((by < 1) | (by > 7)) {
	fprintf(stderr, "In FBN_bitshift_right by=%d - quitting\n", by);
	die("Value of by is outside valid range (1-7)");
    }
    carry = 0;
    for (i = 0; i < 256; i++) {

	tmp = x[i];
	x[i] >>= by;
	tmp = tmp - (x[i] << by);
	x[i] += carry;
	carry = tmp << (8 - by);
    }

    return 0;
}
/*
 * ############################################ 
 */
int FBN_byteshift(unsigned char *x, int by, char *dirn)
{

    int dog, tail = 0;

    if (0 != strcmp("left", dirn) && 0 != strcmp("right", dirn)) {
	die("invalid direction in FBN_byteshift");
    }
    tail = by;

    if ((tail < 1) || (tail > 255)) {
	fprintf(stderr, "Illegal byte shift of %d (line=%d)\n",
		by, __LINE__);
	die("error involving bignumber shifts");
    }
    for (dog = 0; tail < 256; dog++, tail++) {
	if (0 == strcmp("left", dirn)) {
	    x[dog] = x[tail];
	    x[tail] = 0;
	} else {
	    if ((256 - tail) > by)
		x[256 - dog] = x[256 - tail];
	    else
		x[256 - tail] = 0;
	}
    }

    return 0;
}
/*
 * ############################################ 
 */
int FBN_shift(unsigned char *x, int by, char *dirn)
{
    if (by >= 8)
	FBN_byteshift(x, by / 8, dirn);
    if (0 == by % 8)
	return 0;

    if (0 == strcmp("right", dirn))
	return FBN_bitshift_right(x, by % 8);
    if (0 == strcmp("left", dirn))
	return FBN_bitshift_left(x, by % 8);

    fprintf(stderr, "DIRN=%s\n", dirn);
    return die("FBN_shift direction not left or right");
}
/*
 * ############################################ 
 */
int FBN_get_inverse(unsigned char *u,
		    unsigned char *v,
		    unsigned char *inverse)
{
    int i;
    unsigned char a[256], b[256], gcd[256];

    FBNExtBinEuclid(u, v, a, b, gcd);
    memset(inverse, 0, 256);

/*
 * check gcd is 1 
 */
    if (!FBN_isone(gcd)) return -1;

    memcpy(inverse, u, 256);
    FBN_sbc(inverse, b);
    return 0;
}
/*
 * ########################################## 
 */
int FBN_long_mult(unsigned char *left, unsigned char *right)
{
    int i;
    unsigned char running_total[256], this_pass[256];

    memset(running_total, 0, sizeof(running_total));

    for (i = 255; i >= 0; i--) {
	if (0 == right[i])
	    continue;

	memcpy(this_pass, left, sizeof(this_pass));
	FBN_short_mult(this_pass, right[i]);

	if (255 != i)
	    FBN_byteshift(this_pass, 255 - i, "left");
	FBN_add(running_total, this_pass);
    }

    memcpy(left, running_total, sizeof(running_total));
    return 0;
}
/*
 * ########################################## 
 */
int FBN_modreduce(unsigned char *vary, unsigned char *mod)
{
    int i = 0;

    if (FBN_iszero(mod))
	die("mod reductions with mod=0 are daft");
    if (FBN_isone(mod))
	die("mod reductions with mod=1 are daft");

    if (-1 == FBN_cmp(vary, mod))
	return 0;

/*
 * scale mod until .ge. vary 
 */
    for (i = 0; 1 == FBN_cmp(vary, mod); i++) {
	FBN_bitshift_left(mod, 1);
    }
    if (0 == FBN_cmp(vary, mod)) {
	memset(vary, 0, 256);
	return 0;
    }
/*
 * mod .gt. vary for the first time 
 */
    for (; i > 0; i--) {
	FBN_bitshift_right(mod, 1);
	if (1 == FBN_cmp(vary, mod))
	    FBN_sbc(vary, mod);
    }

    return 0;
}
/*
 * ########################################## 
 */
int FBN_powmod(unsigned char *base,
	       unsigned char *exp,
	       unsigned char *modulus,
	       unsigned char *result)
{

/*
 * based on Schneier's code in AC2, page 244 
 */
    unsigned char s[256], t[256], u[256];

    memset(s, 0, 255);
    s[255] = 1;
    memcpy(t, base, 256);
    memcpy(u, exp, 256);

    while (!FBN_iszero(u)) {
	if (isOdd(u[255])) {
	    FBN_long_mult(s, t);
	    if (NULL != modulus)
	        if (!FBN_iszero(modulus))
		    FBN_modreduce(s, modulus);
	}
	FBN_bitshift_right(u, 1);
	FBN_long_mult(t, t);
	if (!FBN_iszero(modulus))
	    FBN_modreduce(t, modulus);
    }

    memcpy(result, s, 256);
    return 0;
}