/src/libgmp/mpz/prodlimbs.c
Line | Count | Source |
1 | | /* mpz_prodlimbs(RESULT, V, LEN) -- Set RESULT to V[0]*V[1]*...*V[LEN-1]. |
2 | | |
3 | | Contributed to the GNU project by Marco Bodrato. |
4 | | |
5 | | THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. |
6 | | IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. |
7 | | IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR |
8 | | DISAPPEAR IN A FUTURE GNU MP RELEASE. |
9 | | |
10 | | Copyright 2010-2012 Free Software Foundation, Inc. |
11 | | |
12 | | This file is part of the GNU MP Library. |
13 | | |
14 | | The GNU MP Library is free software; you can redistribute it and/or modify |
15 | | it under the terms of either: |
16 | | |
17 | | * the GNU Lesser General Public License as published by the Free |
18 | | Software Foundation; either version 3 of the License, or (at your |
19 | | option) any later version. |
20 | | |
21 | | or |
22 | | |
23 | | * the GNU General Public License as published by the Free Software |
24 | | Foundation; either version 2 of the License, or (at your option) any |
25 | | later version. |
26 | | |
27 | | or both in parallel, as here. |
28 | | |
29 | | The GNU MP Library is distributed in the hope that it will be useful, but |
30 | | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
31 | | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
32 | | for more details. |
33 | | |
34 | | You should have received copies of the GNU General Public License and the |
35 | | GNU Lesser General Public License along with the GNU MP Library. If not, |
36 | | see https://www.gnu.org/licenses/. */ |
37 | | |
38 | | #include "gmp-impl.h" |
39 | | |
40 | | /*********************************************************/ |
41 | | /* Section list-prod: product of a list -> mpz_t */ |
42 | | /*********************************************************/ |
43 | | |
44 | | /* FIXME: should be tuned */ |
45 | | #ifndef RECURSIVE_PROD_THRESHOLD |
46 | | #define RECURSIVE_PROD_THRESHOLD (MUL_TOOM22_THRESHOLD) |
47 | | #endif |
48 | | |
49 | | /* Computes the product of the j>1 limbs pointed by factors, puts the |
50 | | * result in x. It assumes that all limbs are non-zero. Above |
51 | | * Karatsuba's threshold it uses a binary splitting strategy, to gain |
52 | | * speed by the asymptotically fast multiplication algorithms. |
53 | | * |
54 | | * The list in {factors, j} is overwritten. |
55 | | * Returns the size of the result |
56 | | */ |
57 | | |
58 | | mp_size_t |
59 | | mpz_prodlimbs (mpz_ptr x, mp_ptr factors, mp_size_t j) |
60 | 1.88k | { |
61 | 1.88k | mp_limb_t cy; |
62 | 1.88k | mp_size_t size, i; |
63 | 1.88k | mp_ptr prod; |
64 | | |
65 | 1.88k | ASSERT (j > 1); |
66 | 1.88k | ASSERT (RECURSIVE_PROD_THRESHOLD > 3); |
67 | | |
68 | 1.88k | if (BELOW_THRESHOLD (j, RECURSIVE_PROD_THRESHOLD)) { |
69 | 986 | j--; |
70 | 986 | size = 1; |
71 | | |
72 | 16.8k | for (i = 1; i < j; i++) |
73 | 15.8k | { |
74 | 15.8k | cy = mpn_mul_1 (factors, factors, size, factors[i]); |
75 | 15.8k | factors[size] = cy; |
76 | 15.8k | size += cy != 0; |
77 | 15.8k | }; |
78 | | |
79 | 986 | prod = MPZ_NEWALLOC (x, size + 1); |
80 | | |
81 | 986 | cy = mpn_mul_1 (prod, factors, size, factors[i]); |
82 | 986 | prod[size] = cy; |
83 | 986 | return SIZ (x) = size + (cy != 0); |
84 | 986 | } else { |
85 | 899 | mpz_t x1, x2; |
86 | 899 | TMP_DECL; |
87 | | |
88 | 899 | i = j >> 1; |
89 | 899 | j -= i; |
90 | 899 | TMP_MARK; |
91 | | |
92 | 899 | MPZ_TMP_INIT (x2, j); |
93 | | |
94 | 899 | PTR (x1) = factors + i; |
95 | 899 | ALLOC (x1) = j; |
96 | 899 | j = mpz_prodlimbs (x2, factors + i, j); |
97 | 899 | i = mpz_prodlimbs (x1, factors, i); |
98 | 899 | size = i + j; |
99 | 899 | prod = MPZ_NEWALLOC (x, size); |
100 | 899 | if (i >= j) |
101 | 420 | cy = mpn_mul (prod, PTR(x1), i, PTR(x2), j); |
102 | 479 | else |
103 | 479 | cy = mpn_mul (prod, PTR(x2), j, PTR(x1), i); |
104 | 899 | TMP_FREE; |
105 | | |
106 | 899 | return SIZ (x) = size - (cy == 0); |
107 | 899 | } |
108 | 1.88k | } |