Home | History | Annotate | Line # | Download | only in generic
remove.c revision 1.1.1.1.2.1 
      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 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 the GNU Lesser General Public License as published by
     16 the Free Software Foundation; either version 3 of the License, or (at your
     17 option) any later version.
     18 
     19 The GNU MP Library is distributed in the hope that it will be useful, but
     20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
     21 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
     22 License for more details.
     23 
     24 You should have received a copy of the GNU Lesser General Public License
     25 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
     26 
     27 #include "gmp.h"
     28 #include "gmp-impl.h"
     29 
     30 #if GMP_LIMB_BITS > 50
     31 #define LOG 50
     32 #else
     33 #define LOG GMP_LIMB_BITS
     34 #endif
     35 
     36 
     37 /* Input: U = {up,un}, V = {vp,vn} must be odd, cap
     38    Ouput  W = {wp,*wn} allocation need is exactly *wn
     39 
     40    Set W = U / V^k, where k is the largest integer <= cap such that the
     41    division yields an integer.
     42 
     43    FIXME: We currently allow any operand overlap.  This is quite non mpn-ish
     44    and might be changed, since it cost significant temporary space.
     45    * If we require W to have space for un limbs, we could save qp or qp2 (but
     46      we will still need to copy things into wp 50% of the time).
     47    * If we allow ourselves to clobber U, we could save the other of qp and qp2.
     48 */
     49 
     50 /* FIXME: We need to wrap mpn_bdiv_qr due to the itch interface.  This need
     51    indicates a flaw in the current itch mechanism: Which operands not greater
     52    than un,un will incur the worst itch?  We need a parallel foo_maxitch set
     53    of functions.  */
     54 static void
     55 mpn_bdiv_qr_wrap (mp_ptr qp, mp_ptr rp,
     56 		  mp_srcptr np, mp_size_t nn,
     57 		  mp_srcptr dp, mp_size_t dn)
     58 {
     59   mp_ptr scratch_out;
     60   TMP_DECL;
     61 
     62   TMP_MARK;
     63   scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (nn, dn));
     64   mpn_bdiv_qr (qp, rp, np, nn, dp, dn, scratch_out);
     65 
     66   TMP_FREE;
     67 }
     68 
     69 mp_bitcnt_t
     70 mpn_remove (mp_ptr wp, mp_size_t *wn,
     71 	    mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn,
     72 	    mp_bitcnt_t cap)
     73 {
     74   mp_ptr    pwpsp[LOG];
     75   mp_size_t pwpsn[LOG];
     76   mp_size_t npowers;
     77   mp_ptr tp, qp, np, pp, qp2;
     78   mp_size_t pn, nn, qn, i;
     79   mp_bitcnt_t pwr;
     80   TMP_DECL;
     81 
     82   ASSERT (un > 0);
     83   ASSERT (vn > 0);
     84   ASSERT (vp[0] % 2 != 0);	/* 2-adic division wants odd numbers */
     85   ASSERT (vn > 1 || vp[0] > 1);	/* else we would loop indefinitely */
     86 
     87   TMP_MARK;
     88 
     89   tp = TMP_ALLOC_LIMBS ((un + 1 + vn) / 2); /* remainder */
     90   qp = TMP_ALLOC_LIMBS (un + 1);	/* quotient, alternating */
     91   qp2 = TMP_ALLOC_LIMBS (un + 1);	/* quotient, alternating */
     92   np = TMP_ALLOC_LIMBS (un + LOG);	/* powers of V */
     93   pp = vp;
     94   pn = vn;
     95 
     96   MPN_COPY (qp, up, un);
     97   qn = un;
     98 
     99   npowers = 0;
    100   while (qn >= pn)
    101     {
    102       qp[qn] = 0;
    103       mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pp, pn);
    104       if (!mpn_zero_p (tp, pn))
    105 	break;			/* could not divide by V^npowers */
    106 
    107       MP_PTR_SWAP (qp, qp2);
    108       qn = qn - pn;
    109       qn += qp[qn] != 0;
    110 
    111       pwpsp[npowers] = pp;
    112       pwpsn[npowers] = pn;
    113       npowers++;
    114 
    115       if (((mp_bitcnt_t) 2 << npowers) - 1 > cap)
    116 	break;
    117 
    118       nn = 2 * pn - 1;		/* next power will be at least this large */
    119       if (nn > qn)
    120 	break;			/* next power would be overlarge */
    121 
    122       mpn_sqr (np, pp, pn);
    123       nn += np[nn] != 0;
    124       pp = np;
    125       pn = nn;
    126       np += nn;
    127     }
    128 
    129   pwr = ((mp_bitcnt_t) 1 << npowers) - 1;
    130 
    131   for (i = npowers - 1; i >= 0; i--)
    132     {
    133       pp = pwpsp[i];
    134       pn = pwpsn[i];
    135       if (qn < pn)
    136 	continue;
    137 
    138       if (pwr + ((mp_bitcnt_t) 1 << i) > cap)
    139 	continue;		/* V^i would bring us past cap */
    140 
    141       qp[qn] = 0;
    142       mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pp, pn);
    143       if (!mpn_zero_p (tp, pn))
    144 	continue;		/* could not divide by V^i */
    145 
    146       MP_PTR_SWAP (qp, qp2);
    147       qn = qn - pn;
    148       qn += qp[qn] != 0;
    149 
    150       pwr += (mp_bitcnt_t) 1 << i;
    151     }
    152 
    153   MPN_COPY (wp, qp, qn);
    154   *wn = qn;
    155 
    156   TMP_FREE;
    157 
    158   return pwr;
    159 }
    160