1 /* Test mpn_inv_divappr_q_n.
2
3 Copyright 2002 Free Software Foundation, Inc.
4 Copyright 2009, 2010 William Hart
5
6 This file is part of the MPIR Library.
7
8 The MPIR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 2.1 of the License, or (at your
11 option) any later version.
12
13 The MPIR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the MPIR Library; see the file COPYING.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "mpir.h"
27 #include "gmp-impl.h"
28 #include "longlong.h"
29 #include "tests.h"
30
31 #if defined(_MSC_VER) || defined(__MINGW32__)
32 #define random rand
33 #endif
34
35 #define MAX_LIMBS 2000
36 #define ITERS 500
37
38 /* Check divide and conquer division routine. */
39 void
check_inv_divappr_q_n(void)40 check_inv_divappr_q_n (void)
41 {
42 mp_limb_t np[2*MAX_LIMBS];
43 mp_limb_t np2[2*MAX_LIMBS];
44 mp_limb_t rp[2*MAX_LIMBS];
45 mp_limb_t dp[MAX_LIMBS];
46 mp_limb_t qp[MAX_LIMBS];
47 mp_limb_t dip[MAX_LIMBS];
48
49 mp_size_t nn, rn, dn, qn;
50
51 gmp_randstate_t rands;
52
53 int i, j, s;
54 gmp_randinit_default(rands);
55
56 for (i = 0; i < ITERS; i++)
57 {
58 dn = (random() % (MAX_LIMBS - 6)) + 6;
59 nn = 2*dn;
60
61 mpn_rrandom (np, rands, nn);
62 mpn_rrandom (dp, rands, dn);
63 dp[dn-1] |= GMP_LIMB_HIGHBIT;
64
65 MPN_COPY(np2, np, nn);
66
67 mpn_invert(dip, dp, dn);
68
69 qn = nn - dn + 1;
70
71 qp[qn - 1] = mpn_inv_divappr_q_n(qp, np, dp, dn, dip);
72
73 MPN_NORMALIZE(qp, qn);
74
75 if (qn)
76 {
77 if (qn >= dn) mpn_mul(rp, qp, qn, dp, dn);
78 else mpn_mul(rp, dp, dn, qp, qn);
79
80 rn = dn + qn;
81 MPN_NORMALIZE(rp, rn);
82
83 s = (rn < nn) ? -1 : (rn > nn) ? 1 : mpn_cmp(rp, np2, nn);
84 if (s <= 0)
85 {
86 mpn_sub(rp, np2, nn, rp, rn);
87 rn = nn;
88 MPN_NORMALIZE(rp, rn);
89 } else
90 {
91 mpn_sub(rp, rp, rn, np2, nn);
92 MPN_NORMALIZE(rp, rn);
93 }
94 } else
95 {
96 rn = nn;
97 MPN_COPY(rp, np, nn);
98 }
99
100 s = (rn < dn) ? -1 : (rn > dn) ? 1 : mpn_cmp(rp, dp, dn);
101 if (s >= 0)
102 {
103 printf ("failed:\n");
104 printf ("nn = %lu, dn = %lu, qn = %lu, rn = %lu\n\n", nn, dn, qn, rn);
105 gmp_printf (" np: %Nx\n\n", np2, nn);
106 gmp_printf (" dp: %Nx\n\n", dp, dn);
107 gmp_printf (" qp: %Nx\n\n", qp, qn);
108 gmp_printf (" rp: %Nx\n\n", rp, rn);
109 abort ();
110 }
111 }
112
113 gmp_randclear(rands);
114 }
115
116 int
main(void)117 main (void)
118 {
119 tests_start ();
120
121 check_inv_divappr_q_n ();
122
123 tests_end ();
124 exit (0);
125 }
126