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

   1 /* mpn_remove -- divide out all multiples of odd mpn number from another mpn
   2    number.
   3
   4    Contributed to the GNU project by Torbjorn Granlund.
   5
   6    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
   7    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   8    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
   9
  10 Copyright 2009, 2012-2014 Free Software Foundation, Inc.
  11
  12 This file is part of the GNU MP Library.
  13
  14 The GNU MP Library is free software; you can redistribute it and/or modify
  15 it under the terms of either:
  16
  17   * the GNU Lesser General Public License as published by the Free
  18     Software Foundation; either version 3 of the License, or (at your
  19     option) any later version.
  20
  21 or
  22
  23   * the GNU General Public License as published by the Free Software
  24     Foundation; either version 2 of the License, or (at your option) any
  25     later version.
  26
  27 or both in parallel, as here.
  28
  29 The GNU MP Library is distributed in the hope that it will be useful, but
  30 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  31 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  32 for more details.
  33
  34 You should have received copies of the GNU General Public License and the
  35 GNU Lesser General Public License along with the GNU MP Library.  If not,
  36 see https://www.gnu.org/licenses/.  */
  37
  38 #include "gmp.h"
  39 #include "gmp-impl.h"
  40
  41 #if GMP_LIMB_BITS > 50
  42 #define LOG 50
  43 #else
  44 #define LOG GMP_LIMB_BITS
  45 #endif
  46
  47
  48 /* Input: U = {up,un}, V = {vp,vn} must be odd, cap
  49    Ouput  W = {wp,*wn} allocation need is exactly *wn
  50
  51    Set W = U / V^k, where k is the largest integer <= cap such that the
  52    division yields an integer.
  53
  54    FIXME: We currently allow any operand overlap.  This is quite non mpn-ish
  55    and might be changed, since it cost significant temporary space.
  56    * If we require W to have space for un + 1 limbs, we could save qp or qp2
  57      (but we will still need to copy things into wp 50% of the time).
  58    * If we allow ourselves to clobber U, we could save the other of qp and qp2,
  59      and the initial COPY (but also here we would need un + 1 limbs).
  60 */
  61
  62 /* FIXME: We need to wrap mpn_bdiv_qr due to the itch interface.  This need
  63    indicates a flaw in the current itch mechanism: Which operands not greater
  64    than un,un will incur the worst itch?  We need a parallel foo_maxitch set
  65    of functions.  */
  66 static void
  67 mpn_bdiv_qr_wrap (mp_ptr qp, mp_ptr rp,
  68                   mp_srcptr np, mp_size_t nn,
  69                   mp_srcptr dp, mp_size_t dn)
  70 {
  71   mp_ptr scratch_out;
  72   TMP_DECL;
  73
  74   TMP_MARK;
  75   scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (nn, dn));
  76   mpn_bdiv_qr (qp, rp, np, nn, dp, dn, scratch_out);
  77
  78   TMP_FREE;
  79 }
  80
  81 mp_bitcnt_t
  82 mpn_remove (mp_ptr wp, mp_size_t *wn,
  83             mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn,
  84             mp_bitcnt_t cap)
  85 {
  86   mp_srcptr pwpsp[LOG];
  87   mp_size_t pwpsn[LOG];
  88   mp_size_t npowers;
  89   mp_ptr tp, qp, np, qp2;
  90   mp_srcptr pp;
  91   mp_size_t pn, nn, qn, i;
  92   mp_bitcnt_t pwr;
  93   TMP_DECL;
  94
  95   ASSERT (un > 0);
  96   ASSERT (vn > 0);
  97   ASSERT (vp[0] % 2 != 0);      /* 2-adic division wants odd numbers */
  98   ASSERT (vn > 1 || vp[0] > 1); /* else we would loop indefinitely */
  99
 100   TMP_MARK;
 101
 102   TMP_ALLOC_LIMBS_3 (qp, un + 1,        /* quotient, alternating */
 103                      qp2, un + 1,       /* quotient, alternating */
 104                      tp, (un + 1 + vn) / 2); /* remainder */
 105   pp = vp;
 106   pn = vn;
 107
 108   MPN_COPY (qp, up, un);
 109   qn = un;
 110
 111   npowers = 0;
 112   while (qn >= pn)
 113     {
 114       qp[qn] = 0;
 115       mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pp, pn);
 116       if (!mpn_zero_p (tp, pn))
 117         break;                  /* could not divide by V^npowers */
 118
 119       MP_PTR_SWAP (qp, qp2);
 120       qn = qn - pn;
 121       qn += qp[qn] != 0;
 122
 123       pwpsp[npowers] = pp;
 124       pwpsn[npowers] = pn;
 125       ++npowers;
 126
 127       if (((mp_bitcnt_t) 2 << npowers) - 1 > cap)
 128         break;
 129
 130       nn = 2 * pn - 1;          /* next power will be at least this large */
 131       if (nn > qn)
 132         break;                  /* next power would be overlarge */
 133
 134       if (npowers == 1)         /* Alloc once, but only if it's needed */
 135         np = TMP_ALLOC_LIMBS (qn + LOG);        /* powers of V */
 136       else
 137         np += pn;
 138
 139       mpn_sqr (np, pp, pn);
 140       pn = nn + (np[nn] != 0);
 141       pp = np;
 142     }
 143
 144   pwr = ((mp_bitcnt_t) 1 << npowers) - 1;
 145
 146   for (i = npowers; --i >= 0;)
 147     {
 148       pn = pwpsn[i];
 149       if (qn < pn)
 150         continue;
 151
 152       if (pwr + ((mp_bitcnt_t) 1 << i) > cap)
 153         continue;               /* V^i would bring us past cap */
 154
 155       qp[qn] = 0;
 156       mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pwpsp[i], pn);
 157       if (!mpn_zero_p (tp, pn))
 158         continue;               /* could not divide by V^i */
 159
 160       MP_PTR_SWAP (qp, qp2);
 161       qn = qn - pn;
 162       qn += qp[qn] != 0;
 163
 164       pwr += (mp_bitcnt_t) 1 << i;
 165     }
 166
 167   MPN_COPY (wp, qp, qn);
 168   *wn = qn;
 169
 170   TMP_FREE;
 171
 172   return pwr;
 173 }