/src/gmp/mpn/mul_basecase.c
Line | Count | Source (jump to first uncovered line) |
1 | | /* mpn_mul_basecase -- Internal routine to multiply two natural numbers |
2 | | of length m and n. |
3 | | |
4 | | THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE. IT IS ONLY |
5 | | SAFE TO REACH THIS FUNCTION THROUGH DOCUMENTED INTERFACES. |
6 | | |
7 | | Copyright 1991-1994, 1996, 1997, 2000-2002 Free Software Foundation, Inc. |
8 | | |
9 | | This file is part of the GNU MP Library. |
10 | | |
11 | | The GNU MP Library is free software; you can redistribute it and/or modify |
12 | | it under the terms of either: |
13 | | |
14 | | * the GNU Lesser General Public License as published by the Free |
15 | | Software Foundation; either version 3 of the License, or (at your |
16 | | option) any later version. |
17 | | |
18 | | or |
19 | | |
20 | | * the GNU General Public License as published by the Free Software |
21 | | Foundation; either version 2 of the License, or (at your option) any |
22 | | later version. |
23 | | |
24 | | or both in parallel, as here. |
25 | | |
26 | | The GNU MP Library is distributed in the hope that it will be useful, but |
27 | | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
28 | | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
29 | | for more details. |
30 | | |
31 | | You should have received copies of the GNU General Public License and the |
32 | | GNU Lesser General Public License along with the GNU MP Library. If not, |
33 | | see https://www.gnu.org/licenses/. */ |
34 | | |
35 | | #include "gmp-impl.h" |
36 | | |
37 | | |
38 | | /* Multiply {up,usize} by {vp,vsize} and write the result to |
39 | | {prodp,usize+vsize}. Must have usize>=vsize. |
40 | | |
41 | | Note that prodp gets usize+vsize limbs stored, even if the actual result |
42 | | only needs usize+vsize-1. |
43 | | |
44 | | There's no good reason to call here with vsize>=MUL_TOOM22_THRESHOLD. |
45 | | Currently this is allowed, but it might not be in the future. |
46 | | |
47 | | This is the most critical code for 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. */ |
51 | | |
52 | | void |
53 | | mpn_mul_basecase (mp_ptr rp, |
54 | | mp_srcptr up, mp_size_t un, |
55 | | mp_srcptr vp, mp_size_t vn) |
56 | 231M | { |
57 | 231M | ASSERT (un >= vn); |
58 | 231M | ASSERT (vn >= 1); |
59 | 231M | ASSERT (! MPN_OVERLAP_P (rp, un+vn, up, un)); |
60 | 231M | ASSERT (! MPN_OVERLAP_P (rp, un+vn, vp, vn)); |
61 | | |
62 | | /* We first multiply by the low order limb (or depending on optional function |
63 | | availability, limbs). This result can be stored, not added, to rp. We |
64 | | also avoid a loop for zeroing this way. */ |
65 | | |
66 | | #if HAVE_NATIVE_mpn_mul_2 |
67 | | if (vn >= 2) |
68 | | { |
69 | | rp[un + 1] = mpn_mul_2 (rp, up, un, vp); |
70 | | rp += 2, vp += 2, vn -= 2; |
71 | | } |
72 | | else |
73 | | { |
74 | | rp[un] = mpn_mul_1 (rp, up, un, vp[0]); |
75 | | return; |
76 | | } |
77 | | #else |
78 | 231M | rp[un] = mpn_mul_1 (rp, up, un, vp[0]); |
79 | 231M | rp += 1, vp += 1, vn -= 1; |
80 | 231M | #endif |
81 | | |
82 | | /* Now accumulate the product of up[] and the next higher limb (or depending |
83 | | on optional function availability, limbs) from vp[]. */ |
84 | | |
85 | 2.74G | #define MAX_LEFT MP_SIZE_T_MAX /* Used to simplify loops into if statements */ |
86 | | |
87 | | |
88 | | #if HAVE_NATIVE_mpn_addmul_6 |
89 | | while (vn >= 6) |
90 | | { |
91 | | rp[un + 6 - 1] = mpn_addmul_6 (rp, up, un, vp); |
92 | | if (MAX_LEFT == 6) |
93 | | return; |
94 | | rp += 6, vp += 6, vn -= 6; |
95 | | if (MAX_LEFT < 2 * 6) |
96 | | break; |
97 | | } |
98 | | #undef MAX_LEFT |
99 | | #define MAX_LEFT (6 - 1) |
100 | | #endif |
101 | | |
102 | | #if HAVE_NATIVE_mpn_addmul_5 |
103 | | while (vn >= 5) |
104 | | { |
105 | | rp[un + 5 - 1] = mpn_addmul_5 (rp, up, un, vp); |
106 | | if (MAX_LEFT == 5) |
107 | | return; |
108 | | rp += 5, vp += 5, vn -= 5; |
109 | | if (MAX_LEFT < 2 * 5) |
110 | | break; |
111 | | } |
112 | | #undef MAX_LEFT |
113 | | #define MAX_LEFT (5 - 1) |
114 | | #endif |
115 | | |
116 | | #if HAVE_NATIVE_mpn_addmul_4 |
117 | | while (vn >= 4) |
118 | | { |
119 | | rp[un + 4 - 1] = mpn_addmul_4 (rp, up, un, vp); |
120 | | if (MAX_LEFT == 4) |
121 | | return; |
122 | | rp += 4, vp += 4, vn -= 4; |
123 | | if (MAX_LEFT < 2 * 4) |
124 | | break; |
125 | | } |
126 | | #undef MAX_LEFT |
127 | | #define MAX_LEFT (4 - 1) |
128 | | #endif |
129 | | |
130 | | #if HAVE_NATIVE_mpn_addmul_3 |
131 | | while (vn >= 3) |
132 | | { |
133 | | rp[un + 3 - 1] = mpn_addmul_3 (rp, up, un, vp); |
134 | | if (MAX_LEFT == 3) |
135 | | return; |
136 | | rp += 3, vp += 3, vn -= 3; |
137 | | if (MAX_LEFT < 2 * 3) |
138 | | break; |
139 | | } |
140 | | #undef MAX_LEFT |
141 | | #define MAX_LEFT (3 - 1) |
142 | | #endif |
143 | | |
144 | | #if HAVE_NATIVE_mpn_addmul_2 |
145 | | while (vn >= 2) |
146 | | { |
147 | | rp[un + 2 - 1] = mpn_addmul_2 (rp, up, un, vp); |
148 | | if (MAX_LEFT == 2) |
149 | | return; |
150 | | rp += 2, vp += 2, vn -= 2; |
151 | | if (MAX_LEFT < 2 * 2) |
152 | | break; |
153 | | } |
154 | | #undef MAX_LEFT |
155 | | #define MAX_LEFT (2 - 1) |
156 | | #endif |
157 | | |
158 | 2.97G | while (vn >= 1) |
159 | 2.74G | { |
160 | 2.74G | rp[un] = mpn_addmul_1 (rp, up, un, vp[0]); |
161 | 2.74G | if (MAX_LEFT == 1) |
162 | 0 | return; |
163 | 2.74G | rp += 1, vp += 1, vn -= 1; |
164 | 2.74G | } |
165 | 231M | } |