Return to aorsmul.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/aorsmul.c, Revision 1.1.1.1

1.1       ohara       1: /* mpz_addmul, mpz_submul -- add or subtract multiple.
                      2:
                      3: Copyright 2001 Free Software Foundation, Inc.
                      4:
                      5: This file is part of the GNU MP Library.
                      6:
                      7: The GNU MP Library is free software; you can redistribute it and/or modify
                      8: it under the terms of the GNU Lesser General Public License as published by
                      9: the Free Software Foundation; either version 2.1 of the License, or (at your
                     10: option) any later version.
                     11:
                     12: The GNU MP Library is distributed in the hope that it will be useful, but
                     13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
                     15: License for more details.
                     16:
                     17: You should have received a copy of the GNU Lesser General Public License
                     18: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     20: MA 02111-1307, USA. */
                     21:
                     22: #include "gmp.h"
                     23: #include "gmp-impl.h"
                     24:
                     25:
                     26: /* expecting x and y both with non-zero high limbs */
                     27: #define mpn_cmp_twosizes_lt(xp,xsize, yp,ysize)                 \
                     28:   ((xsize) < (ysize)                                            \
                     29:    || ((xsize) == (ysize) && mpn_cmp (xp, yp, xsize) < 0))
                     30:
                     31:
                     32: /* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y.
                     33:
                     34:    The signs of w, x and y are fully accounted for by each flipping "sub".
                     35:
                     36:    The sign of w is retained for the result, unless the absolute value
                     37:    submul underflows, in which case it flips.  */
                     38:
                     39: static void __gmpz_aorsmul _PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub))) REGPARM_ATTR (1);
                     40: #define mpz_aorsmul(w,x,y,sub)  __gmpz_aorsmul (REGPARM_3_1 (w, x, y, sub))
                     41:
                     42: static void
                     43: mpz_aorsmul (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub)
                     44: {
                     45:   mp_size_t  xsize, ysize, tsize, wsize, wsize_signed;
                     46:   mp_ptr     wp, tp;
                     47:   mp_limb_t  c, high;
                     48:   TMP_DECL (marker);
                     49:
                     50:   /* w unaffected if x==0 or y==0 */
                     51:   xsize = SIZ(x);
                     52:   ysize = SIZ(y);
                     53:   if (xsize == 0 || ysize == 0)
                     54:     return;
                     55:
                     56:   /* make x the bigger of the two */
                     57:   if (ABS(ysize) > ABS(xsize))
                     58:     {
                     59:       MPZ_SRCPTR_SWAP (x, y);
                     60:       MP_SIZE_T_SWAP (xsize, ysize);
                     61:     }
                     62:
                     63:   sub ^= ysize;
                     64:   ysize = ABS(ysize);
                     65:
                     66:   /* use mpn_addmul_1/mpn_submul_1 if possible */
                     67:   if (ysize == 1)
                     68:     {
                     69:       mpz_aorsmul_1 (w, x, PTR(y)[0], sub);
                     70:       return;
                     71:     }
                     72:
                     73:   sub ^= xsize;
                     74:   xsize = ABS(xsize);
                     75:
                     76:   wsize_signed = SIZ(w);
                     77:   sub ^= wsize_signed;
                     78:   wsize = ABS(wsize_signed);
                     79:
                     80:   tsize = xsize + ysize;
                     81:   MPZ_REALLOC (w, MAX (wsize, tsize) + 1);
                     82:   wp = PTR(w);
                     83:
                     84:   if (wsize_signed == 0)
                     85:     {
                     86:       /* Nothing to add to, just set w=x*y.  No w==x or w==y overlap here,
                     87:          since we know x,y!=0 but w==0.  */
                     88:       high = mpn_mul (wp, PTR(x),xsize, PTR(y),ysize);
                     89:       tsize -= (high == 0);
                     90:       SIZ(w) = (sub >= 0 ? tsize : -tsize);
                     91:       return;
                     92:     }
                     93:
                     94:   TMP_MARK (marker);
                     95:   tp = TMP_ALLOC_LIMBS (tsize);
                     96:
                     97:   high = mpn_mul (tp, PTR(x),xsize, PTR(y),ysize);
                     98:   tsize -= (high == 0);
                     99:   ASSERT (tp[tsize-1] != 0);
                    100:   if (sub >= 0)
                    101:     {
                    102:       mp_srcptr up    = wp;
                    103:       mp_size_t usize = wsize;
                    104:
                    105:       if (usize < tsize)
                    106:         {
                    107:           up    = tp;
                    108:           usize = tsize;
                    109:           tp    = wp;
                    110:           tsize = wsize;
                    111:
                    112:           wsize = usize;
                    113:         }
                    114:
                    115:       c = mpn_add (wp, up,usize, tp,tsize);
                    116:       wp[wsize] = c;
                    117:       wsize += (c != 0);
                    118:     }
                    119:   else
                    120:     {
                    121:       mp_srcptr up    = wp;
                    122:       mp_size_t usize = wsize;
                    123:
                    124:       if (mpn_cmp_twosizes_lt (up,usize, tp,tsize))
                    125:         {
                    126:           up    = tp;
                    127:           usize = tsize;
                    128:           tp    = wp;
                    129:           tsize = wsize;
                    130:
                    131:           wsize = usize;
                    132:           wsize_signed = -wsize_signed;
                    133:         }
                    134:
                    135:       ASSERT_NOCARRY (mpn_sub (wp, up,usize, tp,tsize));
                    136:       wsize = usize;
                    137:       MPN_NORMALIZE (wp, wsize);
                    138:     }
                    139:
                    140:   SIZ(w) = (wsize_signed >= 0 ? wsize : -wsize);
                    141:
                    142:   TMP_FREE (marker);
                    143: }
                    144:
                    145:
                    146: void
                    147: mpz_addmul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v)
                    148: {
                    149:   mpz_aorsmul (w, u, v, (mp_size_t) 0);
                    150: }
                    151:
                    152: void
                    153: mpz_submul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v)
                    154: {
                    155:   mpz_aorsmul (w, u, v, (mp_size_t) -1);
                    156: }