1/* mpz_divexact_gcd -- exact division optimized for GCDs.
2
3   THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE AND ARE ALMOST CERTAIN TO
4   BE SUBJECT TO INCOMPATIBLE CHANGES IN FUTURE GNU MP RELEASES.
5
6Copyright 2000, 2005, 2011, 2012 Free Software Foundation, Inc.
7
8This file is part of the GNU MP Library.
9
10The GNU MP Library is free software; you can redistribute it and/or modify
11it under the terms of either:
12
13  * the GNU Lesser General Public License as published by the Free
14    Software Foundation; either version 3 of the License, or (at your
15    option) any later version.
16
17or
18
19  * the GNU General Public License as published by the Free Software
20    Foundation; either version 2 of the License, or (at your option) any
21    later version.
22
23or both in parallel, as here.
24
25The GNU MP Library is distributed in the hope that it will be useful, but
26WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
28for more details.
29
30You should have received copies of the GNU General Public License and the
31GNU Lesser General Public License along with the GNU MP Library.  If not,
32see https://www.gnu.org/licenses/.  */
33
34#include "gmp-impl.h"
35#include "longlong.h"
36
37
38/* Set q to a/d, expecting d to be from a GCD and therefore usually small.
39
40   The distribution of GCDs of random numbers can be found in Knuth volume 2
41   section 4.5.2 theorem D.
42
43            GCD     chance
44             1       60.8%
45            2^k      20.2%     (1<=k<32)
46           3*2^k      9.0%     (1<=k<32)
47           other     10.1%
48
49   Only the low limb is examined for optimizations, since GCDs bigger than
50   2^32 (or 2^64) will occur very infrequently.
51
52   Future: This could change to an mpn_divexact_gcd, possibly partly
53   inlined, if/when the relevant mpq functions change to an mpn based
54   implementation.  */
55
56
57#if GMP_NUMB_BITS % 2 == 0
58static void
59mpz_divexact_by3 (mpz_ptr q, mpz_srcptr a)
60{
61  mp_size_t  size = SIZ(a);
62  mp_size_t  abs_size = ABS(size);
63  mp_ptr     qp;
64
65  qp = MPZ_REALLOC (q, abs_size);
66
67  mpn_bdiv_dbm1 (qp, PTR(a), abs_size, GMP_NUMB_MASK / 3);
68
69  abs_size -= (qp[abs_size-1] == 0);
70  SIZ(q) = (size>0 ? abs_size : -abs_size);
71}
72#endif
73
74#if GMP_NUMB_BITS % 4 == 0
75static void
76mpz_divexact_by5 (mpz_ptr q, mpz_srcptr a)
77{
78  mp_size_t  size = SIZ(a);
79  mp_size_t  abs_size = ABS(size);
80  mp_ptr     qp;
81
82  qp = MPZ_REALLOC (q, abs_size);
83
84  mpn_bdiv_dbm1 (qp, PTR(a), abs_size, GMP_NUMB_MASK / 5);
85
86  abs_size -= (qp[abs_size-1] == 0);
87  SIZ(q) = (size>0 ? abs_size : -abs_size);
88}
89#endif
90
91static void
92mpz_divexact_limb (mpz_ptr q, mpz_srcptr a, mp_limb_t d)
93{
94  mp_size_t  size = SIZ(a);
95  mp_size_t  abs_size = ABS(size);
96  mp_ptr     qp;
97
98  qp = MPZ_REALLOC (q, abs_size);
99
100  MPN_DIVREM_OR_DIVEXACT_1 (qp, PTR(a), abs_size, d);
101
102  abs_size -= (qp[abs_size-1] == 0);
103  SIZ(q) = (size > 0 ? abs_size : -abs_size);
104}
105
106void
107mpz_divexact_gcd (mpz_ptr q, mpz_srcptr a, mpz_srcptr d)
108{
109  ASSERT (mpz_sgn (d) > 0);
110
111  if (SIZ(a) == 0)
112    {
113      SIZ(q) = 0;
114      return;
115    }
116
117  if (SIZ(d) == 1)
118    {
119      mp_limb_t  dl = PTR(d)[0];
120      int        twos;
121
122      if ((dl & 1) == 0)
123	{
124	  count_trailing_zeros (twos, dl);
125	  dl >>= twos;
126	  mpz_tdiv_q_2exp (q, a, twos);
127	  a = q;
128	}
129
130      if (dl == 1)
131	{
132	  if (q != a)
133	    mpz_set (q, a);
134	  return;
135	}
136#if GMP_NUMB_BITS % 2 == 0
137      if (dl == 3)
138	{
139	  mpz_divexact_by3 (q, a);
140	  return;
141	}
142#endif
143#if GMP_NUMB_BITS % 4 == 0
144      if (dl == 5)
145	{
146	  mpz_divexact_by5 (q, a);
147	  return;
148	}
149#endif
150
151      mpz_divexact_limb (q, a, dl);
152      return;
153    }
154
155  mpz_divexact (q, a, d);
156}
157