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

   1 /* hgcd_matrix.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, 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 /* For input of size n, matrix elements are of size at most ceil(n/2)
  40    - 1, but we need two limbs extra. */
  41 void
  42 mpn_hgcd_matrix_init (struct hgcd_matrix *M, mp_size_t n, mp_ptr p)
  43 {
  44   mp_size_t s = (n+1)/2 + 1;
  45   M->alloc = s;
  46   M->n = 1;
  47   MPN_ZERO (p, 4 * s);
  48   M->p[0][0] = p;
  49   M->p[0][1] = p + s;
  50   M->p[1][0] = p + 2 * s;
  51   M->p[1][1] = p + 3 * s;
  52
  53   M->p[0][0][0] = M->p[1][1][0] = 1;
  54 }
  55
  56 /* Update column COL, adding in Q * column (1-COL). Temporary storage:
  57  * qn + n <= M->alloc, where n is the size of the largest element in
  58  * column 1 - COL. */
  59 void
  60 mpn_hgcd_matrix_update_q (struct hgcd_matrix *M, mp_srcptr qp, mp_size_t qn,
  61                           unsigned col, mp_ptr tp)
  62 {
  63   ASSERT (col < 2);
  64
  65   if (qn == 1)
  66     {
  67       mp_limb_t q = qp[0];
  68       mp_limb_t c0, c1;
  69
  70       c0 = mpn_addmul_1 (M->p[0][col], M->p[0][1-col], M->n, q);
  71       c1 = mpn_addmul_1 (M->p[1][col], M->p[1][1-col], M->n, q);
  72
  73       M->p[0][col][M->n] = c0;
  74       M->p[1][col][M->n] = c1;
  75
  76       M->n += (c0 | c1) != 0;
  77     }
  78   else
  79     {
  80       unsigned row;
  81
  82       /* Carries for the unlikely case that we get both high words
  83          from the multiplication and carries from the addition. */
  84       mp_limb_t c[2];
  85       mp_size_t n;
  86
  87       /* The matrix will not necessarily grow in size by qn, so we
  88          need normalization in order not to overflow M. */
  89
  90       for (n = M->n; n + qn > M->n; n--)
  91         {
  92           ASSERT (n > 0);
  93           if (M->p[0][1-col][n-1] > 0 || M->p[1][1-col][n-1] > 0)
  94             break;
  95         }
  96
  97       ASSERT (qn + n <= M->alloc);
  98
  99       for (row = 0; row < 2; row++)
 100         {
 101           if (qn <= n)
 102             mpn_mul (tp, M->p[row][1-col], n, qp, qn);
 103           else
 104             mpn_mul (tp, qp, qn, M->p[row][1-col], n);
 105
 106           ASSERT (n + qn >= M->n);
 107           c[row] = mpn_add (M->p[row][col], tp, n + qn, M->p[row][col], M->n);
 108         }
 109
 110       n += qn;
 111
 112       if (c[0] | c[1])
 113         {
 114           M->p[0][col][n] = c[0];
 115           M->p[1][col][n] = c[1];
 116           n++;
 117         }
 118       else
 119         {
 120           n -= (M->p[0][col][n-1] | M->p[1][col][n-1]) == 0;
 121           ASSERT (n >= M->n);
 122         }
 123       M->n = n;
 124     }
 125
 126   ASSERT (M->n < M->alloc);
 127 }
 128
 129 /* Multiply M by M1 from the right. Since the M1 elements fit in
 130    GMP_NUMB_BITS - 1 bits, M grows by at most one limb. Needs
 131    temporary space M->n */
 132 void
 133 mpn_hgcd_matrix_mul_1 (struct hgcd_matrix *M, const struct hgcd_matrix1 *M1,
 134                        mp_ptr tp)
 135 {
 136   mp_size_t n0, n1;
 137
 138   /* Could avoid copy by some swapping of pointers. */
 139   MPN_COPY (tp, M->p[0][0], M->n);
 140   n0 = mpn_hgcd_mul_matrix1_vector (M1, M->p[0][0], tp, M->p[0][1], M->n);
 141   MPN_COPY (tp, M->p[1][0], M->n);
 142   n1 = mpn_hgcd_mul_matrix1_vector (M1, M->p[1][0], tp, M->p[1][1], M->n);
 143
 144   /* Depends on zero initialization */
 145   M->n = MAX(n0, n1);
 146   ASSERT (M->n < M->alloc);
 147 }
 148
 149 /* Multiply M by M1 from the right. Needs 3*(M->n + M1->n) + 5 limbs
 150    of temporary storage (see mpn_matrix22_mul_itch). */
 151 void
 152 mpn_hgcd_matrix_mul (struct hgcd_matrix *M, const struct hgcd_matrix *M1,
 153                      mp_ptr tp)
 154 {
 155   mp_size_t n;
 156
 157   /* About the new size of M:s elements. Since M1's diagonal elements
 158      are > 0, no element can decrease. The new elements are of size
 159      M->n + M1->n, one limb more or less. The computation of the
 160      matrix product produces elements of size M->n + M1->n + 1. But
 161      the true size, after normalization, may be three limbs smaller.
 162
 163      The reason that the product has normalized size >= M->n + M1->n -
 164      2 is subtle. It depends on the fact that M and M1 can be factored
 165      as products of (1,1; 0,1) and (1,0; 1,1), and that we can't have
 166      M ending with a large power and M1 starting with a large power of
 167      the same matrix. */
 168
 169   /* FIXME: Strassen multiplication gives only a small speedup. In FFT
 170      multiplication range, this function could be sped up quite a lot
 171      using invariance. */
 172   ASSERT (M->n + M1->n < M->alloc);
 173
 174   ASSERT ((M->p[0][0][M->n-1] | M->p[0][1][M->n-1]
 175            | M->p[1][0][M->n-1] | M->p[1][1][M->n-1]) > 0);
 176
 177   ASSERT ((M1->p[0][0][M1->n-1] | M1->p[0][1][M1->n-1]
 178            | M1->p[1][0][M1->n-1] | M1->p[1][1][M1->n-1]) > 0);
 179
 180   mpn_matrix22_mul (M->p[0][0], M->p[0][1],
 181                     M->p[1][0], M->p[1][1], M->n,
 182                     M1->p[0][0], M1->p[0][1],
 183                     M1->p[1][0], M1->p[1][1], M1->n, tp);
 184
 185   /* Index of last potentially non-zero limb, size is one greater. */
 186   n = M->n + M1->n;
 187
 188   n -= ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) == 0);
 189   n -= ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) == 0);
 190   n -= ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) == 0);
 191
 192   ASSERT ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) > 0);
 193
 194   M->n = n + 1;
 195 }
 196
 197 /* Multiplies the least significant p limbs of (a;b) by M^-1.
 198    Temporary space needed: 2 * (p + M->n)*/
 199 mp_size_t
 200 mpn_hgcd_matrix_adjust (const struct hgcd_matrix *M,
 201                         mp_size_t n, mp_ptr ap, mp_ptr bp,
 202                         mp_size_t p, mp_ptr tp)
 203 {
 204   /* M^-1 (a;b) = (r11, -r01; -r10, r00) (a ; b)
 205      = (r11 a - r01 b; - r10 a + r00 b */
 206
 207   mp_ptr t0 = tp;
 208   mp_ptr t1 = tp + p + M->n;
 209   mp_limb_t ah, bh;
 210   mp_limb_t cy;
 211
 212   ASSERT (p + M->n  < n);
 213
 214   /* First compute the two values depending on a, before overwriting a */
 215
 216   if (M->n >= p)
 217     {
 218       mpn_mul (t0, M->p[1][1], M->n, ap, p);
 219       mpn_mul (t1, M->p[1][0], M->n, ap, p);
 220     }
 221   else
 222     {
 223       mpn_mul (t0, ap, p, M->p[1][1], M->n);
 224       mpn_mul (t1, ap, p, M->p[1][0], M->n);
 225     }
 226
 227   /* Update a */
 228   MPN_COPY (ap, t0, p);
 229   ah = mpn_add (ap + p, ap + p, n - p, t0 + p, M->n);
 230
 231   if (M->n >= p)
 232     mpn_mul (t0, M->p[0][1], M->n, bp, p);
 233   else
 234     mpn_mul (t0, bp, p, M->p[0][1], M->n);
 235
 236   cy = mpn_sub (ap, ap, n, t0, p + M->n);
 237   ASSERT (cy <= ah);
 238   ah -= cy;
 239
 240   /* Update b */
 241   if (M->n >= p)
 242     mpn_mul (t0, M->p[0][0], M->n, bp, p);
 243   else
 244     mpn_mul (t0, bp, p, M->p[0][0], M->n);
 245
 246   MPN_COPY (bp, t0, p);
 247   bh = mpn_add (bp + p, bp + p, n - p, t0 + p, M->n);
 248   cy = mpn_sub (bp, bp, n, t1, p + M->n);
 249   ASSERT (cy <= bh);
 250   bh -= cy;
 251
 252   if (ah > 0 || bh > 0)
 253     {
 254       ap[n] = ah;
 255       bp[n] = bh;
 256       n++;
 257     }
 258   else
 259     {
 260       /* The subtraction can reduce the size by at most one limb. */
 261       if (ap[n-1] == 0 && bp[n-1] == 0)
 262         n--;
 263     }
 264   ASSERT (ap[n-1] > 0 || bp[n-1] > 0);
 265   return n;
 266 }