index-calculus.gp

p=3023;
g=Mod(5,p);
y=Mod(555,p);
l=14;

base=primes(l);
maxprime=base[l];

if(type(rels)=="t_POL",rels=[]);

makerels=() -> {
    my(r,f,m,i,j,rel);
    rel=[];
    for(i=1,100,
        r=random(p-1);
        f=factor(lift(g^r));
	for(j=1,matsize(f)[1],
            if(f[j,1] > maxprime,next(2)));
        rel=concat(rel,[[r,f]]));
    rel;
};

findxrel=() -> {
    my(j,r,f);
    until(0,
          r=random(p-1);
          f=factor(lift(y*g^r));
	  for(j=1,matsize(f)[1],
              if(f[j,1] > maxprime,next(2)));
          return([r,f]);)
};

factoredtosparse=(u) -> {
    my(i,r,j);
    r=vector(l);
    j=1;
    for(i=1,l,if(u[j,1] == base[i],r[i]=u[j,2];j=j+1;if(j>matsize(u)[1],break)));
    r;
};

makematrix=(rels) -> {
    my(A,y,t);
    A = [];
    y = vector(length(rels));
    for(i=1,length(rels),
	y[i] = rels[i][1];
        A = concat(A, [ factoredtosparse(rels[i][2])~ ]););
    [Mat(A)~,y~];
};

dosieve=() -> {
    until((matrank(Aq)==l) && (matrank(A2)==l),
        print("One");
        rels=concat(rels, makerels());
        print("Two");
        R = makematrix(rels);
        print("Three");
        Aq = Mod(R[1],(p-1)/2);
        yq = Mod(R[2],(p-1)/2);
        A2 = Mod(R[1],2);
        y2 = Mod(R[2],2););
};

dolinalg=() -> {
    rq = matindexrank(Aq)[1];
    AAq = vecextract(Aq~,rq)~;
    yyq = vecextract(yq,rq);
    uq = matsolve(AAq,yyq);
    r2 = matindexrank(A2)[1];
    AA2 = vecextract(A2~,r2)~;
    yy2 = vecextract(y2,r2);
    u2 = matsolve(AA2,yy2);
    u=chinese(uq,u2)
};

dofinal=() -> {
    r1 = findxrel();
    r2 = factoredtosparse(r1[2]);
    dlog = lift(Mod(-r1[1],p-1) + sum(i=1,l,r2[i]*u[i]));
};