📄 schoof2.cpp
字号:
// Schoof's Original Algorithm - GF(2^m) Version!
// Mike Scott March 2000 mike@compapp.dcu.ie
//
// Counts points on GF(2^m) Elliptic Curve, y^2+xy = x^3+Ax^2+B a prerequisite
// for implementation of Elliptic Curve Cryptography
//
// The self-contained Windows Command Prompt executable for this program may
// obtained from ftp://ftp.computing.dcu.ie/pub/crypto/schoof2.exe
//
// Basic algorithm is due to Schoof
// "Elliptic Curves Over Finite Fields and the Computation of Square Roots
// mod p", Rene Schoof, Math. Comp., Vol. 44 pp 483-494
// Another very useful reference particularly for GF(2^m) curves is
// "Elliptic Curve Public Key Cryptosystems", Menezes,
// Kluwer Academic Publishers, Chapter 7
// Thanks are due to Richard Crandall for the tip about using prime powers
//
// **
// This program implements Schoof's original algorithm, augmented by
// the use of prime powers. By finding the Number of Points mod the product
// of many small primes and large prime powers, the final search for NP is
// greatly speeded up.
// Final phase search uses Pollard Lambda ("kangaroo") algorithm
// This final phase effectively stretches the range of Schoof's
// algorithm by about 70-80 bits.
// This approach is only feasible due to the use of fast multiplication
// methods for large degree polynomial multiplication
// **
//
// Ref "Monte Carlo Methods for Index Computation"
// by J.M. Pollard in Math. Comp. Vol. 32 1978 pp 918-924
//
// An "ideal" curve is defined as one with with
// (i) 2*(a prime) number of points if trace(A)=1, or
// (ii) 4*(a prime) number of points if trace(A)=0.
//
// Note that the number of points on the curve is invariant if we substitute
// B by B^2 (this is trivial, but I haven't seen it anywhere).
//
// Proof: if (x,y) is a point on
// y^2+xy=x^3+Ax^2+B
// Then clearly (X,Y) = (x^2,y^2) is a point on
// y^4+x^2.y^2=x^6+Ax^4+B^2 (square both sides)
//
// so Y^2+XY=X^3+AX^2+B^2
//
// When using the "schoof2" program, the -s option is particularly useful
// and allows automatic search for an "ideal" curve. If a curve order is
// exactly divisible by a small prime >2, that curve is immediately abandoned,
// and the program moves on to the next, incrementing the B parameter of
// the curve. This is a fairly arbitrary but simple way of moving on to
// the "next" curve (but note the point made above about B -> B^2).
//
// NOTE: The output file can be used directly with for example the ECDSA
// programs IF AND ONLY IF an ideal curve is found. If you wish to use
// a less-than-ideal curve, you will first have to factor NP completely, and
// find a random point of large prime order.
//
// This implementation is free. No terms, no conditions. It requires
// version 4.30 or greater of the MIRACL library (a Shareware, Commercial
// product, but free for non-profit making use),
// available from ftp://ftp.computing.dcu.ie/pub/crypto/miracl.zip
//
// However this program may be used (unmodified) to generate curves for
// commercial use.
//
// 32-bit build only
//
// Revision history
//
// Timings for test curve Y^2+XY=X^3+X^2+52 mod 2^191
// Pentium III 450MHz
//
// Rev. 0 - 45 minutes
// Rev. 1 - 41 minutes
// Rev. 2 - 36 minutes
//
// Note that a small speed-up can be obtained by using an integer-only
// build of MIRACL. See mirdef.hio
//
// Pentium III 450MHz
//
// GF(2^191) - 36 minutes
// GF(2^223) - 1.9 hours
// GF(2^251) - 4 hours
// GF(2^283) - 12.3 hours
//
// Copyright Shamus Software Ltd. 2000
//
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstring>
#include "ec2.h" // GF(2^m) Elliptic Curve Class
#include "crt.h" // Chinese Remainder Theorem Class
// poly2.h implements polynomial arithmetic. Karatsuba multiplication is
// used for maximum speed, as the polynomials get very big. For example when
// searching for the curve cardinality mod the prime 31, the polynomials are
// of degree (31*31-1)/2 = 480. But all that gruesome detail is hidden away.
//
// poly2mod.h implements polynomial arithmetic wrt to a preset polynomial
// modulus. This looks a bit neater. Function setmod() sets the modulus
// to be used.
#include "poly2.h"
#include "poly2mod.h"
using namespace std;
Miracl precision=12; // max. 12x32 bits per big number
Poly2Mod MFX,XX;
Big B;
Big D;
Big A;
void elliptic_dup(Poly2Mod& X,Poly2Mod& Y,Poly2Mod& Yy,Poly2Mod& Z)
{ // point doubling
Poly2Mod W1,W2,W3,W4,W5y,W5;
W1=(Z*Z); // Z^2
W2=X*X; // X^2
W3=W2*W2; // X^4
W4=X*W1; // =Z
X=(X+D*W1);
X=X*X;
X=X*X;
W5=W4+W2+Y*Z;
W5y=Yy*Z;
Z=W4;
Y=W3*Z+W5*X;
Yy=W5y*X;
}
//
// This is addition formula for two distinct points on an elliptic curve
// Works with projective coordinates which are automatically reduced wrt a
// polynomial modulus
// Remember that the expression for the Y coordinate of each point
// is of the form A(x)+B(x).y
// We know Y^2=X^3+AX^2+B+XY, but we don't have an explicit expression for Y
// So if Y^2 ever crops up - substitute for it.
void elliptic_add(Poly2Mod& XT,Poly2Mod& XTy,Poly2Mod& YT,Poly2Mod& YTy,Poly2Mod& ZT,Poly2Mod& X,Poly2Mod& Y,Poly2Mod& Yy)
{ // add (X,Y,Z) to (XT,YT,ZT) on an elliptic curve
// The point Y is of the form A(x)+B(x).y
Poly2Mod ZT2,ZT3,W1,W2,W3,W4,W5,W5y,W6,W6y,W7,W8,W8y,W9,W9y;
ZT2=(ZT*ZT);
W3=XT+X*ZT2;
ZT3=ZT2*ZT;
W5=ZT3*Y;
W5y=ZT3*Yy;
W6=YT+W5;
W6y=YTy+W5y;
if (iszero(W3))
{
if (iszero(W6) && iszero(W6y))
{ // should have doubled!
elliptic_dup(XT,YT,YTy,ZT);
XTy=0;
return;
}
ZT.clear();
return;
}
ZT*=W3;
W8=W6*X+ZT*Y;
W8y=W6y*X+ZT*Yy;
W9=W6+ZT;
W1=ZT*ZT;
W2=W6y*W6y;
XT=A*W1+W3*W3*W3+W6*W9+W2*MFX;
XTy=W9*W6y+W6y*W6+W2*XX;
W2=W6y*XTy;
YT=W9*XT+W2*MFX+W8*W1;
YTy=XT*W6y+W9*XTy+W2*XX+W1*W8y;
}
//
// Program to compute the order of a point on an elliptic curve
// using Pollard's lambda method for catching kangaroos.
//
// As a way of counting points on an elliptic curve, this
// has complexity O(p^(1/4))
//
// However Schoof puts a spring in the step of the kangaroos
// allowing them to make bigger jumps, and lowering overall complexity
// to O(p^(1/4)/sqrt(L)) where L is the product of the Schoof primes
//
// See "Monte Carlo Methods for Index Computation"
// by J.M. Pollard in Math. Comp. Vol. 32 1978 pp 918-924
//
// This code has been considerably speeded up using ideas from
// "Parallel Collision Search with Cryptographic Applications", J. Crypto.,
// Vol. 12, 1-28, 1999
//
#define STORE 80
#define HERD 5
EC2 wild[STORE],tame[STORE];
Big wdist[STORE],tdist[STORE];
int wname[STORE],tname[STORE];
Big kangaroo(Big p,Big order,Big ordermod,int TR,BOOL &found)
{
EC2 ZERO,K[2*HERD],TE[2*HERD],X,P,G,table[50],trap;
Big start[2*HERD],txc,wxc,mean,leaps,upper,lower,middle,a,b,x,y,n,w,t,nrp;
int i,jj,j,m,sp,nw,nt,cw,ct,k,distinguished;
Big D[2*HERD],s,distance[50],real_order;
BOOL bad,collision,abort;
found=FALSE;
forever
{
// find a random point on the curve
do
{
x=rand(p);
} while (!P.set(x,x));
lower=p+1-2*sqrt(p)-3; // lower limit of search
upper=p+1+2*sqrt(p)+3; // upper limit of search
w=1+(upper-lower)/ordermod;
leaps=sqrt(w);
mean=HERD*leaps/2; // ideal mean for set S=1/2*w^(0.5)
distinguished=1<<(bits(leaps/16));
distinguished++; // powers of 2 create "resonances" with 2^m
// so bump it up to odd.
for (s=1,m=1;;m++)
{ /* find table size */
distance[m-1]=s*ordermod;
s*=2;
if ((2*s/m)>mean) break;
}
table[0]=ordermod*P;
for (i=1;i<m;i++)
{ // double last entry
table[i]=table[i-1];
table[i]+=table[i-1];
}
middle=(upper+lower)/2;
if (ordermod>1)
middle+=(ordermod+order-middle%ordermod);
for (i=0;i<HERD;i++) start[i]=middle+13*ordermod*i;
for (i=0;i<HERD;i++) start[HERD+i]=13*ordermod*i;
for (i=0;i<2*HERD;i++)
{
K[i]=start[i]*P; // on your marks
D[i]=0; // distance counter
}
cout << "Releasing " << HERD << " Tame and " << HERD << " Wild Kangaroos\n";
nt=0; nw=0; cw=0; ct=0;
collision=FALSE; abort=FALSE;
forever
{
for (jj=0;jj<HERD;jj++)
{
K[jj].get(txc);
i=txc%m; /* random function */
if (txc%distinguished==0)
{
if (nt>=STORE)
{
abort=TRUE;
break;
}
cout << "." << flush;
tame[nt]=K[jj];
tdist[nt]=D[jj];
tname[nt]=jj;
for (k=0;k<nw;k++)
{
if (wild[k]==tame[nt])
{
ct=nt; cw=k;
collision=TRUE;
break;
}
}
if (collision) break;
nt++;
}
D[jj]+=distance[i];
TE[jj]=table[i];
}
if (collision || abort) break;
for (jj=HERD;jj<2*HERD;jj++)
{
K[jj].get(wxc);
j=wxc%m;
if (wxc%distinguished==0)
{
if (nw>=STORE)
{
abort=TRUE;
break;
}
cout << "." << flush;
wild[nw]=K[jj];
wdist[nw]=D[jj];
wname[nw]=jj;
for (k=0;k<nt;k++)
{
if (tame[k]==wild[nw])
{
ct=k; cw=nw;
collision=TRUE;
break;
}
}
if (collision) break;
nw++;
}
D[jj]+=distance[j];
TE[jj]=table[j];
}
if (collision || abort) break;
multi_add(2*HERD,TE,K); // jump together - its faster
}
cout << endl;
if (abort)
{
cout << "Failed - this should not happen! - trying again" << endl;
continue;
}
nrp=start[tname[ct]]-start[wname[cw]]+tdist[ct]-wdist[cw];
G=P;
G*=nrp;
if (G!=ZERO || nrp%2==1)
{
cout << "Sanity Check Failed. Please report to mike@compapp.dcu.ie" << endl;
exit(0);
}
if (prime(nrp/2))
{
cout << "NP is 2*Prime!" << endl;
cout << "NP= 2*" << nrp/2 << endl;
found=TRUE;
break;
}
if (nrp%4==0)
{
if (prime(nrp/4) && TR==0)
{
cout << "NP is 4*Prime!" << endl;
cout << "NP= 4*" << nrp/4 << endl;
found=TRUE;
break;
}
}
// final checks....
real_order=nrp; i=0;
forever
{
sp=get_mip()->PRIMES[i];
if (sp==0) break;
if (real_order%sp==0)
{
G=P;
G*=(real_order/sp);
if (G==ZERO)
{
real_order/=sp;
continue;
}
}
i++;
}
if (real_order <= 4*sqrt(p))
{
cout << "Low Order point used - trying again" << endl;
continue;
}
real_order=nrp;
for (i=0;(sp=get_mip()->PRIMES[i])!=0;i++)
while (real_order%sp==0) real_order/=sp;
if (real_order==1)
{ // all factors of nrp were considered
cout << "NP= " << nrp << endl;
break;
}
if (prime(real_order))
{ // all factors of NP except for one last one....
G=P;
G*=(nrp/real_order);
if (G==ZERO)
{
cout << "Failed - trying again" << endl;
continue;
}
else
{
cout << "NP= " << nrp << endl;
break;
}
}
// Couldn't be bothered factoring nrp completely
// Probably not an interesting curve for Cryptographic purposes anyway.....
// But if 20 random points are all "killed" by nrp, its almost
// certain to be the true NP, and not a multiple of a small order.
bad=FALSE;
for (i=0;i<20;i++)
{
do
{
x=rand(p);
} while (!P.set(x,x));
G=P;
G*=nrp;
if (G!=ZERO)
{
bad=TRUE;
break;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -