Return to dive_1.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/dive_1.c, Revision 1.1

1.1     ! ohara       1: /* mpn_divexact_1 -- mpn by limb exact division.
        !             2:
        !             3:    THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
        !             4:    CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
        !             5:    FUTURE GNU MP RELEASES.
        !             6:
        !             7: Copyright 2000, 2001, 2002 Free Software Foundation, Inc.
        !             8:
        !             9: This file is part of the GNU MP Library.
        !            10:
        !            11: The GNU MP Library is free software; you can redistribute it and/or modify
        !            12: it under the terms of the GNU Lesser General Public License as published by
        !            13: the Free Software Foundation; either version 2.1 of the License, or (at your
        !            14: option) any later version.
        !            15:
        !            16: The GNU MP Library is distributed in the hope that it will be useful, but
        !            17: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
        !            18: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
        !            19: License for more details.
        !            20:
        !            21: You should have received a copy of the GNU Lesser General Public License
        !            22: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
        !            23: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
        !            24: MA 02111-1307, USA. */
        !            25:
        !            26: #include "gmp.h"
        !            27: #include "gmp-impl.h"
        !            28: #include "longlong.h"
        !            29:
        !            30:
        !            31:
        !            32: /* Divide a={src,size} by d=divisor and store the quotient in q={dst,size}.
        !            33:    q will only be correct if d divides a exactly.
        !            34:
        !            35:    A separate loop is used for shift==0 because n<<BITS_PER_MP_LIMB doesn't
        !            36:    give zero on all CPUs (for instance it doesn't on the x86s).  This
        !            37:    separate loop might run faster too, helping odd divisors.
        !            38:
        !            39:    Possibilities:
        !            40:
        !            41:    mpn_divexact_1c could be created, accepting and returning c.  This would
        !            42:    let a long calculation be done piece by piece.  Currently there's no
        !            43:    particular need for that, and not returning c means that a final umul can
        !            44:    be skipped.
        !            45:
        !            46:    Another use for returning c would be letting the caller know whether the
        !            47:    division was in fact exact.  It would work just to return the carry bit
        !            48:    "c=(l>s)" and let the caller do a final umul if interested.
        !            49:
        !            50:    When the divisor is even, the factors of two could be handled with a
        !            51:    separate mpn_rshift, instead of shifting on the fly.  That might be
        !            52:    faster on some CPUs and would mean just the shift==0 style loop would be
        !            53:    needed.
        !            54:
        !            55:    If n<<BITS_PER_MP_LIMB gives zero on a particular CPU then the separate
        !            56:    shift==0 loop is unnecessary, and could be eliminated if there's no great
        !            57:    speed difference.
        !            58:
        !            59:    It's not clear whether "/" is the best way to handle size==1.  Alpha gcc
        !            60:    2.95 for instance has a poor "/" and might prefer the modular method.
        !            61:    Perhaps a tuned parameter should control this.
        !            62:
        !            63:    If src[size-1] < divisor then dst[size-1] will be zero, and one divide
        !            64:    step could be skipped.  A test at last step for s<divisor (or ls in the
        !            65:    even case) might be a good way to do that.  But if this code is often
        !            66:    used with small divisors then it might not be worth bothering  */
        !            67:
        !            68: void
        !            69: mpn_divexact_1 (mp_ptr dst, mp_srcptr src, mp_size_t size, mp_limb_t divisor)
        !            70: {
        !            71:   mp_size_t  i;
        !            72:   mp_limb_t  c, l, ls, s, s_next, inverse, dummy;
        !            73:   unsigned   shift;
        !            74:
        !            75:   ASSERT (size >= 1);
        !            76:   ASSERT (divisor != 0);
        !            77:   ASSERT (MPN_SAME_OR_SEPARATE_P (dst, src, size));
        !            78:   ASSERT_MPN (src, size);
        !            79:   ASSERT_LIMB (divisor);
        !            80:
        !            81:   if (size == 1)
        !            82:     {
        !            83:       dst[0] = src[0] / divisor;
        !            84:       return;
        !            85:     }
        !            86:
        !            87:   if ((divisor & 1) == 0)
        !            88:     {
        !            89:       count_trailing_zeros (shift, divisor);
        !            90:       divisor >>= shift;
        !            91:     }
        !            92:   else
        !            93:     shift = 0;
        !            94:
        !            95:   modlimb_invert (inverse, divisor);
        !            96:   divisor <<= GMP_NAIL_BITS;
        !            97:
        !            98:   if (shift != 0)
        !            99:     {
        !           100:       s = src[0];
        !           101:       c = 0;
        !           102:       i = 0;
        !           103:       size--;
        !           104:       goto even_entry;
        !           105:
        !           106:       do
        !           107:        {
        !           108:          umul_ppmm (l, dummy, l, divisor);
        !           109:          c += l;
        !           110:
        !           111:        even_entry:
        !           112:          s_next = src[i+1];
        !           113:          ls = ((s >> shift) | (s_next << (GMP_NUMB_BITS-shift))) & GMP_NUMB_MASK;
        !           114:          s = s_next;
        !           115:
        !           116:           SUBC_LIMB (c, l, ls, c);
        !           117:
        !           118:          l = (l * inverse) & GMP_NUMB_MASK;
        !           119:          dst[i] = l;
        !           120:          i++;
        !           121:        }
        !           122:       while (i < size);
        !           123:
        !           124:       umul_ppmm (l, dummy, l, divisor);
        !           125:       c += l;
        !           126:       ls = s >> shift;
        !           127:       l = ls - c;
        !           128:       l = (l * inverse) & GMP_NUMB_MASK;
        !           129:       dst[i] = l;
        !           130:     }
        !           131:   else
        !           132:     {
        !           133:       l = src[0];
        !           134:       l = (l * inverse) & GMP_NUMB_MASK;
        !           135:       dst[0] = l;
        !           136:       i = 1;
        !           137:       c = 0;
        !           138:
        !           139:       do
        !           140:        {
        !           141:          umul_ppmm (l, dummy, l, divisor);
        !           142:          c += l;
        !           143:
        !           144:          s = src[i];
        !           145:           SUBC_LIMB (c, l, s, c);
        !           146:
        !           147:           l = (l * inverse) & GMP_NUMB_MASK;
        !           148:          dst[i] = l;
        !           149:          i++;
        !           150:        }
        !           151:       while (i < size);
        !           152:     }
        !           153: }