contrib/gmp/mpn/generic/sb_div_qr.c

   1 /* mpn_sb_div_qr -- schoolbook division with 2-limb sloppy non-greater
   2    precomputed inverse, returning quotient and remainder.
   3
   4    Contributed to the GNU project by Torbjörn Granlund.
   5
   6    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH A MUTABLE INTERFACE.  IT IS
   7    ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS
   8    ALMOST GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP
   9    RELEASE.
  10
  11 Copyright 2006, 2007 Free Software Foundation, Inc.
  12
  13 This file is part of the GNU MP Library.
  14
  15 The GNU MP Library is free software; you can redistribute it and/or modify
  16 it under the terms of the GNU Lesser General Public License as published by
  17 the Free Software Foundation; either version 3 of the License, or (at your
  18 option) any later version.
  19
  20 The GNU MP Library is distributed in the hope that it will be useful, but
  21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  22 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  23 License for more details.
  24
  25 You should have received a copy of the GNU Lesser General Public License
  26 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  27
  28 #include "gmp.h"
  29 #include "gmp-impl.h"
  30 #include "longlong.h"
  31
  32 /*
  33   CAVEATS:
  34   1. Should it demand normalized operands like now, or normalize on-the-fly?
  35   2. Overwrites {np,nn} instead of writing remainder to a designated area.
  36   3. Uses mpn_submul_1.  It would be nice to somehow make it use mpn_addmul_1
  37      instead.  (That would open for mpn_addmul_2 straightforwardly.)
  38 */
  39
  40 mp_limb_t
  41 mpn_sb_div_qr (mp_ptr qp,
  42                mp_ptr np, mp_size_t nn,
  43                mp_srcptr dp, mp_size_t dn,
  44                mp_srcptr dip)
  45 {
  46   mp_limb_t q, q10, q01a, q00a, q01b, q00b;
  47   mp_limb_t cy;
  48   mp_size_t i;
  49   mp_limb_t qh;
  50   mp_limb_t di1, di0;
  51
  52   ASSERT (dn > 0);
  53   ASSERT (nn >= dn);
  54   ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
  55   ASSERT (! MPN_OVERLAP_P (np, nn, dp, dn));
  56   ASSERT (! MPN_OVERLAP_P (qp, nn-dn, dp, dn));
  57   ASSERT (! MPN_OVERLAP_P (qp, nn-dn, np, nn) || qp+dn >= np);
  58   ASSERT_MPN (np, nn);
  59   ASSERT_MPN (dp, dn);
  60
  61   np += nn;
  62
  63   qh = mpn_cmp (np - dn, dp, dn) >= 0;
  64   if (qh != 0)
  65     mpn_sub_n (np - dn, np - dn, dp, dn);
  66
  67   qp += nn - dn;
  68   di1 = dip[1]; di0 = dip[0];
  69   for (i = nn - dn; i > 0; i--)
  70     {
  71       np--;
  72       umul_ppmm (q, q10, np[0], di1);
  73       umul_ppmm (q01a, q00a, np[-1], di1);
  74       add_ssaaaa (q, q10, q, q10, np[0], q01a);
  75       umul_ppmm (q01b, q00b, np[0], di0);
  76       add_ssaaaa (q, q10, q, q10, 0, q01b);
  77       add_ssaaaa (q, q10, q, q10, 0, np[-1]);
  78
  79       cy = mpn_submul_1 (np - dn, dp, dn, q);
  80
  81       if (UNLIKELY (np[0] > cy || mpn_cmp (np - dn, dp, dn) >= 0))
  82         {
  83           q = q + 1;
  84           mpn_sub_n (np - dn, np - dn, dp, dn);
  85         }
  86
  87       *--qp = q;
  88     }
  89
  90   return qh;
  91 }