contrib/gmp/mpn/generic/gcdext_subdiv_step.c

   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, 2009 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       MPN_COPY (gp, ap, an);
  78       *gn = an;
  79
  80       MPN_NORMALIZE (u1, un);
  81       MPN_COPY (up, u1, un);
  82
  83       *usizep = swapped ? -un : un;
  84
  85       return 0;
  86     }
  87
  88   /* Arrange so that a > b, subtract an -= bn, and maintain
  89      normalization. */
  90   if (an < bn)
  91     {
  92       MPN_PTR_SWAP (ap, an, bp, bn);
  93       MP_PTR_SWAP (u0, u1);
  94       swapped ^= 1;
  95     }
  96   else if (an == bn)
  97     {
  98       int c;
  99       MPN_CMP (c, ap, bp, an);
 100       if (UNLIKELY (c == 0))
 101         {
 102           MPN_COPY (gp, ap, an);
 103           *gn = an;
 104
 105           /* Must return the smallest cofactor, +u1 or -u0 */
 106           MPN_CMP (c, u0, u1, un);
 107           ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1));
 108
 109           if (c < 0)
 110             {
 111               MPN_NORMALIZE (u0, un);
 112               MPN_COPY (up, u0, un);
 113               swapped ^= 1;
 114             }
 115           else
 116             {
 117               MPN_NORMALIZE_NOT_ZERO (u1, un);
 118               MPN_COPY (up, u1, un);
 119             }
 120
 121           *usizep = swapped ? -un : un;
 122           return 0;
 123         }
 124       else if (c < 0)
 125         {
 126           MP_PTR_SWAP (ap, bp);
 127           MP_PTR_SWAP (u0, u1);
 128           swapped ^= 1;
 129         }
 130     }
 131   /* Reduce a -= b, u1 += u0 */
 132   ASSERT_NOCARRY (mpn_sub (ap, ap, an, bp, bn));
 133   MPN_NORMALIZE (ap, an);
 134   ASSERT (an > 0);
 135
 136   u1[un] = mpn_add_n (u1, u1, u0, un);
 137   un += (u1[un] > 0);
 138
 139   /* Arrange so that a > b, and divide a = q b + r */
 140   if (an < bn)
 141     {
 142       MPN_PTR_SWAP (ap, an, bp, bn);
 143       MP_PTR_SWAP (u0, u1);
 144       swapped ^= 1;
 145     }
 146   else if (an == bn)
 147     {
 148       int c;
 149       MPN_CMP (c, ap, bp, an);
 150       if (UNLIKELY (c == 0))
 151         goto return_b;
 152       else if (c < 0)
 153         {
 154           MP_PTR_SWAP (ap, bp);
 155           MP_PTR_SWAP (u0, u1);
 156           swapped ^= 1;
 157         }
 158     }
 159
 160   /* Reduce a -= q b, u1 += q u0 */
 161   qn = an - bn + 1;
 162   mpn_tdiv_qr (qp, ap, 0, ap, an, bp, bn);
 163
 164   if (mpn_zero_p (ap, bn))
 165     goto return_b;
 166
 167   n = bn;
 168
 169   /* Update u1 += q u0 */
 170   u0n = un;
 171   MPN_NORMALIZE (u0, u0n);
 172
 173   if (u0n > 0)
 174     {
 175       qn -= (qp[qn - 1] == 0);
 176
 177       if (qn > u0n)
 178         mpn_mul (tp, qp, qn, u0, u0n);
 179       else
 180         mpn_mul (tp, u0, u0n, qp, qn);
 181
 182       if (qn + u0n > un)
 183         {
 184           mp_size_t u1n = un;
 185           un = qn + u0n;
 186           un -= (tp[un-1] == 0);
 187           u1[un] = mpn_add (u1, tp, un, u1, u1n);
 188         }
 189       else
 190         {
 191           u1[un] = mpn_add (u1, u1, un, tp, qn + u0n);
 192         }
 193
 194       un += (u1[un] > 0);
 195     }
 196
 197   *unp = un;
 198   return n;
 199 }