home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bbs / gnu / gmp-1.3.2-src.lha / src / amiga / gmp-1.3.2 / mpn_dm_1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-03  |  5.2 KB  |  186 lines
  1. /* mpn_divmod_1(quot_ptr, dividend_ptr, dividend_size, divisor_limb) --
  2.    Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
  3.    Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
  4.    Return the single-limb remainder.
  5.    There are no constraints on the value of the divisor.
  6.  
  7.    QUOT_PTR and DIVIDEND_PTR might point to the same limb.
  8.  
  9. Copyright (C) 1991, 1993 Free Software Foundation, Inc.
  10.  
  11. This file is part of the GNU MP Library.
  12.  
  13. The GNU MP Library is free software; you can redistribute it and/or modify
  14. it under the terms of the GNU General Public License as published by
  15. the Free Software Foundation; either version 2, or (at your option)
  16. any later version.
  17.  
  18. The GNU MP Library is distributed in the hope that it will be useful,
  19. but WITHOUT ANY WARRANTY; without even the implied warranty of
  20. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21. GNU General Public License for more details.
  22.  
  23. You should have received a copy of the GNU General Public License
  24. along with the GNU MP Library; see the file COPYING.  If not, write to
  25. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  26.  
  27. #include "gmp.h"
  28. #include "gmp-impl.h"
  29. #include "longlong.h"
  30.  
  31. #ifndef UMUL_TIME
  32. #define UMUL_TIME 1
  33. #endif
  34.  
  35. #ifndef UDIV_TIME
  36. #define UDIV_TIME UMUL_TIME
  37. #endif
  38.  
  39. #if UDIV_TIME > 2 * UMUL_TIME
  40. #undef UDIV_NEEDS_NORMALIZATION
  41. #define UDIV_NEEDS_NORMALIZATION 1
  42. #endif
  43.  
  44. #define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
  45.   do {                                    \
  46.     unsigned long int _q, _ql, _r;                    \
  47.     unsigned long int _xh, _xl;                        \
  48.     umul_ppmm (_q, _ql, (nh), (di));                    \
  49.     _q += (nh);            /* DI is 2**BITS_PER_MP_LIMB too small.  */\
  50.     umul_ppmm (_xh, _xl, _q, (d));                    \
  51.     sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl);                \
  52.     if (_xh != 0)                            \
  53.       {                                    \
  54.     sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                \
  55.     _q += 1;                            \
  56.     if (_xh != 0)                            \
  57.       {                                \
  58.         sub_ddmmss (_xh, _r, _xh, _r, 0, (d));            \
  59.         _q += 1;                            \
  60.       }                                \
  61.       }                                    \
  62.     if (_r >= (d))                            \
  63.       {                                    \
  64.     _r -= (d);                            \
  65.     _q += 1;                            \
  66.       }                                    \
  67.     (r) = _r;                                \
  68.     (q) = _q;                                \
  69.   } while (0)
  70.  
  71. mp_limb
  72. #ifdef __STDC__
  73. mpn_divmod_1 (mp_ptr quot_ptr,
  74.           mp_srcptr dividend_ptr, mp_size dividend_size,
  75.           unsigned long int divisor_limb)
  76. #else
  77. mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
  78.      mp_ptr quot_ptr;
  79.      mp_srcptr dividend_ptr;
  80.      mp_size dividend_size;
  81.      unsigned long int divisor_limb;
  82. #endif
  83. {
  84.   mp_size i;
  85.   mp_limb n1, n0, r;
  86.  
  87.   /* Botch: Should this be handled at all?  Rely on callers?  */
  88.   if (dividend_size == 0)
  89.     return 0;
  90.  
  91.   if (UDIV_NEEDS_NORMALIZATION)
  92.     {
  93.       int normalization_steps;
  94.  
  95.       count_leading_zeros (normalization_steps, divisor_limb);
  96.       if (normalization_steps != 0)
  97.     {
  98.       divisor_limb <<= normalization_steps;
  99.  
  100.       n1 = dividend_ptr[dividend_size - 1];
  101.       r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
  102.  
  103.       /* Possible optimization:
  104.          if (r == 0
  105.          && divisor_limb > ((n1 << normalization_steps)
  106.                  | (dividend_ptr[dividend_size - 2] >> ...)))
  107.          ...one division less...
  108.          [Don't forget to zero most sign. quotient limb!]  */
  109.  
  110.       /* If multiplication is much faster than division, and the
  111.          dividend is large, pre-invert the divisor, and use
  112.          only multiplications in the inner loop.  */
  113.       if (UDIV_TIME > 2 * UMUL_TIME && dividend_size >= 4)
  114.         {
  115.           mp_limb divisor_limb_inverted;
  116.           int dummy;
  117.           /* Compute (2**64 - 2**32 * DIVISOR_LIMB) / DIVISOR_LIMB.
  118.          The result is an 33-bit approximation to 1/DIVISOR_LIMB,
  119.          with the most significant bit (weight 2**32) implicit.  */
  120.  
  121.           /* Special case for DIVISOR_LIMB == 100...000.  */
  122.           if (divisor_limb << 1 == 0)
  123.         divisor_limb_inverted = ~0;
  124.           else
  125.         udiv_qrnnd (divisor_limb_inverted, dummy,
  126.                 -divisor_limb, 0, divisor_limb);
  127.  
  128.           for (i = dividend_size - 2; i >= 0; i--)
  129.         {
  130.           n0 = dividend_ptr[i];
  131.           udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
  132.                      ((n1 << normalization_steps)
  133.                       | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
  134.                      divisor_limb, divisor_limb_inverted);
  135.           n1 = n0;
  136.         }
  137.           udiv_qrnnd_preinv (quot_ptr[0], r, r,
  138.               n1 << normalization_steps,
  139.               divisor_limb, divisor_limb_inverted);
  140.           return r >> normalization_steps;
  141.         }
  142.       else
  143.         {
  144.           for (i = dividend_size - 2; i >= 0; i--)
  145.         {
  146.           n0 = dividend_ptr[i];
  147.           udiv_qrnnd (quot_ptr[i + 1], r, r,
  148.                   ((n1 << normalization_steps)
  149.                    | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
  150.                   divisor_limb);
  151.           n1 = n0;
  152.         }
  153.           udiv_qrnnd (quot_ptr[0], r, r,
  154.               n1 << normalization_steps,
  155.               divisor_limb);
  156.           return r >> normalization_steps;
  157.         }
  158.     }
  159.     }
  160.  
  161.   /* No normalization needed, either because udiv_qrnnd doesn't require
  162.      it, or because DIVISOR_LIMB is already normalized.  */
  163.  
  164.   i = dividend_size - 1;
  165.   r = dividend_ptr[i];
  166.  
  167.   if (r >= divisor_limb)
  168.     {
  169.       r = 0;
  170.     }
  171.   else
  172.     {
  173.       /* Callers expect the quotient to be DIVIDEND_SIZE limbs.  Store
  174.      a leading zero to make that expectation come true.  */
  175.       quot_ptr[i] = 0;
  176.       i--;
  177.     }
  178.  
  179.   for (; i >= 0; i--)
  180.     {
  181.       n0 = dividend_ptr[i];
  182.       udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
  183.     }
  184.   return r;
  185. }
  186.