//
// rsa.c demonstrate rsa encryption/decryption/cracking
//    by W.J.Hengeveld  aug1994
//      itsme@xs4all.nl
//
//    (made to compile under linux feb2001)
//
//
//  publickey = (a, p*q)   secretkey = (x)
//      a*x+b*fi(p*q)=1
//  encrypt :  Cypher = Message^a (mod p*q)
//  decrypt :  Cypher^x (mod p*q)= Message^(a*x) (mod p*q)=
//               Message^(1-b*fi(p*q)) =
//               Message * (Message^fi(p*q))^-b =  Message
//
//  genkey : choose { p, q, a} such that   a*x+b*fi=1 (mod p*q)
//           -> (x)=secret key,  (a, p*q)=public key
//
//  breakkey: compute {x} such that  a*x+b*fi=1 (mod p*q)
//
//
//  { fermat : a^fi(n) (mod n)= 1 }
//
//  rsa [-nMIN] [-mMAX] -g [p] [q] [a]
//     generate a key of MIN-MAX bits [take p, q and a when available
//
//  rsa -b pq a
//     try to break public key  [a], [pq]
//
//  rsa -e pq a msg
//
/* example:
 * generate a key with p and q below 8 bits

$ ./rsa -m8 -g
Generating primes
l(p)=0 l(q)=7
l(p)=7 l(q)=8
p=127  q=251  fi=31500 ... 
choose <a> such that it has no factors of fi:
2^2 * 3^2 * 5^3 * 7
a = 11
(privkey)x=8591, y=-3, (pubkey)pq=31877

 * encrypt secret message '999' (with public key [31877,11])

$ ./rsa -e 31877 11 999
crypted: 15722

 * 'crack' rsa:

$ ./rsa -k 31877 11
a=11  pq=31877 ... fi=31500 ... (privkey)x=8591 y=-3 g=1

 * decrypt encrypted message:

$ ./rsa -e 31877 8591 15722
crypted: 999

*/

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <stdarg.h>
#include <malloc.h>

#if defined(_MSC_)
#include <memory.h>
#elif defined(__TURBOC__)
#include <mem.h>
#else
#include <string.h>
#endif

// generate <x> such that   a*x+fi(pq)*y=1,
//long generatekey(int n, long p, long q)
//{
//
//}

// efficiency of using this method of stepping through primes:
//   l=1   :  100%
//   l=2   :   50%
//   l=6   :   33%
//   l=30  :   27%
//   l=210 :   20%     2*3*5*7
//   l=2310:   15%     2*3*5*7*11
//   ....
//
// 2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67
//  +1 +2 +2 +4 +2 +4 +2 +4 +6 +2 +6 +4 +2 +4 +6 +6 +2 +6
//  start = primes upto p
//  total span of step =  product of primes upto (excluding) p
long start[]={2, 3, 5, 7};
// long start2[]={4, 9, 25, 49};
long step[]={ 4 ,2 ,4 ,2 ,4 ,6 ,2 ,6 };
// int step2[]={ 16 ,4 ,16 ,4 ,16 ,36 ,4 ,36 };

#define STEPLEN sizeof(step)/sizeof(long)
#define STARTLEN sizeof(start)/sizeof(long)

#define MAXFACTORS 64
#define GENSIEVE   16384L
#define MAXTRIES   100

// storage used when factoring
long factors[MAXFACTORS];
int powers[MAXFACTORS];

char *sieve=NULL;

// returns # 1 bits in number
int bitcount(long n)
{
	int count=0;
	while (n)
	{
		count++;
		n&=n-1;
	}
	return count;
}

int bitlen(long n)
{
	int count=0;
	if (n<0) return 32;
	while (n)
	{
		count++;
		n>>=1;
	}
	return count;
}

void print1factor(long f, int k)
{
	printf("%ld", f);
	if (k>1)
		printf("^%d", k);
}

void printfact(long *f, int *k, int n)
{
	int i;
	print1factor(f[0],k[0]);
	for (i=1 ; i<n ; i++)
	{
		printf(" * ");
		print1factor(f[i],k[i]);
	}
	putchar('\n');
}

//  (p+s)^2 = p^2 + 2*p*s + s^2
//  ->  p2:=p2+ 2*step * p + step2
//  ->  p:=p+step
// factor n into  prod(p[i]^k[i], i=0..nfact-1)
//
int factor(long n, long *primefactor, int *powers)
{
	long p;     // prime currently testing
	int  nf=0;  // # of factor currently working on
	int  ex;    // exponent of last prime division
	int  stp;   // steps through start & step array
	int  in_start_up=1;
	stp=0;
	p=start[stp++];
	while (p*p<=n)
	{
		ex=0;
		while (n%p==0)  // determine if it is a factor
		{
			n/=p;
			ex++;
		}
		if (ex)      // if it is a factor store it
		{
			if (nf==MAXFACTORS)
			{
				printf("Too many factors, left over: ");
				print1factor(p, ex);
				if (n>1)
					print1factor(n,1);
				putchar('\n');

				return nf;
			}
			primefactor[nf]=p;
			powers[nf]=ex;
			nf++;
		}

		if (in_start_up)
		{
			p=start[stp++];      // loop through startup primes
			if (stp==STARTLEN)
			{
				in_start_up=0;
				stp=0;
			}
		}
		else
		{
			p+=step[stp++];      // step through numbers using step array
			if (stp==STEPLEN) stp=0;
		}
	}

	if (n!=1)
	{
		if (nf==MAXFACTORS)
		{
			printf("Too many factors, left over: "); print1factor(n, 1); putchar('\n');
			return nf;
		}
		primefactor[nf]=n;
		powers[nf]=1;
		nf++;
	}
	return nf;
}
/*
//-----------------------------------------------------------------------
//          
// for a number n, try to find numbers x and y such that
//  eq1:  n*k=x^2-y^2

// start by finding numbers a and b of form a^2-n

// eq2:  (a^2-n)*(b^2-n)=z^2  ->
// eq3:  (a*b)^2-z^2 = n*(a^2+b^2-n)
//     than (comparing to eq1):
//        x = a*b
//        y = z
//        k = a^2+b^2-n
//
// than a factor of n is the gcd(n, x+y) and gcd(n, x-y)

#define QS_DEPTH 50
int qfactor(long n, long *primefactor, int *powers)
{
	int qflist[QS_DEPTH];  // factors of a^2-n
	int qplist[QS_DEPTH];  // powers of factors of a^2-n
	int sqrt_n;

	sqrt_n=sqr
	for (i=sqrt(n)+1 ; i<...)
		... not finished yet ...

}
*/
// computes  n^k
long power(long n, int k)
{
	long z=n;
	long p=1;
	while(k)
	{
		if (k&1)  p*=z;
		z*=z;
		k>>=1;
	}
	return p;
}

// computes  m^a mod pq
long expmod(long m, long a, long pq)
{
	long z=m;
	long p=1;
	while(a)
	{
		if (a&1)  { p*=z; p%=pq;  }
		z*=z; z%=pq;
		a>>=1;
	}
	return p;
}

// computes eulers totient function: (# of elements in Gn - symmetric group n)
//   the number of numbers in [1..n] relatively prime to n
//   for a prime p :  fi(p^(n))=p^(n)-p^(n-1)
//                    fi(p)=p-1
//  fi(m) = m * product(1-1/p_i)
//     where m = product(p_i^n_i)
//
//  two numbers are relatively prime when ggd(a,b)=1
//
long euler(long n)
{
	int  nfact;
	long  fi=1; // end result
	int   i;

	nfact=factor(n, factors, powers);
//  printf("%ld = ", n); printfact(factors, powers, nfact);

	fi=n;
	for (i=0 ; i<nfact ; i++)
	{
		fi/=factors[i];
		fi*=factors[i]-1;
	}
	return fi;
}

// calculate x & y such that  a*x+b*y=ggd(a,b) { == 1 }
// returns # iterations
int chinese(long a, long b, long *x, long *y, long *g)
{
	long m[2][3];
	int n,i;
	long q;
	m[0][0]=1; m[0][1]=0; m[0][2]=a;
	m[1][0]=0; m[1][1]=1; m[1][2]=b;
	i=0; n=0;
	while (m[1-i][2])
	{
		q=m[i][2]/m[1-i][2];
		m[i][2]-=q*m[1-i][2];
		m[i][1]-=q*m[1-i][1];
		m[i][0]-=q*m[1-i][0];
		i=1-i;
		n++;
	}
	if (m[i][0]<0) { m[i][0]+=b; m[i][1]-=a; }
	*x=m[i][0]; *y=m[i][1];  *g=m[i][2];

	return n;
}

char ormask[]={0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80};
char andmask[]={0xfe,0xfd,0xfb,0xf7,0xef,0xdf,0xbf,0x7f};

int bittest(char *sieve, long bit)
{
	return sieve[(int)(bit/8)]&ormask[(int)(bit&7)];
}

void bitclear(char *sieve, long bit)
{
	sieve[(int)(bit/8)]&=andmask[(int)(bit&7)];
}

void erostatenes(char *sieve, long n)
{
	long i,j;
	memset(sieve, 0xaa, (int)(n/8+1));   // start with 2
//  1 0 1 0 1 1 0 0
//  7 6 5 4 3 2 1 0   (2 is prime,  1 & 0 are not)
	sieve[0]=0xac;
	i=3;
	while (i*i<=n)
	{
		if (bittest(sieve,i))
		{
			for (j=i*i ; j<=n ; j+=i)
				bitclear(sieve, j);
		}
		i++;
	}
}

long printsieve(char *sieve, long n)
{
	long i=0;
	long nprimes=0;

	while (i<=n)
	{
		if (bittest(sieve,i))
		{
			nprimes++;
			printf(" %ld", i);
		}
		i++;
	}
	return nprimes;
}

// generates a prime unequal to <p>
long generateprime(long p, int n, int m)
{
	long q;
	int tries=0;
	if (sieve==NULL)
	{
		sieve=(char*)malloc((int)GENSIEVE);
		if (sieve==NULL)
		{
			printf("Out of memory\n");
			exit(1);
		}
		printf("Generating primes\n");
		erostatenes(sieve, GENSIEVE*8);
		srand(time(NULL));
	}
	if (m>31-bitlen(p))
		m=31-bitlen(p);

	do {
		q=rand()%(1<<(m-3));
		while (q<GENSIEVE && sieve[(int)q]==0)
			q++;
		if (q==GENSIEVE) continue;

		q*=8;
		while (q<GENSIEVE*8 && bittest(sieve, q)==0)
			q++;
		if (q==GENSIEVE*8) continue;
		tries++;
	} while ((q==p || bitlen(q)<n || bitlen(q)>m)
				&& tries<MAXTRIES);
	if (tries==MAXTRIES)
		printf("FAILED to generate prime of %d-%d bits, unequal to %ld\n",
				n, m, p);
	printf("l(p)=%d l(q)=%d\n", bitlen(p), bitlen(q));
	return q;
}


void error(char *s, ...)
{
	va_list ap;
	va_start(ap,s);
	vprintf(s,ap);
	va_end(ap);
	printf("Usage: rsa [-f N] [-nMIN] [-mMAX] [-g [p [q [a]]]] [-k pq a] [-e pq a msg]\n");
	printf("   -f : factor number\n");
	printf("   -g : generate key of MIN-MAX bits(default 8-12 bits)\n");
	printf("   -k : factor key\n");
	printf("   -p : generate primes\n");
	printf("   -e : encrypt\n");
	if (sieve) free(sieve);
	exit(1);
}

long readnumber(void)
{
	char line[80];
	fgets(line, 80, stdin);
	return strtoul(line, 0, 0);
}

int main(int argc, char **argv)
{
	int n=4;          // minimum # bits in primes
	int m=12;         // maximum # bits in primes
	long p=0;
	long q=0;

	long cyp, msg;   // cyphertext / plaintext

	long a=0;        // modular inverse of fi
	long fi=0;       // euler phi of n
	long g;          // contains ggd of a & fi
	long x,y;        // result from chinese()
	long pq=0;       // used

	int  nf;         // # factors

	char *sieve=NULL;// used for prime number generation
	char *arg;       // commandline argument

	while (--argc)
	{
		arg=*++argv;
		if (*arg=='-')
		{
			switch(*++arg)
			{
				case 'n':  n=atoi(arg+1); break;  // minimum # bits
				case 'm':  m=atoi(arg+1); break;  // maximum # bits
				case 'g':  // compute secret key from p, q, a
					if (n>m) error("min > max bits!\n");
					if (argc>1) { argc--; p=atol(*++argv); }
					if (argc>1) { argc--; q=atol(*++argv); }
					if (p==0) p=generateprime(0, n, m);
					if (q==0) q=generateprime(p, n, m);
					fi=(p-1)*(q-1);
					printf("p=%ld  q=%ld  fi=%ld ... ", p, q, fi);

					if (argc>1) { argc--; a=atol(*++argv); }
					if (a==0)
					{
						printf("\nchoose <a> such that it has no factors of fi:\n");
						nf=factor(fi, factors, powers);
						printfact(factors, powers, nf);
						printf("a = "); a=readnumber();
					}
					chinese(a,fi, &x, &y, &g);
					printf("(privkey)x=%ld, y=%ld, (pubkey)pq=%ld\n", x, y, p*q);
					if (g!=1)
						printf("ggd( fi(pq), a ) = %ld !!!! invalid a\n", g);
					break;

				case 'k':  // compute secret key from pq, a
					if (argc>1) { argc--; pq=atol(*++argv); }
					else error("no publickey\n");
					if (argc>1) { argc--; a=atol(*++argv); }
					else error("no publickey\n");
					printf("a=%ld  pq=%ld ... ", a, pq);
					fi=euler(pq);
					printf("fi=%ld ... ", fi);
					chinese(a,fi, &x, &y, &g);
					printf("(privkey)x=%ld y=%ld g=%ld\n", x, y, g);
					break;

				case 'f':  // factor number
					if (argc>1) { argc--; p=atol(*++argv); }
					nf=factor(p, factors, powers);
					printfact(factors, powers, nf);
					break;
/*
				case 'q': // factor number using quadratic sieve
					if (argc>1) { argc--; p=atol(*++argv); }
					nf=qfactor(p, factors, powers);
					printfact(factors, powers, nf);
					break;
*/
				case 'p':  // generate primes
					if (argc>1) { argc--; x=atol(*++argv); } else { x=1000; }
					sieve=(char*)malloc((int)(x/8+1));
					if (sieve==NULL)
						error("Out of memory\n");
					erostatenes(sieve, x);
					y=printsieve(sieve, x);
					printf("\n\ntotal %ld primes upto %ld :\n", y, x);
					free(sieve); sieve=NULL;
					break;
				case 'e':  // encrypt
					if (argc>1) { argc--; pq=atol(*++argv); }
					else error("no key\n");
					if (argc>1) { argc--; a=atol(*++argv); }
					else error("no key\n");
					if (argc>1) { argc--; msg=atol(*++argv); }
					else { error("no data\n");  msg=123456789; }
					cyp=expmod(msg, a, pq);
					printf("crypted: %ld\n", cyp);
					break;
				default:
					error("Invalid option: %s\n", arg);
			}
		}
		else
			error("Invalid argument: %s\n", arg);
	}
	if (sieve)
		free(sieve);
	return 0;
}