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

   1 /* hgcd_reduce.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 2011, 2012 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 "gmp.h"
  36 #include "gmp-impl.h"
  37 #include "longlong.h"
  38
  39 /* Computes R -= A * B. Result must be non-negative. Normalized down
  40    to size an, and resulting size is returned. */
  41 static mp_size_t
  42 submul (mp_ptr rp, mp_size_t rn,
  43         mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
  44 {
  45   mp_ptr tp;
  46   TMP_DECL;
  47
  48   ASSERT (bn > 0);
  49   ASSERT (an >= bn);
  50   ASSERT (rn >= an);
  51   ASSERT (an + bn <= rn + 1);
  52
  53   TMP_MARK;
  54   tp = TMP_ALLOC_LIMBS (an + bn);
  55
  56   mpn_mul (tp, ap, an, bp, bn);
  57   ASSERT ((an + bn <= rn) || (tp[rn] == 0));
  58   ASSERT_NOCARRY (mpn_sub (rp, rp, rn, tp, an + bn - (an + bn > rn)));
  59   TMP_FREE;
  60
  61   while (rn > an && (rp[rn-1] == 0))
  62     rn--;
  63
  64   return rn;
  65 }
  66
  67 /* Computes (a, b)  <--  M^{-1} (a; b) */
  68 /* FIXME:
  69     x Take scratch parameter, and figure out scratch need.
  70
  71     x Use some fallback for small M->n?
  72 */
  73 static mp_size_t
  74 hgcd_matrix_apply (const struct hgcd_matrix *M,
  75                    mp_ptr ap, mp_ptr bp,
  76                    mp_size_t n)
  77 {
  78   mp_size_t an, bn, un, vn, nn;
  79   mp_size_t mn[2][2];
  80   mp_size_t modn;
  81   mp_ptr tp, sp, scratch;
  82   mp_limb_t cy;
  83   unsigned i, j;
  84
  85   TMP_DECL;
  86
  87   ASSERT ( (ap[n-1] | bp[n-1]) > 0);
  88
  89   an = n;
  90   MPN_NORMALIZE (ap, an);
  91   bn = n;
  92   MPN_NORMALIZE (bp, bn);
  93
  94   for (i = 0; i < 2; i++)
  95     for (j = 0; j < 2; j++)
  96       {
  97         mp_size_t k;
  98         k = M->n;
  99         MPN_NORMALIZE (M->p[i][j], k);
 100         mn[i][j] = k;
 101       }
 102
 103   ASSERT (mn[0][0] > 0);
 104   ASSERT (mn[1][1] > 0);
 105   ASSERT ( (mn[0][1] | mn[1][0]) > 0);
 106
 107   TMP_MARK;
 108
 109   if (mn[0][1] == 0)
 110     {
 111       /* A unchanged, M = (1, 0; q, 1) */
 112       ASSERT (mn[0][0] == 1);
 113       ASSERT (M->p[0][0][0] == 1);
 114       ASSERT (mn[1][1] == 1);
 115       ASSERT (M->p[1][1][0] == 1);
 116
 117       /* Put B <-- B - q A */
 118       nn = submul (bp, bn, ap, an, M->p[1][0], mn[1][0]);
 119     }
 120   else if (mn[1][0] == 0)
 121     {
 122       /* B unchanged, M = (1, q; 0, 1) */
 123       ASSERT (mn[0][0] == 1);
 124       ASSERT (M->p[0][0][0] == 1);
 125       ASSERT (mn[1][1] == 1);
 126       ASSERT (M->p[1][1][0] == 1);
 127
 128       /* Put A  <-- A - q * B */
 129       nn = submul (ap, an, bp, bn, M->p[0][1], mn[0][1]);
 130     }
 131   else
 132     {
 133       /* A = m00 a + m01 b  ==> a <= A / m00, b <= A / m01.
 134          B = m10 a + m11 b  ==> a <= B / m10, b <= B / m11. */
 135       un = MIN (an - mn[0][0], bn - mn[1][0]) + 1;
 136       vn = MIN (an - mn[0][1], bn - mn[1][1]) + 1;
 137
 138       nn = MAX (un, vn);
 139       /* In the range of interest, mulmod_bnm1 should always beat mullo. */
 140       modn = mpn_mulmod_bnm1_next_size (nn + 1);
 141
 142       TMP_ALLOC_LIMBS_3 (tp, modn,
 143                          sp, modn,
 144                          scratch, mpn_mulmod_bnm1_itch (modn, modn, M->n));
 145
 146       ASSERT (n <= 2*modn);
 147
 148       if (n > modn)
 149         {
 150           cy = mpn_add (ap, ap, modn, ap + modn, n - modn);
 151           MPN_INCR_U (ap, modn, cy);
 152
 153           cy = mpn_add (bp, bp, modn, bp + modn, n - modn);
 154           MPN_INCR_U (bp, modn, cy);
 155
 156           n = modn;
 157         }
 158
 159       mpn_mulmod_bnm1 (tp, modn, ap, n, M->p[1][1], mn[1][1], scratch);
 160       mpn_mulmod_bnm1 (sp, modn, bp, n, M->p[0][1], mn[0][1], scratch);
 161
 162       /* FIXME: Handle the small n case in some better way. */
 163       if (n + mn[1][1] < modn)
 164         MPN_ZERO (tp + n + mn[1][1], modn - n - mn[1][1]);
 165       if (n + mn[0][1] < modn)
 166         MPN_ZERO (sp + n + mn[0][1], modn - n - mn[0][1]);
 167
 168       cy = mpn_sub_n (tp, tp, sp, modn);
 169       MPN_DECR_U (tp, modn, cy);
 170
 171       ASSERT (mpn_zero_p (tp + nn, modn - nn));
 172
 173       mpn_mulmod_bnm1 (sp, modn, ap, n, M->p[1][0], mn[1][0], scratch);
 174       MPN_COPY (ap, tp, nn);
 175       mpn_mulmod_bnm1 (tp, modn, bp, n, M->p[0][0], mn[0][0], scratch);
 176
 177       if (n + mn[1][0] < modn)
 178         MPN_ZERO (sp + n + mn[1][0], modn - n - mn[1][0]);
 179       if (n + mn[0][0] < modn)
 180         MPN_ZERO (tp + n + mn[0][0], modn - n - mn[0][0]);
 181
 182       cy = mpn_sub_n (tp, tp, sp, modn);
 183       MPN_DECR_U (tp, modn, cy);
 184
 185       ASSERT (mpn_zero_p (tp + nn, modn - nn));
 186       MPN_COPY (bp, tp, nn);
 187
 188       while ( (ap[nn-1] | bp[nn-1]) == 0)
 189         {
 190           nn--;
 191           ASSERT (nn > 0);
 192         }
 193     }
 194   TMP_FREE;
 195
 196   return nn;
 197 }
 198
 199 mp_size_t
 200 mpn_hgcd_reduce_itch (mp_size_t n, mp_size_t p)
 201 {
 202   mp_size_t itch;
 203   if (BELOW_THRESHOLD (n, HGCD_REDUCE_THRESHOLD))
 204     {
 205       itch = mpn_hgcd_itch (n-p);
 206
 207       /* For arbitrary p, the storage for _adjust is 2*(p + M->n) = 2 *
 208          (p + ceil((n-p)/2) - 1 <= n + p - 1 */
 209       if (itch < n + p - 1)
 210         itch = n + p - 1;
 211     }
 212   else
 213     {
 214       itch = 2*(n-p) + mpn_hgcd_itch (n-p);
 215       /* Currently, hgcd_matrix_apply allocates its own storage. */
 216     }
 217   return itch;
 218 }
 219
 220 /* FIXME: Document storage need. */
 221 mp_size_t
 222 mpn_hgcd_reduce (struct hgcd_matrix *M,
 223                  mp_ptr ap, mp_ptr bp, mp_size_t n, mp_size_t p,
 224                  mp_ptr tp)
 225 {
 226   mp_size_t nn;
 227   if (BELOW_THRESHOLD (n, HGCD_REDUCE_THRESHOLD))
 228     {
 229       nn = mpn_hgcd (ap + p, bp + p, n - p, M, tp);
 230       if (nn > 0)
 231         /* Needs 2*(p + M->n) <= 2*(floor(n/2) + ceil(n/2) - 1)
 232            = 2 (n - 1) */
 233         return mpn_hgcd_matrix_adjust (M, p + nn, ap, bp, p, tp);
 234     }
 235   else
 236     {
 237       MPN_COPY (tp, ap + p, n - p);
 238       MPN_COPY (tp + n - p, bp + p, n - p);
 239       if (mpn_hgcd_appr (tp, tp + n - p, n - p, M, tp + 2*(n-p)))
 240         return hgcd_matrix_apply (M, ap, bp, n);
 241     }
 242   return 0;
 243 }