Some modules if required..

[Kaushik Chakraborty]

Hi,
I was trying to implement quadratic-sieve using gmp(people generally 
use freelip). I found that the following three functions are absolutely 
necessary which are not present in gmp-3.0 or gmp-3.1.1...
These are..
1)calculating log of arbitary precision integers..which i named as
mpz_ln()..
2.)calculating log a base b..which i called mpz_log()
3.)finding out the square-root mod a large prime (this is very much
necessary for the sieving stage)..which i called mpz_sqrtmod()...

The first two are adopted from freelip 1.1 & (i guess) requires Lenstra's
permission..
The last one is from Handbook of Applied Cryptography..Menezes..

I successfully implemented quadratic-sieve (with all it's variants)
with the newly added functions..these are tested for various inputs &
working fine...
the log functions can easily extended for mpf,..etc.
I am inserting the codes of the functions as well as a test-pgm. with its 
outputs for various test-cases..
for some cases there are some decisions to be taken..say, ln(a -ve no.)
should be printed as pari or as freelip..etc. ...i put comments on any
decisions so that they can be reverted easily if necessary..

also i have debuglevel enabled which is no longer necessary..the
debug-printf's are useful if someone want's to cross-check all parts of
the code...

HERE A BIG ATTACHMENT FOLLOWS
*******************************

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include "gmp-impl.h"

#define DEBUGGING 1

main(int argc, char **argv)
{

mpz_t n;
mpz_t p,sqrtmod;

/*FUNCTION PROTOTYPES*/
int mpz_sqrtmod(mpz_t sqrtmod, mpz_t n, mpz_t p);
double mpz_ln(mpz_t n);
double mpz_log(mpz_t a, mpz_t b);

mpz_init_set_str(n, argv[1], 10);
mpz_init_set_str(p, argv[2], 10);

printf("legendre= %d\n",mpz_legendre(n, p));/*TO CHECK SQRTMOD EXISTENCE*/
mpz_init(sqrtmod);
printf("sqrtmod-return= %d\n",mpz_sqrtmod(sqrtmod,n, p));/*check return*/
printf("sqrtmod-val= ");/*check value*/
mpz_out_str(stdout , 10, sqrtmod);
printf("\n");

printf("Natural log= %f\n",mpz_ln(n));/*check mpz_ln */
printf("a log (base)b= %f\n",mpz_log(n, p));/*check mpz_log*/

}

/*for checking the individual functions comment out the others' checks*/

/*Functions */

/*mpz_ln*/

/*translated from Lenstra's freelip..requires permission*/
/*for -ve numbers , logarithms is same as positive counterparts (freelip
behaves like this)..pari gives
log(abs)+I(pi)..i.e in the complex form..here we adopt the former
one..printing +I(pi) is trivial but then
we don't have a "complex" datatype in GMP..so doesn't make sense...*/
double
#if __STDC__
mpz_ln(mpz_srcptr n)
#else
mpz_ln(n)
mpz_srcptr n;*/
#endif
{
        extern double log();
        static double log_2 = -1.0;
        register long sn;
        


        if ((PTR(n)[0]==0) || (ALLOC(n) < 0))
        {
                fprintf(stderr,"zero argument in mpz_ln\n");
                fflush(stderr);
                exit(-1);
                /*considered as a fatal error*/
                /*otherwise can return -1 */
                
        }
        
        sn=ALLOC(n);
        if (sn == 1)
        {       
                if (PTR(n)[0]==0)
                {
                fprintf(stderr,"zero argument in mpz_ln\n");
                fflush(stderr);
                /*return (0.0);*/
                exit(-1);
                }
        return (log((double) (PTR(n)[0])));
        }
        
        if (log_2 < 0)
        {
                log_2 = log((double) 2);
        }
        return ((sn - 2) * BITS_PER_MP_LIMB * log_2 +
                log(((double)(1<<BITS_PER_MP_LIMB)) * PTR(n)[sn-1] +
PTR(n)[sn - 2]));
}


/*mpz_log*/

/*translated from Lenstra's freelip..requires permission*/
double 
#if __STDC__
mpz_log(mpz_srcptr a, mpz_srcptr b)
#else
mpz_log(a, b)
        mpz_srcptr a;
        mpz_srcptr b;
#endif

{

        return (mpz_ln(a) / mpz_ln(b));
}


int
#if __STDC__
mpz_sqrtmod(mpz_ptr m, mpz_srcptr n, mpz_srcptr p)
#else
mpz_sqrtmod(m, n, p)
        mpz_ptr m;
        mpz_srcptr n;
        mpz_srcptr p;

#endif
/*compute squareroot mod n prime p result is x*/
/*returning -1 means failure or problem, in that case it does not touches
the existing value of sqrtmod passed as argument*/
/*no side effect*/
/*It is very highlevel*/
/*Tested with various inputs*//*Works correctly*/
/*It always returns the smaller quadratic residue..
the other residue (if not same) is trivially (p-return value)*/
/*returning the smaller one is useful in siutations like quadratic sieve
to avoid both-sides sieving..other-wise we can return (r1, r2)..minor
modification is needed*/
/*I put some debugging info ...no longer necessary*/
/*Source:Handbook of applied cryptography-Menezes pp 100-101*/
{

        mpz_t ntemp,x,xtemp, d,p_minus_one,t,largetemp,verylargetemp,
temp,b,c,ninv,exp;
        unsigned long i,ptempui, ntempui, s;
        mp_size_t sizep,sizen;
        

        
        if ((SIZ(n) <0 ) || (SIZ(p) <0 ))
        {
        #if DEBUGGING
        printf("Sqrtmod:Positive integers expected\n");
        #endif  
        return( -1);
        }
        
        
        
        if ((SIZ(n)== 0)|| (mpz_cmp_ui(p,2)<= 0))
        {
        #if DEBUGGING
        printf("Sqrtmod:Composite should be >= 1 & p an odd prime\n");
        #endif  
        return( -1);
        }
        
        
        
        if (mpz_probab_prime_p(p, 30)== 0)/*composite p*/
        {
        #if DEBUGGING
        printf("Sqrtmod:Composite entered as prime input\n");
        #endif  
        
        return(-1);
        }
        
        sizen=ABSIZ(n);
        sizep=ABSIZ(p);
        
        MPZ_TMP_INIT(ntemp, sizep);
        
        mpz_mod(ntemp, n, p);/*put n in the range 1<=n<=p-1*//*to avoid
side effect*/
        
        if (mpz_cmp_ui(ntemp, 0)==0)
        {
        #if DEBUGGING
        printf("Sqrtmod: prime divides composite\n");
        #endif  
        return(0);
        
        }
        
        if (mpz_legendre(ntemp, p)==-1)
        {
        #if DEBUGGING
        printf("Sqrtmod:Not a valid quadratic residue on this prime\n");
        #endif  
        return(-1);
        }
        
        
        if (mpz_cmp_ui(ntemp, 1)==0)
                {
                 mpz_set_ui(m, 1);
                return(1);
                }
        /* if number is small then search for soultion */       
        
        
        if (mpz_cmp_ui(p, 50)<=0 )
                {
                ptempui=mpz_get_ui(p);/*ptempui is  the (unsigned) prime*/
                ntempui=mpz_get_ui(ntemp);/*ntempui  is the unsigned
reduced composite*/
                for (i=1; i<=ptempui; i++)
                        if (i*i % ptempui == ntempui)
                                {
                                #if DEBUGGING
                                printf("Sqrtmod:Found from small prime
inspection\n");         
                                #endif  
                                mpz_set_ui(m, i);
                                return(1);
                                }
                }

        MPZ_TMP_INIT(x,sizep);
        MPZ_TMP_INIT(temp,sizep+1);
                
        
        if (mpz_mod_ui(temp, p, 4)==3)
                {   /* p == 3 mod 4 *//*special case*/
                        mpz_add_ui(temp, p, 1);
                        mpz_tdiv_q_2exp(temp,temp,2);
                        mpz_powm(x, ntemp,temp, p);
                        mpz_tdiv_q_2exp(temp,p,1);
                        if (mpz_cmp(x, temp) <0)
                        {
                        #if DEBUGGING
                        printf("Sqrtmod:Found from special case p=3 mod
4\n");  
                        #endif  
                        mpz_set(m, x);
                        return(1);
                        }
                                
                        else
                        {
                        #if DEBUGGING
                        printf("Sqrtmod:Found from special case p=3 mod
4\n");  
                        #endif  
                        mpz_sub(m,p, x);
                        return(1);
                        }
                }

          mpz_init(largetemp);/*this will hold different values of
different sizes*/
          /*if (p % 8==5)*/
          if (mpz_mod_ui(temp, p, 8)==5)
                {   /* p == 5 mod 8 *//*special case*/
                        
                        MPZ_TMP_INIT(d, sizep);
                        mpz_sub_ui(temp, p, 1);
                        mpz_tdiv_q_2exp(temp,temp,2);
                        mpz_powm(d, ntemp,temp, p);
                        mpz_abs(d, d);
                        
                        if (mpz_cmp_ui(d, 1)==0)
                        {
                         mpz_add_ui(temp, p, 3);
                         mpz_tdiv_q_2exp(temp,temp,3);
                         mpz_powm(x, ntemp,temp, p);
                        
                        }
                        mpz_sub_ui(temp, p, 1);
                if (mpz_cmp(d, temp)==0)
                        {
                                                                
                        mpz_mul_ui(largetemp,ntemp,4);
                        mpz_sub_ui(temp, p, 5);
                        mpz_tdiv_q_2exp(temp,temp,3);
                        mpz_powm(x, largetemp,temp, p);
                        mpz_div_ui(largetemp,ntemp,2);
                        mpz_mul(largetemp,x,largetemp);
                        mpz_mod(x,largetemp,p);                 
                        }
                
                mpz_tdiv_q_2exp(temp,p,1);
                
                if (mpz_cmp(x, temp) <0)
                        {
                        #if DEBUGGING
                        printf("Sqrtmod:Found from special case p=5 mod
8\n");  
                        #endif  
                        mpz_set(m, x);
                        mpz_clear(largetemp);
                        return(1);
                        }
                        
                else
                        {
                        #if DEBUGGING
                        printf("Sqrtmod:Found from special case p=5 mod
8\n");  
                        #endif  
                        mpz_sub(m,p, x);
                        mpz_clear(largetemp);
                        return(1);
                        }
                        
                }
                
        /*general case*/
        /*these allocations are correct as far as i tested..if something
goes wrong..may be here*/
        MPZ_TMP_INIT(b, sizep);/*this is always mod p*/
        MPZ_TMP_INIT(c, sizep);/*this is always mod p*/
        MPZ_TMP_INIT(d, sizep);/*this is always mod p*/
        MPZ_TMP_INIT(p_minus_one, sizep+1); /*one extra due to
carry/borrow..as in mpz_sub_ui*/
        MPZ_TMP_INIT(ninv, MAX(sizen, sizep)+1);/*as in mpz_invert*/
        MPZ_TMP_INIT(t, sizep+1);/*may be of the same order of
p_minus_one*/
        MPZ_TMP_INIT(exp, sizep); /*mod p*/
        
        
        mpz_sub_ui(p_minus_one, p, 1);
        
        for (mpz_set_ui(b,1);mpz_cmp(b, p_minus_one) <=0;mpz_add_ui(b,b,
1))
                {
                if (mpz_legendre(b, p)==-1)
                break;
                }
        /*we exit with b as a qnr mod p*/
        mpz_set_ui(t,2);
        s=mpz_remove (t,p_minus_one , t);
        if (mpz_invert(ninv, ntemp, p)==0)
        return( -1);
        mpz_powm(c,b,t,p);
        mpz_add_ui(temp, t, 1);
        mpz_tdiv_q_2exp(t,temp,1);
        mpz_powm(x,ntemp,t,p);
        
        for (i=1; i<=s-1;i++)
        {
                mpz_ui_pow_ui(exp,2, s-i-1);
                
                mpz_mul(largetemp, x, x);
                mpz_mul(largetemp, largetemp,ninv);
                mpz_powm(d, largetemp,exp,p);
                mpz_add_ui(temp, d, 1);
                mpz_mod(temp, temp, p);
                if (mpz_cmp_ui(temp, 0)==0)/*d=-1(p)*/
                        {
                        mpz_mul(largetemp,x,c);
                        mpz_mod(x,largetemp,p);
                        
                        }
                mpz_set_ui(temp, 2);/*this is not needed if the doubt
below is answered*/
                /*mpz_powm_ui(c,c,2,p);*/ /*this should be used but
strangely it gives a segmenation fault..it occurs is the last free of
mpz_powm_ui*/
                /*what is the actual limbspace reqd. for it to avoid it??i
checked that psize+2 works but why?*/
                mpz_powm(c,c,temp,p);
        }
        
        mpz_div_ui(temp, p, 1);

                if (mpz_cmp(x, temp) <0)
        
                        {
                        mpz_set(m, x);
                        mpz_clear(largetemp);
                        return(1);
                        }
                else
                        {
                        mpz_sub(m,p, x);
                        mpz_clear(largetemp);
                        return(1);
                        }       

}       



/*****************************TEST CASES****************************/

THESE ARE THE TEST CASES FOR MPZ_SQRTMOD()
**********************************************************************************************
./test -2 3
legendre= 1
Sqrtmod:Positive integers expected
sqrtmod-val= -1
sqrtmod-val= 0
**********************************************
./test 2 -3
legendre= 1
Sqrtmod:Positive integers expected
sqrtmod-val= -1
sqrtmod-val= 0
**********************************************
./test -2 -3
legendre= 0
Sqrtmod:Positive integers expected
sqrtmod-val= -1
sqrtmod-val= 0
**********************************************
./test 0 4
legendre= 1
Sqrtmod:Composite should be >= 1 & p an odd prime
sqrtmod-val= -1
sqrtmod-val= 0
*********************************************
./test 1 4
legendre= 1
Sqrtmod:Composite entered as prime input
sqrtmod-val= -1
sqrtmod-val= 0
*********************************************
./test 15 5
legendre= 0
Sqrtmod: prime divides composite
sqrtmod-val= 0
sqrtmod-val= 0
*********************************************
./test 1 59
legendre= 1
sqrtmod-val= 1
sqrtmod-val= 1
*********************************************
./test 199999999999999999677777777 47
legendre= -1
Sqrtmod:Not a valid quadratic residue on this prime
sqrtmod-val= -1
sqrtmod-val= 0
*********************************************
./test 1999999999999999996777777789 47
legendre= 1
Sqrtmod:Found from small prime inspection
sqrtmod-return= 1
sqrtmod-val= 2
*********************************************
./test
11111111111111222222222222222222222222233333333333333333333333334444444444555555555
610692533270508750441931226384209856405876657993997547171387
legendre= -1
Sqrtmod:Not a valid quadratic residue on this prime
sqrtmod-return= -1
sqrtmod-val= 0
*********************************************
 ./test
11111111111111222222222222222222222222233333333333333333333333334444444444555555556
610692533270508750441931226384209856405876657993997547171387
legendre= 1
Sqrtmod:Found from special case p=3 mod 4
sqrtmod-return= 1
sqrtmod-val= 275516653713458684009505497117074918097093024280973251321004
*********************************************
./test
11111111111111222222222222222222222222233333333333333333333333334444444444555555559
622288097498926496141095869268883999563096063592498055290461
legendre= 1
Sqrtmod:Found from special case p=5 mod 8
sqrtmod-return= 1
sqrtmod-val= 277943860475301908330701701857472838915699845323397799563304
*********************************************
./test
379999999999665499999999996666666667999999999999999999934555555555555555551233333333333333333333
669483106578092405936560831017556154622901950048903016651289
legendre= 1
sqrtmod-return= 1
sqrtmod-val= 615576939515209023669026264567910187110944474760769079588920
*********************************************
./test
379999999999665499999999996666666667999999999999999999934555555555555555551233333333333333333333999999999999999999999944444444444445555555555555567777777777777777777777898
669483106578092405936560831017556154622901950048903016651289
sqrtmod-return= 1
sqrtmod-val= 284741855208397272680783312598950029308734489804234081484029
*********************************************

THESE ARE THE TEST CASES FOR MPZ_LN()
**********************************************************************************************
 ./test 2
Natural log= 0.693147
*********************************************
 ./test -2
Natural log= 0.693147
*********************************************
./test 0
non-positive argument in mpz_ln
*********************************************
./test
9999999999999999999999678888888888888888888888888888888845555555555555566
Natural log= 168.088710
*********************************************
./test
1113333333333444444444556666666666666666666667789999999990000077777777777777888988888898785
Natural log= 207.339726
*********************************************
 ./test
111333333333344444444455666666666666666666666778999999999000007777777777777788898888889800000000099999999999999666666666666666677777777777777777445555555555555567
Natural log= 370.823558
*********************************************

THESE ARE THE TEST CASES FOR MPZ_LOG()
**********************************************************************************************
 ./test 0 0
zero argument in mpz_ln
*********************************************
./test -1 0
zero argument in mpz_ln
*********************************************
./test 0 -999999999
zero argument in mpz_ln
*********************************************
./test
11111222222222222222223344444444444444444444444444444444456666666666666666666666666667888888
2
a log (base)b= 302.447471
*********************************************
./test
8999999999778845222222334444477884444444444444444444444456666666666666666666666666667888888
88777555555555555555555555666666666666666666666666666666666666666666663444444444444444444444444444444444444444444444444555555555555555555567777777777
a log (base)b= 0.610643
*********************************************



Regards,
 


                                                Kaushik Chakraborty
                                                M.Tech Computational Science;
                                                SERC IISC BANGALORE;
                                                E-Mail:address@hidden