//
   //
   //  Magic Flight: A Public Key Exchange based on a
   //  Lossy Commutative Mixer
   //
   //
   //    by Jim Steuert
   //    August 2, 2003
   //
   // (c) Copyright, 2003 by Jim Steuert
   //
   //
   // Warning:
   // This algorithm for public key is experimental.
   // It is based on feedback from the sci.crypt newsgroup,
   // including Scott Fluhrer, Bryan Olson, and Francois Grieu.
   //
   //
   // Abstract
   //
   //
   //  Magic Flight is an algorithm for public key
   //  exchange which is based on a lossy commutative
   //  mixer constructed from products of primes.
   //
   //  The prior version of Magic Flight used a
   //  commutative and distributive semigroup algebra 
   //  based on a number-line analog using an offset and modulus.
   //  
   //  This was broken by Bryan Olson using a recursive attack.
   //
   //  This specific version uses a ring based on addition
   //  which is distributive over the max operation.
   //  I found this as part of the "max-plus" algebraic
   //  theories (references 10 and 11).
   //
   //  The basic operation involves an underlying
   //  ring-like algebra which does not have
   //  multiplicative inverses, and may not have
   //  additive inverses. The elimination of additive
   //  inverses (cancellation) is intended to prevent
   //  standard algebraic methods of solution.
   //
   //  A practical variant of the original program exists
   //  which uses Karatsuba multiplication method
   //  to implement the commutative mixer.
   //  A matrix ring operation with (a,b,kb,a) type
   //  2x2 matrix operations is used to assure that
   //  some fraction (1/4 for 4-bit a and b values)
   //  of the values have no multiplicative inverses.
   //
   //
   //
   //  Add distributes over max, so max is used as
   //  multiply would normally be used, and add is
   //  used as would multiply.
   // 
   //  a simple working program. It is hoped that
   //  this meets the goals of preventing Gaussian
   //  Elimination or other algebraic solution.
   //
   //  Related work (4) by  Chris Monico, or (3) by
   //  Maze, Monico, and Rosenthal, use semirings but
   //  not this specific mixer operation. It does not
   //  implement this particular semiring rng operation.
   //
   //
   //
   //
   //
   // Introduction.
   //
   //
   // The basic operation of this commutative mixer algorithm
   // is representing the secret and digest in a
   // "product of primes" form as shown:
   //
   //  secret = a0 + a1*p0 + a2*p1 + a3*p0p1
   //    and
   //  mixer  = m0 + m1*p0 + m2*p1 + m3*p0p1
   //
   //  where p1 and p2 are not primes, but rather
   //  the square roots of distinct primes.
   //
   //  So "multiplying" the secret times the mixer
   //  gives the output convolution product form:
   //
   //  x0 + x1*p0 + x2*p1 + x3*p0p1
   //
   //  where 
   //    x0 = a0*m0 + p0^2*a1*m1 + p1^2*a2*m2 + p1^2*p2^2*a3*m3
   //    x1 = a0*m1 + a1*m0 + p1^2*a2*m3 + p1^2*a3*m2
   //    x2 = a0*m2 + a2*m0 + p0^2*a1*m3 + p0^2*a3*m1
   //    x3 = a0*m3 + a3*m0 + a1*m2 + a2*m1
   //
   // Note that this is easily generalized to multiple 
   // primes.
   //
   // The only requirement is that the underlying
   // addition and multiplication operations be commutative
   // and associative, and that the multiplication be
   // distributive over addition.
   //
   // There is NO requirement for multiplicative or
   // additive inverses. However, conventional
   // Gaussian Elimination (and other algebraic solutions)
   // require additive inverses or cancellation. Bryan
   // Olson showed how multiplicative inverses are not
   // always needed in Gaussian Elimination pivoting.
   // But additive inverses clearly are, as cancellation
   // is the fundamental operation.
   //
   //
   //
   // 1) Lossy Non-Invertible Multiplication
   // 
   //
   //  The smaller modulo 2^k widths make for situations
   //  in which there is no multiplicative inverse.
   //  That is what I mean by lossy. This makes it
   //  less likely that x can be divided-off from a public key.
   //  Making the mask too small (3 bits) reduces the
   //  entropy of the output dramatically.
   //
   //  A practical variant of the original algorithm, with
   //  some of the ring values not having multiplicative inverses,
   //  involves using 2^20 terms. Gaussian Elimination, as
   //  suggested by Francois Grieu with modifications by
   //  Bryan Olson, would then take n^3 or 2^60 work.
   //  Rather than taking n^2 work to compute this, I have
   //  demonstrated that the mixer can be optimized like
   //  any multiplication, using Karatsuba, Toom-Cook 3-Way, or
   //  FFT optimizations. Karatsuba multiplication takes
   //  n^1.585 rather than n^2 operations to compute,
   //  so a version with 2^20 terms is practical to compute,
   //  but still requires 2^60 to break.
   //
   //  One practical program: matrix.c prevents 
   //  multiplicative inverse for 1/4 of its values, when
   //  the a or b is 4-bits wide. It is based on commutative
   //  matrices of the following form:
   //
   //   | a  b |     |  c d |
   //   |      | +/* |      |
   //   |kb  a |     | kd c |
   //
   //
   //
   //  Details of Karatsuba implementation:
   //
   //  If one looks at the "product of primes" construction as
   //  a number system, then it looks like a vector multiply.
   //
   //    x-> = (a3*p1p0 + a2*p1 + a1*p0 + a0 )
   //    y-> = (b3*p1p0 + b2*p1 + b1*p0 + b0 )
   //
   //  where p1,p2 really means the square root of p1,p2, respectively.
   //  This can also be represented as the following:
   //
   //    x-> =  (a3*p0+a2)*p1 + (a1*p0+a0)
   //    y-> =  (b3*p0+b2)*p1 + (b1*p0+a0)
   //
   //  and so on recursively. The (a1*p0 + a0) can itself be such
   //  that a1 and a0 are themselves of the form of a1 = (m1*pk + m0),
   //  where pk is another prime square root, and so on recursively.
   //
   //    x-> =  u1*p1 * u0
   //    y-> =  v1*p1 + v0
   //
   //  Karatsuba can be applied, or Toom-Cook 3-way, or FFT,
   //
   //  For example, for Karatsuba:
   //
   //    x->*y-> = (u1v1)*p1^2 + ((u1+v1)*(u0+v0)-u1v1-u0v0)*p1 + u0v0
   //    x->*y-> = ((u1+v1)*(u0+v0)-u1v1-u0v0)*p1 + (u0v0+(u1v1*p1^2))
   //
   //  Note that the p1^2 term causes another multiply
   //  (since p1 was a square root) , but it is linear (O(n)) in
   //  the number of the included terms.
   //  
   //  Note that there are only three full-size multiplies (which
   //  are each of half the size (number of terms), not 4, thus 3*O((n/2)^2),
   //  and 2 adds, 2 subtracts (but those are O(n)).
   //  So the result is O(1.585) for Karatsuba.
   //  By expanding the multiplicand terms (points) into 8 or 16 points
   //  instead of two, Toom-Cook 3-way or FFT can be used
   //  for increased performance (at greatly increased complexity).
   //
   //  Karatsuba works only because it works recursively.
   //  And it depends on there being a subtract, or additive 
   //  inverse.
   //
   //  However, the goal of Magic Flight requires that
   //  Gaussian Elimination be prevented completely, and
   //  that cancellation be eliminated, meaning no additive inverses.
   //
   //  That is why the following was implemented.
   //
   //
   //
   // 2) Lossy Non-Invertible Addition
   //
   //
   //  The underlying semigroup rng (ring without an inverse)
   //  has commutative and associative addition without an
   //  inverse, and has multiplication which is comutative,
   //  associative, and distributive over addition.
   //
   //  After extensive research into semigroup algebras,
   //  which is a vast and fascinating subject, some
   //  implementation considerations were arrived at.
   //
   //
   //    a) Distributive Lattices.
   //
   //     The simplest version of this is a boolean algebra,
   //     where addition is "or", and multiplication is "and".
   //     The simplest implementation is bitwise and and or,
   //     which, although simple, would lend itself to
   //     one-per-bit sets of equations. It is not clear that
   //     this is breakable, but there is a lot of literature
   //     on solving sets of boolean equations. Rivest has shown
   //     how to embed boolean algebra into integer linear matrix
   //     equations.
   //
   //     Other distribute lattices based on subsets are possible,
   //     but are complicated to generate , and those based
   //     on subsets can be embedded in the powerset version.
   //     In general, distributive lattices have join (+) and meet
   //     or intersection (*) operations which do not have
   //     additive inverses.
   //
   //
   //
   //    b) "Stateful" Commutative/Associative Rng Algebras.
   //       - Broken by Bryan Olson
   //
   //     Addition in this algebra starts with counting from
   //     0 to the start of the cycle, but then repeats at the
   //     cycle. Multiplication has both the 0 and a 1. Only
   //     when used in a ring with non-invertible multiplication
   //     can multiplicative inverses also be prevented.
   //
   //     As an example, addition could be based on the sequence
   //
   //       0 1 2 3 4 5 6 7 8 9 10 11 12 13 14  5  ...
   //
   //     where the cycle starts at 6, for modulo 9 cycle
   //      (6 7 8 9 10 11 12 13 14). So the offset is 6 and the
   //     modulo is 9.
   //
   //     This was broken by Bryan Olson using a recursive
   //     attack.
   //
   //
   //    c) Semigroup algebras.
   //
   //     There are many methods of implementing semigroups, 
   //     including finite-state-machines, but I could find
   //     scant references for simply generating semigroup
   //     algebras, and none on generating algebras without
   //     additive inverse. The best is "Handbook of Finite
   //     Semigroup Programs" by John A. Hildebrant. This
   //     101-page paper includes fortran programs for 
   //     generating semigroups and semigroup algebras.
   //     However, the semigroup algebras are based on
   //     semigroups over Zp, which has an additive inverse.
   //
   //     I have converted a program (see ref 2) from fortran to c
   //     for generating semigroups, but it only generates
   //     semigroups of orders up to 7 in reasonable time.
   //     This program uses a clever method to eliminate
   //     generation of isomorphic tables, by creating a
   //     table of n-factorial permutations, and lexicographically
   //     generating the semigroups. I have posted it,
   //     named gensemi.c, at <http://users.rcn.com/pjsteuert/gensemi.c>,
   //     I enhanced the program to also generate those
   //     commutative semigroups which are distributive
   //     over other semigroups of the same order. It is still
   //     performance-limited to small (order 6) semigroups.
   //
   //     It could be used for generating semigroup algebras which are
   //     direct products of the smaller ones, but that is
   //     equivalent to using both in parallel, and doesn't
   //     seem to enhance security.
   //
   //     Other methods of generating semigroup algebras of
   //     higher order are being investigated.
   //
   //     However, I stumbled across and interesting
   //     semigroup which is distributive across another
   //     one.
   //
   //
   //
   //   d) Max-plus system.
   //   
   //     The max-plus system is based on work by Stephane Gaubert (10),
   //     with reference by Jeremy Gunanadera (11).
   //
   //     In the max-plus system, add is distributive over
   //     max. So in algebraic notation:
   //
   //         a + ( b max c )  ==   (a+b) max (a+c)
   //  
   //     In this specific usage, the max of the product of
   //     primes (mod 2^8) and the data bytes is then added,
   //     to form the terms. So the "multiply" is really 
   //     mod 2^8 addition, and "add" is really max.
   //
   //     Preliminary analysis shows that this should resist
   //     all linear attacks, as there is no cancellative
   //     inverse to max.
   //
   //
   //
   // 3) Implementation with 2^n terms.
   //
   //   
   //    To make this commutative, the multi-part digest is represented
   //    by the  n terms of the form:
   //
   //     a0 + a1*sqr(r) + a2*sqr(s) + a3*sqr(t) + a4*sqr(rs) + ...
   //
   //    Of course collecting terms gets tedious very fast.
   //    In this scheme, each part of the vector is represented
   //    by an index. For example, index 1001 would
   //    represent the term associated with primes p0p3,
   //    while 1110 would be the term for primes p0p1p2.
   //
   //    For example,  the  1001 term is  a[9]*sqr(p0p3)
   //             and  the  1110 term is  x[14]*sqr(p0p1p2).
   //
   //    Multiplication is a nested loop over the terms,
   //       with xor to index the result term, and "and"
   //       to find the scale factor.
   //     This example has a[9]*x[14]*p0*sqr(p1p2p3).
   //     The xor is 0111 or p1p2p3, and the and is 1000 or p0.
   //
   //
   //  Bibliography:
   //
   //  1. Feedback by Scott Fluhrer, Bryan Olson, and Francois Grieu
   //     on the sci.crypt newsgroup.
   //
   //  2. Handbook of Finite Semigroup Programs
   //      by John A Hildebrant
   //
   //  3. Public Key Cryptography Based on Simple Modules over
   //     Simple Rings, June 20, 2002
   //       by Gerard Maze, Chris Monico, and Joachim Rosenthal
   //       Department of Mathematics, University of Notre Dame
   //
   //  4. Semirings and Semigroup Actions in Public Key Cryptography
   //     A Dissertation Submitted to Graduate School at University
   //     of Notre Dame, April 2002
   //       by Christopher J. Monico
   //
   //  5. Sets, Lattices, and Boolean Algebras, by James C. Abbott
   //
   //  6. The Structure of Abstract Algebra, by Ralph B. Crouch and
   //     David N. Beckman
   //
   //  7. An Introduction to Algebraic Structures, by Azriel Rosenfeld
   //
   //  8. Introduction to Algorithms, by Cormen, Leiserson, and Rivest
   //     First Edition, section 31.3 (page 745 section on Quasirings)
   //
   //  9. Powers of Complex Associative Functions by Jim Steuert
   //      sci.crypt newsgroup Feb 25, 2001
   //
   //  10. Une introductire a l'a;gebre (max,+) by Stephane Gaubert
   //      Seminaire IRISA 19 Sep 1997
   //
   //  11. From max-plus algebra to nonexpansive mappings: a nonlinear
   //      theory for discrete event simulation, by Jeremy Gunanadera
   //      HP Labs, Bristol UK March 7, 1999
   //
   //
#define NUMPRIMES  10
unsigned long prime[30] = {7,11,13,19,17,23,29,37,43,59,67,71,83,89,97};
#define NUMTERMS  1024
#define NUMTERMS2 32
#define INPUTMASK  0xff
//
// capacity should be 4096*256 or 2^20
// loss of entropy is due to max, not masking mod 2^n
//
#define FIELDMASK  0xfffff

extern unsigned long xseed;


unsigned long primeprods[NUMTERMS];
unsigned long digest[ NUMTERMS ], result[ NUMTERMS ];
unsigned long common_global[ NUMTERMS ];
unsigned long alices_secret[ NUMTERMS ], bobs_secret[ NUMTERMS ];
unsigned long alices_public[ NUMTERMS ], bobs_public[ NUMTERMS ];
unsigned long alices_shared[ NUMTERMS ], bobs_shared[ NUMTERMS ];

/*================================================*/


unsigned long add( unsigned long x, unsigned long y)
{
unsigned long rslt;

x = x & 0xfffff;
y = y & 0xfffff;

if ( x > y )
  rslt = x;
else
  rslt = y;

return( rslt );
}

/*================================================*/



unsigned long multiply( unsigned long x, unsigned long y)
{
unsigned long rslt;

x = x & 0xfffff;
y = y & 0xfffff;

rslt = x+y;

rslt = rslt & 0xfffff;

return( rslt );
}


/*================================================*/


void init_rings()
{
unsigned long prod,e,f;
int i,j;

  //
  // build table of products of primes
  //
  for ( i = 0; i < NUMTERMS; i++ )
  {
    prod =  0; // unit 1 in multiply for maxplus system

    for ( j = 0; j < NUMPRIMES; j++)
    {
      if ( ( (0x1<<j) & i) != 0 )
	  {
         prod = multiply( prod, prime[j] );
	  }
    }
    primeprods[i] = prod & 0xfffff;
    //printf("primeprods[%08x] = 0x%04x\n",i,primeprods[i]);
  }
}


/*================================================*/

void domixer( char *title, unsigned long *start,
   unsigned long *secret, unsigned long *result )
{
unsigned long a,x,term,and,pfac,rslt,mask,mrslt,e,f;

mask = FIELDMASK;

for ( x = 0; x < NUMTERMS; x++)
  digest[ x ] = start[ x ];

for ( x = 0; x < NUMTERMS; x++)
  result[ x ] = 0; // 0 is additive unity here



for ( a = 0; a < NUMTERMS; a++)
{
  for ( x = 0; x < NUMTERMS; x++)
  {
      term = (a ^ x) & FIELDMASK;
      and = (a & x) & FIELDMASK;
      pfac = primeprods[and] & FIELDMASK;
      e = secret[a] & FIELDMASK;
      f = digest[x] & FIELDMASK;
	  mrslt = multiply( e, f );
      mrslt = multiply( mrslt, pfac );
	  rslt = result[ term ];
      rslt = add( rslt , mrslt );
      rslt &= FIELDMASK;

      result[ term ] = rslt;
  }  // x
}  // a

for ( a = 0; a < NUMTERMS; a++)
   digest[ a ] = result[ a ];

printf("%s:\n",title);
for ( a = 0; a < NUMTERMS2; a++)
  printf("%03x ",result[a]);
printf("\n\n\n");

return;
}
/*================================================*/
int main( int argc, char **argv )
{
int i,j,fld;
if ( argc < 2 )
{
   printf("Usage: maxplus <seed1>\n");  
   exit(1);
}
xseed = strtoul(argv[1],0,0);
r250_init( xseed );

init_rings();

for ( fld=0; fld < NUMTERMS; fld++)
   common_global[ fld ] = r250n(60000) & INPUTMASK;

for ( fld=0; fld < NUMTERMS; fld++)
   alices_secret[ fld ] = r250n(60000) & INPUTMASK;

for ( fld=0; fld < NUMTERMS; fld++)
   bobs_secret[ fld ] = r250n(60000) & INPUTMASK;

 domixer( "alice's public key", common_global, alices_secret, alices_public );
 domixer( "bob's public key", common_global, bobs_secret, bobs_public );
 domixer( "alice's shared key", alices_public, bobs_secret, bobs_shared );
 domixer( "bob's shared key", bobs_public, alices_secret, alices_shared );
}