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

     1  /* mpn_mod_1s_2p (ap, n, b, cps)
     2     Divide (ap,,n) by b.  Return the single-limb remainder.
     3     Requires that b < B / 2.
     4  
     5     Contributed to the GNU project by Torbjorn Granlund.
     6     Based on a suggestion by Peter L. Montgomery.
     7  
     8     THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
     9     SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
    10     GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
    11  
    12  Copyright 2008-2010 Free Software Foundation, Inc.
    13  
    14  This file is part of the GNU MP Library.
    15  
    16  The GNU MP Library is free software; you can redistribute it and/or modify
    17  it under the terms of either:
    18  
    19    * the GNU Lesser General Public License as published by the Free
    20      Software Foundation; either version 3 of the License, or (at your
    21      option) any later version.
    22  
    23  or
    24  
    25    * the GNU General Public License as published by the Free Software
    26      Foundation; either version 2 of the License, or (at your option) any
    27      later version.
    28  
    29  or both in parallel, as here.
    30  
    31  The GNU MP Library is distributed in the hope that it will be useful, but
    32  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    33  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    34  for more details.
    35  
    36  You should have received copies of the GNU General Public License and the
    37  GNU Lesser General Public License along with the GNU MP Library.  If not,
    38  see https://www.gnu.org/licenses/.  */
    39  
    40  #include "gmp.h"
    41  #include "gmp-impl.h"
    42  #include "longlong.h"
    43  
    44  void
    45  mpn_mod_1s_2p_cps (mp_limb_t cps[5], mp_limb_t b)
    46  {
    47    mp_limb_t bi;
    48    mp_limb_t B1modb, B2modb, B3modb;
    49    int cnt;
    50  
    51    ASSERT (b <= (~(mp_limb_t) 0) / 2);
    52  
    53    count_leading_zeros (cnt, b);
    54  
    55    b <<= cnt;
    56    invert_limb (bi, b);
    57  
    58    cps[0] = bi;
    59    cps[1] = cnt;
    60  
    61    B1modb = -b * ((bi >> (GMP_LIMB_BITS-cnt)) | (CNST_LIMB(1) << cnt));
    62    ASSERT (B1modb <= b);		/* NB: not fully reduced mod b */
    63    cps[2] = B1modb >> cnt;
    64  
    65    udiv_rnnd_preinv (B2modb, B1modb, CNST_LIMB(0), b, bi);
    66    cps[3] = B2modb >> cnt;
    67  
    68    udiv_rnnd_preinv (B3modb, B2modb, CNST_LIMB(0), b, bi);
    69    cps[4] = B3modb >> cnt;
    70  
    71  #if WANT_ASSERT
    72    {
    73      int i;
    74      b = cps[2];
    75      for (i = 3; i <= 4; i++)
    76        {
    77  	b += cps[i];
    78  	ASSERT (b >= cps[i]);
    79        }
    80    }
    81  #endif
    82  }
    83  
    84  mp_limb_t
    85  mpn_mod_1s_2p (mp_srcptr ap, mp_size_t n, mp_limb_t b, const mp_limb_t cps[5])
    86  {
    87    mp_limb_t rh, rl, bi, ph, pl, ch, cl, r;
    88    mp_limb_t B1modb, B2modb, B3modb;
    89    mp_size_t i;
    90    int cnt;
    91  
    92    ASSERT (n >= 1);
    93  
    94    B1modb = cps[2];
    95    B2modb = cps[3];
    96    B3modb = cps[4];
    97  
    98    if ((n & 1) != 0)
    99      {
   100        if (n == 1)
   101  	{
   102  	  rl = ap[n - 1];
   103  	  bi = cps[0];
   104  	  cnt = cps[1];
   105  	  udiv_rnnd_preinv (r, rl >> (GMP_LIMB_BITS - cnt),
   106  			     rl << cnt, b, bi);
   107  	  return r >> cnt;
   108  	}
   109  
   110        umul_ppmm (ph, pl, ap[n - 2], B1modb);
   111        add_ssaaaa (ph, pl, ph, pl, CNST_LIMB(0), ap[n - 3]);
   112        umul_ppmm (rh, rl, ap[n - 1], B2modb);
   113        add_ssaaaa (rh, rl, rh, rl, ph, pl);
   114        n--;
   115      }
   116    else
   117      {
   118        rh = ap[n - 1];
   119        rl = ap[n - 2];
   120      }
   121  
   122    for (i = n - 4; i >= 0; i -= 2)
   123      {
   124        /* rr = ap[i]				< B
   125  	    + ap[i+1] * (B mod b)		<= (B-1)(b-1)
   126  	    + LO(rr)  * (B^2 mod b)		<= (B-1)(b-1)
   127  	    + HI(rr)  * (B^3 mod b)		<= (B-1)(b-1)
   128        */
   129        umul_ppmm (ph, pl, ap[i + 1], B1modb);
   130        add_ssaaaa (ph, pl, ph, pl, CNST_LIMB(0), ap[i + 0]);
   131  
   132        umul_ppmm (ch, cl, rl, B2modb);
   133        add_ssaaaa (ph, pl, ph, pl, ch, cl);
   134  
   135        umul_ppmm (rh, rl, rh, B3modb);
   136        add_ssaaaa (rh, rl, rh, rl, ph, pl);
   137      }
   138  
   139    umul_ppmm (rh, cl, rh, B1modb);
   140    add_ssaaaa (rh, rl, rh, rl, CNST_LIMB(0), cl);
   141  
   142    cnt = cps[1];
   143    bi = cps[0];
   144  
   145    r = (rh << cnt) | (rl >> (GMP_LIMB_BITS - cnt));
   146    udiv_rnnd_preinv (r, r, rl << cnt, b, bi);
   147  
   148    return r >> cnt;
   149  }