/***********************************************************
FILE : pohlig_hellman.c
Implementation of the Pohlig-Hellman Algorithm to solve
the Discrete Logarithm Problem in Zp, in other words:
Find the unique integer a such that y=a^result (mod p)
where p is prime and a is a primitive element.
The program accepts as inputs the big variables y, a,
p, and the prime factorization of p-1. Its output is the
big variable result.
PROGRAMMERS : Pedro Soria-Rodriguez and Jorge Guajardo
LAST UPDATE : 12/6/95
************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "mirdef.h"
#include "miracl.h"
typedef struct pp { /*Defines a structure which will be used to */
big value1; /*to store the factors of p-1, the partial */
big value2; /*results of Shank's Algorithm,and the partial */
struct pp *next; /*results of the Pohlig Hellman Algorithm before */
} pair; /*applying the Chinese Reminder Theorem */
/*
* Function Prototypes
*/
void shanks(big, big, big, big, big );
int insert(pair *, big, big);
void find_elem_list(int ,int , pair *, big , big );
int search(pair * ,big , big , big );
void main(void)
{
int number_of_factors,k,i,comp;
big y, a, p, result, dummy_base, dummy_exp, comparator, gamma, pminus1, pi, dummy,zero;
big zj, delta, j, temp, bee, uno, *x_base, *moduli_base, *x_ptr, *moduli_ptr;
pair factorization;
pair *ptemp, bees;
mirsys (155, 10); /* Initialize the Big Number library to use */
/* integers of up to ~500 bits */
dummy_base=mirvar(1);
dummy_exp =mirvar(1);
dummy=mirvar(0);
comparator=mirvar(0);
pminus1=mirvar(1); /* Initialization of all Variables */
gamma=mirvar(0);
pi=mirvar(0);
y=mirvar(1);
a=mirvar(1);
p=mirvar(1);
zj=mirvar(0);
delta=mirvar(0);
temp=mirvar(0);
bee=mirvar(0);
uno=mirvar(1);
result=mirvar(0);
zero=mirvar(0);
factorization.value1=mirvar(0); /* Initialize the two 'big' variables contatined */
factorization.value2=mirvar(0); /* in the 'factorization' and 'bees' structures */
factorization.next=NULL;
bees.value1=mirvar(0);
bees.value2=mirvar(0);
bees.next=NULL;
printf ("This program solves the Discrete Logarithm Problem ");
printf (" using the Pohlig-Hellman Algorithm\n");
printf (" x \n");
printf ("It solves the equation y=a (mod p) for the value x given y and a \n");
printf ("Enter the value of y: "); innum (y,stdin);
printf ("Enter the value of a: "); innum (a,stdin);
printf ("Enter the value of p: "); innum (p,stdin);
decr(p,1,pminus1); /*Create p-1*/
/* Ask the user for all prime factors of p-1 */
number_of_factors=1;
printf ("Enter the factors of p-1 and the powers to which they are raised\n");
printf ("(When you are done entering all the prime factors please enter 0):\n");
while(1)
{
printf ("Enter the prime factor #%d: ", number_of_factors);
innum(dummy_base,stdin);
if (0==compare(dummy_base,comparator))
break;
printf ("Enter the exponent for prime factor #%d: ", number_of_factors );
innum(dummy_exp,stdin);
insert(&factorization, dummy_exp, dummy_base);
number_of_factors++;
}
/* Create two pointers to big numbers that will contain the address of what is
to be the first two elements of two arrays */
x_base = (big *)calloc( (number_of_factors-1), sizeof(big) );
moduli_base = (big *)calloc( (number_of_factors-1), sizeof(big) );
x_ptr=x_base;
moduli_ptr=moduli_base;
for(k=0;k<(number_of_factors-1);k++) /* Initialize all elements in the arrays */
{ /* 'x' and 'moduli' */
(*x_ptr)=mirvar(0);
(*moduli_ptr)=mirvar(0);
x_ptr++;
moduli_ptr++;
}
/* Run the Pohlig-Hellman Algorithm */
/* Find x mod p_i for all 1 <= i < # of factors */
for ( i = 1 , moduli_ptr = moduli_base , x_ptr=x_base ;
i < number_of_factors ;
i++, x_ptr++, moduli_ptr++ )
{
copy(pminus1,dummy);
find_elem_list(i,(number_of_factors-1),&factorization,pi,dummy_exp);
divide(dummy,pi,dummy);
powmod(a,dummy,p,gamma);
copy(y,zj);
comp = -1;
for( j=mirvar(1) ; compare(j,dummy_exp) != (1) ; incr(j,1,j) )
{
powmod(pi,j,pminus1,temp);
copy(pminus1,dummy);
divide(dummy,temp,dummy);
powmod(zj,dummy,p,delta); /* Compute delta to be passed to Shank's alg. */
shanks(gamma,delta,p,pi,bee); /* Given gamma,delta and p, find b using Shank's */
insert(&bees,bee,j); /* Put all b's in a linked list */
decr(j,1,temp);
powmod(pi,temp,pminus1,temp);
multiply(bee,temp,temp);
powmod(a,temp,p,temp);
xgcd(temp,p,temp,temp,temp);
powmod2(temp,uno,zj,uno,p,zj);
}
powmod(pi,dummy_exp,p,dummy); /* Create p_i ^ exp modulus in variable dummy */
bee=mirvar(0);
comp=-1;
for( ptemp=bees.next, copy(zero,j) ; compare(j,dummy_exp) == (-1); incr(j,1,j) )
{
powmod2(ptemp->value1,uno,pi,j,dummy,temp);
add(bee,temp,bee);
ptemp = ptemp->next;
}
copy(bee,*x_ptr);
copy(dummy,*moduli_ptr);
free(&bees); /* Free the memory used in the list of b's */
bees.value1=mirvar(0); /* A new list of b's is to be created. */
bees.value2=mirvar(0);
bees.next=NULL;
}
/* Apply the Chinese Reminder Theorem */
/* 'moduli_base' is an array with the bases of the factors of p-1 */
crt_init((number_of_factors-1),moduli_base);
crt(x_base,result);
crt_end();
printf("Thus, the exponent is "); otnum(result,stdout);
}
/*
* Function find_elem_list:
*
* Given the position number that an element occupies on a linked list
* determined by *factors, the funtion returns the values of the base and
* the exponent corresponding to that element number.
*
* Parameters: elem_num: The element number that is to be searched for.
* maximum: The maximum number of elements in this list.
* *factors: Pointer to the first element in the linked list.
* base : Big number to which the found result is to be
* written.
* exponent: Same as base.
*/
void find_elem_list(int elem_num,int maximum, pair *factors, big base, big exponent)
{
int i;
pair *ptemp;
ptemp=factors;
if (elem_num>maximum)
{
printf ("This is not a possible element\n");
exit(0);
}
i=0;
while (1)
{
if (i==elem_num)
{
copy(ptemp->value1,exponent);
copy(ptemp->value2,base);
return;
}
ptemp=ptemp->next;
i++;
}
}
/*
* Function shanks():
*
* This function uses Shanks algorithm to find a = log_alpha of beta mod p.
* The result is returned through the pointer ret to a big.
*/
void shanks(big alpha, big beta, big p, big pi, big ret)
{
int comp;
big j,i,dummy;
big t;
big m; /* Variables Declaration */
big pminus;
big ein; /* Internal Joke */
pair list1,list2;
t = mirvar(0); /* Variables initialization */
list1.value1 = mirvar(0);
list1.value2 = mirvar(0);
list2.value1 = mirvar(0);
list2.value2 = mirvar(0);
list1.next = NULL;
list2.next = NULL;
m = mirvar(0);
j = mirvar(0);
ein = mirvar(1);
pminus = mirvar(0);
dummy = mirvar(0);
decr(p,1,pminus); /* Create a variable with p-1 */
nroot(pi,2,m); /* Determine m = sqrt(pi) */
incr(m,1,m); /* Increase m by 1 to obtain m = ceil( sqrt(pminus)) */
comp = (-1);
while ( comp == (-1) ) { /* Create 1st list with pairs (j, alpha^(mj) */
comp = compare(j,m);
multiply(m,j,t);
powmod(alpha,t,p,t);
insert(&list1,j,t);
incr(j,1,j);
}
/* Once the 1st list is created, compute each value of the second list */
/* and look for it in the 1st list as it is created. Therefore there is */
/* no such 'real' 2nd list. */
i = mirvar(0);
comp = (-1);
while ( comp == (-1) ) {
comp = compare(i,m);
powmod(alpha,i,p,t);
xgcd(t,p,t,t,t); /* alpha^(-i) */
powmod2(t,ein,beta,ein,p,t);
if (search (&list1, t, m,j)==1)
{
mad(m,j,i,pminus,dummy,t);
copy(t,ret);
return;
}
incr(i,1,i);
}
printf("Matching pairs not found. Error!\n");
}
/*
* Function insert():
*
* The function insert() takes as its argumets:
* pointer to the head of a list: *plist
* first value to be stored: el1
* second value to be stored: el2
*
*/
int insert(pair *plist, big el1, big el2)
{
pair *ptemp = NULL;
pair *pnew = NULL; /* Variable definition */
pair *pprevious = NULL;
pair *pchange = NULL;
/* Check that there is enough memory available to store the new element */
if ( (pnew = (pair *)malloc(sizeof(pair))) == NULL ) {
printf("Error - out of memory.\n");
exit(1);
}
/* If there is enough memory, it will create a new structure, pointed to by pnew */
else {
pnew->value1 = mirvar(0);
pnew->value2 = mirvar(0);
copy(el1,pnew->value1);
copy(el2,pnew->value2);
pnew->next = NULL;
}
/* Positioning the structure within the linked list ordered by the second value */
if (plist->next==NULL) {
plist->next = pnew;
/*Insert here if the list is empty */
}
else {
/* Compare new element with those on the list to determine the position
of the new element */
pprevious = plist;
ptemp = plist;
while (ptemp->next != NULL) {
ptemp=ptemp->next;
if ( (-1) == compare(pnew->value2,ptemp->value2) ) {
pchange = ptemp;
pprevious->next = pnew;
pnew->next = pchange;
return(0);
}
pprevious=ptemp;
}
ptemp->next = pnew;
}
return(1);
} /* End of the 'insert' function. */
/*
* Function search:
*
* This function performs a binary search in a linked list. The
* parameters to the funcion are:
* *list2 : a pointer to the first element of the list
* el : a big number to search for.
* m : the number of elements in the list
* ret : A big number used to return the value found.
*
* Return values:
* 0 : search unsuccessful
* 1 : search successful
*/
int search(pair *list2,big el, big m, big ret)
{
big top,bottom,middle;
pair *temp2;
big i,neg1;
big temp;
top = mirvar(0);
bottom = mirvar(0);
middle = mirvar(0);
neg1 = mirvar(-1);
temp=mirvar(0);
i=mirvar(0);
decr(m,1,top);
while( compare(top,bottom) != (-1) )
{
add(top,bottom,temp);
subdiv(temp,2,middle);
for (i=mirvar(-1), temp2 = list2->next; (compare(i,middle) == (-1)) ;
temp2 = temp2->next, incr(i,1,i));
if ( compare(el,temp2->value2) == 0 ) {
copy(temp2->value1,ret);
return(1);
}
else if ( compare(el,temp2->value2) == 1 )
incr(middle,1,bottom);
else
decr(middle,1,top);
}
return(0);
} /* End of search */