A Basic ECC Crypto C Library

$ gcc -Wall -O2 -o crypto.out crypto.c 
$ ./crypto.out pgen 

x: [0][7][8]=[4880449572248532107701327692677150037572104965674444437385142643512]
y^2: [0][8][8]=[55380869580354224246510320945966028555317502525493429017926966133843051295792]
y: [0][8][8]=[29139577269574826820926935201858599433659429201914524387051377879407947248399]

$ ./crypto.out ecdh 

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    17101372528228046186917539708449150102286202153232625914396359604563
      * P(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)

      = Q(x, y) = (38616894501258673852970036866543659072635089662921568068340331697101576554270, 11643169480892671629858370512120070259223559167488419079368716073871895962191)


  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    12761708460950637051849206961623904377595162559992784492118523127724
      * P(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)

      = Q(x, y) = (32477568348650104756356831653582245276448206004887493791495532411741897328458, 28063399834915231969701720976238207473037393607389857712493721633327966471074)


  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    17101372528228046186917539708449150102286202153232625914396359604563
      * P(x, y) = (32477568348650104756356831653582245276448206004887493791495532411741897328458, 28063399834915231969701720976238207473037393607389857712493721633327966471074)

      = Q(x, y) = (50351522368134241142819091458530332577142187603512687882689541070260459723459, 50221873001104206223375144395698504985316658966829066883833880358017679689691)


  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    12761708460950637051849206961623904377595162559992784492118523127724
      * P(x, y) = (38616894501258673852970036866543659072635089662921568068340331697101576554270, 11643169480892671629858370512120070259223559167488419079368716073871895962191)

      = Q(x, y) = (50351522368134241142819091458530332577142187603512687882689541070260459723459, 50221873001104206223375144395698504985316658966829066883833880358017679689691)

$ ./crypto.out psig 1234
Sign:

d=[0][7][8]=[6538729008379423383518670268808646995288081785172776771199052697947]
P=(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)
dP=(x, y) = (14803911054618354585053684993263408691271520722671720792800356230205975217771, 11341018286853926408444620055947798723342444873304700039478772171855673884697)

h=[0][1][1]=[1234]

k=[0][7][8]=[8232207578348459001023088138872857936131335029042440438075457201646]
kP=(x, y) = (31369530948108575720266905846722071920323191212477028674668598541699791483098, 9139751042718096616566254849289158387164951389944576741553517798024823610265)
r=[0][15][40]=[149879290832663482850783820956010355668942537394521944880131051841396204961042266046319521062426903751537961621638691535228251607511236059222514]

Verify:

P=(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)
dP=(x, y) = (14803911054618354585053684993263408691271520722671720792800356230205975217771, 11341018286853926408444620055947798723342444873304700039478772171855673884697)

h=[0][1][1]=[1234]

kP=(x, y) = (31369530948108575720266905846722071920323191212477028674668598541699791483098, 9139751042718096616566254849289158387164951389944576741553517798024823610265)
r=[0][15][40]=[149879290832663482850783820956010355668942537394521944880131051841396204961042266046319521062426903751537961621638691535228251607511236059222514]

(ek+hd)P=(x, y) = (51608966471038425598161188547671772160458281462630354621077601757763167308540, 17341325372123419326180270058693292718501038313000707228169935748980907943680)
==
rP=(x, y) = (51608966471038425598161188547671772160458281462630354621077601757763167308540, 17341325372123419326180270058693292718501038313000707228169935748980907943680)

[GOOD]
$ ./crypto.out penc 1234
Encrypt:

P=(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)
dP=(x, y) = (19547865285467707344496388019666506767309887017163137556078134582670804301851, 36180672281078333342162157694938953906019362565570590334768364807054263309922)

h=[0][1][1]=[1234]

k=[0][7][8]=[3873062938517225834483697070677896532844396051944377142545351178872]
kP=(x, y) = (31307487414183777373465427652384958249613501222520858769326720798351611112341, 42327343332069595141617190350930625663948508466032477124028143798481079201841)
kdP=(x, y) = (17546014782981223520910010062776248455503090940217295620750099274078977893385, 40213720617506594547897783150731613585768378209283692064354760696398013829358)

Decrypt:

P=(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)
dP=(x, y) = (19547865285467707344496388019666506767309887017163137556078134582670804301851, 36180672281078333342162157694938953906019362565570590334768364807054263309922)

d=[0][7][8]=[22873181922607864299317581582398892665871088710637471622228302193764]
dkP=(x, y) = (17546014782981223520910010062776248455503090940217295620750099274078977892151, 40213720617506594547897783150731613585768378209283692064354760696398013829358)

t=[0][1][8]=[1234]

$ 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>

#include "../lib/sha256.c"
#include "../lib/aes256.c"
#include "../lib/ec.c"

bnum *bnrnd(int size)
{
	FILE *f;
	bnum *r;
	
	r = bninit(size);
	f = fopen("/dev/urandom", "r");
	fread(r->nums, sizeof(unsigned int), size - 1, f);
	fclose(f);
	r->leng = size;
	while ((r->leng > 1) && (r->nums[r->leng - 1] == 0)) { r->leng -= 1; }
	
	return r;
}

void eccp(ecc *e, char *xs)
{
	int psiz = (e->p)->size;
	bnum *x = bninit(psiz), *y = bninit(psiz), *t;
	
	bnum *xa = bninit(psiz * 4), *xx = bninit(psiz * 4);
	bnum *yy = bninit(psiz * 4), *vv = bninit(psiz * 4);
	
	if (xs != NULL) { t = bndec(xs); }
	else { t = bnrnd(psiz); }
	bncopy(t, x);
	
	while (1)
	{
		// yy = ((x * (x * x)) + (a * (x * x)) + x) % m
		bnzero(xx); bnmul(x, x, xx);
		bnzero(xa); bnmul(xx, e->a, xa);
		bnzero(vv); bnmul(x, xx, vv);
		
		bnzero(yy); bnadd(vv, xa, yy, 1); bnadd(yy, x, vv, 1);
		bndiv(vv, e->p, xx, t);
		// todo: divide t by b?
		
		printf("\n");
		bnout("x: ", x, "\n");
		bnout("y^2: ", t, "\n");
		
		if (sqrtmod(t, e->p, y) == 0) { break; }
		
		if (x != NULL) { bnfree(x); }
		x = bnrnd(psiz);
	}
	
	bnout("y: ", y, "\n\n");
	// note: another y point == p - y
	
	bnzero(yy); bncopy(t, yy);
	bnzero(vv); bnmul(y, y, vv);
	bnzero(xx); bndiv(vv, e->p, xa, xx);
	
	bnfree(xa); bnfree(vv); bnfree(t);
	
	if (bncmp(xx, yy) != 0)
	{
		bnout("y^2: ",xx," != "); bnout("y^2: ",yy,"\n");
	}
	
	bnfree(xx); bnfree(yy);
	bnfree(x); bnfree(y);
}

/*
Curve25519: a=486662, b=1, p=2^255 - 19
ECC DH: o * (n * P) == onP == noP == n * (o * P)
*/

char *ecdh(ecc *e, char *n, int o)
{
	char *r = n;
	bnum *m;
	ecc *f;
	
	if (r == NULL)
	{
		m = bnrnd((e->p)->size);
		r = bnstr(m);
	}
	
	else
	{
		m = bndec(n);
	}
	
	f = ecdup(e);
	pmul(m, e, f);
	
	if (o == 1)
	{
		char v[2050];
		bzero(v, 2050 * sizeof(char));
		strncat(v, "    "       , max(0, 2048 - strlen(v)));
		strncat(v, r            , max(0, 2048 - strlen(v)));
		strncat(v, "\n      * P", max(0, 2048 - strlen(v)));
		printf("\n");
		ecout(1,v,e,"\n\n");
		ecout(0,"      = Q",f,"\n");
		printf("\n");
	}
	
	bnfree(e->x); e->x = bndup(f->x);
	bnfree(e->y); e->y = bndup(f->y);
	bnfree(m);
	ecfree(f);
	
	return r;
}

/*
ECC - Private Sign / Public Verify

* initialize
- generate secret key          : d
- publish                      : P, dP

* private sign
- hash message                 : h = SHA256(m)
- secret unique multiplier     : k
- calculate public multipliers : kP
                               : e = hP(x)
                               : r = ((e * k) + (h * d))
- publish                      : kP, r

* public verify
                               : e = hP(x)
-                              : (e * kP) + (h * dP) == rP
*/

int ecsig(ecc *e, bnum *d, ecc *dp, char *hh, ecc *kp, bnum *r, int o)
{
	int a = 0, psiz = (e->p)->size;
	bnum *k, *h;
	bnum *t = bninit(psiz * 5), *u = bninit(psiz * 5);
	ecc *hp = ecdup(e), *ekp = ecdup(e), *hdp = ecdup(e), *ekhdp = ecdup(e), *rp = ecdup(e);
	ect *tadd = etinit(e);
	
	if (hh != NULL)
	{
		if (((kp->x)->leng == 1) && ((kp->x)->nums[0] == 0))
		{
			if ((d->leng == 1) && (d->nums[0] == 0))
			{
				k = bnrnd(psiz); bncopy(k, d); bnfree(k);
				pmul(d, e, dp);
				
				if (o == 1) { bnout("d=", d, "\n"); }
			}
		}
		
		if (o == 1)
		{
			ecout(0, "P=", e, "\n");
			ecout(0, "dP=", dp, "\n\n");
		}
		
		h = bndec(hh);
		pmul(h, e, hp);
		
		if (o == 1) { bnout("h=", h, "\n\n"); }
		
		if (((kp->x)->leng == 1) && ((kp->x)->nums[0] == 0))
		{
			k = bnrnd(psiz);
			pmul(k, e, kp);
			
			if (o == 1) { bnout("k=", k, "\n"); }
			
			bnzero(t); bnmul(hp->x, k, t);
			bnzero(u); bnmul(h, d, u);
			bnadd(t, u, r, 1);
			
			bnfree(k);
		}
		
		if (o == 1)
		{
			ecout(0, "kP=", kp, "\n");
			bnout("r=", r, "\n\n");
		}
		
		pmul(hp->x, kp, ekp);
		pmul(h, dp, hdp);
		padd(ekp, hdp, ekhdp, tadd);
		pmul(r, e, rp);
		
		if ((d->leng == 1) && (d->nums[0] == 0))
		{
			if (o == 1)
			{
				ecout(0, "(ek+hd)P=", ekhdp, "\n==\n");
				ecout(0, "rP=", rp, "\n\n");
			}
		}
		
		if (bncmp(ekhdp->x, rp->x) == 0)
		{
			if (bncmp(ekhdp->y, rp->y) == 0)
			{
				a = 1;
			}
		}
		
		bnfree(h);
	}
	
	bnfree(t); bnfree(u);
	ecfree(hp); ecfree(ekp); ecfree(hdp); ecfree(ekhdp); ecfree(rp);
	etfree(tadd);
	
	return a;
}

/*
ECC ElGamal - Public Encrypt / Private Decrypt

* initialize
- generate secret integer  : d
- publish                  : P, dP

* public encrypt
- generate secret key      : t = (SHA256(pwd) || tmpkey)
- secret unique multiplier : k = rand()
- public key encrypt       : r = k * P        = kP
                           : s = t + (k * dP) = t + kdP = kdP(x + t, y + t)
- encrypt optional message : i = rand()
                           : e = AES256CBC(m, i, t)
- publish                  : [ (r, s) , (i, e) ]

* private decrypt
- private key decrypt      : u = s - (d * r)
                               = (t + kdP) - (d * kP)
                               = kdP(x + t, y + t) - dkP(x, y)
                               = kdP(x - x + t)
                               = t
*/

void ecenc(ecc *e, bnum *d, ecc *dp, char *hh, ecc *kp, ecc *kdp, int o)
{
	int psiz = (e->p)->size;
	bnum *k, *h, *t;
	ecc *dkp;
	
	if (hh != NULL)
	{
		if ((d->leng == 1) && (d->nums[0] == 0))
		{
			k = bnrnd(psiz); bncopy(k, d); bnfree(k); k = NULL;
			pmul(d, e, dp);
		}
		
		if (o == 1)
		{
			ecout(0, "P=", e, "\n");
			ecout(0, "dP=", dp, "\n\n");
		}
		
		h = bndec(hh);
		
		if (((kp->x)->leng == 1) && ((kp->x)->nums[0] == 0))
		{
			if (o == 1) { bnout("h=", h, "\n\n"); }
			
			k = bnrnd(psiz);
			pmul(k, e, kp);
			
			pmul(k, dp, kdp);
			t = kdp->x; kdp->x = bninit(psiz + 1); bncopy(t, kdp->x); bnfree(t);
			bnadd(kdp->x, h, kdp->x, 1);
			
			if (o == 1)
			{
				bnout("k=", k, "\n");
				ecout(0, "kP=", kp, "\n");
				ecout(0, "kdP=", kdp, "\n\n");
			}
			
			bnfree(k);
		}
		
		else
		{
			dkp = ecdup(e);
			pmul(d, kp, dkp);
			t = bninit(psiz); bnsub(kdp->x, dkp->x, t, 1);
			
			if (o == 1)
			{
				bnout("d=", d, "\n");
				ecout(0, "dkP=", dkp, "\n\n");
				bnout("t=", t, "\n\n");
			}
			
			bnfree(t);
			ecfree(dkp);
		}
		
		bnfree(h);
	}
}

int main(int argc, char **argv)
{
	char *a = "486662", *b = "1", *p = "57896044618658097711785492504343953926634992332820282019728792003956564819949";
	char *x = "9", *y = "43114425171068552920764898935933967039370386198203806730763910166200978582548";
	char *n, *m, *z = NULL;
	char hash[256];
	unsigned char msg[16], key[256];
	bnum *d, *r, *s;
	bnum *t, *u, *v, *w;
	ecc *e, *f, *dp, *kp, *kdp;
	
	if (argc > 2) { z = argv[2]; }
	
	t = bndec(p);
	u = bninit(t->size); w = bndec(x); bncopy(w, u); bnfree(w);
	v = bninit(t->size); w = bndec(y); bncopy(w, v); bnfree(w);
	e = ecinit(bndec(a), bndec(b), t, u, v);
	
	if (strcmp(argv[1], "pgen") == 0)
	{
		eccp(e, z);
	}
	
	if (strcmp(argv[1], "pmul") == 0)
	{
		ecdh(e, argv[2], 1);
	}
	
	if (strcmp(argv[1], "ecdh") == 0)
	{
		f = ecdup(e);
		
		n = ecdh(e, NULL, 1);
		m = ecdh(f, NULL, 1);
		
		t = e->x; u = e->y;
		e->x = f->x; e->y = f->y;
		f->x = t; f->y = u;
		
		ecdh(e, n, 1);
		ecdh(f, m, 1);
		
		free(n); free(m);
		ecfree(f);
	}
	
	if (strcmp(argv[1], "psig") == 0)
	{
		d = bninit((e->p)->size);
		dp = ecdup(e); kp = ecdup(e);
		r = bninit((e->p)->size * 5);
		s = bninit((e->p)->size * 5);
		
		printf("Sign:\n\n");
		(kp->x)->leng = 1; (kp->x)->nums[0] = 0;
		(kp->y)->leng = 1; (kp->y)->nums[0] = 0;
		ecsig(e, d, dp, z, kp, r, 1);
		
		printf("Verify:\n\n");
		bnzero(d);
		if (ecsig(e, d, dp, z, kp, r, 1) != 0) { printf("[GOOD]\n"); }
		else { printf("\n[FAILED]\n"); }
		
		bnfree(d); bnfree(r); bnfree(s);
		ecfree(dp); ecfree(kp);
	}
	
	if (strcmp(argv[1], "penc") == 0)
	{
		d = bninit((e->p)->size);
		dp = ecdup(e); kp = ecdup(e); kdp = ecdup(e);
		
		printf("Encrypt:\n\n");
		(kp->x)->leng = 1; (kp->x)->nums[0] = 0;
		(kp->y)->leng = 1; (kp->y)->nums[0] = 0;
		ecenc(e, d, dp, z, kp, kdp, 1);
		
		printf("Decrypt:\n\n");
		ecenc(e, d, dp, z, kp, kdp, 1);
		
		bnfree(d);
		ecfree(dp); ecfree(kp); ecfree(kdp);
	}
	
	if (strcmp(argv[1], "hash") == 0)
	{
		sha256 hobj;
		sha256init(&hobj);
		sha256update(&hobj, (unsigned char *)argv[2], strlen(argv[2]));
		sha256final(&hobj, hash);
		printf("%s\n", hash);
	}
	
	if (strcmp(argv[1], "menc") == 0)
	{
		bzero(key, 256);
		strcpy((char *)key, argv[3]);
		aes256keys(key);
		
		bzero(msg, 16);
		strcpy((char *)msg, argv[2]);
		aes256core(msg, key, 0);
		
		int i;
		for (i = 0; i < 16; ++i) { printf("%02x", msg[i]); }
		printf("\n");
	}
	
	ecfree(e);
	
	return 0;
}

A Basic ECC Crypto C Library

Curve25519 ECDH (Using bn.c)

$ a="486662"; b="1"; p="57896044618658097711785492504343953926634992332820282019728792003956564819949"
$ x="9"; y="43114425171068552920764898935933967039370386198203806730763910166200978582548"
$ ./ec.out "$a" "$b" "$p" "$x" "$y"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    25401056036235045220436215630227531782477232918721757282083167532085
      * P(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)

      = Q(x, y) = (2152181765955508802144811366391548358206596269149408838626058763387933371482, 11578758651574146923540708489424215527254544219548915855368449621649032022464)

$ ./ec.out "$a" "$b" "$p" "$x" "$y"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    14923829837058818803507170043209691418601607100463140519943405743260
      * P(x, y) = (9, 43114425171068552920764898935933967039370386198203806730763910166200978582548)

      = Q(x, y) = (10616717952671290844210554362351196131260739681824015303968067548318407877443, 44474405861832920486283648404346377112839096552743000067399554821542048806382)

$ x="10616717952671290844210554362351196131260739681824015303968067548318407877443"; y="44474405861832920486283648404346377112839096552743000067399554821542048806382"
$ ./ec.out "$a" "$b" "$p" "$x" "$y" "25401056036235045220436215630227531782477232918721757282083167532085"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    25401056036235045220436215630227531782477232918721757282083167532085
      * P(x, y) = (10616717952671290844210554362351196131260739681824015303968067548318407877443, 44474405861832920486283648404346377112839096552743000067399554821542048806382)

      = Q(x, y) = (5796465362698668629448044005371479448781495853974848652041497206752313461872, 30198337259834952114045673142150509214022566212693868468740722556827195886413)

$ x="2152181765955508802144811366391548358206596269149408838626058763387933371482"; y="11578758651574146923540708489424215527254544219548915855368449621649032022464"
$ ./ec.out "$a" "$b" "$p" "$x" "$y" "14923829837058818803507170043209691418601607100463140519943405743260"

  1*(y^2) = x^3 + 486662*(x^2) + x (mod 57896044618658097711785492504343953926634992332820282019728792003956564819949)
    14923829837058818803507170043209691418601607100463140519943405743260
      * P(x, y) = (2152181765955508802144811366391548358206596269149408838626058763387933371482, 11578758651574146923540708489424215527254544219548915855368449621649032022464)

      = Q(x, y) = (5796465362698668629448044005371479448781495853974848652041497206752313461872, 30198337259834952114045673142150509214022566212693868468740722556827195886413)

$ 
#include "bn.c"

struct ecurve {
	bnum *a, *b, *p, *x, *y;
};

#define ecc struct ecurve

struct ectemp {
	bnum *i, *s, *xr, *yr;
	bnum *t, *u, *v;
	bnum *w, *h, *g;
};

#define ect struct ectemp

ecc *ecinit(bnum *a, bnum *b, bnum *p, bnum *x, bnum *y)
{
	ecc *r = malloc(1 * sizeof(ecc));
	r->a = a; r->b = b; r->p = p;
	r->x = x; r->y = y;
	return r;
}

ecc *ecdup(ecc *e)
{
	ecc *r = malloc(1 * sizeof(ecc));
	r->a = bndup(e->a); r->b = bndup(e->b); r->p = bndup(e->p);
	r->x = bndup(e->x); r->y = bndup(e->y);
	return r;
}

void ecfree(ecc *e)
{
	bnfree(e->a); bnfree(e->b); bnfree(e->p);
	bnfree(e->x); bnfree(e->y);
	free(e);
}

void ecout(int d, char *s, ecc *e, char *t)
{
	char *a = bnstr(e->a), *b = bnstr(e->b), *p = bnstr(e->p);
	char *x = bnstr(e->x), *y = bnstr(e->y);
	char as[2], bs[2];
	as[0] = '+'; as[1] = '\0';
	bs[0] = '\0'; bs[1] = '\0';
	if ((e->a)->sign == 1) { as[0] = '-'; }
	if ((e->b)->sign == 1) { bs[0] = '-'; }
	if (d == 1) { printf("  %s%s*(y^2) = x^3 %s %s*(x^2) + x (mod %s)\n", bs, b, as, a, p); }
	printf("%s", s);
	printf("(x, y) = (%s, %s)", x, y);
	printf("%s", t);
	free(a); free(b); free(p);
	free(x); free(y);
}

ect *etinit(ecc *e)
{
	ect *t = malloc(1 * sizeof(ect));
	int ss = max(1, (e->b)->size);
	ss = max(((e->a)->size * 2) + 2, ((e->p)->size * 2) + 2);
	ss = max(((e->x)->size * 2) + 2, ((e->y)->size * 2) + 2);
	t->i = bninit(ss); t->s = bninit(ss); t->xr = bninit(ss); t->yr = bninit(ss);
	t->t = bninit(ss); t->u = bninit(ss); t->v = bninit(ss);
	int tt = ((ss * 2) + 2);
	t->w = bninit(tt); t->h = bninit(tt); t->g = bninit(tt + 4);
	return t;
}

void etfree(ect *t)
{
	bnfree(t->i); bnfree(t->s); bnfree(t->xr); bnfree(t->yr);
	bnfree(t->t); bnfree(t->u); bnfree(t->v);
	bnfree(t->w); bnfree(t->h); bnfree(t->g);
	free(t);
}

// modular multiplicative inverse

void egcd(bnum *a, bnum *b, bnum *g)
{
	int size = ((a->size + b->size) * 3);
	// s = 0; news = 1
	bnum *s = bninit(size);
	bnum *news = bninit(size); news->nums[0] = 1;
	// r = b; newr = a
	bnum *r = bninit(size); bncopy(b, r);
	bnum *newr = bninit(size); bncopy(a, newr);
	// init some temp vars
	bnum *prev = bninit(size), *quot = bninit(size), *temp = bninit(size);
	while ((r->leng > 1) || (r->nums[0] > 0))
	{
		// quot = (newr / r)
		if ((r->leng == 1) && (r->nums[0] < 3))
		{
			bncopy(newr, quot);
			if (r->nums[0] > 1) { bnrshift(quot, 1); }
		}
		else { bndiv(newr, r, quot, temp); }
		// prev = s
		bncopy(s, prev);
		// s = (news - (quot * prev))
		bnzero(temp); bnmul(quot, prev, temp);
		bnsub(news, temp, s, 0);
		// news = prev
		bncopy(prev, news);
		// prev = r
		bncopy(r, prev);
		// r = (newr - (quot * prev))
		bnzero(temp); bnmul(quot, prev, temp);
		bnsub(newr, temp, r, 0);
		// newr = prev
		bncopy(prev, newr);
	}
	if (news->sign == 1)
	{
		// news = news + b
		bnadd(news, b, news, 0);
	}
	bncopy(news, g);
	bnfree(s); bnfree(news);
	bnfree(r); bnfree(newr);
	bnfree(prev); bnfree(quot); bnfree(temp);
}

// modular square root

int sqrtmod(bnum *a, bnum *p, bnum *r)
{
	bnum *o = bninit(1);
	o->nums[0] = 1; o->leng = 1; o->sign = 0;
	
	// legendre symbol
	// define if a is a quadratic residue modulo odd prime
	
	// g = (p - 1) / 2
	
	// p - 1
	bnum *qq = bndup(p);
	bnsub(qq, o, qq, 1);
	// (p - 1) / 2
	bnum *g = bndup(qq);
	bnrshift(g, 1);
	
	// l = pow(a, g, p)
	
	// pow(a, g, p)
	bnum *l = bninit(max(max(a->size, g->size), p->size) * 3);
	bnpowmod(a, g, p, l);
	if (bncmp(l, qq) == 0)
	{
		bnfree(o); bnfree(qq); bnfree(g); bnfree(l);
		return -1;
	}
	
	// factor p - 1 on the form q * (2 ^ s) (with Q odd)
	// q = p - 1; s = 0
	bnum *q = bndup(qq);
	bnum *s = bninit(p->size);
	while ((q->nums[0] % 2) == 0)
	{
		// s += 1; q /= 2
		bnadd(s, o, s, 1);
		bnrshift(q, 1);
	}
	
	// select a z which is a quadratic non resudue modulo p
	// z = 1
	bnum *z = bninit(p->size);
	z->nums[0] = 1; z->leng = 1; z->sign = 0;
	while (1)
	{
		// while (lsym(z, p) != -1)
		bnpowmod(z, g, p, l);
		if (bncmp(l, qq) == 0) { break; }
		// z += 1
		bnadd(z, o, z, 1);
	}
	// c = pow(z, q, p)
	bnum *c = bninit(max(max(z->size, q->size), p->size) * 3);
	bnpowmod(z, q, p, c);
	
	// search for a solution
	// f = ((q + 1) / 2)
	bnum *f = bndup(q);
	bnadd(f, o, f, 1); bnrshift(f, 1);
	// x = pow(a, f, p)
	bnpowmod(a, f, p, r);
	// t = pow(a, q, p)
	bnum *t = bninit(max(max(a->size, q->size), p->size) * 3);
	bnpowmod(a, q, p, t);
	// m = s
	bnum *m = bninit(p->size), *i = bninit(p->size), *e = bninit(p->size);
	bncopy(s, m);
	// u = 2
	bnum *u = bninit(1);
	u->nums[0] = 2; u->leng = 1; u->sign = 0;
	bnum *b = bninit(p->size * 4), *v = bninit(p->size * 4), *w = bninit(p->size * 4);
	while ((t->leng > 1) || (t->nums[0] != 1))
	{
		// find the lowest i such that t ^ (2 ^ i) = 1
		// i = 1; e = 2
		i->nums[0] = 1; i->leng = 1; i->sign = 0;
		e->nums[0] = 2; e->leng = 1; e->sign = 0;
		while (bncmp(i, m) < 0)
		{
			bnpowmod(t, e, p, l);
			if ((l->leng == 1) && (l->nums[0] == 1)) { break; }
			bnlshift(e, 1);
			bnadd(i, o, i, 1);
		}
		// update next value to iterate
		// (m - i - 1)
		bnsub(m, i, v, 0);
		bnsub(v, o, v, 0);
		// 2 ^ (m - i - 1)
		bnpowmod(u, v, p, l);
		// b = (c ^ (2 ^ (m - i - 1))) % p
		bnpowmod(c, l, p, b);
		// x = ((x * b) % p)
		bnzero(v); bnmul(r, b, v);
		bndiv(v, p, w, r);
		// b = (b * b) % p
		bnzero(v); bnmul(b, b, v);
		bndiv(v, p, w, b);
		// t = ((t * b) % p)
		bnzero(v); bnmul(t, b, v);
		bndiv(v, p, w, t);
		// c = b; m = i
		bncopy(b, c);
		bncopy(i, m);
	}
	
	bnfree(o); bnfree(qq); bnfree(g); bnfree(l);
	bnfree(q); bnfree(s); bnfree(z); bnfree(c);
	bnfree(f); bnfree(t); bnfree(m);
	bnfree(i); bnfree(e); bnfree(b);
	bnfree(u); bnfree(v); bnfree(w);
	
	// r = [x, p - x]
	
	return 0;
}

// montgomery curve arithmetic

void nmod(bnum *a, bnum *b)
{
	if (a->sign == 1) { bnadd(b, a, a, 0); }
}

void pdub(ecc *p, ecc *r, ect *t)
{
	// printf("2P=\n");
	
	// l = 3*x^2 + 2*a*x + 1 / 2*b*y
	
	// x^2
	bnzero(t->w); bnmul(p->x, p->x, t->w); bndiv(t->w, p->p, t->t, t->v);
	// 3*x^2
	bnadd(t->v, t->v, t->w, 1); bnadd(t->w, t->v, t->w, 1);
	// 2*a*x
	bnzero(t->h); bnmul(p->a, p->x, t->h); bnlshift(t->h, 1);
	bndiv(t->h, p->p, t->t, t->u); nmod(t->u, p->p);
	// 3*x^2 + 2*a*x + 1
	bnadd(t->w, t->u, t->g, 1);
	int x, o = 1;
	for (x = 0; (x < (t->g)->leng) && (o == 1); ++x)
	{
		o = 0; if ((t->g)->nums[x] == 0xffffffff) { o = 1; } (t->g)->nums[x] += 1;
	}
	if (o == 1) { (t->g)->nums[x] = 1; (t->g)->leng += 1; }
	bndiv(t->g, p->p, t->t, t->yr);
	// 1 / 2*b*y
	bnzero(t->w); bnmul(p->b, p->y, t->w); bnlshift(t->w, 1);
	bndiv(t->w, p->p, t->t, t->u); nmod(t->u, p->p); egcd(t->u, p->p, t->i);
	// 3*x^2 + 2*a*x + 1 / 2*b*y
	bnzero(t->w); bnmul(t->yr, t->i, t->w); bndiv(t->w, p->p, t->t, t->s);
	
	// xr = b*l^2 - a - 2*x
	
	// l^2
	bnzero(t->g); bnmul(t->s, t->s, t->g);
	// b*l^2 - a
	bnzero(t->w); bnmul(p->b, t->g, t->w);
	bnsub(t->w, p->a, t->w, 0);
	// 2*x
	bnzero(t->h); bncopy(p->x, t->h); bnlshift(t->h, 1);
	// b*l^2 - a - 2*x
	bnsub(t->w, t->h, t->w, 0);
	bndiv(t->w, p->p, t->t, t->xr); nmod(t->xr, p->p);
	
	// yr = ((3*x + a) * l) - b*l^3 - y
	
	// (3*x + a) * l
	bnadd(t->h, p->x, t->w, 1); bnadd(t->w, p->a, t->w, 0);
	bndiv(t->w, p->p, t->t, t->u); nmod(t->u, p->p);
	bnzero(t->w); bnmul(t->u, t->s, t->w);
	// l^3
	bncopy(t->g, t->h);
	bnzero(t->g); bnmul(t->h, t->s, t->g); bndiv(t->g, p->p, t->t, t->u);
	// b*l^3
	bnzero(t->h); bnmul(p->b, t->u, t->h);
	bndiv(t->h, p->p, t->t, t->u); nmod(t->u, p->p);
	// ((3*x + a) * l) - b*l^3 - y
	bnsub(t->w, t->u, t->w, 0); bnsub(t->w, p->y, t->w, 0);
	bndiv(t->w, p->p, t->t, t->yr); nmod(t->yr, p->p);
	
	(t->xr)->leng = (p->p)->size; bncopy(t->xr, r->x);
	(t->yr)->leng = (p->p)->size; bncopy(t->yr, r->y);
}

void padd(ecc *p, ecc *q, ecc *r, ect *t)
{
	// printf("P+Q=\n");
	
	// l = (Qy - Py) / (Qx - Px)
	
	// Qy - Py
	bnsub(q->y, p->y, t->yr, 0);
	// Qx - Px
	bnsub(q->x, p->x, t->xr, 0);
	bndiv(t->xr, p->p, t->t, t->u); nmod(t->u, p->p);
	// 1 / (Qx - Px)
	egcd(t->u, p->p, t->i);
	// (Qy - Py) / (Qx - Px)
	bnzero(t->w); bnmul(t->yr, t->i, t->w);
	bndiv(t->w, p->p, t->t, t->s);
	
	// xr = b*l^2 - a - Px - Qx
	
	// b*l^2 - a - Px - Qx
	bnzero(t->w); bnmul(t->s, t->s, t->w);
	bnzero(t->g); bnmul(p->b, t->w, t->g);
	bnsub(t->g, p->a, t->g, 0);
	bnsub(t->g, p->x, t->g, 0);
	bnsub(t->g, q->x, t->g, 0);
	bndiv(t->g, p->p, t->t, t->xr); nmod(t->xr, p->p);
	
	// yr = ((2*Px + Qx + a) * l) - b*l^3 - Py
	
	// 2*Px + Qx + a
	bnadd(p->x, p->x, t->t, 1);
	bnadd(t->t, q->x, t->t, 1);
	bnadd(t->t, p->a, t->t, 0);
	bndiv(t->t, p->p, t->u, t->v); nmod(t->v, p->p);
	// (2*Px + Qx + a) * l
	bnzero(t->w); bnmul(t->v, t->s, t->w); bndiv(t->w, p->p, t->t, t->u);
	// b*l^3
	bnzero(t->w); bnmul(t->s, t->s, t->w);
	bnzero(t->g); bnmul(t->w, t->s, t->g); bndiv(t->g, p->p, t->t, t->v);
	bnzero(t->w); bnmul(t->v, p->b, t->w);
	// ((2*Px + Qx + a) * l) - b*l^3 - Py
	bnsub(t->u, t->w, t->v, 0);
	bnsub(t->v, p->y, t->v, 0);
	bndiv(t->v, p->p, t->t, t->yr); nmod(t->yr, p->p);
	
	(t->xr)->leng = (p->p)->size; bncopy(t->xr, r->x);
	(t->yr)->leng = (p->p)->size; bncopy(t->yr, r->y);
}

void pmul(bnum *m, ecc *p, ecc *r)
{
	int init = 0;
	bnum *mul = bndup(m);
	ecc *b = ecdup(p);
	ect *t = etinit(p);
	while ((mul->leng > 1) || (mul->nums[0] > 0))
	{
		if ((mul->nums[0] % 2) == 1)
		{
			if (init == 0)
			{
				bnfree(r->x); r->x = bndup(b->x);
				bnfree(r->y); r->y = bndup(b->y);
			}
			else
			{
				padd(r, b, r, t);
			}
			init = 1;
		}
		pdub(b, b, t);
		bnrshift(mul, 1);
	}
	bnfree(mul);
	ecfree(b);
	etfree(t);
}

Curve25519 ECDH (Using bn.c)

bn.c – A Basic Big Number/Integer Calculator

$ time ./bn.out
35052111338673026690212423937053328511880760811579981620642802346685810623109850235943049080973386241113784040794704193978215378499765413083646438784740952306932534945195080183861574225226218879827232453912820596886440377536082465681750074417459151485407445862511023472235560823053497791518928820272257787786
^
89489425009274444368228545921773093919669586065884257445497854456487674839629818390934941973262879616797970608917283679875499331574161113854088813275488110588247193077582527278437906504015680623423550067240042466665654232383502922215493623289472138866445818789127946123407807725702626644091036502372545139713
%
145906768007583323230186939349070635292401872375357164399581871019873438799005358938369571402670149802121818086292467422828157022922076746906543401224889672472407926969987100581290103199317858753663710862357656510507883714297115637342788911463535102712032765166518411726859837988672111837205085526346618740053
=
1976620216402300889624482718775150

real	0m0.783s
user	0m0.774s
sys	0m0.005s
struct bigint {
	int sign, leng, size;
	unsigned int *nums;
};

#define bnum struct bigint

int min(int a, int b)
{
	if (a < b) { return a; }
	return b;
}

int max(int a, int b)
{
	if (a > b) { return a; }
	return b;
}

int oshift(char **a, int b, char c)
{
	// right shift string a ("234") and insert overflow char c ('1')
	int x, n = (b + 1);
	char *s = *a;
	s = realloc(s, n * sizeof(char));
	for (x = (n - 2); x > 0; --x) { s[x] = s[x - 1]; }
	s[0] = c; s[n - 1] = '\0';
	*a = s;
	return n;
}

char *bnstr(bnum *bint)
{
	int x, y, w, z, o, n;
	int r = 2, t = 2;
	char *result = malloc(r * sizeof(char));
	char *twos = malloc(t * sizeof(char));
	// result = 0
	result[0] = '0'; result[1] = '\0';
	// twos = 1
	twos[0] = '1'; twos[1] = '\0';
	// loop through the bn binary
	for (x = 0; x < bint->leng; ++x)
	{
		for (y = 0; y < 32; ++y)
		{
			// if bn binary[n] == 1 then add the 2^n to result
			if (((bint->nums[x] >> y) & 0x1) == 1)
			{
				o = 0; w = (t - 2); z = (r - 2);
				// if the length of twos is bigger than length of result then realloc result
				if (t > r)
				{
					result = realloc(result, t * sizeof(char));
					result[t - 1] = '\0';
					r = t;
				}
				// add twos to result
				while ((w > -1) || (z > -1))
				{
					int a = 0; if (w > -1) { a = (twos[w] - '0'); }
					int b = 0; if (z > -1) { b = (result[z] - '0'); }
					n = (a + b + o);
					result[w] = ('0' + (n % 10));
					o = (n / 10);
					--w; --z;
				}
				if (o > 0) { r = oshift(&result, r, '0' + o); }
			}
			// multiply twos * 2
			o = 0; w = (t - 2);
			while (w > -1)
			{
				n = (o + ((twos[w] - '0') * 2));
				twos[w] = ('0' + (n % 10));
				o = (n / 10);
				--w;
			}
			if (o > 0) { t = oshift(&twos, t, '0' + o); }
		}
	}
	free(twos);
	return result;
}

void bnout(char *s, bnum *p, char *t)
{
	char *o = bnstr(p);
	printf("%s[%d][%d][%d]=[%s]%s", s, p->sign, p->leng, p->size, o, t);
	free(o);
}

void bnzero(bnum *a)
{
	a->sign = 0;
	bzero(a->nums, a->size * sizeof(unsigned int));
	a->leng = 1;
}

void bncopy(bnum *src, bnum *dst)
{
	dst->sign = src->sign;
	bcopy(src->nums, dst->nums, src->leng * sizeof(unsigned int));
	dst->leng = src->leng;
}

void bnfree(bnum *bint)
{
	//todo:zero out all values
	free(bint->nums);
	free(bint);
}

void bndfree(bnum **bint)
{
	bnfree(bint[0]);
	bnfree(bint[1]);
	free(bint);
}

int bncmp(bnum *a, bnum *b)
{
	int x = (a->leng - 1), y = (b->leng - 1);
	while ((x > -1) && (a->nums[x] < 1)) { --x; }
	while ((y > -1) && (b->nums[y] < 1)) { --y; }
	if (x > y) { return 1; }
	if (y > x) { return -1; }
	while (x > -1)
	{
		if (a->nums[x] > b->nums[x]) { return 1; }
		if (b->nums[x] > a->nums[x]) { return -1; }
		--x;
	}
	return 0;
}

int bnhigh(bnum *a)
{
	int x, y, f = 0, hib = -1;
	for (x = (a->leng - 1); (x > -1) && (f == 0); --x)
	{
		for (y = 31; (y > -1) && (f == 0); --y)
		{
			if (((a->nums[x] >> y) & 0x1) == 1)
			{
				hib = ((x * 32) + y);
				f = 1;
			}
		}
	}
	return hib;
}

void bnrshift(bnum *a, int s)
{
	int x, leng = 0, i = (s / 32);
	s = (s % 32);
	for (x = 0; x < (a->leng - i); ++x)
	{
		if (i > 0)
		{
			a->nums[x] = a->nums[x + i];
		}
		if (s > 0)
		{
			if (x > 0)
			{
				a->nums[x - 1] = ((a->nums[x] << (32 - s)) | a->nums[x - 1]);
				if (a->nums[x - 1] > 0) { leng = x; }
			}
			a->nums[x] = (a->nums[x] >> s);
		}
		if (a->nums[x] > 0) { leng = (x + 1); }
	}
	if (leng < 1)
	{
		a->nums[0] = 0;
		leng = 1;
	}
	a->leng = leng;
}

void bnlshift(bnum *a, int s)
{
	int x, o = 0, i = (s / 32);
	unsigned int over = (a->nums[a->leng - 1] >> (32 - (s % 32)));
	s = (s % 32);
	// check to see if the left shift will have any carry over bits
	if ((s > 0) && (over > 0)) { o = 1; }
	// pre adjust the size of the bn to hold the left shift bits
	if (a->size < (a->leng + o + i))
	{
		a->size = (a->leng + o + i);
		a->nums = realloc(a->nums, a->size * sizeof(unsigned int));
	}
	a->leng += (o + i);
	// shift the bits from left to right (last to first) minus the overflow spot
	for (x = (a->leng - o - 1); x > -1; --x)
	{
		if (i > 0)
		{
			if ((x - i) > -1) { a->nums[x] = a->nums[x - i]; }
			else { a->nums[x] = 0; }
		}
		if (s > 0)
		{
			if ((x + 1) < (a->leng - o)) { a->nums[x + 1] = (a->nums[x + 1] | (a->nums[x] >> (32 - s))); }
			a->nums[x] = (a->nums[x] << s);
		}
	}
	// if we have some overflow bits then set them last now
	if (o == 1) { a->nums[a->leng - 1] = over; }
}

bnum *bninit(int ssiz)
{
	int x;
	bnum *a = malloc(1 * sizeof(bnum));
	a->sign = 0; a->leng = 1, a->size = ssiz;
	a->nums = malloc(ssiz * sizeof(unsigned int));
	for (x = 0; x < ssiz; ++x) { a->nums[x] = 0; }
	return a;
}

bnum *bndup(bnum *a)
{
	int x;
	bnum *r = malloc(1 * sizeof(bnum));
	r->sign = a->sign; r->leng = a->leng; r->size = a->size;
	r->nums = malloc(a->size * sizeof(unsigned int));
	for (x = 0; x < a->size; ++x) { r->nums[x] = a->nums[x]; }
	return r;
}

bnum *bndec(char *decstr)
{
	int x, y, r, n;
	int i = 0;
	unsigned long m, l = strlen(decstr);
	char numstr[l + 1], outstr[l + 1];
	// result = 0
	bnum *result = bninit(1);
	// copy the input decimal string into a temp num str
	strncpy(numstr, decstr, l);
	numstr[l] = '\0';
	m = l;
	// while there are numbers left to be divided by 2
	while ((m > 1) || (numstr[0] > '0'))
	{
		// if we are shifting above 32 bits then append a new bit block
		if (i > 31)
		{
			result->leng += 1; result->size += 1;
			result->nums = realloc(result->nums, result->size * sizeof(unsigned int));
			result->nums[result->leng - 1] = 0;
			i = 0;
		}
		// divide the number by 2 and store the result in a temp string
		r = 0;
		outstr[0] = '\0';
		for (x = 0, y = 0; (x < m) && (y < l); ++x, ++y)
		{
			// get the last remainder and the current digit
			n = ((r * 10) + (numstr[x] - '0'));
			// if the digit num is less than the divider 2 then add on the next digit
			if ((n < 2) && ((x + 1) < m) && ((y + 1) < l))
			{
				n = ((n * 10) + (numstr[x + 1] - '0'));
				if (y > 0) { outstr[y] = '0'; ++y; }
				++x;
			}
			// store the division and the remainder
			outstr[y] = ('0' + (n / 2));
			outstr[y + 1] = '\0';
			r = (n % 2);
		}
		// save the binary result
		result->nums[result->leng - 1] = ((r << i) | result->nums[result->leng - 1]);
		++i;
		// copy the divided result into the temp num str
		strncpy(numstr, outstr, l);
		numstr[y] = '\0';
		m = strlen(numstr);
	}
	return result;
}

int bnsub(bnum *, bnum *, bnum *, int);

int bnadd(bnum *a, bnum *b, bnum *r, int s)
{
	int x = 0, y = 0, z = 0, g = 0;
	unsigned int over = 0;
	if (s == 0)
	{
		if ((a->sign == 0) && (b->sign == 1))
		{
			b->sign = 0;
			bnsub(a, b, r, 0); g = r->sign;
			b->sign = 1;
			r->sign = g;
			return 0;
		}
		else if ((a->sign == 1) && (b->sign == 0))
		{
			a->sign = 0;
			bnsub(b, a, r, 0); g = r->sign;
			a->sign = 1;
			r->sign = g;
			return 0;
		}
	}
	while ((x < a->leng) || (y < b->leng))
	{
		unsigned int c = 0; if (x < a->leng) { c = a->nums[x]; }
		unsigned int d = 0; if (y < b->leng) { d = b->nums[y]; }
		unsigned int lo = ((c & 0xffff) + (d & 0xffff) + (over & 0xffff));
		unsigned int hi = (((c >> 16) & 0xffff) + ((d >> 16) & 0xffff) + ((over >> 16) & 0xffff) + ((lo >> 16) & 0xffff));
		over = ((hi >> 16) & 0xffff);
		r->nums[z] = (((hi << 16) & 0xffff0000) | (lo & 0xffff));
		++x; ++y; ++z;
	}
	if (over > 0) { r->nums[z] = over; ++z; }
	r->sign = (a->sign | b->sign);
	r->leng = z;
	return 1;
}

int bnsub(bnum *a, bnum *b, bnum *r, int s)
{
	int x = 0, y = 0;
	int n = 1, g = 0;
	bnum *t = a, *m = b, *temp;
	if (s == 0)
	{
		if ((a->sign == 0) && (b->sign == 1))
		{
			b->sign = 0;
			bnadd(a, b, r, 0); g = r->sign;
			b->sign = 1;
			r->sign = g;
			return 0;
		}
		else if ((a->sign == 1) && (b->sign == 0))
		{
			a->sign = 0;
			bnadd(a, b, r, 0);
			a->sign = 1;
			r->sign = 1;
			return 0;
		}
		else if ((a->sign == 1) && (b->sign == 1))
		{
			t = b; m = a;
		}
		s = bncmp(t, m);
		if (s < 0)
		{
			g = 1;
			temp = t;
			t = m;
			m = temp;
		}
		else if (s == 0)
		{
			r->sign = 0; r->leng = 1;
			r->nums[0] = 0;
			return 2;
		}
	}
	int indx = -1;
	unsigned int over = 0;
	while ((x < t->leng) || (y < m->leng))
	{
		unsigned int c = 0; if (x < t->leng) { c = t->nums[x]; }
		unsigned int d = 0; if (y < m->leng) { d = m->nums[y]; }
		if (x == indx) { c = over; indx = -1; }
		if (c < d)
		{
			int i = x, z;
			for (z = (x + 1); z < t->leng; ++z)
			{
				if (t->nums[z] > 0)
				{
					unsigned int hic = ((1 << 16) + (c >> 16)), loc = (c & 0xffff);
					unsigned int hid = (d >> 16)              , lod = (d & 0xffff);
					if (loc < lod) { hic -= 1; loc += (1 << 16); }
					r->nums[i] = (((hic - hid) << 16) + (loc - lod));
					indx = z; over = (t->nums[z] - 1);
					if (r->nums[i] > 0) { n = max(n, (i + 1)); }
					break;
				}
				++x; ++y;
				if (y < m->leng) { r->nums[z] = (0xffffffff - m->nums[y]); }
				else { r->nums[z] = 0xffffffff; }
				if (r->nums[z] > 0) { n = max(n, (z + 1)); }
			}
		}
		else { r->nums[x] = (c - d); }
		if (r->nums[x] > 0) { n = max(n, (x + 1)); }
		++x; ++y;
	}
	r->sign = g;
	r->leng = n;
	return 1;
}

void bnmul(bnum *a, bnum *b, bnum *r)
{
	int x, y;
	bnum *dub = bninit(r->size);
	bncopy(b, dub);
	for (x = 0; x < a->leng; ++x)
	{
		for (y = 0; y < 32; ++y)
		{
			if (((a->nums[x] >> y) & 0x1) == 1)
			{
				bnadd(r, dub, r, 1);
			}
			bnadd(dub, dub, dub, 1);
		}
	}
	r->sign = (a->sign ^ b->sign);
	bnfree(dub);
}

void bndiv(bnum *a, bnum *b, bnum *r, bnum *m)
{
	int s = 0, d = 0, c = bncmp(a, b);
	bnum *t = bninit(max(a->size, b->size) * 2);
	bncopy(a, t);
	r->nums[0] = 0; r->leng = 1;
	m->nums[0] = 0; m->leng = 1;
	if (c < 0)
	{
		//printf("return <\n");
		r->nums[0] = 0; r->leng = 1;
		bncopy(a, m);
	}
	else if (c == 0)
	{
		//printf("return ==\n");
		r->nums[0] = 1; r->leng = 1;
		m->nums[0] = 0; m->leng = 1;
	}
	else if ((b->leng == 1) && (b->nums[0] == 0))
	{
		//printf("return /0!\n");
	}
	else if (c > 0)
	{
		int hir = 0;
		
		//printf("> right shift the number to match the divider\n");
		int hia = bnhigh(a), hib = bnhigh(b), hid = (hia - hib - 1);
		bnrshift(t, hid + 1);
		c = bncmp(t, b);
		if (c < 0) { s = 1; }
		
		while (1)
		{
			if (s == 1)
			{
				//printf("equal bit length but number < divider\n");
				int abyte = (hid / 32), abit = (hid % 32);
				int bitv = ((a->nums[abyte] & (1 << abit)) >> abit);
				bnlshift(t, 1); t->nums[0] |= bitv;
				--hid;
				if (d > 0) { bnlshift(r, 1); ++hir; }
			}
			s = 0;
			
			//printf("substract num - div and shift 1 into result\n");
			bnsub(t, b, t, 1);
			bnlshift(r, 1); r->nums[0] |= 1; ++hir;
			
			int hit = bnhigh(t);
			d = 0;
			while ((hit < hib) && (hid > -1))
			{
				int abyte = (hid / 32), abit = (hid % 32);
				int diff = min(hib - hit, 32);
				if ((hid + 1) >= diff)
				{
					//printf("large: shift another bit into the number to match the divider\n");
					hid -= (diff - 1);
					int bbyte = (hid / 32), bbit = (hid % 32);
					--hid;
					bnlshift(t, diff);
					t->nums[0] = (t->nums[0] | ((a->nums[abyte] << (31 - abit)) >> (32 - diff)));
					if (abyte != bbyte) { t->nums[0] = (t->nums[0] | (a->nums[bbyte] >> bbit)); }
					d += diff;
				}
				else
				{
					//printf("small: shift another bit into the number to match the divider\n");
					int bitv = ((a->nums[abyte] & (1 << abit)) >> abit);
					bnlshift(t, 1); t->nums[0] |= bitv;
					--hid;
					++d;
				}
				hit = bnhigh(t);
			}
			if (d > 1) { bnlshift(r, d - 1); hir += (d - 1); }
			
			//printf("check resulting number and break if done dividing else repeat\n");
			c = bncmp(t, b);
			if (c < 0)
			{
				if (hid < 0)
				{
					//printf("exiting <\n");
					if (d > 0) { bnlshift(r, 1); ++hir; }
					bncopy(t, m);
					break;
				}
				s = 1;
			}
			else if (c == 0)
			{
				if (hid < 0)
				{
					//printf("exiting ==\n");
					bnlshift(r, 1); r->nums[0] |= 1; ++hir;
					break;
				}
			}
		}
	}
	r->sign = (a->sign ^ b->sign);
	m->sign = r->sign;
	bnfree(t);
}

void bnpowmod(bnum *b, bnum *e, bnum *m, bnum *r)
{
	int s = max(max(b->size, e->size), max(m->size, r->size));
	bnum *t = bninit(s * 3), *u = bninit(s * 3);
	bnum *base = bninit(s * 3), *exp = bninit(s * 3);
	bncopy(b, base); bncopy(e, exp);
	r->nums[0] = 1; r->leng = 1;
	while ((exp->leng > 1) || (exp->nums[0] > 0))
	{
		if ((exp->nums[0] % 2) == 1)
		{
			// r = r * base
			bnzero(t); bnmul(r, base, t);
			// r = r % m
			bndiv(t, m, u, r);
		}
		// exp = exp / 2
		bnrshift(exp, 1);
		// b = b * b
		bnzero(t); bnmul(base, base, t);
		// b = b % m
		bndiv(t, m, u, base);
	}
	// free the rest
	bnfree(base); bnfree(exp);
	bnfree(t); bnfree(u);
}

bn.c – A Basic Big Number/Integer Calculator

Helping people study – my recurring job in life (Python, Classes, LinkedLists, QuickSort, BinarySearch)

$ python study.py 
Unsorted:[
  {'last': 'Chiappetta', 'user': 'johnny', 'key': 'abcxyz', 'first': 'Jonathan'}
  {'last': 'Dada', 'user': 'zdada', 'key': 'iowa', 'first': 'Zain'}
  {'last': 'Xray', 'user': 'yx', 'key': 'xy', 'first': 'Yvan'}
  {'last': 'Bob', 'user': 'abob', 'key': 'zyxcba', 'first': 'Alex'}
  {'last': 'Bob', 'user': 'bchar', 'key': 'hmmm', 'first': 'Charlie'}
]
Sorted:[
  {'last': 'Bob', 'user': 'abob', 'key': 'zyxcba', 'first': 'Alex'}
  {'last': 'Bob', 'user': 'bchar', 'key': 'hmmm', 'first': 'Charlie'}
  {'last': 'Chiappetta', 'user': 'johnny', 'key': 'abcxyz', 'first': 'Jonathan'}
  {'last': 'Xray', 'user': 'yx', 'key': 'xy', 'first': 'Yvan'}
  {'last': 'Dada', 'user': 'zdada', 'key': 'iowa', 'first': 'Zain'}
]
('Search=Dave:', None)
('Search=Jon:', {'last': 'Chiappetta', 'user': 'johnny', 'key': 'abcxyz', 'first': 'Jonathan'})
('Search=J Chi:', {'last': 'Chiappetta', 'user': 'johnny', 'key': 'abcxyz', 'first': 'Jonathan'})
('Search=J Chip:', None)
class contact:
	def __init__(self, firstName, lastName, username, pubKey):
		self.data = {"first":firstName, "last":lastName, "user":username, "key":pubKey}
		self.next = None
	
	def getData(self):
		return self.data
	
	def setData(self, dataContact):
		self.data = dataContact
	
	def getNext(self):
		return self.next
	
	def setNext(self, nextContact):
		self.next = nextContact

class contactList:
	def __init__(self):
		self.head = None
	
	def append(self, contact):
		if (self.head == None):
			self.head = contact
		else:
			temp = self.head
			while (temp.getNext() != None):
				temp = temp.getNext()
			temp.setNext(contact)
	
	def get(self):
		contacts = ""
		temp = self.head
		while (temp != None):
			contacts += ("  " + str(temp.getData()) + "\n")
			temp = temp.getNext()
		return ("[\n" + contacts + "]")
	
	def leng(self):
		indx = 0
		temp = self.head
		while (temp != None):
			indx += 1
			temp = temp.getNext()
		return indx
	
	def getin(self, gind):
		indx = 0
		temp = self.head
		while (temp != None):
			if (indx == gind):
				return temp
			temp = temp.getNext()
			indx += 1
		return None
	
	def swap(self, a, b):
		ac = self.getin(a)
		bc = self.getin(b)
		temp = ac.getData()
		ac.setData(bc.getData())
		bc.setData(temp)
	
	def sort(self, beg=-1, end=-1):
		if (beg < 0):
			beg = 0
		if (end < 0):
			end = self.leng()
		if ((end - beg) < 2):
			return 0
		# todo: use head && next instead of index lookup
		pivot = beg
		p = self.getin(pivot)
		z = p.getData()
		indx = (pivot + 1)
		hi = -1
		while (indx < end):
			a = self.getin(indx)
			b = a.getData()
			if ((b["first"] + b["last"]) < (z["first"] + z["last"])):
				if (hi > -1):
					self.swap(hi, indx)
					hi += 1
			else:
				if (hi < 0):
					hi = indx
			indx += 1
		if (hi < 0):
			hi = (end - 1)
		elif (hi > 0):
			hi = (hi - 1)
		self.swap(pivot, hi)
		self.sort(beg, hi)
		self.sort(hi + 1, end)
	
	def search(self, fullname, beg=-1, end=-1):
		if (beg < 0):
			beg = 0
		if (end < 0):
			end = self.leng()
		if ((end - beg) < 1):
			return None
		fulllist = fullname.split(" ")
		fulllist.append("")
		first = fulllist[0]; flen = len(first)
		last = fulllist[1]; llen = len(last)
		mid = (((end - beg) / 2) + beg)
		mc = self.getin(mid)
		md = mc.getData()
		fname = md["first"][:flen]
		lname = md["last"][:llen]
		if ((fname == first) and (lname == last)):
			return mc
		if ((fname > first) or (lname > last)):
			end = (mid - 1)
			mid = -1
		elif ((fname < first) or (lname < last)):
			beg = (mid + 1)
			mid = -1
		if (mid != -1):
			return None
		return self.search(fullname, beg=beg, end=end)

def main():
	contactBook = contactList()
	
	newContact = contact("Jonathan", "Chiappetta", "johnny", "abcxyz")
	contactBook.append(newContact)
	
	newContact = contact("Zain", "Dada", "zdada", "iowa")
	contactBook.append(newContact)
	
	newContact = contact("Yvan", "Xray", "yx", "xy")
	contactBook.append(newContact)
	
	newContact = contact("Alex", "Bob", "abob", "zyxcba")
	contactBook.append(newContact)
	
	newContact = contact("Charlie", "Bob", "bchar", "hmmm")
	contactBook.append(newContact)
	
	print("Unsorted:" + contactBook.get())
	
	contactBook.sort()
	print("Sorted:" + contactBook.get())
	
	print("Search=Dave:", contactBook.search("Charlie"))
	print("Search=Jon:", contactBook.search("Jon").getData())
	print("Search=J Chi:", contactBook.search("J Chi").getData())
	print("Search=J Chip:", contactBook.search("J Chip"))

if (__name__ == "__main__"):
	main()

Helping people study – my recurring job in life (Python, Classes, LinkedLists, QuickSort, BinarySearch)

AES 256 CTR Mode Hash Function

Proof of concept, do not use, just for theory!

With the recent news of SHA1 everybody should be switching to SHA2 or SHA3 by now. I tried turning AES 256 in CTR mode into a 256 bit hash function by XORing the encrypted outputs together. For example, splitting your message into 32 byte blocks and using it as the keys (m0, m1, …, mN):

hash(m) = [AES256(nonce || counter0, m0) || AES256(nonce || counter1, m0)]
XOR [AES256(nonce || counter2, m1) || AES256(nonce || counter3, m1)]
...
XOR [AES256(nonce || counterN, mN) || AES256(nonce || counterO, mN)]
('test', '8ea2b7ca516745bfeafc49904b496089')
-
('', '09975b45dc8ecebd1519328dbec1c54d66b7fa7a2a4b762368497f81d864819b')
('a', 'c458ace20dc6458d6e589dbf0e560cbcad5b482e5934a86b6e98202a6cd5a6d5')
('abc', '304dfabb945406b40d53160080d078a1a0e5f8330263f578725ecfbccbb2c0ad')
('cba', '76cdfb5e58e06b5bf191656986aa35e405d3f582f307ea3847fe2a48136bf34e')
('the quick brown fox jumps over the lazy dog', '205639b3097f25a1c3a7bee4a3a66c4c3786129dd92c57fdddafc81568ab66f5')
('the quick brown fox jumps over the lazy eog', '3a7a025b5eb99f93caffd09331786cf6a15e1cb5bf41ae68342ea6b2208cf539')
import sys

def subbytes(matrix):
	sbox = [0x63, 0x7C, 0x77, 0x7B, 0xF2, 0x6B, 0x6F, 0xC5, 0x30, 0x01, 0x67, 0x2B, 0xFE, 0xD7, 0xAB, 0x76,
		0xCA, 0x82, 0xC9, 0x7D, 0xFA, 0x59, 0x47, 0xF0, 0xAD, 0xD4, 0xA2, 0xAF, 0x9C, 0xA4, 0x72, 0xC0,
		0xB7, 0xFD, 0x93, 0x26, 0x36, 0x3F, 0xF7, 0xCC, 0x34, 0xA5, 0xE5, 0xF1, 0x71, 0xD8, 0x31, 0x15,
		0x04, 0xC7, 0x23, 0xC3, 0x18, 0x96, 0x05, 0x9A, 0x07, 0x12, 0x80, 0xE2, 0xEB, 0x27, 0xB2, 0x75,
		0x09, 0x83, 0x2C, 0x1A, 0x1B, 0x6E, 0x5A, 0xA0, 0x52, 0x3B, 0xD6, 0xB3, 0x29, 0xE3, 0x2F, 0x84,
		0x53, 0xD1, 0x00, 0xED, 0x20, 0xFC, 0xB1, 0x5B, 0x6A, 0xCB, 0xBE, 0x39, 0x4A, 0x4C, 0x58, 0xCF,
		0xD0, 0xEF, 0xAA, 0xFB, 0x43, 0x4D, 0x33, 0x85, 0x45, 0xF9, 0x02, 0x7F, 0x50, 0x3C, 0x9F, 0xA8,
		0x51, 0xA3, 0x40, 0x8F, 0x92, 0x9D, 0x38, 0xF5, 0xBC, 0xB6, 0xDA, 0x21, 0x10, 0xFF, 0xF3, 0xD2,
		0xCD, 0x0C, 0x13, 0xEC, 0x5F, 0x97, 0x44, 0x17, 0xC4, 0xA7, 0x7E, 0x3D, 0x64, 0x5D, 0x19, 0x73,
		0x60, 0x81, 0x4F, 0xDC, 0x22, 0x2A, 0x90, 0x88, 0x46, 0xEE, 0xB8, 0x14, 0xDE, 0x5E, 0x0B, 0xDB,
		0xE0, 0x32, 0x3A, 0x0A, 0x49, 0x06, 0x24, 0x5C, 0xC2, 0xD3, 0xAC, 0x62, 0x91, 0x95, 0xE4, 0x79,
		0xE7, 0xC8, 0x37, 0x6D, 0x8D, 0xD5, 0x4E, 0xA9, 0x6C, 0x56, 0xF4, 0xEA, 0x65, 0x7A, 0xAE, 0x08,
		0xBA, 0x78, 0x25, 0x2E, 0x1C, 0xA6, 0xB4, 0xC6, 0xE8, 0xDD, 0x74, 0x1F, 0x4B, 0xBD, 0x8B, 0x8A,
		0x70, 0x3E, 0xB5, 0x66, 0x48, 0x03, 0xF6, 0x0E, 0x61, 0x35, 0x57, 0xB9, 0x86, 0xC1, 0x1D, 0x9E,
		0xE1, 0xF8, 0x98, 0x11, 0x69, 0xD9, 0x8E, 0x94, 0x9B, 0x1E, 0x87, 0xE9, 0xCE, 0x55, 0x28, 0xDF,
		0x8C, 0xA1, 0x89, 0x0D, 0xBF, 0xE6, 0x42, 0x68, 0x41, 0x99, 0x2D, 0x0F, 0xB0, 0x54, 0xBB, 0x16
	]
	for i in range(0, len(matrix)):
		matrix[i] = sbox[matrix[i]]
	return matrix

def transform(matrix):
	t = []
	y = 0
	for x in range(0, 16):
		z = (x % 4)
		if ((x > 0) and (z == 0)):
			y += 1
		t.append(matrix[(z * 4) + y])
	return t

def shiftrows(matrix):
	val = matrix.pop(4)
	matrix.insert(7, val)
	# end row 1
	val = matrix.pop(8)
	matrix.insert(11, val)
	val = matrix.pop(8)
	matrix.insert(11, val)
	# end row 2
	val = matrix.pop(15)
	matrix.insert(12, val)
	# end row 3
	return matrix

def gmul(a, b):
	p = 0
	for c in range(0, 8):
		if ((b & 1) != 0):
			p = (p ^ a)
		hi_bit_set = (a & 0x80)
		a = ((a << 1) & 0xff)
		if (hi_bit_set != 0):
			a = (a ^ 0x1b)
		b = (b >> 1)
	return p

def mixcols(s):
	t = []
	for c in range(0, 16):
		t.append(0)
	for c in range(0, 4):
		t[c +  0] = (gmul(0x02, s[c + 0]) ^ gmul(0x03, s[c + 4]) ^ s[c + 8] ^ s[c + 12]);
		t[c +  4] = (s[c + 0] ^ gmul(0x02, s[c + 4]) ^ gmul(0x03, s[c + 8]) ^ s[c + 12]);
		t[c +  8] = (s[c + 0] ^ s[c + 4] ^ gmul(0x02, s[c + 8]) ^ gmul(0x03, s[c + 12]));
		t[c + 12] = (gmul(0x03, s[c + 0]) ^ s[c + 4] ^ s[c + 8] ^ gmul(0x02, s[c + 12]));
	return t

def addkey(matrix, keyval):
	for i in range(0, 16):
		matrix[i] = (matrix[i] ^ keyval[i])
	return matrix

def rcon(ind):
	c = 1
	if (ind == 0):
		return 0
	while (ind != 1):
		c = gmul(c, 2)
		ind -= 1
	return c

def keycore(subkey, i):
	val = subkey.pop(0)
	subkey.append(val)
	subkey = subbytes(subkey)
	subkey[0] = (subkey[0] ^ rcon(i))
	return subkey

def keyexp(matrix):
	c = 32
	i = 1
	t = [0, 0, 0, 0]
	while (c < 240):
		for a in range(0, 4):
			t[a] = matrix[a + c - 4]
		if ((c % 32) == 0):
			t = keycore(t, i)
			i += 1
		if ((c % 32) == 16):
			t = subbytes(t)
		for a in range(0, 4):
			if (c >= len(matrix)):
				matrix.append(0)
			matrix[c] = (matrix[c - 32] ^ t[a])
			c += 1
	return matrix

def hexout(matrix):
	o = ""
	for d in matrix:
		h = hex(d)
		h = str(h)
		h = h[2:]
		if (len(h) < 2):
			h = ("0" + h)
		o += h
	return o

def aescoree(msg, key):
	out = []
	
	t = []
	i = 0
	l = len(key)
	for x in range(0, 32):
		t.append(0)
		if (i < l):
			t[x] = ord(key[i])
			i += 1
	keyval = keyexp(t)
	
	i = 0
	l = len(msg)
	while (i < l):
		input = []
		for x in range(0, 16):
			input.append(0)
		for x in range(0, 16):
			if (i < l):
				input[x] = ord(msg[i])
				i += 1
		
		for r in range(0, 15):
			if (r == 0):
				input = addkey(input, keyval[r*16:])
			
			if (r > 0):
				input = subbytes(input)
				input = transform(input)
				input = shiftrows(input)
				
				if (r < 14):
					input = mixcols(input)
				
				input = transform(input)
			
			if ((r > 0) and (r < 15)):
				input = addkey(input, keyval[r*16:])
		
		for o in input:
			out.append(o)
	
	return out

def aesctr_hash(message):
	x = 0
	l = len(message)
	n = 0
	h = []
	c = "aesctrhash31337!"
	d = 0
	for e in c:
		d = ((d << 8) + ord(e))
	while ((x == 0) or (x < l)):
		u = ((d + n) & 0xffffffffffffffffffffffffffffffff); n += 1
		v = ((d + n) & 0xffffffffffffffffffffffffffffffff); n += 1
		r = ""
		s = ""
		while (u > 0):
			r = (chr(u & 0xff) + r)
			u = (u >> 8)
		while (v > 0):
			s = (chr(v & 0xff) + s)
			v = (v >> 8)
		m = message[x:x+32]
		a = aescoree(r, m)
		b = aescoree(s, m)
		t = (a + b)
		if (x == 0):
			h = t
		else:
			for i in range(0, 32):
				h[i] = (h[i] ^ t[i])
		x += 32
	print(message, hexout(h))
	return h

print("test", hexout(aescoree("\x00\x11\x22\x33\x44\x55\x66\x77\x88\x99\xaa\xbb\xcc\xdd\xee\xff", "\x00\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10\x11\x12\x13\x14\x15\x16\x17\x18\x19\x1a\x1b\x1c\x1d\x1e\x1f")))
print("-")
aesctr_hash("")
aesctr_hash("a")
aesctr_hash("abc")
aesctr_hash("cba")
aesctr_hash("the quick brown fox jumps over the lazy dog")
aesctr_hash("the quick brown fox jumps over the lazy eog")

AES 256 CTR Mode Hash Function

The Monty Hall Problem

I was watching MythBusters and they visited this problem which at Seneca we discussed briefly as well (I might have done this one before but thought I’d repost it anyway):

  1. There are 3 doors, 2 are losers and 1 is a winner
  2. You pick a door (2/3 chance it’s a loser)
  3. They show you the other loser door (leaving the winner door [1/3 chance] and the loser door [2/3 chance])
  4. You decide to switch your door (since it was 2/3 likely you originally chose a losing door and they took the other losing door away, switching increases your odds since the other door is likely to be the winner door)
  5. You win! (still a small chance of losing though)

$ python monty.py 0 100 1

('Doors:', {1: 'win', 2: 'lose', 3: 'lose'})
('Pick:', 1, 'win')
('Show:', 2, 'lose')
('Final:', 1, 'win')

('Doors:', {1: 'lose', 2: 'lose', 3: 'win'})
('Pick:', 2, 'lose')
('Show:', 1, 'lose')
('Final:', 2, 'lose')

('Doors:', {1: 'lose', 2: 'lose', 3: 'win'})
('Pick:', 1, 'lose')
('Show:', 2, 'lose')
('Final:', 1, 'lose')


('Wins:', 32, '/', 100, '=', 0.32)

$ python monty.py 1 100 1

('Doors:', {1: 'lose', 2: 'win', 3: 'lose'})
('Pick:', 2, 'win')
('Show:', 1, 'lose')
('Switching:', 2, 'to', 3)
('Final:', 3, 'lose')

('Doors:', {1: 'lose', 2: 'win', 3: 'lose'})
('Pick:', 3, 'lose')
('Show:', 1, 'lose')
('Switching:', 3, 'to', 2)
('Final:', 2, 'win')

('Doors:', {1: 'lose', 2: 'win', 3: 'lose'})
('Pick:', 1, 'lose')
('Show:', 3, 'lose')
('Switching:', 1, 'to', 2)
('Final:', 2, 'win')


('Wins:', 74, '/', 100, '=', 0.74)
import random
import sys
try:
	switch = int(sys.argv[1])
except:
	switch = 0
try:
	rounds = int(sys.argv[2])
except:
	rounds = 0
try:
	debug = int(sys.argv[3])
except:
	debug = 0
wins = 0
for x in range(0, rounds):
	doors = ["lose", "lose", "win"]
	random.shuffle(doors)
	doors = {1:doors[0], 2:doors[1], 3:doors[2]}
	if (debug != 0):
		print("Doors:", doors)
	pick = random.randint(1, 3)
	if (debug != 0):
		print("Pick:", pick, doors[pick])
	show = 1
	for k in doors.keys():
		if ((k != pick) and (doors[k] == "lose")):
			show = k
			break
	if (debug != 0):
		print("Show:", show, doors[show])
	if (switch != 0):
		for k in doors.keys():
			if ((k != pick) and (k != show)):
				other = k
				break
		print("Switching:", pick, "to", other)
		pick = other
	print("Final:", pick, doors[pick])
	if (doors[pick] != "lose"):
		wins += 1
	if (debug != 0):
		print("")
print("")
print("Wins:", wins, "/", rounds, "=", float(wins) / float(rounds))

The Monty Hall Problem

Numerical divider function

def div(num, den):
	s = str(num); o = ""
	x = 0; l = len(s)
	f = 0; t = []
	while (x < l):
		d = int(s[x]); p = o
		while (d < den):
			if ((x + 1) < l):
				d = ((d * 10) + int(s[x + 1]))
				x += 1
			else:
				if (p == ""):
					p += "0"
				if (not "." in p):
					p += "."
					f = len(p)
				d = (d * 10)
				if (d < den):
					p += "0"
		i = (d / den)
		r = (i * den)
		if (f > 0):
			if (d in t):
				p = o[f+t.index(d):]
				break
			else:
				t.append(d)
		o = (p + str(i)); p = ""
		if (r != d):
			s = str(d - r)
			x = 0; l = len(s)
		else:
			x += 1
	if (p == ""):
		p = "0"
	print(str(num) + " / " + str(den) + " = " + o + " (" + p + ") ")
div(2, 4)
div(7, 22)
div(1, 999)
div(22, 7)
div(1, 3)

2 / 4 = 0.5 (0) 
7 / 22 = 0.318 (18) 
1 / 999 = 0.001 (001) 
22 / 7 = 3.142857 (142857) 
1 / 3 = 0.3 (3) 
Numerical divider function

[boredness] ARC4-drop[8192] cipher based on ROT(256) instead of XOR – don’t ask why! :)

arc4-drop8192-rot256.py
import sys
def ksa(key):
	s = []
	for i in range(0, 256):
		s.append(i)
	j = 0
	for i in range(0, 256):
		j = ((j + s[i] + ord(key[i % len(key)])) % 256)
		t = s[i]
		s[i] = s[j]
		s[j] = t
	return s
def rot(inp, num):
	t = (inp + num)
	while (t > 255):
		t = (t % 256)
	while (t < 0):
		t = (t + 256)
	return t
def cipher(inp, key, dir):
	i = 0
	j = 0
	x = 0
	s = ksa(key)
	o = ""
	while (x < 8192):
		i = ((i + 1) % 256)
		j = ((j + s[i]) % 256)
		t = s[i]
		s[i] = s[j]
		s[j] = t
		x += 1
	x = 0
	while (x < len(inp)):
		i = ((i + 1) % 256)
		j = ((j + s[i]) % 256)
		t = s[i]
		s[i] = s[j]
		s[j] = t
		k = s[(s[i] + s[j]) % 256]
		o = (o + chr(rot(ord(inp[x]), dir * k)))
		x += 1
	return o
sys.stdout.write(cipher(sys.stdin.read(), sys.argv[1], int(sys.argv[2])))

$ echo "secret messageAAAA" | python arc4-drop8192-rot256.py "secret kez" "1" | hexdump -C
00000000  0c 06 1d 6c 9f 24 cb 6e  10 5f 5d 87 10 2e 2b 9a  |...l.$.n._]...+.|
00000010  90 1c b2                                          |...|
00000013
$ echo "secret messageAAAA" | python arc4-drop8192-rot256.py "secret kez" "1" | python arc4-drop8192-rot256.py "secret kez" "-1"
secret messageAAAA
[boredness] ARC4-drop[8192] cipher based on ROT(256) instead of XOR – don’t ask why! :)

A small simple SSH log watcher && iptables blocker in Python

TRYTIME = 8#s
TRYNUMB = 10#x
DROPTIME = 60#s
import os
import re
import signal
import subprocess
import sys
import time
os.system("iptables -F INPUT")
p = subprocess.Popen(["tail", "-n", "0", "-f", sys.argv[1]], stdout=subprocess.PIPE)
h = {}
def chek(i, t):
	s = int(time.time())
	d = pow(2, i[2] - 1)
	if (i[3] > 1):
		return 2
	if ((i[1] > 0) and (i[2] > 0) and (i[3] > 0) and ((s - i[0]) >= (t * d))):
		return 1
	return 0
while (1):
	q = os.fork()
	if (q == 0):
		while (1):
			s = int(time.time())
			for a in h.keys():
				if (chek(h[a], DROPTIME) == 1):
					c = ("iptables -D INPUT -s '%s' -j DROP > /dev/null 2>&1 ; #for %ds" % (a, DROPTIME * pow(2, h[a][2] - 1)))
					print("del",s,c)
					os.system(c)
					del h[a]
			time.sleep(1)
		sys.exit(0)
	else:
		s = int(time.time())
		for a in h.keys():
			if (chek(h[a], DROPTIME) == 1):
				h[a][3] = 2
				h[a][0] = s
				print("~",s,a,h[a])
			if (h[a][3] != 1):
				if ((s - h[a][0]) > TRYTIME):
					h[a][1] = max(0, h[a][1] - int((s - h[a][0]) / TRYTIME))
					h[a][0] = s
				if (h[a][1] < 1):
					print("x",s,a,h[a])
					del h[a]
		l = p.stdout.readline().strip()
		r = re.match("^.*sshd.*fail.*[^0-9]([0-9]+\.[0-9]+\.[0-9]+\.[0-9]+).*$", l, re.I)
		if (r):
			a = str(r.group(1))
			s = int(time.time())
			if (not a in h.keys()):
				h[a] = [s, 0, 0, 0];#address-string->[last-time,try-number,block-number,state-flag]
			if ((s - h[a][0]) <= TRYTIME):
				h[a][1] += 1
				print("+",s,a,h[a])
			if ((s - h[a][0]) > TRYTIME):
				h[a][1] = max(1, h[a][1] - int((s - h[a][0]) / TRYTIME))
				print("-",s,a,h[a])
			if (h[a][1] >= TRYNUMB):
				h[a][1] = (int(TRYNUMB / 2) + 1)
				h[a][2] += 1
				h[a][3] = 1
				c = ("iptables -A INPUT -s '%s' -j DROP ; #for %ds" % (a, DROPTIME * pow(2, h[a][2] - 1)))
				print("add",s,c)
				os.system(c)
			h[a][0] = s
		try:
			os.kill(q, signal.SIGKILL)
		except:
			pass
		try:
			os.waitpid(q, 0)
		except:
			pass

python sshipt.py /var/log/auth.log

('+', 1412665540, '1.93.29.79', [1412665540, 1, 0, 0])
('+', 1412665542, '1.93.29.79', [1412665540, 2, 0, 0])
('+', 1412665544, '1.93.29.79', [1412665542, 3, 0, 0])
('+', 1412665545, '1.93.29.79', [1412665544, 4, 0, 0])
('+', 1412665546, '1.93.29.79', [1412665545, 5, 0, 0])
('+', 1412665547, '1.93.29.79', [1412665546, 6, 0, 0])
('+', 1412665548, '1.93.29.79', [1412665547, 7, 0, 0])
('+', 1412665549, '1.93.29.79', [1412665548, 8, 0, 0])
('+', 1412665551, '1.93.29.79', [1412665549, 9, 0, 0])
('+', 1412665552, '1.93.29.79', [1412665551, 10, 0, 0])
('add', 1412665552, "iptables -A INPUT -s '1.93.29.79' -j DROP ; #for 60s")
A small simple SSH log watcher && iptables blocker in Python

New Buffalo DD-WRT Config && nix box && dnsmasq-dhcpd-tftpd-pxe

I got a new router (for use with Comcast) and I found a Buffalo one at Frys which runs full DD-WRT. Here is some sample configurations I’m using on it at the moment:

  • Home Network Goals: (w/ magic help from AP isolation mode)
    • Auths to FreeRADIUS with EAP-TTLS-MSCHAPv2 WPA2-CCMP-AES
    • Blocks ARP replies not from the correct modem/server/wifi
    • Blocks DHCP replies coming in from the wifi lan
    • Maps/learns/watches/monitors ARP replies from the correct clients
    • Supports IPv6 wan with the same requirements as above for IPv4

 

#!/usr/bin/python
import os
import re
import sys
import subprocess
import time
os.system("rm -fv /root/icmp*")
os.system("tcpdump -lnni br0 icmp6 -s65535 -C 1 -W 3 -w /root/icmp6 &")
hist = []
while (1):
	pids = []
	for x in range(0, 10):
		pobj = subprocess.Popen(["tcpdump", "-lnnr", "/root/icmp6"+str(x)], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
		plis = pobj.stdout.readlines()
		for line in plis:
			line = line.strip()
			indx = line.find("who has 2601:9:3400:aa9:1337:")
			if ((indx > -1) and (not line in hist)):
				addr = re.sub("[^0-9A-Fa-f:]+.*$", "", line[indx+8:])
				neih = subprocess.check_output(["ip", "-6", "neigh", "show"])
				if (not addr in neih):
					print("address",addr)
					pidn = os.fork()
					if (pidn == 0):
						os.system("ip -6 addr add '"+addr+"/128' dev br0 ; sleep 3 ; ip -6 addr del '"+addr+"/128' dev br0")
						sys.exit()
					else:
						pids.append(pidn)
				hist.append(line)
	while (len(hist) > 1000000):
		hist.pop(0)
	#print("sleeping...")
	time.sleep(5)
	for pidn in pids:
		try:
			os.waitpid(pidn, 0)
		except:
			pass
#!/bin/bash
while true
do
	cat /var/lib/misc/dnsmasq.leases | while read line
	do
		m=`echo "$line" | awk '{ print $2 }'`
		i=`echo "$line" | awk '{ print $3 }'`
		c=`echo "$i" | grep -i '[0-9][0-9]*\.[0-9][0-9]*\.[0-9][0-9]*\.[0-9][0-9]*'`
		if [ "$c" == "" ]
		then
			continue
		fi
		arp -s "$i" "$m"
	done
	sleep 5
done
#!/bin/bash


echo > /etc/resolv.conf
echo 'nameserver 4.2.2.1' >> /etc/resolv.conf
echo 'nameserver 8.8.8.8' >> /etc/resolv.conf


brctl addbr br0
brctl addif br0 eth0 eth2

ip link set dev br0 up
ip link set dev eth0 up
ip link set dev eth1 up
ip link set dev eth2 up


iptables -F ; iptables -X
iptables -F -t nat ; iptables -X -t nat

ip address add 10.0.0.10/24 dev br0
ip route add 0.0.0.0/0 via 10.0.0.1

iptables -t nat -A POSTROUTING -o br0 -s 10.10.10.0/24 -j SNAT --to 10.0.0.10

echo 1 > /proc/sys/net/ipv4/ip_forward


ip6tables -F ; ip6tables -X
ip6tables -F -t nat ; ip6tables -X -t nat

ip -6 address add 2601:9:3400:aa9::1337/80 dev br0
ip -6 route add ::/0 via 2601:9:3400:aa9::1

ip6tables -t nat -A POSTROUTING -o br0 -s 1337:1337:1337:1337::/64 -j SNAT --to 1337:1337:1337:1337::1337

echo 1 > /proc/sys/net/ipv6/conf/all/forwarding


ip address add 10.10.10.10/24 dev eth1
ip -6 address add 2601:9:3400:aa9:1337::1337/80 dev eth1

cat > /etc/dns.cfg << EOF
interface=eth1
listen-address=10.10.10.10
port=0
bind-interfaces
dhcp-range=10.10.10.20,10.10.10.90,24,10.10.10.255,1h
dhcp-range=2601:9:3400:aa9:1337:0000:0000:aaaa,2601:9:3400:aa9:1337:ffff:ffff:cccc,80,1h
dhcp-option=3,10.10.10.10
dhcp-option=6,4.2.2.1,8.8.8.8
enable-ra
EOF

killall tcpdump ; killall tcpdump
killall python ; killall python

/usr/sbin/dnsmasq -C /etc/dns.cfg
/bin/bash /root/sarp.sh &
/usr/bin/python /root/ipvs.py &
#firewall

echo > /etc/resolv.conf
echo 'nameserver 4.2.2.1' >> /etc/resolv.conf
echo 'nameserver 8.8.8.8' >> /etc/resolv.conf

ifconfig br0 2601:9:3400:aa9:1337::7331/80 up
ifconfig vlan4 2601:9:3400:aa9::7331/80 up

route -A inet6 add ::/0 gw 2601:9:3400:aa9:1337::1337

ebtables -F ; ebtables -X
ebtables -t nat -F ; ebtables -t nat -X

iptables -F ; iptables -X
iptables -t nat -F ; iptables -t nat -X
#pxe server install

curl -sL 'http://www.thekelleys.org.uk/dnsmasq/dnsmasq-2.71.tar.gz' > dns.tgz
tar -xzvf dns.tgz
cd dnsmasq*
make ; ( echo 'dhcp-leasefile=/tmp/dnsmasq.leases' ; echo 'dhcp-range=10.0.0.20,10.0.0.30,255.255.255.0,1h' ; echo 'dhcp-option=3,10.0.0.10' ; echo 'dhcp-boot=pxelinux.0' ; echo 'enable-tftp' ; echo 'tftp-root=/tmp/tftpd' ) > dns.cfg
mkdir -p /tmp/tftpd

curl -sL 'http://ftp.openbsd.org/pub/OpenBSD/5.5/amd64/pxeboot' > /tmp/tftpd/pxeboot ; curl -sL 'http://ftp.openbsd.org/pub/OpenBSD/5.5/amd64/bsd.rd' > /tmp/tftpd/bsd.rd
cp /tmp/tftpd/bsd.rd /tmp/tftpd/bsd

curl -sL 'http://ftp.nl.debian.org/debian/dists/wheezy/main/installer-amd64/current/images/netboot/netboot.tar.gz' > /tmp/tftpd/netboot.tar.gz
tar -xzvf /tmp/tftpd/netboot.tar.gz -C /tmp/tftpd/ ; cp -frv /tmp/tftpd/debian-installer/amd64/* /tmp/tftpd/

#sudo ifconfig en4 inet 10.0.0.10 netmask 255.255.255.0 up
sudo killall dnsmasq
sudo ./src/dnsmasq -C ./dns.cfg
New Buffalo DD-WRT Config && nix box && dnsmasq-dhcpd-tftpd-pxe