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

     1  /* mpn_dcpi1_bdiv_q -- divide-and-conquer Hensel division with precomputed
     2     inverse, returning quotient.
     3  
     4     Contributed to the GNU project by Niels Möller and Torbjorn Granlund.
     5  
     6     THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
     7     SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
     8     GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
     9  
    10  Copyright 2006, 2007, 2009-2011 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  
    42  mp_size_t
    43  mpn_dcpi1_bdiv_q_n_itch (mp_size_t n)
    44  {
    45    /* NOTE: Depends on mullo_n interface */
    46    return n;
    47  }
    48  
    49  /* Computes Q = N / D mod B^n, destroys N.
    50  
    51     N = {np,n}
    52     D = {dp,n}
    53  */
    54  
    55  void
    56  mpn_dcpi1_bdiv_q_n (mp_ptr qp,
    57  		    mp_ptr np, mp_srcptr dp, mp_size_t n,
    58  		    mp_limb_t dinv, mp_ptr tp)
    59  {
    60    while (ABOVE_THRESHOLD (n, DC_BDIV_Q_THRESHOLD))
    61      {
    62        mp_size_t lo, hi;
    63        mp_limb_t cy;
    64  
    65        lo = n >> 1;			/* floor(n/2) */
    66        hi = n - lo;			/* ceil(n/2) */
    67  
    68        cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, lo, dinv, tp);
    69  
    70        mpn_mullo_n (tp, qp, dp + hi, lo);
    71        mpn_sub_n (np + hi, np + hi, tp, lo);
    72  
    73        if (lo < hi)
    74  	{
    75  	  cy += mpn_submul_1 (np + lo, qp, lo, dp[lo]);
    76  	  np[n - 1] -= cy;
    77  	}
    78        qp += lo;
    79        np += lo;
    80        n -= lo;
    81      }
    82    mpn_sbpi1_bdiv_q (qp, np, n, dp, n, dinv);
    83  }
    84  
    85  /* Computes Q = N / D mod B^nn, destroys N.
    86  
    87     N = {np,nn}
    88     D = {dp,dn}
    89  */
    90  
    91  void
    92  mpn_dcpi1_bdiv_q (mp_ptr qp,
    93  		  mp_ptr np, mp_size_t nn,
    94  		  mp_srcptr dp, mp_size_t dn,
    95  		  mp_limb_t dinv)
    96  {
    97    mp_size_t qn;
    98    mp_limb_t cy;
    99    mp_ptr tp;
   100    TMP_DECL;
   101  
   102    TMP_MARK;
   103  
   104    ASSERT (dn >= 2);
   105    ASSERT (nn - dn >= 0);
   106    ASSERT (dp[0] & 1);
   107  
   108    tp = TMP_SALLOC_LIMBS (dn);
   109  
   110    qn = nn;
   111  
   112    if (qn > dn)
   113      {
   114        /* Reduce qn mod dn in a super-efficient manner.  */
   115        do
   116  	qn -= dn;
   117        while (qn > dn);
   118  
   119        /* Perform the typically smaller block first.  */
   120        if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD))
   121  	cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv);
   122        else
   123  	cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp);
   124  
   125        if (qn != dn)
   126  	{
   127  	  if (qn > dn - qn)
   128  	    mpn_mul (tp, qp, qn, dp + qn, dn - qn);
   129  	  else
   130  	    mpn_mul (tp, dp + qn, dn - qn, qp, qn);
   131  	  mpn_incr_u (tp + qn, cy);
   132  
   133  	  mpn_sub (np + qn, np + qn, nn - qn, tp, dn);
   134  	  cy = 0;
   135  	}
   136  
   137        np += qn;
   138        qp += qn;
   139  
   140        qn = nn - qn;
   141        while (qn > dn)
   142  	{
   143  	  mpn_sub_1 (np + dn, np + dn, qn - dn, cy);
   144  	  cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, dn, dinv, tp);
   145  	  qp += dn;
   146  	  np += dn;
   147  	  qn -= dn;
   148  	}
   149        mpn_dcpi1_bdiv_q_n (qp, np, dp, dn, dinv, tp);
   150      }
   151    else
   152      {
   153        if (BELOW_THRESHOLD (qn, DC_BDIV_Q_THRESHOLD))
   154  	mpn_sbpi1_bdiv_q (qp, np, qn, dp, qn, dinv);
   155        else
   156  	mpn_dcpi1_bdiv_q_n (qp, np, dp, qn, dinv, tp);
   157      }
   158  
   159    TMP_FREE;
   160  }