Blame mpz/prodlimbs.c

Packit 5c3484
/* mpz_prodlimbs(RESULT, V, LEN) -- Set RESULT to V[0]*V[1]*...*V[LEN-1].
Packit 5c3484
Packit 5c3484
Contributed to the GNU project by Marco Bodrato.
Packit 5c3484
Packit 5c3484
THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.
Packit 5c3484
IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.
Packit 5c3484
IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR
Packit 5c3484
DISAPPEAR IN A FUTURE GNU MP RELEASE.
Packit 5c3484
Packit 5c3484
Copyright 2010-2012 Free Software Foundation, Inc.
Packit 5c3484
Packit 5c3484
This file is part of the GNU MP Library.
Packit 5c3484
Packit 5c3484
The GNU MP Library is free software; you can redistribute it and/or modify
Packit 5c3484
it under the terms of either:
Packit 5c3484
Packit 5c3484
  * the GNU Lesser General Public License as published by the Free
Packit 5c3484
    Software Foundation; either version 3 of the License, or (at your
Packit 5c3484
    option) any later version.
Packit 5c3484
Packit 5c3484
or
Packit 5c3484
Packit 5c3484
  * the GNU General Public License as published by the Free Software
Packit 5c3484
    Foundation; either version 2 of the License, or (at your option) any
Packit 5c3484
    later version.
Packit 5c3484
Packit 5c3484
or both in parallel, as here.
Packit 5c3484
Packit 5c3484
The GNU MP Library is distributed in the hope that it will be useful, but
Packit 5c3484
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
Packit 5c3484
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
Packit 5c3484
for more details.
Packit 5c3484
Packit 5c3484
You should have received copies of the GNU General Public License and the
Packit 5c3484
GNU Lesser General Public License along with the GNU MP Library.  If not,
Packit 5c3484
see https://www.gnu.org/licenses/.  */
Packit 5c3484
Packit 5c3484
#include "gmp.h"
Packit 5c3484
#include "gmp-impl.h"
Packit 5c3484
Packit 5c3484
/*********************************************************/
Packit 5c3484
/* Section list-prod: product of a list -> mpz_t         */
Packit 5c3484
/*********************************************************/
Packit 5c3484
Packit 5c3484
/* FIXME: should be tuned */
Packit 5c3484
#ifndef RECURSIVE_PROD_THRESHOLD
Packit 5c3484
#define RECURSIVE_PROD_THRESHOLD (MUL_TOOM22_THRESHOLD)
Packit 5c3484
#endif
Packit 5c3484
Packit 5c3484
/* Computes the product of the j>1 limbs pointed by factors, puts the
Packit 5c3484
 * result in x. It assumes that all limbs are non-zero. Above
Packit 5c3484
 * Karatsuba's threshold it uses a binary splitting strategy, to gain
Packit 5c3484
 * speed by the asymptotically fast multiplication algorithms.
Packit 5c3484
 *
Packit 5c3484
 * The list in  {factors, j} is overwritten.
Packit 5c3484
 * Returns the size of the result
Packit 5c3484
 */
Packit 5c3484
Packit 5c3484
mp_size_t
Packit 5c3484
mpz_prodlimbs (mpz_ptr x, mp_ptr factors, mp_size_t j)
Packit 5c3484
{
Packit 5c3484
  mp_limb_t cy;
Packit 5c3484
  mp_size_t size, i;
Packit 5c3484
  mp_ptr    prod;
Packit 5c3484
Packit 5c3484
  ASSERT (j > 1);
Packit 5c3484
  ASSERT (RECURSIVE_PROD_THRESHOLD > 3);
Packit 5c3484
Packit 5c3484
  if (BELOW_THRESHOLD (j, RECURSIVE_PROD_THRESHOLD)) {
Packit 5c3484
    j--;
Packit 5c3484
    size = 1;
Packit 5c3484
Packit 5c3484
    for (i = 1; i < j; i++)
Packit 5c3484
      {
Packit 5c3484
	cy = mpn_mul_1 (factors, factors, size, factors[i]);
Packit 5c3484
	factors[size] = cy;
Packit 5c3484
	size += cy != 0;
Packit 5c3484
      };
Packit 5c3484
Packit 5c3484
    prod = MPZ_NEWALLOC (x, size + 1);
Packit 5c3484
Packit 5c3484
    cy = mpn_mul_1 (prod, factors, size, factors[i]);
Packit 5c3484
    prod[size] = cy;
Packit 5c3484
    return SIZ (x) = size + (cy != 0);
Packit 5c3484
  } else {
Packit 5c3484
    mpz_t x1, x2;
Packit 5c3484
    TMP_DECL;
Packit 5c3484
Packit 5c3484
    i = j >> 1;
Packit 5c3484
    j -= i;
Packit 5c3484
    TMP_MARK;
Packit 5c3484
Packit 5c3484
    MPZ_TMP_INIT (x2, j);
Packit 5c3484
Packit 5c3484
    PTR (x1) = factors + i;
Packit 5c3484
    ALLOC (x1) = j;
Packit 5c3484
    j = mpz_prodlimbs (x2, factors + i, j);
Packit 5c3484
    i = mpz_prodlimbs (x1, factors, i);
Packit 5c3484
    size = i + j;
Packit 5c3484
    prod = MPZ_NEWALLOC (x, size);
Packit 5c3484
    if (i >= j)
Packit 5c3484
      cy = mpn_mul (prod, PTR(x1), i, PTR(x2), j);
Packit 5c3484
    else
Packit 5c3484
      cy = mpn_mul (prod, PTR(x2), j, PTR(x1), i);
Packit 5c3484
    TMP_FREE;
Packit 5c3484
Packit 5c3484
    return SIZ (x) = size - (cy == 0);
Packit 5c3484
  }
Packit 5c3484
}