1/*
2
3Copyright 2012, 2016 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library test suite.
6
7The GNU MP Library test suite is free software; you can redistribute it
8and/or modify it under the terms of the GNU General Public License as
9published by the Free Software Foundation; either version 3 of the License,
10or (at your option) any later version.
11
12The GNU MP Library test suite is distributed in the hope that it will be
13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15Public License for more details.
16
17You should have received a copy of the GNU General Public License along with
18the 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
31static void
32test_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
66static void
67test_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
102void
103testmain (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