1 /* Test mpn_gcdext.
2
3 Copyright 2010 William Hart
4
5 This file is part of the MPIR Library.
6
7 The MPIR Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 2.1 of the License, or (at your
10 option) any later version.
11
12 The MPIR Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with the MPIR Library; see the file COPYING.LIB. If not, write to
19 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
20 MA 02110-1301, USA. */
21
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include "mpir.h"
26 #include "gmp-impl.h"
27 #include "longlong.h"
28 #include "tests.h"
29
30 #if defined(_MSC_VER) || defined(__MINGW32__)
31 #define random rand
32 #endif
33
34 #define MAX_LIMBS 2000
35 #define ITERS 1000
36
check_normalisation(const mpz_t G,const mpz_t S,const mpz_t U,const mpz_t V,mpz_t T,mpz_t T1,mpz_t T2)37 void check_normalisation(const mpz_t G, const mpz_t S, const mpz_t U, const mpz_t V,
38 mpz_t T, mpz_t T1, mpz_t T2)
39 {
40 /* compute T */
41 mpz_mul(T1, U, S);
42 mpz_sub(T, G, T1);
43 mpz_divexact(T, T, V);
44
45 /* T2 = V/(2G) */
46 mpz_divexact(T2, V, G);
47 mpz_cdiv_q_ui(T2, T2, 2);
48
49 /* T1 = abs(S) */
50 mpz_abs(T1, S);
51
52 if ((mpz_cmp(T1, T2) >= 0) && (mpz_cmp_ui(S, 1) != 0))
53 {
54 printf("FAIL\n");
55 gmp_printf("U = %Zd\n", U);
56 gmp_printf("V = %Zd\n", V);
57 gmp_printf("G = %Zd\n", G);
58 gmp_printf("S = %Zd\n", S);
59 gmp_printf("T = %Zd\n", T);
60 abort();
61 }
62 }
63
64 /* Check extended gcd routine. */
check_gcdext(void)65 void check_gcdext(void)
66 {
67 mp_limb_t up[MAX_LIMBS + 1];
68 mp_limb_t vp[MAX_LIMBS + 1];
69 mp_limb_t up2[MAX_LIMBS];
70 mp_limb_t vp2[MAX_LIMBS];
71 mp_limb_t sp[MAX_LIMBS + 1];
72 mp_limb_t tp[MAX_LIMBS];
73 mp_limb_t gp[MAX_LIMBS + 1];
74
75 mp_size_t un, vn, sn, tn, bits, gn;
76 mp_bitcnt_t u_bits, v_bits, g_bits;
77 long i, j;
78 mpz_t U, V, S, T, G, T1, T2;
79
80 gmp_randstate_t rands;
81 gmp_randinit_default(rands);
82
83 mpz_init(T);
84 mpz_init(T1);
85 mpz_init(T2);
86
87 for (i = 0; i < ITERS; i++)
88 {
89 u_bits = (random()%(GMP_LIMB_BITS*MAX_LIMBS)) + 1;
90 v_bits = (random()%u_bits) + 1;
91 g_bits = (random()%v_bits) + 1;
92
93 gn = (g_bits + GMP_LIMB_BITS - 1)/GMP_LIMB_BITS;
94 sn = (u_bits - g_bits + GMP_LIMB_BITS - 1)/GMP_LIMB_BITS;
95 tn = (v_bits - g_bits + GMP_LIMB_BITS - 1)/GMP_LIMB_BITS;
96
97 /* generate quasi gcd gp*/
98 do{
99 mpn_urandomb(gp, rands, g_bits);
100 gn = (g_bits + GMP_LIMB_BITS - 1)/GMP_LIMB_BITS;
101 MPN_NORMALIZE(gp, gn);
102 } while (gn == 0);
103
104 /*
105 generate random {sp, sn}
106 set {up, un} = {gp, gn}*{sp, sn}
107 */
108 if (u_bits > g_bits)
109 {
110 mpn_urandomb(sp, rands, u_bits - g_bits);
111 up[0] |= (mp_limb_t) 1;
112
113 if (gn >= sn)
114 mpn_mul(up, gp, gn, sp, sn);
115 else
116 mpn_mul(up, sp, sn, gp, gn);
117 un = gn + sn;
118 MPN_NORMALIZE(up, un);
119 } else
120 {
121 MPN_COPY(up, gp, gn);
122 un = gn;
123 MPN_NORMALIZE(up, un);
124 }
125
126 do
127 {
128 /*
129 generate random odd {tp, tn}
130 set {vp, vn} = {gp, gn}*{tp, tn}
131 */
132 if (v_bits > g_bits)
133 {
134 mpn_urandomb(tp, rands, v_bits - g_bits);
135 tp[0] |= (mp_limb_t) 1;
136
137 if (gn >= tn)
138 mpn_mul(vp, gp, gn, tp, tn);
139 else
140 mpn_mul(vp, tp, tn, gp, gn);
141 vn = gn + tn;
142 MPN_NORMALIZE(vp, vn);
143 } else
144 {
145 MPN_COPY(vp, gp, gn);
146 vn = gn;
147 MPN_NORMALIZE(vp, vn);
148 }
149 } while ((un == vn) && (mpn_cmp(up, vp, un) < 0));
150
151 /* Save a copy of up and vp */
152 MPN_COPY(up2, up, un);
153 MPN_COPY(vp2, vp, vn);
154 if(un<vn)continue;
155 gn = mpn_gcdext(gp, sp, &sn, up, un, vp, vn);
156
157 U->_mp_d = up2;
158 U->_mp_size = un;
159 V->_mp_d = vp2;
160 V->_mp_size = vn;
161 S->_mp_d = sp;
162 S->_mp_size = sn;
163 G->_mp_d = gp;
164 G->_mp_size = gn;
165
166 check_normalisation(G, S, U, V, T, T1, T2);
167 }
168
169 mpz_clear(T);
170 mpz_clear(T1);
171 mpz_clear(T2);
172 gmp_randclear(rands);
173 }
174
main(void)175 int main(void)
176 {
177 tests_start();
178 check_gcdext();
179 tests_end();
180
181 return 0;
182 }
183