[BACK]Return to mul_n.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2-ssh-2 / mpn / generic

Annotation of OpenXM/src/kan96xx/gmp-2.0.2-ssh-2/mpn/generic/mul_n.c, Revision 1.1.1.1

1.1       takayama    1: /* mpn_mul_n -- Multiply two natural numbers of length n.
                      2:
                      3: Copyright (C) 1991, 1992, 1993, 1994, 1996 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 Library General Public License as published by
                      9: the Free Software Foundation; either version 2 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 Library General Public
                     15: License for more details.
                     16:
                     17: You should have received a copy of the GNU Library 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: /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
                     26:    both with SIZE limbs, and store the result at PRODP.  2 * SIZE limbs are
                     27:    always stored.  Return the most significant limb.
                     28:
                     29:    Argument constraints:
                     30:    1. PRODP != UP and PRODP != VP, i.e. the destination
                     31:       must be distinct from the multiplier and the multiplicand.  */
                     32:
                     33: /* If KARATSUBA_THRESHOLD is not already defined, define it to a
                     34:    value which is good on most machines.  */
                     35: #ifndef KARATSUBA_THRESHOLD
                     36: #define KARATSUBA_THRESHOLD 32
                     37: #endif
                     38:
                     39: /* The code can't handle KARATSUBA_THRESHOLD smaller than 2.  */
                     40: #if KARATSUBA_THRESHOLD < 2
                     41: #undef KARATSUBA_THRESHOLD
                     42: #define KARATSUBA_THRESHOLD 2
                     43: #endif
                     44:
                     45: /* Handle simple cases with traditional multiplication.
                     46:
                     47:    This is the most critical code of multiplication.  All multiplies rely
                     48:    on this, both small and huge.  Small ones arrive here immediately.  Huge
                     49:    ones arrive here as this is the base case for Karatsuba's recursive
                     50:    algorithm below.  */
                     51:
                     52: void
                     53: #if __STDC__
                     54: impn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
                     55: #else
                     56: impn_mul_n_basecase (prodp, up, vp, size)
                     57:      mp_ptr prodp;
                     58:      mp_srcptr up;
                     59:      mp_srcptr vp;
                     60:      mp_size_t size;
                     61: #endif
                     62: {
                     63:   mp_size_t i;
                     64:   mp_limb_t cy_limb;
                     65:   mp_limb_t v_limb;
                     66:
                     67:   /* Multiply by the first limb in V separately, as the result can be
                     68:      stored (not added) to PROD.  We also avoid a loop for zeroing.  */
                     69:   v_limb = vp[0];
                     70:   if (v_limb <= 1)
                     71:     {
                     72:       if (v_limb == 1)
                     73:        MPN_COPY (prodp, up, size);
                     74:       else
                     75:        MPN_ZERO (prodp, size);
                     76:       cy_limb = 0;
                     77:     }
                     78:   else
                     79:     cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
                     80:
                     81:   prodp[size] = cy_limb;
                     82:   prodp++;
                     83:
                     84:   /* For each iteration in the outer loop, multiply one limb from
                     85:      U with one limb from V, and add it to PROD.  */
                     86:   for (i = 1; i < size; i++)
                     87:     {
                     88:       v_limb = vp[i];
                     89:       if (v_limb <= 1)
                     90:        {
                     91:          cy_limb = 0;
                     92:          if (v_limb == 1)
                     93:            cy_limb = mpn_add_n (prodp, prodp, up, size);
                     94:        }
                     95:       else
                     96:        cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
                     97:
                     98:       prodp[size] = cy_limb;
                     99:       prodp++;
                    100:     }
                    101: }
                    102:
                    103: void
                    104: #if __STDC__
                    105: impn_mul_n (mp_ptr prodp,
                    106:             mp_srcptr up, mp_srcptr vp, mp_size_t size, mp_ptr tspace)
                    107: #else
                    108: impn_mul_n (prodp, up, vp, size, tspace)
                    109:      mp_ptr prodp;
                    110:      mp_srcptr up;
                    111:      mp_srcptr vp;
                    112:      mp_size_t size;
                    113:      mp_ptr tspace;
                    114: #endif
                    115: {
                    116:   if ((size & 1) != 0)
                    117:     {
                    118:       /* The size is odd, the code code below doesn't handle that.
                    119:         Multiply the least significant (size - 1) limbs with a recursive
                    120:         call, and handle the most significant limb of S1 and S2
                    121:         separately.  */
                    122:       /* A slightly faster way to do this would be to make the Karatsuba
                    123:         code below behave as if the size were even, and let it check for
                    124:         odd size in the end.  I.e., in essence move this code to the end.
                    125:         Doing so would save us a recursive call, and potentially make the
                    126:         stack grow a lot less.  */
                    127:
                    128:       mp_size_t esize = size - 1;      /* even size */
                    129:       mp_limb_t cy_limb;
                    130:
                    131:       MPN_MUL_N_RECURSE (prodp, up, vp, esize, tspace);
                    132:       cy_limb = mpn_addmul_1 (prodp + esize, up, esize, vp[esize]);
                    133:       prodp[esize + esize] = cy_limb;
                    134:       cy_limb = mpn_addmul_1 (prodp + esize, vp, size, up[esize]);
                    135:
                    136:       prodp[esize + size] = cy_limb;
                    137:     }
                    138:   else
                    139:     {
                    140:       /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
                    141:
                    142:         Split U in two pieces, U1 and U0, such that
                    143:         U = U0 + U1*(B**n),
                    144:         and V in V1 and V0, such that
                    145:         V = V0 + V1*(B**n).
                    146:
                    147:         UV is then computed recursively using the identity
                    148:
                    149:                2n   n          n                     n
                    150:         UV = (B  + B )U V  +  B (U -U )(V -V )  +  (B + 1)U V
                    151:                        1 1        1  0   0  1              0 0
                    152:
                    153:         Where B = 2**BITS_PER_MP_LIMB.  */
                    154:
                    155:       mp_size_t hsize = size >> 1;
                    156:       mp_limb_t cy;
                    157:       int negflg;
                    158:
                    159:       /*** Product H.   ________________  ________________
                    160:                        |_____U1 x V1____||____U0 x V0_____|  */
                    161:       /* Put result in upper part of PROD and pass low part of TSPACE
                    162:         as new TSPACE.  */
                    163:       MPN_MUL_N_RECURSE (prodp + size, up + hsize, vp + hsize, hsize, tspace);
                    164:
                    165:       /*** Product M.   ________________
                    166:                        |_(U1-U0)(V0-V1)_|  */
                    167:       if (mpn_cmp (up + hsize, up, hsize) >= 0)
                    168:        {
                    169:          mpn_sub_n (prodp, up + hsize, up, hsize);
                    170:          negflg = 0;
                    171:        }
                    172:       else
                    173:        {
                    174:          mpn_sub_n (prodp, up, up + hsize, hsize);
                    175:          negflg = 1;
                    176:        }
                    177:       if (mpn_cmp (vp + hsize, vp, hsize) >= 0)
                    178:        {
                    179:          mpn_sub_n (prodp + hsize, vp + hsize, vp, hsize);
                    180:          negflg ^= 1;
                    181:        }
                    182:       else
                    183:        {
                    184:          mpn_sub_n (prodp + hsize, vp, vp + hsize, hsize);
                    185:          /* No change of NEGFLG.  */
                    186:        }
                    187:       /* Read temporary operands from low part of PROD.
                    188:         Put result in low part of TSPACE using upper part of TSPACE
                    189:         as new TSPACE.  */
                    190:       MPN_MUL_N_RECURSE (tspace, prodp, prodp + hsize, hsize, tspace + size);
                    191:
                    192:       /*** Add/copy product H.  */
                    193:       MPN_COPY (prodp + hsize, prodp + size, hsize);
                    194:       cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
                    195:
                    196:       /*** Add product M (if NEGFLG M is a negative number).  */
                    197:       if (negflg)
                    198:        cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
                    199:       else
                    200:        cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
                    201:
                    202:       /*** Product L.   ________________  ________________
                    203:                        |________________||____U0 x V0_____|  */
                    204:       /* Read temporary operands from low part of PROD.
                    205:         Put result in low part of TSPACE using upper part of TSPACE
                    206:         as new TSPACE.  */
                    207:       MPN_MUL_N_RECURSE (tspace, up, vp, hsize, tspace + size);
                    208:
                    209:       /*** Add/copy Product L (twice).  */
                    210:
                    211:       cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
                    212:       if (cy)
                    213:        mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
                    214:
                    215:       MPN_COPY (prodp, tspace, hsize);
                    216:       cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
                    217:       if (cy)
                    218:        mpn_add_1 (prodp + size, prodp + size, size, 1);
                    219:     }
                    220: }
                    221:
                    222: void
                    223: #if __STDC__
                    224: impn_sqr_n_basecase (mp_ptr prodp, mp_srcptr up, mp_size_t size)
                    225: #else
                    226: impn_sqr_n_basecase (prodp, up, size)
                    227:      mp_ptr prodp;
                    228:      mp_srcptr up;
                    229:      mp_size_t size;
                    230: #endif
                    231: {
                    232:   mp_size_t i;
                    233:   mp_limb_t cy_limb;
                    234:   mp_limb_t v_limb;
                    235:
                    236:   /* Multiply by the first limb in V separately, as the result can be
                    237:      stored (not added) to PROD.  We also avoid a loop for zeroing.  */
                    238:   v_limb = up[0];
                    239:   if (v_limb <= 1)
                    240:     {
                    241:       if (v_limb == 1)
                    242:        MPN_COPY (prodp, up, size);
                    243:       else
                    244:        MPN_ZERO (prodp, size);
                    245:       cy_limb = 0;
                    246:     }
                    247:   else
                    248:     cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
                    249:
                    250:   prodp[size] = cy_limb;
                    251:   prodp++;
                    252:
                    253:   /* For each iteration in the outer loop, multiply one limb from
                    254:      U with one limb from V, and add it to PROD.  */
                    255:   for (i = 1; i < size; i++)
                    256:     {
                    257:       v_limb = up[i];
                    258:       if (v_limb <= 1)
                    259:        {
                    260:          cy_limb = 0;
                    261:          if (v_limb == 1)
                    262:            cy_limb = mpn_add_n (prodp, prodp, up, size);
                    263:        }
                    264:       else
                    265:        cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
                    266:
                    267:       prodp[size] = cy_limb;
                    268:       prodp++;
                    269:     }
                    270: }
                    271:
                    272: void
                    273: #if __STDC__
                    274: impn_sqr_n (mp_ptr prodp,
                    275:             mp_srcptr up, mp_size_t size, mp_ptr tspace)
                    276: #else
                    277: impn_sqr_n (prodp, up, size, tspace)
                    278:      mp_ptr prodp;
                    279:      mp_srcptr up;
                    280:      mp_size_t size;
                    281:      mp_ptr tspace;
                    282: #endif
                    283: {
                    284:   if ((size & 1) != 0)
                    285:     {
                    286:       /* The size is odd, the code code below doesn't handle that.
                    287:         Multiply the least significant (size - 1) limbs with a recursive
                    288:         call, and handle the most significant limb of S1 and S2
                    289:         separately.  */
                    290:       /* A slightly faster way to do this would be to make the Karatsuba
                    291:         code below behave as if the size were even, and let it check for
                    292:         odd size in the end.  I.e., in essence move this code to the end.
                    293:         Doing so would save us a recursive call, and potentially make the
                    294:         stack grow a lot less.  */
                    295:
                    296:       mp_size_t esize = size - 1;      /* even size */
                    297:       mp_limb_t cy_limb;
                    298:
                    299:       MPN_SQR_N_RECURSE (prodp, up, esize, tspace);
                    300:       cy_limb = mpn_addmul_1 (prodp + esize, up, esize, up[esize]);
                    301:       prodp[esize + esize] = cy_limb;
                    302:       cy_limb = mpn_addmul_1 (prodp + esize, up, size, up[esize]);
                    303:
                    304:       prodp[esize + size] = cy_limb;
                    305:     }
                    306:   else
                    307:     {
                    308:       mp_size_t hsize = size >> 1;
                    309:       mp_limb_t cy;
                    310:
                    311:       /*** Product H.   ________________  ________________
                    312:                        |_____U1 x U1____||____U0 x U0_____|  */
                    313:       /* Put result in upper part of PROD and pass low part of TSPACE
                    314:         as new TSPACE.  */
                    315:       MPN_SQR_N_RECURSE (prodp + size, up + hsize, hsize, tspace);
                    316:
                    317:       /*** Product M.   ________________
                    318:                        |_(U1-U0)(U0-U1)_|  */
                    319:       if (mpn_cmp (up + hsize, up, hsize) >= 0)
                    320:        {
                    321:          mpn_sub_n (prodp, up + hsize, up, hsize);
                    322:        }
                    323:       else
                    324:        {
                    325:          mpn_sub_n (prodp, up, up + hsize, hsize);
                    326:        }
                    327:
                    328:       /* Read temporary operands from low part of PROD.
                    329:         Put result in low part of TSPACE using upper part of TSPACE
                    330:         as new TSPACE.  */
                    331:       MPN_SQR_N_RECURSE (tspace, prodp, hsize, tspace + size);
                    332:
                    333:       /*** Add/copy product H.  */
                    334:       MPN_COPY (prodp + hsize, prodp + size, hsize);
                    335:       cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
                    336:
                    337:       /*** Add product M (if NEGFLG M is a negative number).  */
                    338:       cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
                    339:
                    340:       /*** Product L.   ________________  ________________
                    341:                        |________________||____U0 x U0_____|  */
                    342:       /* Read temporary operands from low part of PROD.
                    343:         Put result in low part of TSPACE using upper part of TSPACE
                    344:         as new TSPACE.  */
                    345:       MPN_SQR_N_RECURSE (tspace, up, hsize, tspace + size);
                    346:
                    347:       /*** Add/copy Product L (twice).  */
                    348:
                    349:       cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
                    350:       if (cy)
                    351:        mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
                    352:
                    353:       MPN_COPY (prodp, tspace, hsize);
                    354:       cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
                    355:       if (cy)
                    356:        mpn_add_1 (prodp + size, prodp + size, size, 1);
                    357:     }
                    358: }
                    359:
                    360: /* This should be made into an inline function in gmp.h.  */
                    361: void
                    362: #if __STDC__
                    363: mpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
                    364: #else
                    365: mpn_mul_n (prodp, up, vp, size)
                    366:      mp_ptr prodp;
                    367:      mp_srcptr up;
                    368:      mp_srcptr vp;
                    369:      mp_size_t size;
                    370: #endif
                    371: {
                    372:   TMP_DECL (marker);
                    373:   TMP_MARK (marker);
                    374:   if (up == vp)
                    375:     {
                    376:       if (size < KARATSUBA_THRESHOLD)
                    377:        {
                    378:          impn_sqr_n_basecase (prodp, up, size);
                    379:        }
                    380:       else
                    381:        {
                    382:          mp_ptr tspace;
                    383:          tspace = (mp_ptr) TMP_ALLOC (2 * size * BYTES_PER_MP_LIMB);
                    384:          impn_sqr_n (prodp, up, size, tspace);
                    385:        }
                    386:     }
                    387:   else
                    388:     {
                    389:       if (size < KARATSUBA_THRESHOLD)
                    390:        {
                    391:          impn_mul_n_basecase (prodp, up, vp, size);
                    392:        }
                    393:       else
                    394:        {
                    395:          mp_ptr tspace;
                    396:          tspace = (mp_ptr) TMP_ALLOC (2 * size * BYTES_PER_MP_LIMB);
                    397:          impn_mul_n (prodp, up, vp, size, tspace);
                    398:        }
                    399:     }
                    400:   TMP_FREE (marker);
                    401: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>