home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / gmp202.zip / mpn / generic / divrem.c < prev    next >
C/C++ Source or Header  |  1996-05-08  |  6KB  |  246 lines
  1. /* mpn_divrem -- Divide natural numbers, producing both remainder and
  2.    quotient.
  3.  
  4. Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
  5.  
  6. This file is part of the GNU MP Library.
  7.  
  8. The GNU MP Library is free software; you can redistribute it and/or modify
  9. it under the terms of the GNU Library General Public License as published by
  10. the Free Software Foundation; either version 2 of the License, or (at your
  11. option) any later version.
  12.  
  13. The GNU MP Library is distributed in the hope that it will be useful, but
  14. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
  16. License for more details.
  17.  
  18. You should have received a copy of the GNU Library General Public License
  19. along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  20. the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  21. MA 02111-1307, USA. */
  22.  
  23. #include "gmp.h"
  24. #include "gmp-impl.h"
  25. #include "longlong.h"
  26.  
  27. /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
  28.    the NSIZE-DSIZE least significant quotient limbs at QP
  29.    and the DSIZE long remainder at NP.  If QEXTRA_LIMBS is
  30.    non-zero, generate that many fraction bits and append them after the
  31.    other quotient limbs.
  32.    Return the most significant limb of the quotient, this is always 0 or 1.
  33.  
  34.    Preconditions:
  35.    0. NSIZE >= DSIZE.
  36.    1. The most significant bit of the divisor must be set.
  37.    2. QP must either not overlap with the input operands at all, or
  38.       QP + DSIZE >= NP must hold true.  (This means that it's
  39.       possible to put the quotient in the high part of NUM, right after the
  40.       remainder in NUM.
  41.    3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero.  */
  42.  
  43. mp_limb_t
  44. #if __STDC__
  45. mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
  46.         mp_ptr np, mp_size_t nsize,
  47.         mp_srcptr dp, mp_size_t dsize)
  48. #else
  49. mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize)
  50.      mp_ptr qp;
  51.      mp_size_t qextra_limbs;
  52.      mp_ptr np;
  53.      mp_size_t nsize;
  54.      mp_srcptr dp;
  55.      mp_size_t dsize;
  56. #endif
  57. {
  58.   mp_limb_t most_significant_q_limb = 0;
  59.  
  60.   switch (dsize)
  61.     {
  62.     case 0:
  63.       /* We are asked to divide by zero, so go ahead and do it!  (To make
  64.      the compiler not remove this statement, return the value.)  */
  65.       return 1 / dsize;
  66.  
  67.     case 1:
  68.       {
  69.     mp_size_t i;
  70.     mp_limb_t n1;
  71.     mp_limb_t d;
  72.  
  73.     d = dp[0];
  74.     n1 = np[nsize - 1];
  75.  
  76.     if (n1 >= d)
  77.       {
  78.         n1 -= d;
  79.         most_significant_q_limb = 1;
  80.       }
  81.  
  82.     qp += qextra_limbs;
  83.     for (i = nsize - 2; i >= 0; i--)
  84.       udiv_qrnnd (qp[i], n1, n1, np[i], d);
  85.     qp -= qextra_limbs;
  86.  
  87.     for (i = qextra_limbs - 1; i >= 0; i--)
  88.       udiv_qrnnd (qp[i], n1, n1, 0, d);
  89.  
  90.     np[0] = n1;
  91.       }
  92.       break;
  93.  
  94.     case 2:
  95.       {
  96.     mp_size_t i;
  97.     mp_limb_t n1, n0, n2;
  98.     mp_limb_t d1, d0;
  99.  
  100.     np += nsize - 2;
  101.     d1 = dp[1];
  102.     d0 = dp[0];
  103.     n1 = np[1];
  104.     n0 = np[0];
  105.  
  106.     if (n1 >= d1 && (n1 > d1 || n0 >= d0))
  107.       {
  108.         sub_ddmmss (n1, n0, n1, n0, d1, d0);
  109.         most_significant_q_limb = 1;
  110.       }
  111.  
  112.     for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
  113.       {
  114.         mp_limb_t q;
  115.         mp_limb_t r;
  116.  
  117.         if (i >= qextra_limbs)
  118.           np--;
  119.         else
  120.           np[0] = 0;
  121.  
  122.         if (n1 == d1)
  123.           {
  124.         /* Q should be either 111..111 or 111..110.  Need special
  125.            treatment of this rare case as normal division would
  126.            give overflow.  */
  127.         q = ~(mp_limb_t) 0;
  128.  
  129.         r = n0 + d1;
  130.         if (r < d1)    /* Carry in the addition? */
  131.           {
  132.             add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
  133.             qp[i] = q;
  134.             continue;
  135.           }
  136.         n1 = d0 - (d0 != 0);
  137.         n0 = -d0;
  138.           }
  139.         else
  140.           {
  141.         udiv_qrnnd (q, r, n1, n0, d1);
  142.         umul_ppmm (n1, n0, d0, q);
  143.           }
  144.  
  145.         n2 = np[0];
  146.       q_test:
  147.         if (n1 > r || (n1 == r && n0 > n2))
  148.           {
  149.         /* The estimated Q was too large.  */
  150.         q--;
  151.  
  152.         sub_ddmmss (n1, n0, n1, n0, 0, d0);
  153.         r += d1;
  154.         if (r >= d1)    /* If not carry, test Q again.  */
  155.           goto q_test;
  156.           }
  157.  
  158.         qp[i] = q;
  159.         sub_ddmmss (n1, n0, r, n2, n1, n0);
  160.       }
  161.     np[1] = n1;
  162.     np[0] = n0;
  163.       }
  164.       break;
  165.  
  166.     default:
  167.       {
  168.     mp_size_t i;
  169.     mp_limb_t dX, d1, n0;
  170.  
  171.     np += nsize - dsize;
  172.     dX = dp[dsize - 1];
  173.     d1 = dp[dsize - 2];
  174.     n0 = np[dsize - 1];
  175.  
  176.     if (n0 >= dX)
  177.       {
  178.         if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
  179.           {
  180.         mpn_sub_n (np, np, dp, dsize);
  181.         n0 = np[dsize - 1];
  182.         most_significant_q_limb = 1;
  183.           }
  184.       }
  185.  
  186.     for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
  187.       {
  188.         mp_limb_t q;
  189.         mp_limb_t n1, n2;
  190.         mp_limb_t cy_limb;
  191.  
  192.         if (i >= qextra_limbs)
  193.           {
  194.         np--;
  195.         n2 = np[dsize];
  196.           }
  197.         else
  198.           {
  199.         n2 = np[dsize - 1];
  200.         MPN_COPY_DECR (np + 1, np, dsize);
  201.         np[0] = 0;
  202.           }
  203.  
  204.         if (n0 == dX)
  205.           /* This might over-estimate q, but it's probably not worth
  206.          the extra code here to find out.  */
  207.           q = ~(mp_limb_t) 0;
  208.         else
  209.           {
  210.         mp_limb_t r;
  211.  
  212.         udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
  213.         umul_ppmm (n1, n0, d1, q);
  214.  
  215.         while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
  216.           {
  217.             q--;
  218.             r += dX;
  219.             if (r < dX)    /* I.e. "carry in previous addition?"  */
  220.               break;
  221.             n1 -= n0 < d1;
  222.             n0 -= d1;
  223.           }
  224.           }
  225.  
  226.         /* Possible optimization: We already have (q * n0) and (1 * n1)
  227.            after the calculation of q.  Taking advantage of that, we
  228.            could make this loop make two iterations less.  */
  229.  
  230.         cy_limb = mpn_submul_1 (np, dp, dsize, q);
  231.  
  232.         if (n2 != cy_limb)
  233.           {
  234.         mpn_add_n (np, np, dp, dsize);
  235.         q--;
  236.           }
  237.  
  238.         qp[i] = q;
  239.         n0 = np[dsize - 1];
  240.       }
  241.       }
  242.     }
  243.  
  244.   return most_significant_q_limb;
  245. }
  246.