1 /*
2
3 Copyright 2012, 2016 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library test suite.
6
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 Public License for more details.
16
17 You should have received a copy of the GNU General Public License along with
18 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
19
20 #include <assert.h>
21 #include <limits.h>
22 #include <stdlib.h>
23 #include <stdio.h>
24
25 #include "testutils.h"
26
27 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
28
29 #define COUNT 10000
30
31 static void
test_2by1(const mpz_t u)32 test_2by1(const mpz_t u)
33 {
34 mpz_t m, p, t;
35
36 mpz_init (m);
37 mpz_init (p);
38 mpz_init (t);
39
40 assert (mpz_size (u) == 1);
41
42 mpz_set_ui (m, mpn_invert_limb (u->_mp_d[0]));
43 mpz_setbit (m, GMP_LIMB_BITS);
44
45 mpz_mul (p, m, u);
46
47 mpz_set_ui (t, 0);
48 mpz_setbit (t, 2* GMP_LIMB_BITS);
49 mpz_sub (t, t, p);
50
51 /* Should have 0 < B^2 - m u <= u */
52 if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0)
53 {
54 fprintf (stderr, "mpn_invert_limb failed:\n");
55 dump ("u", u);
56 dump ("m", m);
57 dump ("p", p);
58 dump ("t", t);
59 abort ();
60 }
61 mpz_clear (m);
62 mpz_clear (p);
63 mpz_clear (t);
64 }
65
66 static void
test_3by2(const mpz_t u)67 test_3by2(const mpz_t u)
68 {
69 mpz_t m, p, t;
70
71 mpz_init (m);
72 mpz_init (p);
73 mpz_init (t);
74
75 assert (mpz_size (u) == 2);
76
77 mpz_set_ui (m, mpn_invert_3by2 (u->_mp_d[1], u[0]._mp_d[0]));
78
79 mpz_setbit (m, GMP_LIMB_BITS);
80
81 mpz_mul (p, m, u);
82
83 mpz_set_ui (t, 0);
84 mpz_setbit (t, 3 * GMP_LIMB_BITS);
85 mpz_sub (t, t, p);
86
87 /* Should have 0 < B^3 - m u <= u */
88 if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0)
89 {
90 fprintf (stderr, "mpn_invert_3by2 failed:\n");
91 dump ("u", u);
92 dump ("m", m);
93 dump ("p", p);
94 dump ("t", t);
95 abort ();
96 }
97 mpz_clear (m);
98 mpz_clear (p);
99 mpz_clear (t);
100 }
101
102 void
testmain(int argc,char ** argv)103 testmain (int argc, char **argv)
104 {
105 unsigned i;
106 mpz_t u, m, p, t;
107
108 mpz_init (u);
109 mpz_init (m);
110 mpz_init (p);
111 mpz_init (t);
112
113 /* These values trigger 32-bit overflow of ql in mpn_invert_3by2. */
114 if (GMP_LIMB_BITS == 64)
115 {
116 mpz_set_str (u, "80007fff3ffe0000", 16);
117 test_2by1 (u);
118 mpz_set_str (u, "80007fff3ffe000000000000000003ff", 16);
119 test_3by2 (u);
120 }
121
122 for (i = 0; i < COUNT; i++)
123 {
124 mini_urandomb (u, GMP_LIMB_BITS);
125 mpz_setbit (u, GMP_LIMB_BITS -1);
126
127 test_2by1 (u);
128 }
129
130 for (i = 0; i < COUNT; i++)
131 {
132 mini_urandomb (u, 2*GMP_LIMB_BITS);
133 mpz_setbit (u, 2*GMP_LIMB_BITS -1);
134
135 test_3by2 (u);
136 }
137
138 mpz_clear (u);
139 }
140