/******    bigint.h   *****/
/****************************************************************
*                                                               *
*       defines and structures for crypto section.              *
*               Jan. 21, 1995                                   *
*                                                               *
****************************************************************/

#define WORDSIZE        32
/*#define       MAXBITS         1024*/
#define MAXBITS         256
#define MAXLONG         (MAXBITS/WORDSIZE)
#define SUBMASK         0x80000000
#define LONGPOS         (MAXLONG-1)
#define MAXSHIFT        (WORDSIZE-1)
/*#define field_prime   1019
#define log2_fp         9
/*  that's the most significant power of 2, or floor(log2(field_prime)) */
#define field_prime   149
#define log2_fp         7
/*#define field_prime     13
#define log2_fp         3 */
#define NUMBITS         (field_prime-1)
#define NUMWORD         (NUMBITS/WORDSIZE)
#define UPRSHIFT        (NUMBITS%WORDSIZE)
#define UPRBIT          (1L<<(UPRSHIFT-1))
#define STRTPOS         (LONGPOS-NUMWORD)
#define UPRMASK         (~(-1L<<UPRSHIFT))

typedef short int       INDEX;
typedef unsigned long   ELEMENT;
typedef struct  {
        ELEMENT b[MAXLONG];
} BIGINT;

/******   eliptic.h   *****/
/****************************************************************
*                                                               *
*       These are structures used to create elliptic curve      *
*  points and parameters.  "form" is a just a fast way to check *
*  if a2 == 0.                                                  *
*               form            equation                        *
*                                                               *
*                0              y^2 + xy = x^3 + a_6            *
*                1              y^2 + xy = x^3 + a_2*x^2 + a_6  *
*                                                               *
****************************************************************/

typedef struct {
        INDEX   form;
        BIGINT  a2;
        BIGINT  a6;
} CURVE;

/*  projective coordinates for a point  */

typedef struct {
        BIGINT  x;
        BIGINT  y;
        BIGINT  z;
} POINT;

/*  started getting tired of writing this */

#define SUMLOOP(i) for(i=STRTPOS; i<MAXLONG; i++)

/******    eliptic_keys.h   *****/
/*  higher level structures to manipulate elliptic curve public keys  */

#define MAX_NAME_SIZE  64
#define MAX_PHRASE_SIZE  256

typedef struct {
  POINT  p;
  POINT  q;
  CURVE  crv;
  char  name[MAX_NAME_SIZE];
  char  address[MAX_NAME_SIZE];
}  PUBKEY;

/******   bigint.c   *****/
/************************************************************************
*                                                                       *
*       Routines to implement optimal normal basis multiplication.      *
*  See Mullin, Onyszchuk, Vanstone, Wilson, "Optimal Normal Bases in    *
*  GF(p^n)", Discrete Applied Math, V22, 1988, p149                     *
*   Ash, Blake, Vanstone, "Low Complexity Normal Bases", Discrete       *
*   Applied Math, V25, 1989, p191                                       *
* Agnew, Mullin, Onyszchuk, Vanstone, "An Implementation for a Fast     *
*  Public-Key Cryptosystem", Jour. Cryptology, 1991, V3, p 63           *
* "Elliptic Curve Public Key Cryptosystems", A. Menezes, Kluwer, 1993   *
*  pages 83-86.                                                         *
*                                                                       *
*                       Jan. 22, 1995                                   *
*                                                                       *
************************************************************************/

#include <stdio.h>
#include "bigint.h"

/*  global structure for all multiply routines.  compute once at begining
to save disk space.  */

INDEX  Lambda[field_prime];   /*  for multiply  */
ELEMENT mask_table[WORDSIZE];   /*  for shift_index  */

/*  shift routines assume bigendian structure.  operate in place */

void shift_left(a)
BIGINT *a;
{
        register INDEX i;
        register ELEMENT bit,temp;

        bit = 0;
        for (i=LONGPOS; i>=STRTPOS; i--) {
           temp = ( a->b[i] << 1) | bit;
           bit = (a->b[i] & SUBMASK) ? 1L : 0L;
           a->b[i] = temp;
        }
        a->b[STRTPOS] &= UPRMASK;
}

void shift_right(a)
BIGINT *a;
{
        register INDEX i;
        register ELEMENT bit,temp;

        bit = (a->b[STRTPOS] & 1) ? SUBMASK : 0;
        for (i=STRTPOS; i< MAXLONG; i++) {
           temp = ( a->b[i] >> 1) | bit;
           bit = (a->b[i] & 1) ? SUBMASK : 0;
           a->b[i] = temp;
        }
}

/* an entirely different way to do the same thing */
void rot_left(a)
BIGINT *a;
{
        register INDEX i;
        register ELEMENT bit,temp;

        bit = (a->b[STRTPOS] & UPRBIT) ? 1L : 0L;
        for (i=LONGPOS; i>=STRTPOS; i--) {
           temp = (a->b[i] & SUBMASK) ? 1L : 0L;
           a->b[i] = ( a->b[i] << 1) | bit;
           bit = temp;
        }
        a->b[STRTPOS] &= UPRMASK;
}

void rot_right(a)
BIGINT *a;
{
        register INDEX i;
        register ELEMENT bit,temp;

        bit = (a->b[LONGPOS] & 1) ? UPRBIT : 0L;
        for (i=STRTPOS; i< MAXLONG; i++) {
           temp = ( a->b[i] >> 1)  | bit;
           bit = (a->b[i] & 1) ? SUBMASK : 0L;
           a->b[i] = temp;
        }
        a->b[STRTPOS] &= UPRMASK;
}

void null(a)
BIGINT *a;
{
        register INDEX i;

        for (i=0; i<MAXLONG; i++)  a->b[i] = 0;
}

void copy (a,b)
BIGINT *a,*b;
{
        register INDEX i;

        for (i=0; i<MAXLONG; i++)  b->b[i] = a->b[i];
}

/* create Lambda [i,j] table.  indexed by j, each entry contains the
value of i which satisfies 2^i + 2^j = 1 || 0 mod field_prime.  There are
two 16 bit entries per index j except for zero.  See references for
details.  Since 2^0 = 1 and 2^2n = 1, 2^n = -1 and the first entry would
be 2^0 + 2^n = 0.  Multiplying both sides by 2, it stays congruent to
zero.  So Half the table is unnecessary since multiplying exponents by
2 is the same as squaring is the same as rotation once.  Lambda[0] stores
n = (field_prime - 1)/2.  The terms congruent to one must be found via
lookup in the log table.  Since every entry for (i,j) also generates an
entry for (j,i), the whole 1D table can be built quickly.
*/

void genlambda()
{
        INDEX i,logof,n,index;
        INDEX log2[field_prime],twoexp;

        for (i=0; i<field_prime; i++) log2[i] = -1;

/*  build antilog table first  */

        twoexp = 1;
        for (i=0; i<field_prime; i++) {
          log2[twoexp] = i;
          twoexp = (twoexp << 1) % field_prime;
        }

/*  compute n for easy reference */

        n = (field_prime - 1)/2;
        Lambda[0] = n;
        Lambda[1] = n;
        Lambda[n] = 1;

/*  loop over result space.  Since we want 2^i + 2^j = 1 mod field_prime
        it's a ton easier to loop on 2^i and look up i then solve the silly
        equations.  Think about it, make a table, and it'll be obvious.  */

        for (i=2; i<=n; i++) {
          index = log2[i];
          logof = log2[field_prime - i + 1];
          Lambda[index] = logof;
          Lambda[logof] = index;
        }
/*  last term, it's the only one which equals itself.  See references.  */

        Lambda[log2[n+1]] = log2[n+1];
}

/*  for use with multiply routine.  indexed rotation uses:
        a is source
        b is destination
        amount is bit position which gets moved to zero bit
  For the moment, starting with rotations if shift amount is same
        size as first or last words.  May want to divide that by 2
        depending on how slow shift is compared with bit chunk mover.
  Bit chunk mover needs mask table, which should be built during
        initialization.
*/

void initmask()
{
        register INDEX i;

        mask_table[0] = -1L;
        for (i=1; i<WORDSIZE; i++) mask_table[i] = (ELEMENT)~(-1L << i);
}

void shift_index(src,dst,amount)
BIGINT *src,*dst;
INDEX amount;
{
        INDEX srcptr, srcstrt, srcwrd, limit;
        INDEX dstptr, dststrt, dstwrd, bitsmvd, strt;
        ELEMENT temp, mask;
        INDEX srcbits, dstbits, shift;

        null(dst);

/*  first determine if simple rotation is worth doing.  This could use
        some tuning, depending on processor, bus time, etc.
*/
        limit = 3;
        if ((amount < limit) || (amount > (NUMBITS - limit))){
           copy(src,dst);
           if (! amount) return;
           if (amount < limit) {
              while (amount) {
                rot_right(dst);
                amount--;
              }
           } else {
              while (amount < NUMBITS) {
                rot_left(dst);
                amount++;
              }
           }
           return;
        }

/*  rather than move every word amount times, try and only read it 2 or 3 and
        put each sub component in the right place.
*/

/*  initialize bit pointers */

        srcptr = amount;
        dstptr = 0;

/*  scan at bit level  */

        while (dstptr < NUMBITS) {

/*  find word offsets into bigint array  */

           srcwrd = LONGPOS - srcptr/WORDSIZE;
           dstwrd = LONGPOS - dstptr/WORDSIZE;

/*  locate starting bit positions within each ELEMENT  */

           srcstrt = srcptr % WORDSIZE;
           dststrt = dstptr % WORDSIZE;

/* smallest value determines maximum bits that can be moved  */

           srcbits = (srcwrd == STRTPOS) ? UPRSHIFT - srcstrt
                                        : WORDSIZE - srcstrt;
           dstbits = (dstwrd == STRTPOS) ? UPRSHIFT - dststrt
                                        : WORDSIZE - dststrt;
           bitsmvd = (srcbits < dstbits) ? srcbits : dstbits;

/* max bits set in zero index of mask_table, then adjust mask over
        result and compute source shift amount.   */

           bitsmvd %= WORDSIZE;
           mask = mask_table[bitsmvd] << dststrt;
           temp = src->b[srcwrd];
           shift = abs(srcstrt - dststrt);

/*  create result, or into place */

           if (dststrt < srcstrt)
             temp = (temp >> shift) & mask;
           else
             temp = (temp << shift) & mask;
           dst->b[dstwrd] |= temp;

/*  bump source cyclicly, destination linearly */

           if (!bitsmvd) bitsmvd = WORDSIZE;
           srcptr = (srcptr + bitsmvd) % NUMBITS;
           dstptr = (dstptr + bitsmvd);
        }  /* until dstptr >= NUMBITS */
}

/*  Normal Basis Multiplication.  Assumes Lambda vector already initialized
        for type 1 normal basis.  See above references for details
                Output = c = a*b over GF(2^NUMBITS)
*/

void opt_mul(a,b,c)
BIGINT *a,*b,*c;
{
        register INDEX i,j;
        BIGINT a1,a2,bcopy;

        null(c);
        copy(b,&bcopy);

/*  for each rotation of B vector there are at most two rotations of A vector
in a type 1 normal basis.  Lambda is lookup table of A rotations from 2^i +
2^j = 1 mod field_prime.  For 2^i + 2^j = 0, need only one master shift and
single rotations thereafter.  The Lambda table is uniformly scrambled and
shows why this is an efficient bit mixing algorithm.
*/

        shift_index(a,&a1,Lambda[0]);
        for (i=STRTPOS; i<MAXLONG; i++)
                c->b[i] = bcopy.b[i] & a1.b[i];

/*  main loop, two lookups for every position */

        for (j=1; j<NUMBITS; j++) {
           rot_right(&bcopy);
           rot_right(&a1);
           shift_index(a,&a2,Lambda[j]);
           for (i=STRTPOS; i<MAXLONG; i++) {
/*              c->b[i] ^= b->b[i] & (a1.b[i] ^ a2.b[i]);
        expression too complex!!??!!    */
              a2.b[i] ^= a1.b[i];
              a2.b[i] &= bcopy.b[i];
              c->b[i] ^= a2.b[i];
           }
        }
}

/*  opt_inv computes the inverse of a normal basis, GF(2^n) "number".
        Enter with pointers to source number, destination storage,
Lambda table and mask_table.  Leaves source alone and puts its inverse
into destination.  The algorithm is based on the explanation given in
Menezes book.  the bit pattern of n-1 determines the number of multiplies.
a sliding window starting from msb determines the power of 2 (shift index)
and the next bit determines the number of multiplies for that power.
Seems pretty cool to this programmer!
*/

void opt_inv(src,dst)
BIGINT *src, *dst;
{
        INDEX m, x, y, bitpos;
        BIGINT ap1,ap2;

        copy(src, &ap1);
        null(dst);
        m = field_prime - 2; /* = NUMBITS - 1 */
        bitpos = log2_fp;

        x = m >> bitpos;
        if (x != 1) {
           printf("Error attempting to compute inverse!!\n");
           printf("x=%d m=%d bitpos=%d\n",x,m,bitpos);
           return;
        }
        while (x < m) {
           bitpos--;
           y = m >> bitpos;
           shift_index( &ap1, &ap2, NUMBITS-x);
           opt_mul( &ap1, &ap2, dst);
           if (y & 1) {
              rot_left(dst);            /* this squares destination */
              opt_mul( src, dst, &ap1);
           } else
              copy(dst, &ap1);
           x = y;
        }
        rot_left( &ap1); /* final squaring */
        copy (&ap1,dst);
}

void init_opt_math()
{

        initmask();             /*  create shift_index mask table  */
        genlambda();            /*  create Lambda pointer table    */
}

/******   eliptic.c   *****/
/************************************************************************
*                                                                       *
*       Elliptic curves over Galois Fields.  These routines find points *
*  on a curve, add and double points for type 1 optimal normal bases.   *
*                                                                       *
*  For succint explanaition, see Menezes, page 92                  *
*                                                                       *
************************************************************************/

#include <stdio.h>
#include "bigint.h"
#include "eliptic.h"

/************************************************************************
*  Note that the following is obvious to mathematicians.  I thought it  *
*  was pretty cool when I discovered it myself, <sigh>.      *
*                                                                       *
*       Routine to solve quadradic equation.  Enter with coeficients    *
*  a and b and it returns solutions y[2]: y^2 + ay + b = 0.             *
*  If Tr(b/a^2) != 0, returns y=0 and error code 1.                     *
*  If Tr(b/a^2) == 0, returns y[2] and error code 0.                    *
*  If solution fails, returns y=0 and error code 2.                     *
*                                                                       *
*      Algorithm used based on normal basis GF math.  Since (a+b)^2 =   *
*  a^2 + b^2 it follows that (a+b)^.5 = a^.5 + b^.5.  Note that squaring*
*  is a left shift and rooting is a right shift in a normal basis.      *
*  Transforming the source equation with y = ax and dividing by a^2     *
*  gives:                                                               *
*               x^2 + x + b/a^2 = 0                                     *
*                                                                       *
*       or      x = x^.5 + (b/a^2)^.5                                   *
*                                                                       *
*  Let k_i = the ith significant bit of (b/a^2)^.5 and                  *
*      x_i = the ith significant bit of x.                              *
*  The above equation is equivelent to a bitwise representation as      *
*                                                                       *
*               x_i = x_(i+1) + k_i                                     *
*       or                                                              *
*               x(i+1) = x_i + k_i.                                     *
*                                                                       *
*  Since both x and x+1 are solutions, and 1 is represented by all      *
*  bits set in a normal basis, we can start with x_0 = 0 or 1 at our    *
*  pleasure and use the recursion relation to discover every bit of x.  *
*  The answer is then ax and ax+a returned in y[0] and y[1] respectively*
*  If the sum of x_(n-1) + k_(n-1) != x_0, returns error code 2 and     *
*  y = 0.                                                               *
*                                                                       *
*       error code                      returns                         *
*          0                    y[0] and y[1] values                    *
*          1                    y[0] = y[1] = 0                         *
*          2                    mathematicly impossible !!!!            *
*                                                                       *
************************************************************************/

extern  ELEMENT  mask_table[WORDSIZE];
extern  opt_inv(), rot_left(), rot_right(), null(), opt_mul();
extern  big_print(),copy();

int gf_quadradic(a, b, y)
BIGINT *a, *b, *y;
{
        INDEX   i, l, bits;
        BIGINT  x, k, a2;
        ELEMENT  r, t, mask;

/*  test for b=0.  Won't work if it is.  */

        l = 1;
        for (i=STRTPOS; i<MAXLONG; i++) {
           if (b->b[i] == 0) continue;
           else {
              l = 0;
              break;
           }
        }
        if (l) {
           null(&y[0]);
           null(&y[1]);
           return(1);
        }

/*  find a^-2  */

        opt_inv( a, &a2);
        rot_left(&a2);

/*  find k=(b/a^2)^.5 */

        opt_mul( b, &a2, &k);
        rot_right(&k);
        r = 0;

/*  check that Tr(k) is zero.  Combine all words first. */

        for (i=STRTPOS; i<MAXLONG; i++)
                r ^= k.b[i];

/*  take trace of word, combining half of all the bits each time */

        for (bits = WORDSIZE/2; bits > 0; bits /= 2)
                r = ((r & mask_table[bits]) ^ (r >> bits));

/*  if not zero, return error code 1.  */

        if (r) {
           null(&y[0]);
           null(&y[1]);
           return(1);
        }

/*  point is valid, proceed with solution.  mask points to bit i,
which is known, in x bits previously found and k (=b/a^2).  */

        null(&x);
        mask = 1;
        for (bits=0; bits < NUMBITS ; bits++) {

/* source long word could be different than destination  */

           i = LONGPOS - bits/WORDSIZE;
           l = LONGPOS - (bits + 1)/WORDSIZE;

/*  use present bits to compute next one */

           r = k.b[i] & mask;
           t = x.b[i] & mask;
           r ^= t;

/*  same word, so just shift result up */

           if ( l == i ) {
              r <<= 1;
              x.b[l] |= r;
              mask <<= 1;
           } else {

/*  different word, reset mask and use a 1 */

              mask = 1;
              if (r) x.b[l] = 1;
           }
        }

/*  test that last bit generates a zero */

        r = k.b[STRTPOS] & UPRBIT;
        t = x.b[STRTPOS] & UPRBIT;
        if ( r^t ) {
           null(&y[0]);
           null(&y[1]);
           return(2);
        }

/*  convert solution back via y = ax */

        opt_mul(a, &x, &y[0]);

/*  and create complementary soultion y = ax + a */

  null (&y[1]);
        for (i = STRTPOS; i < MAXLONG;  i++)
                y[1].b[i] = y[0].b[i] ^ a->b[i];

/*  no errors, bye!  */

        return(0);
}

/*  compute R.H.S. f(x) = x^3 + a2*x^2 + a6  
    curv.form = 0 implies a2 = 0, so no extra multiply.  
    curv.form = 1 is the "twist" curve.
*/

fofx(x, curv, f)
BIGINT *x, *f;
CURVE *curv;
{

        BIGINT x2,x3;
        INDEX i;

        copy(x, &x2);
        rot_left(&x2);
        opt_mul(x, &x2, &x3);
        if (curv->form) opt_mul(&x2, &curv->a2, f);
  else null(f);
        for (i = STRTPOS; i < MAXLONG; i++)
             f->b[i] ^= (x3.b[i] ^ curv->a6.b[i]);
}

/****************************************************************************
*                                                                           *
*   Implement elliptic curve point addition for optimal normal basis form.  *
*   Enter with pointers to p1, p2, p3 and curve data for E.  Computes       *
*   p3 = p1 + p2 in projective coordinates over curve E.  See Menezes       *
*   "Elliptic Curve Public Key Cryptosystems."                              *
*   Had to derive this several times to get it correct.  Menezes form not   *
*   generic for both p1 and p2 being arbitrary projections.                 *
*   Note that this is not optimized to minimize multiplies, but to follow   *
*    the math easier.                  *
****************************************************************************/

esum (p1, p2, p3, curv)
POINT   *p1, *p2, *p3;
CURVE   *curv;
{
    INDEX   i;
    BIGINT  x1z2, y1z2, z1z2, x2z1, y2z1;
    BIGINT  A, B, AplusB, Bsqr, C, D, E;
    BIGINT  temp1, temp2;

/*  start with basic building blocks  */

    opt_mul (&p1->x, &p2->z, &x1z2);
    opt_mul (&p1->y, &p2->z, &y1z2);
    opt_mul( &p1->z, &p2->z, &z1z2);    /* z_1 * z_2 */
    opt_mul (&p2->x, &p1->z, &x2z1);
    opt_mul (&p2->y, &p1->z, &y2z1);

/*  find terms A and B  */

    SUMLOOP(i) {
       B.b[i] = x1z2.b[i] ^ x2z1.b[i];
       A.b[i] = y1z2.b[i] ^ y2z1.b[i];
    }

/*  Next group of terms A+B and B^2  */

    SUMLOOP(i) AplusB.b[i] = A.b[i] ^ B.b[i];
    copy( &B, &Bsqr);
    rot_left( &Bsqr);

/*  set up term C.  No multiply if a2 = 0  */

    copy( &B, &C);
    if (curv->form) {
       opt_mul(&curv->a2, &z1z2, &temp1);
       SUMLOOP(i) C.b[i] ^= temp1.b[i];
    }

/*  main term for x3, also used in y3 (D in Menezes book) */

    opt_mul(&A, &AplusB, &temp1);
    opt_mul(&z1z2, &temp1, &temp2);
    opt_mul(&Bsqr, &C, &D);
    SUMLOOP(i) D.b[i] ^= temp2.b[i];

/*  compute x3  */

    opt_mul( &B, &D, &p3->x);

/*  compute z3  */

    opt_mul( &z1z2, &B, &temp1);
    opt_mul(&Bsqr, &temp1, &p3->z);

/*  next major term for y  */

    opt_mul( &x1z2, &A, &temp1);
    opt_mul( &y1z2, &B, &E);
    SUMLOOP(i) E.b[i] ^= temp1.b[i];

/*  finally compute y3  */

    opt_mul( &AplusB, &D, &temp1);
    opt_mul( &Bsqr, &E, &temp2);
    SUMLOOP(i) p3->y.b[i] = temp1.b[i] ^ temp2.b[i];
}

/*  elliptic curve doubling routine for projective coordinates over normal
    basis.  Enter with p1, p3 as source and destination as well as curv
    to operate on.  Returns p3 = 2*p1.
*/

edbl (p1, p3, curv)
POINT *p1, *p3;
CURVE *curv;
{
    BIGINT  x1z1, x14, termd, temp1, temp2;
    INDEX   i;

    opt_mul( &p1->x, &p1->z, &x1z1);
    copy(&p1->x, &x14);              /*  x_1 ^ 4 */
    rot_left( &x14);
    rot_left( &x14);
    copy(&p1->z, &termd);       /*  z_1 ^ 4 */
    rot_left ( &termd);
    rot_left ( &termd);
    opt_mul ( &termd, &curv->a6, &temp1);
    SUMLOOP(i) termd.b[i] = x14.b[i] ^ temp1.b[i];  /*  D = x1^4 + a6*z1^4 */

/*  now compute results  */

    opt_mul ( &termd, &x1z1, &p3->x);
    copy ( &x1z1, &temp1);
    rot_left( &temp1);
    opt_mul ( &temp1, &x1z1, &p3->z);    /* z3 = A^3 */
    opt_mul ( &x1z1, &x14, &p3->y);      /* x1^5 * z1 */
    copy ( &p1->x, &temp1);
    rot_left( &temp1);          /* x1 ^ 2 */
    opt_mul ( &p1->y, &p1->z, &temp2);  /*  y1 * z1 */
    SUMLOOP(i) {
       temp2.b[i] ^= temp1.b[i];        /*  proto sum for D */
       p3->y.b[i] ^= p3->x.b[i];        /*  first 2 terms of y3 found */
    }
    opt_mul ( &termd, &temp2, &temp1);  /*  do multiply */
    SUMLOOP(i) p3->y.b[i] ^= temp1.b[i];    /* and finish y3 */
}

/*  subtract two points on a curve.  just negates p2 and does a sum.
    Returns p3 = p1 - p2 over curv.
*/

esub (p1, p2, p3, curv)
POINT   *p1, *p2, *p3;
CURVE   *curv;
{
    POINT   negp;
    INDEX   i;

    copy ( &p2->x, &negp.x);
    null (&negp.y);
    SUMLOOP(i) negp.y.b[i] = p2->x.b[i] ^ p2->y.b[i];
    copy ( &p2->z, &negp.z);
    esum (p1, &negp, p3, curv);
}

/*  need to move points around, not just values.  Optimize later.  */

void copy_point (p1, p2)
POINT *p1, *p2;
{
  copy (&p1->x, &p2->x);
  copy (&p1->y, &p2->y);
  copy (&p1->z, &p2->z);
}

/*  Routine to compute kP where k is an integer (base 2, not normal basis)
  and P is a point on an elliptic curve.  This routine assumes that K
  is representable in the same bit field as x, y or z values of P.
  This is for simplicity, larger or smaller fields can be independently 
  implemented.
    Enter with: integer k, source point P, curve to compute over (curv) and 
    Returns with: result point R.

  Reference: Koblitz, "CM-Curves with good Cryptografic Properties", 
  Springer-Verlag LNCS #576, p279 (pg 284 really), 1992
*/
ELEMENT  bit_table[32] = {1, 2, 4, 8, 16, 32, 64, 128, 256,
      0x200, 0x400, 0x800, 0x1000, 0x2000, 0x4000, 0x8000,
      0x10000, 0x20000, 0x40000, 0x80000,
      0x100000, 0x200000, 0x400000, 0x800000,
      0x1000000, 0x2000000, 0x4000000, 0x8000000,
      0x10000000, 0x20000000, 0x40000000, 0x80000000 };

void  elptic_mul(k, p, r, curv)
BIGINT  *k;
POINT  *p, *r;
CURVE  *curv;
{
  char    blncd[NUMBITS+1];
  ELEMENT    wordptr, bit_count, mskindx, sequnce;
  POINT    temp;

/*  scan across k from right to left to expand bits to chars  */

  for (bit_count = 0; bit_count < NUMBITS; bit_count++) {
     wordptr = LONGPOS - bit_count/WORDSIZE;
     mskindx = bit_count % WORDSIZE;
     if (k->b[wordptr] & bit_table[mskindx]) blncd[bit_count] = 1;
     else blncd[bit_count] = 0;
  }
  blncd[NUMBITS] = 0;

/*  Next convert balanced array to ballanced representation  */
/*  sequence flag used for multiple 1's being switched to 0's.  */

  sequnce = 0;
  bit_count = 0;
  while (bit_count <= NUMBITS) {
     if ( blncd[bit_count]) {
        if (sequnce)    /*  in middle of 1's sequence */
    blncd[bit_count] = 0;
        else if (blncd[bit_count+1]) {
    sequnce = 1;    /*  next bit also set, begin sequnce */
    blncd[bit_count] = -1;
        }
        bit_count++;
     } else {  
        if (sequnce) {     /* last bit of sequence */
    blncd[bit_count] = 1;
    sequnce = 0;
        } else
    bit_count++;
     }
  }

/*  now follow ballanced representation and compute kP  */

  bit_count = NUMBITS;
  while (!blncd[bit_count] && bit_count >= 0) bit_count--; /* find first bit */
  if (bit_count < 0) {
     null (&r->x);
     null (&r->y);
     null (&r->z);
     return;
  }
  copy_point(p,r);    /* first bit always set */
  while (bit_count > 0) {
    edbl(r, &temp, curv);
    bit_count--;
    switch (blncd[bit_count]) {
       case 1: esum (p, &temp, r, curv);
      break;
       case -1: esub (&temp, p, r, curv);
      break;
       case 0: copy_point (&temp, r);
     }
  }
}

/*  One is not what it appears to be.  In any normal basis, "1" is the sum of
all powers of the generator.  So this routine puts ones to fill the number size
being used in the address of the BIGINT supplied.  */

void one (place)
BIGINT *place;
{
  INDEX i;

  SUMLOOP(i) place->b[i] = -1L;
  place->b[STRTPOS] &= UPRMASK;
}

/*  This routine "flatens" a projective 3D point back to the 2D plane.  */

void flaten (point)
POINT *point;
{
  POINT invers;

  opt_inv( &point->z, &invers.z);
  opt_mul( &invers.z, &point->x, &invers.x);
  opt_mul( &invers.z, &point->y, &invers.y);
  copy (&invers.x, &point->x);
  copy (&invers.y, &point->y);
  one(&point->z);
}

/******   support.c   *****/
/*  Simple support functions for elliptic curve data manipulation  */

#include <stdio.h>
#include <strings.h>
#include "bigint.h"
#include "eliptic.h"
#include "eliptic_keys.h"

extern void init_opt_math();
extern gf_quadradic( BIGINT*, BIGINT*, BIGINT*);
extern void fofx( BIGINT*, CURVE*, BIGINT*);
extern void copy( BIGINT*, BIGINT*);
extern void null( BIGINT*);
extern void elptic_mul(BIGINT*, POINT*, POINT*, CURVE*);
extern void opt_inv(BIGINT*, BIGINT*);
extern void edbl(POINT*, POINT*, CURVE*);
extern void esum(POINT*, POINT*, POINT*, CURVE*);
extern void esub(POINT*, POINT*, POINT*, CURVE*);
extern void one( BIGINT*);
extern void flaten( POINT*);
extern ELEMENT bit_table[WORDSIZE];

/*  random seed is accessable to everyone, not best way, but functional.  */

unsigned long random_seed;

/*  below is from Mother code, till end of mother.  Above is all my fault.  */

#include <string.h>

static short mother1[10];
static short mother2[10];
static short mStart=1;

#define m16Long 65536L                          /* 2^16 */
#define m16Mask 0xFFFF          /* mask for lower 16 bits */
#define m15Mask 0x7FFF                  /* mask for lower 15 bits */
#define m31Mask 0x7FFFFFFF     /* mask for 31 bits */
#define m32Double  4294967295.0  /* 2^32-1 */

/* Mother **************************************************************
|       George Marsaglia's The mother of all random number generators
|               producing uniformly distributed pseudo random 32 bit values with
|               period about 2^250.
|
|       The arrays mother1 and mother2 store carry values in their
|               first element, and random 16 bit numbers in elements 1 to 8.
|               These random numbers are moved to elements 2 to 9 and a new
|               carry and number are generated and placed in elements 0 and 1.
|       The arrays mother1 and mother2 are filled with random 16 bit values
|               on first call of Mother by another generator.  mStart is the switch.
|
|       Returns:
|       A 32 bit random number is obtained by combining the output of the
|               two generators and returned in *pSeed.  It is also scaled by
|               2^32-1 and returned as a double between 0 and 1
|
|       SEED:
|       The inital value of *pSeed may be any long value
|
|       Bob Wheeler 8/8/94
|
|  removed double return since I don't need it.
*/


void Mother(unsigned long *pSeed)
{
        unsigned long  number,
                       number1,
                       number2;
        short          n,
                       *p;
        unsigned short sNumber;

                /* Initialize motheri with 9 random values the first time */
        if (mStart) {
                sNumber= *pSeed&m16Mask;   /* The low 16 bits */
                number= *pSeed&m31Mask;   /* Only want 31 bits */

                p=mother1;
                for (n=18;n--;) {
                        number=30903*sNumber+(number>>16);   
        /* One line multiply-with-cary */
                        *p++=sNumber=number&m16Mask;
                        if (n==9)
                                p=mother2;
                }
                /* make cary 15 bits */
                mother1[0]&=m15Mask;
                mother2[0]&=m15Mask;
                mStart=0;
        }

                /* Move elements 1 to 8 to 2 to 9 */
        memmove(mother1+2,mother1+1,8*sizeof(short));
        memmove(mother2+2,mother2+1,8*sizeof(short));

                /* Put the carry values in numberi */
        number1=mother1[0];
        number2=mother2[0];

                /* Form the linear combinations */

number1+=1941*mother1[2]+1860*mother1[3]+1812*mother1[4]+1776*mother1[5]+
         1492*mother1[6]+1215*mother1[7]+1066*mother1[8]+12013*mother1[9];

number2+=1111*mother2[2]+2222*mother2[3]+3333*mother2[4]+4444*mother2[5]+
         5555*mother2[6]+6666*mother2[7]+7777*mother2[8]+9272*mother2[9];

                /* Save the high bits of numberi as the new carry */
        mother1[0]=number1/m16Long;
        mother2[0]=number2/m16Long;
                /* Put the low bits of numberi into motheri[1] */
        mother1[1]=m16Mask&number1;
        mother2[1]=m16Mask&number2;

                /* Combine the two 16 bit random numbers into one 32 bit */
        *pSeed=(((long)mother1[1])<<16)+(long)mother2[1];

                /* Return a double value between 0 and 1 
        return ((double)*pSeed)/m32Double;  */
}

/*  save data associated with a point on a curve.  
  Enter with name of file, a curve and a valid point on that curve.
  Returns 0 on success, -1 on failure.
*/

int save_curve (name, curv, point)
char *name;
CURVE *curv;
POINT   *point;
{
  FILE  *save;
  int  err1, err2;

  save = fopen(name, "w");
  if ( !save) return(-1);
  err1 = fwrite(curv, sizeof(CURVE), 1, save);
  err2 = fwrite(point, sizeof(POINT), 1, save);
  if (!(err1 && err2)) return(-1);
  fclose(save);
  return(0);
}

/* get data saved to disk.
  Enter with name of file.
  Returns with curve and point on that curve restored, 0 function value
  or null results and -1 function value.
*/

int get_curve (name, curv, point)
char  *name;
CURVE  *curv;
POINT  *point;
{
  FILE  *getcrv;
  int  err1, err2;

  getcrv = fopen (name, "r");
  if (!getcrv) return(-1);
  err1 = fread(curv, sizeof(CURVE), 1, getcrv);
  err2 = fread(point, sizeof(POINT), 1, getcrv);
  fclose(getcrv);
  if (!(err1 && err2)) return(-1);
  return(0);
}

/*  Random number initialization requires some arbitrary input.  Request
  32 bits from coin tossing.  Feed to Mother of all random number
  generators.  Supposedly this generator has field size ~2^250.  That's
  roughly the size of the universe measured in cubic angstroms, so
  this should be as close to "random" as deterministic can get.
*/

void init_rand()
{
  FILE  *rand;
  INDEX  i;
  unsigned long mask;
  char  z1,cr;

  if (! (rand = fopen("random.seed", "r"))) {
     printf("\n pull out a coin.\n");
     printf(" chose one side as '0'\n");
     printf("  and the other side as '1'\n");
     printf("   flip me some bits please\n");
     random_seed = 0;
     mask = 1;
     i = 0;
     while (i<32) {
        printf("bit %d: ",i);
        scanf("%c%c", &z1, &cr);
        if (z1 == '1') {
    random_seed |= mask;
    i++;
    mask <<= 1;
        } else if (z1 == '0') {
    i++;
    mask <<= 1;
        }
     }
     return;
  }
  fread (&random_seed, sizeof(long), 1, rand);
  fclose(rand);
}

void close_rand()
{
  FILE  *rand;

  if (rand = fopen("random.seed", "w")) {
     fwrite(&random_seed, sizeof(long), 1, rand);
     fclose(rand);
     return;
  }
  printf("\n Catastrophic Failure: random.seed not saved\n");
}

/*  print out a BIGINT with a label, to standard out.  */

void big_print (strng, a)
char *strng;
BIGINT *a;
{
  int i;

  printf("%s",strng);
  SUMLOOP(i) printf("%lx ",a->b[i]);
  printf("\n");
}

/*  Starting to get lazy, This needs to be part of a "point object" ultimately */

void print_point(title, p3)
char *title;
POINT *p3;
{
        printf("\n%s\n",title);
        big_print("x : ",&p3->x);
        big_print("y : ",&p3->y);
        big_print("z : ",&p3->z);
}

/*  generate a random point on a random curve.  Since these values are public anyway
  no need to attempt security measures on random seed or resulting points  */

void rand_curv_pnt( point, curve)
POINT  *point;
CURVE  *curve;
{
  BIGINT  f, y[2];
  INDEX  j;

/*  generate a random regular curve  */

  curve->form = 0;
  SUMLOOP(j) {
     Mother(&random_seed);
     curve->a6.b[j] = random_seed;
  }
  curve->a6.b[STRTPOS] &= UPRMASK;

/*  generate a random point on that curve */

  one (&point->z);
  SUMLOOP(j) {
     Mother(&random_seed);
     point->x.b[j] = random_seed;
  }
  point->x.b[STRTPOS] &= UPRMASK;
  fofx (&point->x, curve, &f);
  while (gf_quadradic(&point->x, &f, &y[0]) > 0) {
     point->x.b[LONGPOS] += 1L;
     fofx(&point->x, curve, &f);
  }
  copy (&y[0], &point->y);
}

/*  This hash function is for educational purposes.  elliptic curves 
  have the property that there are some x's for which y^2 + x*y = f(x)
  has no solution for y.  Further, it takes 30 seconds to perform an
  elliptic multiply for one block of data.  A meg would be slow to hash.
*/

#define  WORDS_NEEDED  (LONGPOS-STRTPOS)

void eliptic_hash(num_words, data_ptr, result)
INDEX  num_words;
ELEMENT *data_ptr;
BIGINT  *result;
{
  static init=0;
  static CURVE hcurv;
  static POINT hpnt;
  BIGINT  nxt_blok;
  INDEX  j, wrd_cnt;
  POINT  hashed, dashed;

/*  initialize hash curve and point only once  */

  if (!init) {
     if (get_curve( "hash.curve", &hcurv, &hpnt) < 0) {
        printf("Error, can't open hash.curve\n");
        exit(0);
     }
     init = -1;
  }

/*  initialize hash output value  */

  SUMLOOP(j) {
     hashed.x.b[j] = hpnt.x.b[j];
     hashed.y.b[j] = hpnt.y.b[j];
     hashed.z.b[j] = hpnt.z.b[j];
  }

/*  grab a block of data that completely fills long words.  Zero fill for
  last unused block < max length.
*/
  wrd_cnt = num_words;
  null(&nxt_blok);
  while (wrd_cnt) {
     if (wrd_cnt >= WORDS_NEEDED) {
        for (j=0; j < WORDS_NEEDED; j++)
    nxt_blok.b[STRTPOS + j] = *data_ptr++;
        wrd_cnt -= WORDS_NEEDED;
     } else {
        null(&nxt_blok);
        while (wrd_cnt) nxt_blok.b[STRTPOS + wrd_cnt--] = *data_ptr++;
     }

/*  use block of data as multiplier to find next point on curve.  */

     printf(".");
     elptic_mul(&nxt_blok, &hashed, &dashed, &hcurv);
     copy_point(&dashed, &hashed);
  }
  copy (&hashed.x, result);
  printf("\n");
}

/*  smash bits around to create a key.  experimental and for fun ok?!
  Take last 4 bits of each pair of input characters to create one byte
  of data block.  Hash data block using LONGPOS-STRTPOS ELEMENTS as
  multipliers of a single point, until whole data block consumed.  More
  than one or 2 factors is probably stupid (unless they are relatively
  prime).  And it would take too long.  The 256 character input restriction
  is much larger than most people use anyway.
*/

void elptic_key_gen( string, key)
char *string;
BIGINT  *key;
{
  char  bit_string[128], *bs_ptr;
  char  byte0, byte1;
  int  byt_cnt;
  INDEX  num_elements;
  INDEX  i;

  byt_cnt = 0;
  bs_ptr = bit_string;
  for (i=0; i<128; i++) bit_string[i] = 0;
  while (*string) {
     byte0 = (*string & 0xf) << 4;
     string++;
     byte1 = *string & 0xf;
    *bs_ptr++ = byte0 | byte1;
     byt_cnt++;
     if (! *string) break;
     else string++;
  }
  num_elements = byt_cnt/(WORDSIZE/8);
  if (! num_elements) {
     printf("key size too small\n");
     return;
  }
  eliptic_hash( num_elements, bit_string, key);
}

/*  gnu complains about gets, build my own. replace with something better, please! */

int get_string(buf, max)
char *buf;
int  max;
{
  int limit;
  char *pointer;
  char ch;

  pointer = buf;
  limit = 0;
  fpurge(stdin);
  while ((ch = getchar()) != '\n') {
    *pointer = ch;
    pointer++;
    limit++;
    if (limit > max) return(max);
  }
  return(limit);
}

/*  generate a public key.  If full = 0, only generates secret key.  
  for full = 1, fills entire public key with new random value.
*/

void public_key_gen( skey, pkey, full)
BIGINT  *skey;
PUBKEY  *pkey;
INDEX  full;
{
  char  pass[MAX_PHRASE_SIZE];
  char  bit_string[MAX_PHRASE_SIZE/2];
  char  byte0, byte1;
  INDEX  byt_cnt, chr_cnt;

/*  NOTE:  this is not done correctly.  turning off echo is platform dependent.
  see routines in PGP (getstring(char *strbuf, unsigned maxlen, int echo) in
  random.c) for much better ways to do this.
*/

  printf("Enter pass phrase:\n");
  get_string(pass, MAX_PHRASE_SIZE);
  printf("\nGenerating secret key.\n");

  byt_cnt = 0;
  chr_cnt = 0;
  while (pass[chr_cnt]) {
     byte0 = (pass[chr_cnt++] & 0xF) << 4;
     byte1 = pass[chr_cnt] & 0xF;
     bit_string[byt_cnt++] = byte0 | byte1;
     if (! pass[chr_cnt]) break;
     else chr_cnt++;
  }
  eliptic_hash(byt_cnt/(WORDSIZE/8), bit_string, skey);
  if (!full) return;

/*  create random point and curve.  for large enough fields this is not too
  dangerous, but cardinality of curve and order of point really ought to
  be checked.
*/

  printf("\nGenerating public key.\n");
  rand_curv_pnt(&pkey->p, &pkey->crv);
  elptic_mul(skey, &pkey->p, &pkey->q, &pkey->crv);
  flaten(&pkey->q);
  printf("\nOK, now what name and address for this key?\n");
  printf("Name: ");
  get_string(&pkey->name, MAX_NAME_SIZE);
  printf("Address: ");
  get_string(&pkey->address, MAX_NAME_SIZE);
}

/*  Save a public key to dsik file in ascii format.  File name taken from field
  of public key, up to first blank or tab.  .PUB added as extension.

  Returns 0 if ok, -1 on failure.
*/

static const char extensn[] = ".PUB";

int save_pub_key (PUBKEY *pub)
{
  FILE  *save;
  char  *cpy,*src;
  BIGINT  px, qx, ax, qbit, qxinv;
  char  filename[MAX_NAME_SIZE+5];
  INDEX  i,j;

/*  create key name file */

  src = pub->name;
  while (*src == ' ' || *src == '\t') src++;
  cpy = filename;
  while (*src != ' ' && *src != '\t') *cpy++ = *src++;
  *cpy = '\0';
  strcat (filename, extensn);
  if (! (save = fopen(filename, "w"))) {
     printf("can't create %s\n",filename);
     return(-1);
  }

/*  use last bit of y to define all of it.  Store in msb of first ELEMENT used.
  So far as I can tell, this is valid for type 1 normal basis (i.e.
  there are no 2^m+1 for m congruent to 5 valid field primes.)
  Point p is always y[0], so no need to encode it.  For q, compute y/x
  to determine which y to use from quadradic solution.  See Menezes, pg 92.
*/

  copy(&pub->p.x, &px);
  copy(&pub->q.x, &qx);
  copy(&pub->crv.a6, &ax);
  opt_inv(&pub->q.x, &qxinv);
  opt_mul(&pub->q.y, &qxinv, &qbit);
  if (1 & qbit.b[LONGPOS]) qx.b[STRTPOS] |= SUBMASK;

/*  may want to add the following some time if a2 ever used, can be attached after
  a6 and full curve takes minimal storage:
  if (pub->crv.form) ax.b[STRTPOS] |= SUBMASK;
*/

  fprintf(save, "%s\n", pub->name);
  fprintf(save, "%s\n", pub->address);
  SUMLOOP(i) fprintf(save, "%lx ", px.b[i]);
  fprintf(save, "\n");
  SUMLOOP(i) fprintf(save, "%lx ", qx.b[i]);
  fprintf(save, "\n");
  SUMLOOP(i) fprintf(save, "%lx ", ax.b[i]);
  fprintf(save, "\n");
/*  output a2 here  */
  fclose(save);
}

/*  Recover a public key from disk file.  Assume ascii format of save_pub_key.
  Enter with file name (with or without extension) and pointer to a
  storage block.  Returns 0 on success, -1 on failure.
*/

int restore_pub_key( char *name, PUBKEY *pub)
{
  FILE  *restore;
  char  filename[MAX_NAME_SIZE+5];
  INDEX  i,j;
  BIGINT  px, qx, ax, y[2], f;

/*  check for extension on file name and open file  */

  if ( !( j = (INDEX)strlen(name))) {
     printf("filename too long: %s\n",name);
     return(-1);
  }
  strcpy(filename,name);
  if (strcmp(&name[j-4], extensn) && j<=MAX_NAME_SIZE) 
     strcat(filename,extensn);
  if ( !(restore = fopen(filename, "r"))) {
     printf("can't open file %s\n");
     return(-1);
  }

/*  read in raw data  */

  null( &px);
  null( &qx);
  null( &ax);
  fgets(pub->name, (size_t) MAX_NAME_SIZE, restore);
  pub->name[strlen( pub->name) - 1] = '\0';
  fgets(pub->address, (size_t)MAX_NAME_SIZE, restore);
  pub->address[strlen( pub->address) - 1] = '\0';
  SUMLOOP(i) fscanf(restore, "%lx", &px.b[i]);
  SUMLOOP(i) fscanf(restore, "%lx", &qx.b[i]);
  SUMLOOP(i) fscanf(restore, "%lx", &ax.b[i]);
/*  read in a2 here  */
  fclose(restore);

/*  create curve parameters.  fix this if form == 1 ever used.  */

  null(&pub->crv.a2);
  pub->crv.form = 0;
  copy(&ax, &pub->crv.a6);

  copy( &px, &pub->p.x);
  fofx( &px, &pub->crv, &f);
  if (gf_quadradic( &px, &f, &y[0])) {
     printf("Key in file %s does not have valid point on given curve.\n",
      name);
     printf("This is a major malfunction.\n");
     return(-1);
  }
  copy (&y[0], &pub->p.y);
  one( &pub->p.z);

/*  get last bit of y/x for subscript into quadradic results  */

  if (qx.b[STRTPOS] & SUBMASK) {
     j = 1;
     qx.b[STRTPOS] &= UPRMASK;
  } else
     j = 0;
  copy( &qx, &pub->q.x);
  fofx( &qx, &pub->crv, &f);
  if (gf_quadradic( &qx, &f, &y[0])) {
     printf("Key in file %s does not have valid point on given curve.\n",
      name);
     printf("This is a major malfunction.\n");
     return(-1);
  }
  copy (&y[j], &pub->q.y);
  one( &pub->q.z);
  return(0);
}

void print_pubkey(PUBKEY *pk)
{

  printf("name: %s\n",pk->name);
  printf("address: %s\n",pk->address);
  printf("public key is:\n");
  big_print("curve: ",&pk->crv.a6);
  print_point("  P:",&pk->p);
  print_point("  Q:",&pk->q);
}
/*main()
{
  char  file[256];
  BIGINT  secret_key;
  PUBKEY  public_key, recovered_key;

  init_opt_math();
  init_rand();

  public_key_gen(&secret_key, &public_key, 1);
  big_print("secret key: ",&secret_key);
  print_pubkey(&public_key);
  save_pub_key(&public_key); 
  printf("input file name for previous key: ");
  scanf("%s",&file);
  restore_pub_key(file, &recovered_key);
  print_pubkey(&recovered_key);
  close_rand();
}*/
/****** krypto_knot.c *****/
/********************************************************************************
*                    *
*  The purpose of this file is to tie together all the subroutines from  *
*  the elliptic curve optimal normal basis suite to create a key hiding system.  *
*  The actual compression and encoding should be performed by standard     *
*  algorithms such as arith-n (compression book) and blowfish (Dr. Dobbs's  *
*  Journal or Applied Cryptography.  The complete system would then violate ITAR*
*  and not be accessable, so you'll have to hack that yourself.      *
*                    *
********************************************************************************/

#include <stdio.h>
#include "bigint.h"
#include "eliptic.h"
#include "eliptic_keys.h"

extern void null(BIGINT*);
extern void copy(BIGINT*, BIGINT*);
extern void fofx(BIGINT*, CURVE*, BIGINT*);
extern int gf_quadradic(BIGINT*, BIGINT*, BIGINT*);
extern void one(BIGINT*);
extern void Mother(unsigned long*);
extern void esum(POINT*, POINT*, POINT*, CURVE*);
extern void esub(POINT*, POINT*, POINT*, CURVE*);
extern void elptic_mul(BIGINT*, POINT*, POINT*, CURVE*);
extern void flaten(POINT*);
extern void public_key_gen(BIGINT*, PUBKEY*, INDEX);
extern void restore_pub_key( char*, PUBKEY*);
extern void print_pubkey( PUBKEY*);
extern void big_print(char*, BIGINT*);

extern unsigned long random_seed;

/*   encrypt a session key.  Enter with given session key to hide, public key
to hide it in, and storage block for result.  It is a waste of space to use a
PUBKEY for this, too bad.

session == pointer to key to be hidden
pk == pointer to public key block to use
ek == pointer to resultant key block.  ek->p holds kP, ek->q holds S+kQ
*/

void elptic_encrypt(BIGINT *session, PUBKEY *pk, PUBKEY *ek)
{
  INDEX  i;
  BIGINT  k, f, y[2];
  POINT  s, t;

/*  encode session key onto a random point using this public key.  */

  null(&k);
  Mother(&random_seed);
  copy(session, &k);
  k.b[STRTPOS] = random_seed & UPRMASK;

/*  note that this assumes session key < NUMBITS and that STRTPOS ELEMENT is free
to be clobbered.  For all reasonable encoding schemes this shouldn't be a problem.
*/

  fofx( &k, &pk->crv, &f);
  while (gf_quadradic(&k, &f, y)) {
     k.b[STRTPOS]++;
     fofx(&k, &pk->crv, &f);
  }
  copy( &k, &s.x);
  copy( &y[1],&s.y);  /*  use 1 just to be different, why not eh?  */
  one(&s.z);

/*  next generate a random multiplier k  */

  null(&k);
  SUMLOOP(i) {
     Mother(&random_seed);
     k.b[i] = random_seed;
  }
  k.b[STRTPOS] &= UPRMASK;

/*  do 2 multiplies, kp and kq  */

  elptic_mul(&k, &pk->p, &ek->p, &pk->crv);
  flaten( &ek->p);
  elptic_mul(&k, &pk->q, &t, &pk->crv);

/*  add s to kQ as final step  */

  esum(&s, &t, &ek->q, &pk->crv);
  flaten(&ek->q);
  ek->crv.form = pk->crv.form;
  copy( &pk->crv.a2, &ek->crv.a2);
  copy( &pk->crv.a6, &ek->crv.a6);
}

/*  decrypt session key from public and encrypted key.
  returns 0 if successful, -1 on failure (wrong pass phrase).
*/

int elptic_decrypt(BIGINT *session, PUBKEY *pk, PUBKEY *ek)
{
  INDEX  i;
  BIGINT  skey;
  POINT  check, t, s;

/*  first ensure you can generate secret key.  */

  public_key_gen(&skey, pk, 0);
  elptic_mul(&skey, &pk->p, &check, &pk->crv);
  flaten(&check);
  SUMLOOP(i) {
     if (check.x.b[i] != pk->q.x.b[i]) {
        printf("Invalid pass phrase.\n");
        return(-1);
     }
  }

/*  next compute T = aR and subtract from R' to get S  */

  elptic_mul(&skey, &ek->p, &t, &pk->crv);
  esub(&ek->q, &t, &s, &pk->crv);
  flaten(&s);

/*  clear out encoding garbage and return session key */

  copy( &s.x, session);
  session->b[STRTPOS] = 0;
  return(0);
}

main()
{
        char    file[256];
  BIGINT  session_key, recovered_key, secret_key;
  PUBKEY  public_key, hidden_key;
  INDEX  i;

        init_opt_math();
        init_rand();

/*  for this test, assume you have already generated a public key (and remember
the pass phrase!).  Read in key file from disk.  */

/*        printf("input file name for previous key: ");
        scanf("%s",&file);
        restore_pub_key(file, &public_key);
*/
  public_key_gen(&secret_key, &public_key, 1);
  print_pubkey(&public_key);
  big_print("secret key: ", &secret_key);

/*  generate a session key...  */

  printf("Generating session key and encryping it...\n");
  SUMLOOP(i) {
     Mother(&random_seed);
     session_key.b[i] = random_seed;
  }
  elptic_encrypt(&session_key, &public_key, &hidden_key);
  print_pubkey(&hidden_key);
  big_print("session_key: ",&session_key);
  elptic_decrypt(&recovered_key, &public_key, &hidden_key);
  big_print("recovered key: ",&recovered_key);

  close_rand();
}

#  simple makefile for elliptic curve stuff

test: eliptic.o bigint.o support.o krypto_knot.o
  ld -o eliptic /lib/crt0.o krypto_knot.o support.o eliptic.o bigint.o -lc

krypto_knot.o:  eliptic.h bigint.h eliptic_keys.h krypto_knot.o
  cc -c -g krypto_knot.c

support.o:  eliptic_keys.h eliptic.h bigint.h support.c
#  cc -c -O support.c
  cc -c -g support.c

eliptic.o:  eliptic.h bigint.h eliptic.c
#  cc -c -O eliptic.c
  cc -c -g eliptic.c

bigint.o:  bigint.h bigint.c
  cc -c -O bigint.c

/***** optprime.dat *****/
The following primes are valid for type 1 optimal normal
basis math.  They have the property that 2 is a generator
among other things.

    2     3    11    13    29    53    59    83   107   149
  173   179   227   269   293   317   347   389   467   509
  557   563   587   653   773   797  1019  1109  1187  1229
 1283  1307  1493  1523  1619  1637  1733  1907  1949  1997
 2027  2099  2309  2459  2477  2579  2693  2819  2837  2909
 2957  2963  3203  3413  3467  3533  3677  3779  3803  3947
 3989  4133  4139  4157  4253  4259  4283  4349  4373  4493
 4517  4547  4787  5099  5189  5309  5387  5483  5507  5693
 5717  5813  5939  6173  6197  6269  6317  6389  6653  6659
 6779  6827  6899  7013  7109  7187  7517  7523  7643  7949
 8069  8117  8147  8573  8699  8747  8819  8963  9173  9467
 9533  9587  9749 10163 10589 10667 10709 10733 10853 10883
11003 11069 11213 11483 11549 11699 11813 12107 12149 12197
12203 12227 12269 12347 12437 12539 12653 12659 12899 12917
13037 13043 13163 13229 13523 13829 13877 13997 14243 14387
14549 14699 14867 14957 15077 15083 15173 15299 15413 15629
15683 15773 15803 16139 16187 16229 16547 17027 17093 17189
17387 17483 17789 17939 18059 18077 18269 18413 18443 18587
18917 18947 19037 19157 19259 19379 19949 19973 19997 20123
20477 20507 20627 20717 20789 21059 21179 21227 21419 21467
21563 21773 22013 22109 22229 22277 22613 22637 22643 22739
22787 22973 23099 23117 23357 23603 23627 23813 23819 24029
24083 24203 24317 24533 24659 24917 24989 25307 25349 25373
25469 25589 25643 26003 26099 26189 26309 26459 26627 26693
26717 26813 27107 27299 27803 27827 28019 28109 28163 28229
28277 28307 28499 28517 28643 28949 28979 29123 29243 29333
29339 29483 29573 29669 29837 30029 30197 30203 30293 30323
30347 30467 30539 30557 30677 30803 30869 30893 31013 31139
31259 31469 31517 31547 31973 32003 32069 32213 32237 32507
32603 32717 32843 32933 32987 33053 33107 33149 33317 33347
33413 33773 34157 34253 34589 34667 34757 34877 34949 35117
35339 35363 35573 35963 36083 36269 36299 36467 36629 36683
36749 36923 37277 37397 37547 37589