Borderbasix

reconstruct.hpp
Go to the documentation of this file.
1 void chinese(mpz_t res,
2  const mpz_t r1, const mpz_t m1,
3  const mpz_t r2, const mpz_t m2)
4 {
5  mpz_t tmpres;
6  int swap=1;
7  mpz_t bez1_1,bez1_2;
8  mpz_t bez2_1,bez2_2;
9  mpz_t a,b;
10  mpz_t q;
11  mpz_init(tmpres);
12  mpz_init(q);
13  mpz_init_set(a,m1);
14  mpz_init_set(b,m2);
15  mpz_init_set_si(bez1_1,1);
16  mpz_init_set_si(bez1_2,0);
17  mpz_init_set_si(bez2_1,0);
18  mpz_init_set_si(bez2_2,1);
19 
20  if (mpz_cmp(a,b)<0)
21  {
22  mpz_swap(a,b);
23  mpz_swap(bez1_1,bez2_1);
24  mpz_swap(bez1_2,bez2_2);
25  swap*=-1;
26  }
27  for(;mpz_cmp_si(b,1);)
28  {
29 
30  mpz_tdiv_qr(q,a,a,b);
31  mpz_submul(bez1_1,bez2_1,q);
32  mpz_submul(bez1_2,bez2_2,q);
33  swap*=-1;
34  mpz_swap(a,b);
35  mpz_swap(bez1_1,bez2_1);
36  mpz_swap(bez1_2,bez2_2);
37  }
38  mpz_mul(tmpres,bez2_2,m2);
39  mpz_mul(tmpres,tmpres,r1);
40  mpz_mul(q,bez2_1,m1);
41  mpz_mul(q,q,r2);
42  mpz_add(tmpres,tmpres,q);
43  mpz_mul(q,m1,m2);
44  mpz_mod(tmpres,tmpres,q);
45 
46  mpz_set(res,tmpres);
47  mpz_clear(tmpres);
48  mpz_clear(q);
49  mpz_clear(bez1_1);
50  mpz_clear(bez1_2);
51  mpz_clear(bez2_1);
52  mpz_clear(bez2_2);
53  mpz_clear(a);
54  mpz_clear(b);
55 
56 }
57 
58 int reconstruct(mpq_t res, const mpz_t u, const mpz_t m)
59 {
60  int ret=-1;
61  mpz_t N;
62  mpz_t r0,t0;
63  mpz_t r1,t1;
64  mpz_t q,tmp;
65  mpz_init_set(N,m);
66  mpz_fdiv_q_ui(N,N,2);
67  mpz_sqrt(N,N);
68  mpz_init_set(r0,m);
69  mpz_init_set(r1,u);
70  mpz_mod(r1,u,m);
71  mpz_init(t0);
72  mpz_init_set_ui(t1,1);
73  mpz_init(q);
74  mpz_init(tmp);
75  //printf(" N "); mpz_out_str(stdout,10,N);printf("\n");
76  while(mpz_cmp(r1,N)>0)
77  {
78  // printf("passage\n");
79  mpz_fdiv_q(q,r0,r1);
80  mpz_set(tmp,r1);
81  //r1=r0-q*r1;
82  mpz_mul(r1,r1,q);
83  mpz_neg(r1,r1);
84  mpz_add(r1,r0,r1);
85  //
86  mpz_set(r0,tmp);
87  //
88  //t1=t0-q*t1;
89  mpz_set(tmp,t1);
90  mpz_set(t1,t0);
91  mpz_submul(t1,q,tmp);
92  //
93  //t0=tmp;
94  mpz_set(t0,tmp);
95  //printf("t0 ");mpz_out_str(stdout,10,t0);printf("\n");
96  //printf("r0 ");mpz_out_str(stdout,10,r0);printf("\n");
97  //printf("t1 ");mpz_out_str(stdout,10,t1);printf("\n");
98  //printf("r1 ");mpz_out_str(stdout,10,r1);printf("\n");
99  }
100  if(mpz_cmp_si(t1,0)<0)
101  {
102  //t1=-1*t1;
103  mpz_neg(t1,t1);
104  //r1=-1*r1;
105  mpz_neg(r1,r1);
106  }
107  mpz_gcd(tmp,r1,t1);
108  if((mpz_cmp(t1,N)<0) && (mpz_cmp_si(tmp,1)==0))
109  {
110  // printf("ici \n");
111  mpq_set_num(res,r1);
112  mpq_set_den(res,t1);
113  ret=0;
114  //mpq_out_str(stdout,10,res);printf("\n");
115  }
116  else
117  {
118  mpq_set_si(res,-0,1);
119  ret=1;
120  //printf("je suis pas arrive a reconstruire un nombre\n");
121  //mpz_out_str(stdout,10,m);
122  //printf("\n");
123  }
124  mpz_clear(N);
125  mpz_clear(r0);
126  mpz_clear(r1);
127  mpz_clear(t0);
128  mpz_clear(t1);
129  mpz_clear(q);
130  mpz_clear(tmp);
131  return ret;
132 }
133 
int reconstruct(mpq_t res, const mpz_t u, const mpz_t m)
Definition: reconstruct.hpp:58
void swap(X *a, X *b)
Definition: Swap.H:5
MSKstreamtypee MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t a
Definition: mosek.h:3833
void chinese(mpz_t res, const mpz_t r1, const mpz_t m1, const mpz_t r2, const mpz_t m2)
Definition: reconstruct.hpp:1
Home  |  Download & InstallContributions