github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpn/generic/remove.c (about)

     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  }