1/* hgcd.c.
2
3   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
4   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
5   GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
6
7Copyright 2003-2005, 2008, 2011, 2012 Free Software Foundation, Inc.
8
9This file is part of the GNU MP Library.
10
11The GNU MP Library is free software; you can redistribute it and/or modify
12it under the terms of either:
13
14  * the GNU Lesser General Public License as published by the Free
15    Software Foundation; either version 3 of the License, or (at your
16    option) any later version.
17
18or
19
20  * the GNU General Public License as published by the Free Software
21    Foundation; either version 2 of the License, or (at your option) any
22    later version.
23
24or both in parallel, as here.
25
26The GNU MP Library is distributed in the hope that it will be useful, but
27WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
29for more details.
30
31You should have received copies of the GNU General Public License and the
32GNU Lesser General Public License along with the GNU MP Library.  If not,
33see https://www.gnu.org/licenses/.  */
34
35#include "gmp-impl.h"
36#include "longlong.h"
37
38
39/* Size analysis for hgcd:
40
41   For the recursive calls, we have n1 <= ceil(n / 2). Then the
42   storage need is determined by the storage for the recursive call
43   computing M1, and hgcd_matrix_adjust and hgcd_matrix_mul calls that use M1
44   (after this, the storage needed for M1 can be recycled).
45
46   Let S(r) denote the required storage. For M1 we need 4 * (ceil(n1/2) + 1)
47   = 4 * (ceil(n/4) + 1), for the hgcd_matrix_adjust call, we need n + 2,
48   and for the hgcd_matrix_mul, we may need 3 ceil(n/2) + 8. In total,
49   4 * ceil(n/4) + 3 ceil(n/2) + 12 <= 10 ceil(n/4) + 12.
50
51   For the recursive call, we need S(n1) = S(ceil(n/2)).
52
53   S(n) <= 10*ceil(n/4) + 12 + S(ceil(n/2))
54	<= 10*(ceil(n/4) + ... + ceil(n/2^(1+k))) + 12k + S(ceil(n/2^k))
55	<= 10*(2 ceil(n/4) + k) + 12k + S(ceil(n/2^k))
56	<= 20 ceil(n/4) + 22k + S(ceil(n/2^k))
57*/
58
59mp_size_t
60mpn_hgcd_itch (mp_size_t n)
61{
62  unsigned k;
63  int count;
64  mp_size_t nscaled;
65
66  if (BELOW_THRESHOLD (n, HGCD_THRESHOLD))
67    return n;
68
69  /* Get the recursion depth. */
70  nscaled = (n - 1) / (HGCD_THRESHOLD - 1);
71  count_leading_zeros (count, nscaled);
72  k = GMP_LIMB_BITS - count;
73
74  return 20 * ((n+3) / 4) + 22 * k + HGCD_THRESHOLD;
75}
76
77/* Reduces a,b until |a-b| fits in n/2 + 1 limbs. Constructs matrix M
78   with elements of size at most (n+1)/2 - 1. Returns new size of a,
79   b, or zero if no reduction is possible. */
80
81mp_size_t
82mpn_hgcd (mp_ptr ap, mp_ptr bp, mp_size_t n,
83	  struct hgcd_matrix *M, mp_ptr tp)
84{
85  mp_size_t s = n/2 + 1;
86
87  mp_size_t nn;
88  int success = 0;
89
90  if (n <= s)
91    /* Happens when n <= 2, a fairly uninteresting case but exercised
92       by the random inputs of the testsuite. */
93    return 0;
94
95  ASSERT ((ap[n-1] | bp[n-1]) > 0);
96
97  ASSERT ((n+1)/2 - 1 < M->alloc);
98
99  if (ABOVE_THRESHOLD (n, HGCD_THRESHOLD))
100    {
101      mp_size_t n2 = (3*n)/4 + 1;
102      mp_size_t p = n/2;
103
104      nn = mpn_hgcd_reduce (M, ap, bp, n, p, tp);
105      if (nn)
106	{
107	  n = nn;
108	  success = 1;
109	}
110
111      /* NOTE: It appears this loop never runs more than once (at
112	 least when not recursing to hgcd_appr). */
113      while (n > n2)
114	{
115	  /* Needs n + 1 storage */
116	  nn = mpn_hgcd_step (n, ap, bp, s, M, tp);
117	  if (!nn)
118	    return success ? n : 0;
119
120	  n = nn;
121	  success = 1;
122	}
123
124      if (n > s + 2)
125	{
126	  struct hgcd_matrix M1;
127	  mp_size_t scratch;
128
129	  p = 2*s - n + 1;
130	  scratch = MPN_HGCD_MATRIX_INIT_ITCH (n-p);
131
132	  mpn_hgcd_matrix_init(&M1, n - p, tp);
133
134	  /* FIXME: Should use hgcd_reduce, but that may require more
135	     scratch space, which requires review. */
136
137	  nn = mpn_hgcd (ap + p, bp + p, n - p, &M1, tp + scratch);
138	  if (nn > 0)
139	    {
140	      /* We always have max(M) > 2^{-(GMP_NUMB_BITS + 1)} max(M1) */
141	      ASSERT (M->n + 2 >= M1.n);
142
143	      /* Furthermore, assume M ends with a quotient (1, q; 0, 1),
144		 then either q or q + 1 is a correct quotient, and M1 will
145		 start with either (1, 0; 1, 1) or (2, 1; 1, 1). This
146		 rules out the case that the size of M * M1 is much
147		 smaller than the expected M->n + M1->n. */
148
149	      ASSERT (M->n + M1.n < M->alloc);
150
151	      /* Needs 2 (p + M->n) <= 2 (2*s - n2 + 1 + n2 - s - 1)
152		 = 2*s <= 2*(floor(n/2) + 1) <= n + 2. */
153	      n = mpn_hgcd_matrix_adjust (&M1, p + nn, ap, bp, p, tp + scratch);
154
155	      /* We need a bound for of M->n + M1.n. Let n be the original
156		 input size. Then
157
158		 ceil(n/2) - 1 >= size of product >= M.n + M1.n - 2
159
160		 and it follows that
161
162		 M.n + M1.n <= ceil(n/2) + 1
163
164		 Then 3*(M.n + M1.n) + 5 <= 3 * ceil(n/2) + 8 is the
165		 amount of needed scratch space. */
166	      mpn_hgcd_matrix_mul (M, &M1, tp + scratch);
167	      success = 1;
168	    }
169	}
170    }
171
172  for (;;)
173    {
174      /* Needs s+3 < n */
175      nn = mpn_hgcd_step (n, ap, bp, s, M, tp);
176      if (!nn)
177	return success ? n : 0;
178
179      n = nn;
180      success = 1;
181    }
182}
183