Return to mod_34lsub1.c CVS log

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

Annotation of OpenXM_contrib/gmp/mpn/generic/mod_34lsub1.c, Revision 1.1.1.1

1.1       ohara       1: /* mpn_mod_34lsub1 -- remainder modulo 2^(BITS_PER_MP_LIMB*3/4)-1.
                      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:
                     27: #include "gmp.h"
                     28: #include "gmp-impl.h"
                     29:
                     30:
                     31: /* Calculate a remainder from {p,n} divided by 2^(BITS_PER_MP_LIMB*3/4)-1.
                     32:    The remainder is not fully reduced, it's any limb value congruent to
                     33:    {p,n} modulo that divisor.
                     34:
                     35:    This implementation is only correct when BITS_PER_MP_LIMB is a multiple
                     36:    of 4, but that suffices for all current uses.  */
                     37:
                     38: #if BITS_PER_MP_LIMB % 4 == 0
                     39:
                     40: #define B1  (BITS_PER_MP_LIMB / 4)
                     41: #define B2  (B1 * 2)
                     42: #define B3  (B1 * 3)
                     43:
                     44: #define M1  ((CNST_LIMB(1) << B1) - 1)
                     45: #define M2  ((CNST_LIMB(1) << B2) - 1)
                     46: #define M3  ((CNST_LIMB(1) << B3) - 1)
                     47:
                     48: #define LOW0(n)      ((n) & M3)
                     49: #define HIGH0(n)     ((n) >> B3)
                     50:
                     51: #define LOW1(n)      (((n) & M2) << B1)
                     52: #define HIGH1(n)     ((n) >> B2)
                     53:
                     54: #define LOW2(n)      (((n) & M1) << B2)
                     55: #define HIGH2(n)     ((n) >> B1)
                     56:
                     57: #define PARTS0(n)    (LOW0(n) + HIGH0(n))
                     58: #define PARTS1(n)    (LOW1(n) + HIGH1(n))
                     59: #define PARTS2(n)    (LOW2(n) + HIGH2(n))
                     60:
                     61: #define ADD(c,a,val)    \
                     62:   do {                  \
                     63:     mp_limb_t  l = val; \
                     64:     a += l;             \
                     65:     c += (a < l);       \
                     66:   } while (0)
                     67:
                     68: mp_limb_t
                     69: mpn_mod_34lsub1 (mp_srcptr p, mp_size_t n)
                     70: {
                     71:   mp_limb_t  c0 = 0;
                     72:   mp_limb_t  c1 = 0;
                     73:   mp_limb_t  c2 = 0;
                     74:   mp_limb_t  a0, a1, a2;
                     75:
                     76:   ASSERT (n >= 1);
                     77:
                     78:   a0 = a1 = a2 = 0;
                     79:   c0 = c1 = c2 = 0;
                     80:
                     81:   while ((n -= 3) >= 0)
                     82:     {
                     83:       ADD (c0, a0, p[0]);
                     84:       ADD (c1, a1, p[1]);
                     85:       ADD (c2, a2, p[2]);
                     86:       p += 3;
                     87:     }
                     88:
                     89:   if (n != -3)
                     90:     {
                     91:       ADD (c0, a0, p[0]);
                     92:       if (n != -2)
                     93:        ADD (c1, a1, p[1]);
                     94:     }
                     95:
                     96:   return
                     97:     PARTS0 (a0) + PARTS1 (a1) + PARTS2 (a2)
                     98:     + PARTS1 (c0) + PARTS2 (c1) + PARTS0 (c2);
                     99: }
                    100:
                    101: #endif