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 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 the GNU Lesser General Public License as published by
12the Free Software Foundation; either version 3 of the License, or (at your
13option) any later version.
14
15The GNU MP Library is distributed in the hope that it will be useful, but
16WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
18License for more details.
19
20You should have received a copy of the GNU Lesser General Public License
21along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
22
23#include "gmp.h"
24#include "gmp-impl.h"
25#include "longlong.h"
26
27
28/* Set q to a/d, expecting d to be from a GCD and therefore usually small.
29
30   The distribution of GCDs of random numbers can be found in Knuth volume 2
31   section 4.5.2 theorem D.
32
33            GCD     chance
34             1       60.8%
35            2^k      20.2%     (1<=k<32)
36           3*2^k      9.0%     (1<=k<32)
37           other     10.1%
38
39   Only the low limb is examined for optimizations, since GCDs bigger than
40   2^32 (or 2^64) will occur very infrequently.
41
42   Future: This could change to an mpn_divexact_gcd, possibly partly
43   inlined, if/when the relevant mpq functions change to an mpn based
44   implementation.  */
45
46
47static void
48mpz_divexact_by3 (mpz_ptr q, mpz_srcptr a)
49{
50  mp_size_t  size = SIZ(a);
51  if (size == 0)
52    {
53      SIZ(q) = 0;
54      return;
55    }
56  else
57    {
58      mp_size_t  abs_size = ABS(size);
59      mp_ptr     qp;
60
61      MPZ_REALLOC (q, abs_size);
62
63      qp = PTR(q);
64      mpn_divexact_by3 (qp, PTR(a), abs_size);
65
66      abs_size -= (qp[abs_size-1] == 0);
67      SIZ(q) = (size>0 ? abs_size : -abs_size);
68    }
69}
70
71void
72mpz_divexact_gcd (mpz_ptr q, mpz_srcptr a, mpz_srcptr d)
73{
74  ASSERT (mpz_sgn (d) > 0);
75
76  if (SIZ(d) == 1)
77    {
78      mp_limb_t  dl = PTR(d)[0];
79      int        twos;
80
81      if (dl == 1)
82        {
83          if (q != a)
84            mpz_set (q, a);
85          return;
86        }
87      if (dl == 3)
88        {
89          mpz_divexact_by3 (q, a);
90          return;
91        }
92
93      count_trailing_zeros (twos, dl);
94      dl >>= twos;
95
96      if (dl == 1)
97        {
98          mpz_tdiv_q_2exp (q, a, twos);
99          return;
100        }
101      if (dl == 3)
102        {
103          mpz_tdiv_q_2exp (q, a, twos);
104          mpz_divexact_by3 (q, q);
105          return;
106        }
107    }
108
109  mpz_divexact (q, a, d);
110}
111