working code Elliptic curve cryptography

at last...
now its there. done. It only took like 4 year :)

Go check it out - contains all sourcecode for a working ECC solution! - look for the Master Thesis for source
I've used Menezes Vanstone ECC for cleartext encoding.
Also available: general application for calculation on any curve - any size!

binary expansion in C# - converting BigInteger to binary form

This source code will do the job.
Add appropriate controls and fire the number: 174050332293622031404857552280219410364023488927386650641
and you'll get:
01110001100100101011100101011111111111001000110110100111100001100011000100000001000111101101011010110010010011001101110101010111001111111001011101111010000100011110011110010100100000010001
with a bunch of more noughts in front.
What a happy day. :)

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using System.Numerics;

namespace binaryexpansion
{
    public partial class Form1 : Form
    {   
 
        public Form1()
        {
            InitializeComponent();
        }

        private void button1_Click(object sender, EventArgs e)
        {
            BigInteger t = BigInteger.Parse(inputnr.Text); 
            selResult.Text = calcBinary(t);

            //the line below would work if t was an integer - so nice and easy, but no dice with BigIntegers
            //string t2 = Convert.ToString(t, 2);
             
        }

      
        private string calcBinary(BigInteger tee) {

            BigInteger t = tee;
            //4 times the length because it takes about 3 times the number of digits to write a number i binary form
            int len = 4*(inputnr.Text).Length;
            
            string fullnum= " ";
            for (int i = len; i >= 0; i--) 
            {
                int inum;
                
                BigInteger num =  (t >> i)&1;
                if(num == 1){
                //ugly but easy
                inum = 1;}
                else if(num == 0){
                inum = 0;}
                else{
                inum=666;
                }
                fullnum =  fullnum  + Convert.ToString(inum);
                
            }
            return fullnum;

            }
        }
    }

Huge cryptographically secure random numbers

Just a tiny implementation example on how to use .NET crypto-package to create really huge pseudo random numbers.
Add approprite controls and rock.
(in my example I just kicked out the number on a webpage - feel free to do it where ever you like)

using System;
using System.Collections.Generic;
using System.Linq;
using System.Web;
using System.Web.UI;
using System.Web.UI.WebControls;
using System.Numerics;
using System.Security.Cryptography;

namespace numbers
{
    public partial class WebForm1 : System.Web.UI.Page
    {
        protected void Page_Load(object sender, EventArgs e)
        {

        }

        private static BigInteger GetRandom()
        {
            RandomNumberGenerator ran = RNGCryptoServiceProvider.Create();
            byte[] ranbytes = new Byte[256];
            //ex: Byte[20] gives 8 bit * 20 = 160 bit ~  50 digits
            //ex: Byte[2] gives 8 bit * 2 = 16 bit  = +- 32768
            //ex: Byte[128] gives 8 bit * 128 = 1024 bit ~ 300 digits
    
            //fill the byte array with cryptographically strong random bytes
            ran.GetBytes(ranbytes);

            //put the random bytes into a bigInteger
            BigInteger mynum = new BigInteger(ranbytes);
            return mynum;
        }

        protected void Button1_Click(object sender, EventArgs e)
            {

            //Roll the really really big dice :)
            BigInteger roll = GetRandom();
            result.Text= roll.ToString();
 
        }
    }
}

This will work! P + Q Point addition in C#.NET code

Add appropriate controls and it will work to calculate P + Q, addition of two points on an elliptic curve over a finite field Fp.
And it all works with big Integers.
Be sure that you only use valid numbers (prime for Fp and points that acctually ARE on the curve)
Next version will contain valid curves generated from startpoint and startvalues , plus it will contain an implementation of the double and add algorithm.

Sorry for the horrible buggy primitive spaghetti code without validation and exception handling for now - but I just made it and it works! - tested against a huge amount of numeric examples calculated by hand + toy examples provided by Hoffman, Pipher and Silverman.

It will not handle input of point at infinity. But you already know that P + Ø = P so... here we go.

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using System.Numerics;

namespace WindowsFormsApplication2
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }

        private void Form1_Load(object sender, EventArgs e)
        {

        }

        private void button2_Click(object sender, EventArgs e)
        {
            calcPoints();

        }


        private BigInteger calcLambda1(BigInteger y2, BigInteger y1, BigInteger x2, BigInteger x1, BigInteger modspace) {
            
            BigInteger lam1 = (y2 - y1) * CalcInverse((x2 - x1), modspace);
           
            return lam1;
        }


        private  BigInteger calcLambda2(BigInteger x1, BigInteger y1, BigInteger A, BigInteger modspace)
        {
           
            BigInteger lam2= (3 * x1 * x1 + A)*CalcInverse((2 * y1), modspace);
            return lam2;
        }


        private  BigInteger calcXFinal(BigInteger Lambda, BigInteger x1, BigInteger x2, BigInteger modspace) {
            BigInteger bigX= (Lambda*Lambda - x1 - x2)%modspace;
            if (bigX < 0) { bigX = bigX + modspace; }
            return bigX;
        
        }
        private BigInteger calcYFinal(BigInteger Lambda, BigInteger x1, BigInteger x3, BigInteger y1,  BigInteger modspace)
        {
            BigInteger bigY= (Lambda*(x1-x3) - y1)%modspace;
            if (bigY < 0) { bigY = bigY + modspace; }
            return bigY;
        }


        private  BigInteger CalcInverse(BigInteger n,  BigInteger p)
        {

            BigInteger x = 1;
            BigInteger y = 0;
            BigInteger a = p;
            BigInteger b = n;
            BigInteger q, t;
            BigInteger res;
            while (b != 0)
            {
                t = b;
                //q2 = Math.floor(a/t);
                q = BigInteger.Divide(a, t);
                b = a - q * t;
                a = t;
                t = x;
                x = y - q * t;
                y = t;
            }
            if (y < 0)
            {
                res = y + p;
            }
            else
            {
                res = y;
            }
            return res;

        }

        private void calcPoints() {

            BigInteger inputA = BigInteger.Parse(paramA.Text);
            BigInteger inputB = BigInteger.Parse(paramB.Text);

            BigInteger inputx1 = BigInteger.Parse(paramx1.Text);
            BigInteger inputy1 = BigInteger.Parse(paramy1.Text);

            BigInteger inputx2 = BigInteger.Parse(paramx2.Text);
            BigInteger inputy2 = BigInteger.Parse(paramy2.Text);

            BigInteger modp = BigInteger.Parse(paramp.Text);
        
        
            if ((inputx1 == inputx2) && (inputy1 != inputy2) )
            {
                resultx3.Text = "point at infinity";
                resulty3.Text = "point at infinity";
            }
            else {


                if ((inputx1 == inputx2) && (inputy1 == inputy2))
                {

                    BigInteger l2 = (calcLambda2(inputx1,inputy1,inputA, modp))%modp;
                    dontpanic.Text = l2.ToString();

                    BigInteger resl2x3 = (calcXFinal(l2, inputx1, inputx2, modp))%modp;
                    BigInteger resl2y3 = (calcYFinal(l2, inputx1, resl2x3, inputy1, modp))%modp;


                    resultx3.Text = resl2x3.ToString();
                    resulty3.Text = resl2y3.ToString();


                }
                else {

                    BigInteger  l1 = (calcLambda1(inputy2,inputy1, inputx2,inputx1, modp))%modp;

                    dontpanic.Text = l1.ToString();
                    BigInteger resl1x3 = (calcXFinal(l1, inputx1, inputx2, modp))%modp;
                    BigInteger resl1y3 = (calcYFinal(l1, inputx1, resl1x3, inputy1, modp))%modp;

                    resultx3.Text = resl1x3.ToString();
                    resulty3.Text = resl1y3.ToString();
                
                
                }

        }


       }
}

The BigInteger is a part of the brand new package System.Numerics in .NET 4.0
more about that one here -->
http://msdn.microsoft.com/en-us/library/system.numerics.biginteger.aspx
For some reason the System.Numerics is not automatically added to projects - select menu project/Add Reference../ select tab: .NET and find System. Numerics.

algorithm for modular inverse in Fp

not the most effective maybe - but at least ONE way of doing it:
A binary algorithm:

You need a prime p for Fp and a number a less than p-1 (larger than 0)
The output will be a^-1 mod p

If you need to calculate b/a mod p ( equals ba^-1 mod a) just substitute x1 with b instead of 1 in the pseudo-code below:

1. s=a, t=p
2. x1 = 1, x2 = 0
3. while (s!= 1 && t != 1) do
   3.1. While s even do
           s=s/2
           if (x1 even) x1 = x1/2 else x1 = (x1+p)/2
   3.2. While t even do
           t = t/2
           if (x2 even) x2 = x2/2 else x2 = (x2+p)/2
   3.3. if (s>=t) s = s-t, x1 = x1-x2
        else t= t-s, x2 = x2-x1
4. if (s==1) return x1 mod p else return x2 mod p

twisted curve

Et - read E twist

Et(Fq) : y^2 = x^3+ d^2Ax+ d^3B

Where d is any Quadratic Non Residue mod q

more on Quadratic Non Residue mod q later

Torsion points and bilinear pairing

if you have a point P where mP=O, then you have a point of order m in the group E
If you have a set of points E[m]={P belongs to E:[m]P=O}
these points are points of finite order, they are called torsion points!

E[m] is a subgroup af E

the great "news" ;Bilinear pairing, is that every point P can be written as a linear combination
P = aP1 + bP2 for unique choice of a and b (in Z/mZ)

If m is large its difficult to find a and b. If b=0, then finding a is solving ECDLP!
P=aP1.

Weil - some light please

Weil:

input two points P and Q on E so that Dvs mP= O og mQ=O
output: m'th root of unity
.......................
P,Q with order m, meaning belongs to E[m]
em(P,Q)=(fP(Q+S)/fP(S))/(fQ(P-S)/fQ(-S)), S not being O,P,-P,Q or -Q
fP and fQ are rational functions on E so that
div(fP)= m[P]-m[O] and div(fQ)= m[Q]-m[O]
(the value of em(P,Q) does for some reason not depend on the choice of fP and fQ)

and the nice thing is that em(P,Q)^m = 1
that means that em(P,Q) is an m'th root of unity

more properties:
em(P,P) = 1 which implies that em(P,Q)= em(Q,P)^-1

if em(P,Q) = 1 for all Q then P=O
..............................................................

What does m'th root of unity means

Normally
2nd root of unity means= square x equals 1 x^2 = 1
and cubic for m=3
A primitiv m'th root of unity means is one where there is no solution to
x^k=1 for every k smaller than m

Koblitz curves

Elliptic curve defined over F2 on the form
Y^2 + XY = X^3 + aX^2 +1, with a 1 or 0. Koblitz idea: define the curve in F2, but take points on E with coordinates i F2^k.
An easy way to calculate the order of the curve - i.e. the number of points on E(F2^k) (with a=0) is
#E(F2^k) = 2^k + 1 - ((-1+sqrt(-7))/2)^k - ((-1-sqrt(-7))/2)^k

for k=11
we have a cute tiny one with order 2116

for k=97
it's a little larger:
158456325028528296935114828764

for k=163
the number of points is:
11692013098647223345629473816263631617836683539492

for k=233
the number of points is:
13803492693581127574869511724554051042283763955449008505312348098965372

for k= 283
the number of points is:
15541351137805832567355695254588151253139246935172245297183499990119263\
318817690415492

k=409 gives:
13221119375804971979038306160655420796568093659285624385692975800915228\
45156996764202693033831109832056385466362470925434684

k=571 gives:
7729075046034516689390703781863974688597854659412869997314470502903038\
2845791208490725359140908268473388268512033014058450946998962664692477\
18729686468370014222934741106692

isn't that amazing?

Galois fields - F2^n

Elliptic curves over F2^n:
characteristic: 2 expansion degree: n
Elements in F2^4
0000, 0001, 0010, 0011, 0100, 0101, 0110, 0111, 1000, 1001, 1010, 1011, 1100, 1101, 1110, 1111
characteristic 2, expansion degree 4
Other fields:
Elements in F3^2
00,01,02,10,11,,12,20,21,22
characteristic 3, expansion degree 2
suitable for cryptography:
y^2 + xy = x^3 + ax^2 + b
b cannot be zero, a may be zero.

Supersingular curve not suitable for cryptography:
y^2 + y = x^3 + ax + b
points (x,y) are variables in F2^n and plus operations is the just XOR of x bits with y bits.

...to be tested!

random curve with point

If you just need any short form weierstrass curve at random with a random point the easiest way is to create the point first, then the curve.
suitable for lenstra.
Given a composite you would like to factor.
toy example:
given N= 1289
random point P=(x,y)=(235,156)
random A= 21
calculate B as
B=y^2 -x^3 - A*x mod N
156^2-235^3-21*235 mod 1289 gives B=-157
then you got your random curve:
Y^2=X^3 + 21X-157 with a point P that by construction is on the curve.

the solution for B ofcourse can be plus/minus any multiplication of N. equivalence class: B + Z/ZN
More about lenstra later

Double and Add algorithm

Very much like the fast powering algorithm - binary expansion can be used to calculate Q= nP fast.
ex n=947 = 2^9+2^8+2^7+2^5+2^4+2^1+2^0

which means that 947P = P+2P+16P+32P+128P+256P+512P. Thats 9 doublings and 6 additions, much better than 946 additions....
if we allow subtraction as well (ternary expansion) its even better.

want to calculate nP?

algorithm:

1. input point P from E(Fp) and n>=1
2. set Q=P and R=Ø (point at infinity)
3. loop while n>0{
     if n = 1 (mod2) set R= R+Q
     set Q=2Q and n=n/2 (floor)}
4.return R (which equals nP)

I've now posted working sourcecode available at http://christelbach.com
(source code contained in master thesis)

for instance this:
//double and add algorithm

private BigInteger[] calculatenP(BigInteger x, BigInteger y, BigInteger n, BigInteger A, BigInteger modp)
{
BigInteger[] Q = { x, y };
BigInteger[] R = { 0, 0 };
while (n > 0)
{
   BigInteger t = BigInteger.ModPow(n, 1, 2); 
   if ((t) == 1)
   {
   R = addPoints(Q[0], Q[1], R[0], R[1], modp);
   } 
   Q = doublePoints(Q[0], Q[1], A, modp);
   n = n / 2;
}
BigInteger[] resPoint = { R[0], R[1] };
return resPoint;
}

Look in the pdf file for the methods addPoints() and doublePoints()

square roots modulus

IF a has a square root modulus p,we can solve it like this:

we have x^2= a mod p, then
b= a^((p+1)/4)mod p
will be a solution. It will satisfy b^2=a mod p.
Isn't that a miracle!

Remember ONLY if a has a squareroot mod p...that certainly not allways the case.

more inverse

how could I forget fermat's little..
a^(p-1) = 1 mod p if p does not devide a, and 0 otherwise.
Multiply a^(p-2) with a and we will get:

a^(-1)= a^(p-2) mod p


together with fast powering algorithm this is handy and to some limit efficient!
An alternativ to the euclidian algorithm.

About fast powering algorithm:

toy example

5^11 mod 13
write 11 as a sum of power of 2 (like binary)
11= 1011, which means 11= 5^(2^3+2^1+2^0)= 5^(2^3)*5^2*5. Reduce mod 13 in each step and the numbers will stay low
so the result is 1*12*5= 60 mod 13= 8

In general:
create a table with relevant powers mod p to some limit, The easy part is that each number is just the square of the preceding one.

Multiplicative inverse modulus p

JavaScript snippet that will implement this nicely.

Use approprite html elements that fits the script.

window.onload= initForm;

function initForm(){
document.getElementById("ok").onclick = calcInverse;
}

function calcInverse(){

var n= parseInt(document.getElementById("nn").value);
var p= parseInt(document.getElementById("pp").value);

var x = 1;
var y = 0;
var a=p;
var b=n;
var q,t;
var res;
while (b != 0) {
                t = b;
                q = Math.floor(a/t);
                b = a - q*t;
                a = t;
                t = x;
                x = y - q*t;
                y = t;
        }
if(y<0){
res= y+p;
}  
else{
res=y;
}

document.getElementById("r1").innerHTML="a^-1  er "+res;
return false;

}

The hunt for multiplicative inverses

The inverse of a, a^1=z mod m. Calculation aint easy for large m, this is where all the computing power goes.
a does have an inverse mod m if and only if a*z=1 mod m. Otherwice there is no inverse.
If our finite field is prime, then all members sport this property

We can find the inverse by using euclid. We know that in the Euclidean Algorithm you repeatedly divide the divisor by the remainder until the remainder is 0.
The gcd is then the last non-zero remainder.
We can write
za + mb = 1. This says that za = 1 + (-b)m, which means za=1 (mod m). We dont care about -b in this case of cause. So z is the inverse of a.

The extended version gives us z.
Nice tiny toy example here:
http://www-math.cudenver.edu/~wcherowi/courses/m5410/exeucalg.html

projective space and homogenous coordinates.

Guess I've got the projective space right now. Love the explanation at wiki.
http://en.wikipedia.org/wiki/Projective_space
But its urgent to really understand that and the homogeneous coordinates.
Projective space: really a projection! 3d space projected into a tiny 2d picture. A projection from 3d to 2d.

Quote:'the set of equivalence classes of R3\(0, 0, 0), i.e. 3-space without the origin, where two points P = (x, y, z) and Pˈ = (xˈ, yˈ, zˈ) are equivalent if there is a nonzero real number λ such that P = λ·Pˈ, i.e. x = λxˈ, y = λyˈ, z = λzˈ. The usual way to write an element of the projective plane, i.e. the equivalence class corresponding to an honest point (x, y, z) in R3, is

[x : y : z].

The last formula goes under the name of homogeneous coordinates'

You can homogenize like this: you add z until all parts of the equation have the same degree:

y^2=x^3 + Ax + B becomes y^2z=x^3 + Axz^2+Bz^3

Calculating in projective space gives us a way to work with the point at infinity without just panicking. The point at infinity just becomes a special plane i a 3 dimensional space, where parallel lines meet.
lucky bastards.

change of language

ok, I'll change the main language of this blog.
Today I found an indonesian girl doing exaclty the same as I, trying to focus and trying to study elliptic curves and implement some cryptography using EC. I was so happy she was not writing in whatever language - indonesian probably..

So if you read this blog, you'll have to deal with my not so perfect language skills.

Feel free to ask me translate any of my articles if some of them seems to be relevant for you.

Suitable curves

english summary:
Curves suitable for cryptography according to Rosing:
Elliptic Curves over Galois Fields.
Meaning curves over F2^n, meaning fields with characteristic 2 - and only nonsupersingular.
y^2 + xy = x^3 + a2x^2 + a6.
a6 cannot be zero
a2 can be zero.
.........
Ifølge Michael Rosing er de eneste kurver der er velegnede til kryptografisk brug kurver over Galois legemer. Dvs kurver på F2^n , med andre ord kurver med karakteristikken 2, og udvidelsesgraden n.- og det da kun hvis de ikke er supersingulære. Med andre ord skal de være af formen:

y^2 + xy = x^3 + a2x^2 + a6.
a6 må IKKE være nul.
a2 kan godt være nul.
Jeg tror fortsat polynominal basis er rarere at regne med end optinormal basis.

Hvordan mon denne klasse af anbefalede kurver matcher de i suite B angivne kurver. To be tested.

Punktaddition

Simpelt demonstationseksempel på R, med værdier max 'integer' størrelse:
Der konstrueres passende felter i html dokument der matcher id'erne i denne kode.
Denne meget lille primitive løsning tester IKKE for om punkterne overhovedet ligger på kurven - det er på eget ansvar. Det er kun i tilfældet P1=p2 hvor y1 er forskellig fra nul at kurvens koeficienter direkte indgår i additionsberegningen - og derfor der fejlkilden skal findes.
fgl placeres i eksternt javascript

window.onload= initForm;
function initForm(){
document.getElementById("ok").onclick = calcPoint;
}
function calcPoint(){
var x1= parseInt(document.getElementById("xx1").value);
var y1= parseInt(document.getElementById("yy1").value);

var x2= parseInt(document.getElementById("xx2").value);
var y2= parseInt(document.getElementById("yy2").value);
if((x2-x1) != 0 ){
var m = (y2-y1)/(x2-x1);
var x3 = (m*m)-x1-x2;
var y3= m*(x1-x3)-y1;
document.getElementById("r1").innerHTML="P3, (x3,y3) er "+x3+","+y3;
}
else if (((x2-x1) == 0) && ((y2-y1) != 0)){
alert("det er farligt at dividere med nul, P1 + P2 giver punktet i uendelig");}
else if ((x2-x1)==0 && (y2-y1)==0){
if(y2 != 0){
var A= parseInt(document.getElementById("AA").value);
var m = ((3*x1*x1)+A)/(2*y1);
var x3 = (m*m)-2*x1;
var y3= m*(x1-x3)-y1;
document.getElementById("r1").innerHTML="P3, (x3,y3) er "+x3+","+y3;
}
else{
alert("P1 + P2 giver punktet i uendelig");
document.getElementById("r1").innerHTML="eternal sunshine in a spotless mind";
}
}
else{
alert("noget gik helt galt");
document.getElementById("r1").innerHTML="42";
}

return false;
}