source/libs/gmp/gmp-src/mpn/generic/gcd_subdiv_step.c

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