[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Some modules if required..(please disregard the previous mail)
From: |
Kaushik Chakraborty |
Subject: |
Some modules if required..(please disregard the previous mail) |
Date: |
Sat, 16 Dec 2000 17:55:13 +0530 (IST) |
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
>
>
- Some modules if required.., Kaushik Chakraborty, 2000/12/16
- Some modules if required..(please disregard the previous mail),
Kaushik Chakraborty <=