Some modules if required..(please disregard the previous mail)

[Kaushik Chakraborty]

OOPS!
mistakenly attached a earlier version of code..
sorry for that...
pls ignore the previous mail...
thanx,


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

On Sat, 16 Dec 2000, Kaushik Chakraborty wrote:

> 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,largetemp,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
> 
>