Return to addmul_1.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / cray / ieee

Annotation of OpenXM_contrib/gmp/mpn/cray/ieee/addmul_1.c, Revision 1.1

1.1     ! ohara       1: /* Cray PVP/IEEE mpn_addmul_1 -- multiply a limb vector with a limb and add the
        !             2:    result to a second limb vector.
        !             3:
        !             4: Copyright 2000, 2001, 2002 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 Lesser General Public License as published by
        !            10: the Free Software Foundation; either version 2.1 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 Lesser General Public
        !            16: License for more details.
        !            17:
        !            18: You should have received a copy of the GNU Lesser 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: /* This code runs at just under 9 cycles/limb on a T90.  That is not perfect,
        !            24:    mainly due to vector register shortage in the main loop.  Assembly code
        !            25:    should bring it down to perhaps 7 cycles/limb.  */
        !            26:
        !            27: #include <intrinsics.h>
        !            28: #include "gmp.h"
        !            29: #include "gmp-impl.h"
        !            30:
        !            31: mp_limb_t
        !            32: mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
        !            33: {
        !            34:   mp_limb_t cy[n];
        !            35:   mp_limb_t a, b, r, s0, s1, c0, c1;
        !            36:   mp_size_t i;
        !            37:   int more_carries;
        !            38:
        !            39:   if (up == rp)
        !            40:     {
        !            41:       /* The algorithm used below cannot handle overlap.  Handle it here by
        !            42:         making a temporary copy of the source vector, then call ourselves.  */
        !            43:       mp_limb_t xp[n];
        !            44:       MPN_COPY (xp, up, n);
        !            45:       return mpn_addmul_1 (rp, xp, n, vl);
        !            46:     }
        !            47:
        !            48:   a = up[0] * vl;
        !            49:   r = rp[0];
        !            50:   s0 = a + r;
        !            51:   rp[0] = s0;
        !            52:   c0 = ((a & r) | ((a | r) & ~s0)) >> 63;
        !            53:   cy[0] = c0;
        !            54:
        !            55:   /* Main multiply loop.  Generate a raw accumulated output product in rp[]
        !            56:      and a carry vector in cy[].  */
        !            57: #pragma _CRI ivdep
        !            58:   for (i = 1; i < n; i++)
        !            59:     {
        !            60:       a = up[i] * vl;
        !            61:       b = _int_mult_upper (up[i - 1], vl);
        !            62:       s0 = a + b;
        !            63:       c0 = ((a & b) | ((a | b) & ~s0)) >> 63;
        !            64:       r = rp[i];
        !            65:       s1 = s0 + r;
        !            66:       rp[i] = s1;
        !            67:       c1 = ((s0 & r) | ((s0 | r) & ~s1)) >> 63;
        !            68:       cy[i] = c0 + c1;
        !            69:     }
        !            70:   /* Carry add loop.  Add the carry vector cy[] to the raw result rp[] and
        !            71:      store the new result back to rp[].  */
        !            72:   more_carries = 0;
        !            73: #pragma _CRI ivdep
        !            74:   for (i = 1; i < n; i++)
        !            75:     {
        !            76:       r = rp[i];
        !            77:       c0 = cy[i - 1];
        !            78:       s0 = r + c0;
        !            79:       rp[i] = s0;
        !            80:       c0 = (r & ~s0) >> 63;
        !            81:       more_carries += c0;
        !            82:     }
        !            83:   /* If that second loop generated carry, handle that in scalar loop.  */
        !            84:   if (more_carries)
        !            85:     {
        !            86:       mp_limb_t cyrec = 0;
        !            87:       /* Look for places where rp[k] == 0 and cy[k-1] == 1 or
        !            88:         rp[k] == 1 and cy[k-1] == 2.
        !            89:         These are where we got a recurrency carry.  */
        !            90:       for (i = 1; i < n; i++)
        !            91:        {
        !            92:          r = rp[i];
        !            93:          c0 = r < cy[i - 1];
        !            94:          s0 = r + cyrec;
        !            95:          rp[i] = s0;
        !            96:          c1 = (r & ~s0) >> 63;
        !            97:          cyrec = c0 | c1;
        !            98:        }
        !            99:       return _int_mult_upper (up[n - 1], vl) + cyrec + cy[n - 1];
        !           100:     }
        !           101:
        !           102:   return _int_mult_upper (up[n - 1], vl) + cy[n - 1];
        !           103: }