1 /* gcdext_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 
7 Copyright 2003, 2004, 2005, 2008 Free Software Foundation, Inc.
8 
9 This file is part of the GNU MP Library.
10 
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15 
16 The GNU MP Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
19 License for more details.
20 
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
23 
24 #include "gmp.h"
25 #include "gmp-impl.h"
26 #include "longlong.h"
27 
28 /* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or
29    b is small, or the difference is small. Perform one subtraction
30    followed by one division. If the gcd is found, stores it in gp and
31    *gn, and returns zero. Otherwise, compute the reduced a and b,
32    return the new size, and cofactors. */
33 
34 /* Temporary storage: Needs n limbs for the quotient, at qp. tp must
35    point to an area large enough for the resulting cofactor, plus one
36    limb extra. All in all, 2N + 1 if N is a bound for both inputs and
37    outputs. */
38 mp_size_t
39 mpn_gcdext_subdiv_step (mp_ptr gp, mp_size_t *gn, mp_ptr up, mp_size_t *usizep,
40 			mp_ptr ap, mp_ptr bp, mp_size_t n,
41 			mp_ptr u0, mp_ptr u1, mp_size_t *unp,
42 			mp_ptr qp, mp_ptr tp)
43 {
44   mp_size_t an, bn, un;
45   mp_size_t qn;
46   mp_size_t u0n;
47 
48   int swapped;
49 
50   an = bn = n;
51 
52   ASSERT (an > 0);
53   ASSERT (ap[an-1] > 0 || bp[an-1] > 0);
54 
55   MPN_NORMALIZE (ap, an);
56   MPN_NORMALIZE (bp, bn);
57 
58   un = *unp;
59 
60   swapped = 0;
61 
62   if (UNLIKELY (an == 0))
63     {
64     return_b:
65       MPN_COPY (gp, bp, bn);
66       *gn = bn;
67 
68       MPN_NORMALIZE (u0, un);
69       MPN_COPY (up, u0, un);
70 
71       *usizep = swapped ? un : -un;
72 
73       return 0;
74     }
75   else if (UNLIKELY (bn == 0))
76     {
77     return_a:
78       MPN_COPY (gp, ap, an);
79       *gn = an;
80 
81       MPN_NORMALIZE (u1, un);
82       MPN_COPY (up, u1, un);
83 
84       *usizep = swapped ? -un : un;
85 
86       return 0;
87     }
88 
89   /* Arrange so that a > b, subtract an -= bn, and maintain
90      normalization. */
91   if (an < bn)
92     {
93       MPN_PTR_SWAP (ap, an, bp, bn);
94       MP_PTR_SWAP (u0, u1);
95       swapped ^= 1;
96     }
97   else if (an == bn)
98     {
99       int c;
100       MPN_CMP (c, ap, bp, an);
101       if (UNLIKELY (c == 0))
102 	goto return_a;
103       else if (c < 0)
104 	{
105 	  MP_PTR_SWAP (ap, bp);
106 	  MP_PTR_SWAP (u0, u1);
107 	  swapped ^= 1;
108 	}
109     }
110   /* Reduce a -= b, u1 += u0 */
111   ASSERT_NOCARRY (mpn_sub (ap, ap, an, bp, bn));
112   MPN_NORMALIZE (ap, an);
113   ASSERT (an > 0);
114 
115   u1[un] = mpn_add_n (u1, u1, u0, un);
116   un += (u1[un] > 0);
117 
118   /* Arrange so that a > b, and divide a = q b + r */
119   if (an < bn)
120     {
121       MPN_PTR_SWAP (ap, an, bp, bn);
122       MP_PTR_SWAP (u0, u1);
123       swapped ^= 1;
124     }
125   else if (an == bn)
126     {
127       int c;
128       MPN_CMP (c, ap, bp, an);
129       if (UNLIKELY (c == 0))
130 	goto return_a;
131       else if (c < 0)
132 	{
133 	  MP_PTR_SWAP (ap, bp);
134 	  MP_PTR_SWAP (u0, u1);
135 	  swapped ^= 1;
136 	}
137     }
138 
139   /* Reduce a -= q b, u1 += q u0 */
140   qn = an - bn + 1;
141   mpn_tdiv_qr (qp, ap, 0, ap, an, bp, bn);
142 
143   if (mpn_zero_p (ap, bn))
144     goto return_b;
145 
146   n = bn;
147 
148   /* Update u1 += q u0 */
149   u0n = un;
150   MPN_NORMALIZE (u0, u0n);
151 
152   if (u0n > 0)
153     {
154       qn -= (qp[qn - 1] == 0);
155 
156       if (qn > u0n)
157 	mpn_mul (tp, qp, qn, u0, u0n);
158       else
159 	mpn_mul (tp, u0, u0n, qp, qn);
160 
161       if (qn + u0n > un)
162 	{
163 	  ASSERT_NOCARRY (mpn_add (u1, tp, qn + u0n, u1, un));
164 	  un = qn + u0n;
165 	  un -= (u1[un-1] == 0);
166 	}
167       else
168 	{
169 	  u1[un] = mpn_add (u1, u1, un, tp, qn + u0n);
170 	  un += (u1[un] > 0);
171 	}
172     }
173 
174   *unp = un;
175   return n;
176 }
177