1/* gcd_subdiv_step.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, 2010, 2011 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 <stdlib.h>		/* for NULL */
36
37#include "gmp-impl.h"
38#include "longlong.h"
39
40/* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or
41   b is small, or the difference is small. Perform one subtraction
42   followed by one division. The normal case is to compute the reduced
43   a and b, and return the new size.
44
45   If s == 0 (used for gcd and gcdext), returns zero if the gcd is
46   found.
47
48   If s > 0, don't reduce to size <= s, and return zero if no
49   reduction is possible (if either a, b or |a-b| is of size <= s). */
50
51/* The hook function is called as
52
53     hook(ctx, gp, gn, qp, qn, d)
54
55   in the following cases:
56
57   + If A = B at the start, G is the gcd, Q is NULL, d = -1.
58
59   + If one input is zero at the start, G is the gcd, Q is NULL,
60     d = 0 if A = G and d = 1 if B = G.
61
62   Otherwise, if d = 0 we have just subtracted a multiple of A from B,
63   and if d = 1 we have subtracted a multiple of B from A.
64
65   + If A = B after subtraction, G is the gcd, Q is NULL.
66
67   + If we get a zero remainder after division, G is the gcd, Q is the
68     quotient.
69
70   + Otherwise, G is NULL, Q is the quotient (often 1).
71
72 */
73
74mp_size_t
75mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, mp_size_t s,
76		     gcd_subdiv_step_hook *hook, void *ctx,
77		     mp_ptr tp)
78{
79  static const mp_limb_t one = CNST_LIMB(1);
80  mp_size_t an, bn, qn;
81
82  int swapped;
83
84  ASSERT (n > 0);
85  ASSERT (ap[n-1] > 0 || bp[n-1] > 0);
86
87  an = bn = n;
88  MPN_NORMALIZE (ap, an);
89  MPN_NORMALIZE (bp, bn);
90
91  swapped = 0;
92
93  /* Arrange so that a < b, subtract b -= a, and maintain
94     normalization. */
95  if (an == bn)
96    {
97      int c;
98      MPN_CMP (c, ap, bp, an);
99      if (UNLIKELY (c == 0))
100	{
101	  /* For gcdext, return the smallest of the two cofactors, so
102	     pass d = -1. */
103	  if (s == 0)
104	    hook (ctx, ap, an, NULL, 0, -1);
105	  return 0;
106	}
107      else if (c > 0)
108	{
109	  MP_PTR_SWAP (ap, bp);
110	  swapped ^= 1;
111	}
112    }
113  else
114    {
115      if (an > bn)
116	{
117	  MPN_PTR_SWAP (ap, an, bp, bn);
118	  swapped ^= 1;
119	}
120    }
121  if (an <= s)
122    {
123      if (s == 0)
124	hook (ctx, bp, bn, NULL, 0, swapped ^ 1);
125      return 0;
126    }
127
128  ASSERT_NOCARRY (mpn_sub (bp, bp, bn, ap, an));
129  MPN_NORMALIZE (bp, bn);
130  ASSERT (bn > 0);
131
132  if (bn <= s)
133    {
134      /* Undo subtraction. */
135      mp_limb_t cy = mpn_add (bp, ap, an, bp, bn);
136      if (cy > 0)
137	bp[an] = cy;
138      return 0;
139    }
140
141  /* Arrange so that a < b */
142  if (an == bn)
143    {
144      int c;
145      MPN_CMP (c, ap, bp, an);
146      if (UNLIKELY (c == 0))
147	{
148	  if (s > 0)
149	    /* Just record subtraction and return */
150	    hook (ctx, NULL, 0, &one, 1, swapped);
151	  else
152	    /* Found gcd. */
153	    hook (ctx, bp, bn, NULL, 0, swapped);
154	  return 0;
155	}
156
157      hook (ctx, NULL, 0, &one, 1, swapped);
158
159      if (c > 0)
160	{
161	  MP_PTR_SWAP (ap, bp);
162	  swapped ^= 1;
163	}
164    }
165  else
166    {
167      hook (ctx, NULL, 0, &one, 1, swapped);
168
169      if (an > bn)
170	{
171	  MPN_PTR_SWAP (ap, an, bp, bn);
172	  swapped ^= 1;
173	}
174    }
175
176  mpn_tdiv_qr (tp, bp, 0, bp, bn, ap, an);
177  qn = bn - an + 1;
178  bn = an;
179  MPN_NORMALIZE (bp, bn);
180
181  if (UNLIKELY (bn <= s))
182    {
183      if (s == 0)
184	{
185	  hook (ctx, ap, an, tp, qn, swapped);
186	  return 0;
187	}
188
189      /* Quotient is one too large, so decrement it and add back A. */
190      if (bn > 0)
191	{
192	  mp_limb_t cy = mpn_add (bp, ap, an, bp, bn);
193	  if (cy)
194	    bp[an++] = cy;
195	}
196      else
197	MPN_COPY (bp, ap, an);
198
199      MPN_DECR_U (tp, qn, 1);
200    }
201
202  hook (ctx, NULL, 0, tp, qn, swapped);
203  return an;
204}
205