1 /* mpn_add_n_sub_n -- Add and Subtract two limb vectors of equal, non-zero length.
2
3 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
4 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
5 GUARANTEED THAT IT'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
6
7 Copyright 1999-2001, 2006 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.h"
36 #include "gmp-impl.h"
37
38 #ifndef L1_CACHE_SIZE
39 #define L1_CACHE_SIZE 8192 /* only 68040 has less than this */
40 #endif
41
42 #define PART_SIZE (L1_CACHE_SIZE / GMP_LIMB_BYTES / 6)
43
44
45 /* mpn_add_n_sub_n.
46 r1[] = s1[] + s2[]
47 r2[] = s1[] - s2[]
48 All operands have n limbs.
49 In-place operations allowed. */
50 mp_limb_t
mpn_add_n_sub_n(mp_ptr r1p,mp_ptr r2p,mp_srcptr s1p,mp_srcptr s2p,mp_size_t n)51 mpn_add_n_sub_n (mp_ptr r1p, mp_ptr r2p, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
52 {
53 mp_limb_t acyn, acyo; /* carry for add */
54 mp_limb_t scyn, scyo; /* carry for subtract */
55 mp_size_t off; /* offset in operands */
56 mp_size_t this_n; /* size of current chunk */
57
58 /* We alternatingly add and subtract in chunks that fit into the (L1)
59 cache. Since the chunks are several hundred limbs, the function call
60 overhead is insignificant, but we get much better locality. */
61
62 /* We have three variant of the inner loop, the proper loop is chosen
63 depending on whether r1 or r2 are the same operand as s1 or s2. */
64
65 if (r1p != s1p && r1p != s2p)
66 {
67 /* r1 is not identical to either input operand. We can therefore write
68 to r1 directly, without using temporary storage. */
69 acyo = 0;
70 scyo = 0;
71 for (off = 0; off < n; off += PART_SIZE)
72 {
73 this_n = MIN (n - off, PART_SIZE);
74 #if HAVE_NATIVE_mpn_add_nc
75 acyo = mpn_add_nc (r1p + off, s1p + off, s2p + off, this_n, acyo);
76 #else
77 acyn = mpn_add_n (r1p + off, s1p + off, s2p + off, this_n);
78 acyo = acyn + mpn_add_1 (r1p + off, r1p + off, this_n, acyo);
79 #endif
80 #if HAVE_NATIVE_mpn_sub_nc
81 scyo = mpn_sub_nc (r2p + off, s1p + off, s2p + off, this_n, scyo);
82 #else
83 scyn = mpn_sub_n (r2p + off, s1p + off, s2p + off, this_n);
84 scyo = scyn + mpn_sub_1 (r2p + off, r2p + off, this_n, scyo);
85 #endif
86 }
87 }
88 else if (r2p != s1p && r2p != s2p)
89 {
90 /* r2 is not identical to either input operand. We can therefore write
91 to r2 directly, without using temporary storage. */
92 acyo = 0;
93 scyo = 0;
94 for (off = 0; off < n; off += PART_SIZE)
95 {
96 this_n = MIN (n - off, PART_SIZE);
97 #if HAVE_NATIVE_mpn_sub_nc
98 scyo = mpn_sub_nc (r2p + off, s1p + off, s2p + off, this_n, scyo);
99 #else
100 scyn = mpn_sub_n (r2p + off, s1p + off, s2p + off, this_n);
101 scyo = scyn + mpn_sub_1 (r2p + off, r2p + off, this_n, scyo);
102 #endif
103 #if HAVE_NATIVE_mpn_add_nc
104 acyo = mpn_add_nc (r1p + off, s1p + off, s2p + off, this_n, acyo);
105 #else
106 acyn = mpn_add_n (r1p + off, s1p + off, s2p + off, this_n);
107 acyo = acyn + mpn_add_1 (r1p + off, r1p + off, this_n, acyo);
108 #endif
109 }
110 }
111 else
112 {
113 /* r1 and r2 are identical to s1 and s2 (r1==s1 and r2==s2 or vice versa)
114 Need temporary storage. */
115 mp_limb_t tp[PART_SIZE];
116 acyo = 0;
117 scyo = 0;
118 for (off = 0; off < n; off += PART_SIZE)
119 {
120 this_n = MIN (n - off, PART_SIZE);
121 #if HAVE_NATIVE_mpn_add_nc
122 acyo = mpn_add_nc (tp, s1p + off, s2p + off, this_n, acyo);
123 #else
124 acyn = mpn_add_n (tp, s1p + off, s2p + off, this_n);
125 acyo = acyn + mpn_add_1 (tp, tp, this_n, acyo);
126 #endif
127 #if HAVE_NATIVE_mpn_sub_nc
128 scyo = mpn_sub_nc (r2p + off, s1p + off, s2p + off, this_n, scyo);
129 #else
130 scyn = mpn_sub_n (r2p + off, s1p + off, s2p + off, this_n);
131 scyo = scyn + mpn_sub_1 (r2p + off, r2p + off, this_n, scyo);
132 #endif
133 MPN_COPY (r1p + off, tp, this_n);
134 }
135 }
136
137 return 2 * acyo + scyo;
138 }
139
140 #ifdef MAIN
141 #include <stdlib.h>
142 #include <stdio.h>
143 #include "timing.h"
144
145 long cputime ();
146
147 int
main(int argc,char ** argv)148 main (int argc, char **argv)
149 {
150 mp_ptr r1p, r2p, s1p, s2p;
151 double t;
152 mp_size_t n;
153
154 n = strtol (argv[1], 0, 0);
155
156 r1p = malloc (n * GMP_LIMB_BYTES);
157 r2p = malloc (n * GMP_LIMB_BYTES);
158 s1p = malloc (n * GMP_LIMB_BYTES);
159 s2p = malloc (n * GMP_LIMB_BYTES);
160 TIME (t,(mpn_add_n(r1p,s1p,s2p,n),mpn_sub_n(r1p,s1p,s2p,n)));
161 printf (" separate add and sub: %.3f\n", t);
162 TIME (t,mpn_add_n_sub_n(r1p,r2p,s1p,s2p,n));
163 printf ("combined addsub separate variables: %.3f\n", t);
164 TIME (t,mpn_add_n_sub_n(r1p,r2p,r1p,s2p,n));
165 printf (" combined addsub r1 overlap: %.3f\n", t);
166 TIME (t,mpn_add_n_sub_n(r1p,r2p,r1p,s2p,n));
167 printf (" combined addsub r2 overlap: %.3f\n", t);
168 TIME (t,mpn_add_n_sub_n(r1p,r2p,r1p,r2p,n));
169 printf (" combined addsub in-place: %.3f\n", t);
170
171 return 0;
172 }
173 #endif
174