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, 2000, 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 the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15
16 The GNU MP Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 License for more details.
20
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
23
24 #include "gmp.h"
25 #include "gmp-impl.h"
26
27 #ifndef L1_CACHE_SIZE
28 #define L1_CACHE_SIZE 8192 /* only 68040 has less than this */
29 #endif
30
31 #define PART_SIZE (L1_CACHE_SIZE / BYTES_PER_MP_LIMB / 6)
32
33
34 /* mpn_add_n_sub_n.
35 r1[] = s1[] + s2[]
36 r2[] = s1[] - s2[]
37 All operands have n limbs.
38 In-place operations allowed. */
39 mp_limb_t
mpn_add_n_sub_n(mp_ptr r1p,mp_ptr r2p,mp_srcptr s1p,mp_srcptr s2p,mp_size_t n)40 mpn_add_n_sub_n (mp_ptr r1p, mp_ptr r2p, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
41 {
42 mp_limb_t acyn, acyo; /* carry for add */
43 mp_limb_t scyn, scyo; /* carry for subtract */
44 mp_size_t off; /* offset in operands */
45 mp_size_t this_n; /* size of current chunk */
46
47 /* We alternatingly add and subtract in chunks that fit into the (L1)
48 cache. Since the chunks are several hundred limbs, the function call
49 overhead is insignificant, but we get much better locality. */
50
51 /* We have three variant of the inner loop, the proper loop is chosen
52 depending on whether r1 or r2 are the same operand as s1 or s2. */
53
54 if (r1p != s1p && r1p != s2p)
55 {
56 /* r1 is not identical to either input operand. We can therefore write
57 to r1 directly, without using temporary storage. */
58 acyo = 0;
59 scyo = 0;
60 for (off = 0; off < n; off += PART_SIZE)
61 {
62 this_n = MIN (n - off, PART_SIZE);
63 #if HAVE_NATIVE_mpn_add_nc
64 acyo = mpn_add_nc (r1p + off, s1p + off, s2p + off, this_n, acyo);
65 #else
66 acyn = mpn_add_n (r1p + off, s1p + off, s2p + off, this_n);
67 acyo = acyn + mpn_add_1 (r1p + off, r1p + off, this_n, acyo);
68 #endif
69 #if HAVE_NATIVE_mpn_sub_nc
70 scyo = mpn_sub_nc (r2p + off, s1p + off, s2p + off, this_n, scyo);
71 #else
72 scyn = mpn_sub_n (r2p + off, s1p + off, s2p + off, this_n);
73 scyo = scyn + mpn_sub_1 (r2p + off, r2p + off, this_n, scyo);
74 #endif
75 }
76 }
77 else if (r2p != s1p && r2p != s2p)
78 {
79 /* r2 is not identical to either input operand. We can therefore write
80 to r2 directly, without using temporary storage. */
81 acyo = 0;
82 scyo = 0;
83 for (off = 0; off < n; off += PART_SIZE)
84 {
85 this_n = MIN (n - off, PART_SIZE);
86 #if HAVE_NATIVE_mpn_sub_nc
87 scyo = mpn_sub_nc (r2p + off, s1p + off, s2p + off, this_n, scyo);
88 #else
89 scyn = mpn_sub_n (r2p + off, s1p + off, s2p + off, this_n);
90 scyo = scyn + mpn_sub_1 (r2p + off, r2p + off, this_n, scyo);
91 #endif
92 #if HAVE_NATIVE_mpn_add_nc
93 acyo = mpn_add_nc (r1p + off, s1p + off, s2p + off, this_n, acyo);
94 #else
95 acyn = mpn_add_n (r1p + off, s1p + off, s2p + off, this_n);
96 acyo = acyn + mpn_add_1 (r1p + off, r1p + off, this_n, acyo);
97 #endif
98 }
99 }
100 else
101 {
102 /* r1 and r2 are identical to s1 and s2 (r1==s1 and r2==s2 or vice versa)
103 Need temporary storage. */
104 mp_limb_t tp[PART_SIZE];
105 acyo = 0;
106 scyo = 0;
107 for (off = 0; off < n; off += PART_SIZE)
108 {
109 this_n = MIN (n - off, PART_SIZE);
110 #if HAVE_NATIVE_mpn_add_nc
111 acyo = mpn_add_nc (tp, s1p + off, s2p + off, this_n, acyo);
112 #else
113 acyn = mpn_add_n (tp, s1p + off, s2p + off, this_n);
114 acyo = acyn + mpn_add_1 (tp, tp, this_n, acyo);
115 #endif
116 #if HAVE_NATIVE_mpn_sub_nc
117 scyo = mpn_sub_nc (r2p + off, s1p + off, s2p + off, this_n, scyo);
118 #else
119 scyn = mpn_sub_n (r2p + off, s1p + off, s2p + off, this_n);
120 scyo = scyn + mpn_sub_1 (r2p + off, r2p + off, this_n, scyo);
121 #endif
122 MPN_COPY (r1p + off, tp, this_n);
123 }
124 }
125
126 return 2 * acyo + scyo;
127 }
128
129 #ifdef MAIN
130 #include <stdlib.h>
131 #include <stdio.h>
132 #include "timing.h"
133
134 long cputime ();
135
136 int
main(int argc,char ** argv)137 main (int argc, char **argv)
138 {
139 mp_ptr r1p, r2p, s1p, s2p;
140 double t;
141 mp_size_t n;
142
143 n = strtol (argv[1], 0, 0);
144
145 r1p = malloc (n * BYTES_PER_MP_LIMB);
146 r2p = malloc (n * BYTES_PER_MP_LIMB);
147 s1p = malloc (n * BYTES_PER_MP_LIMB);
148 s2p = malloc (n * BYTES_PER_MP_LIMB);
149 TIME (t,(mpn_add_n(r1p,s1p,s2p,n),mpn_sub_n(r1p,s1p,s2p,n)));
150 printf (" separate add and sub: %.3f\n", t);
151 TIME (t,mpn_add_n_sub_n(r1p,r2p,s1p,s2p,n));
152 printf ("combined addsub separate variables: %.3f\n", t);
153 TIME (t,mpn_add_n_sub_n(r1p,r2p,r1p,s2p,n));
154 printf (" combined addsub r1 overlap: %.3f\n", t);
155 TIME (t,mpn_add_n_sub_n(r1p,r2p,r1p,s2p,n));
156 printf (" combined addsub r2 overlap: %.3f\n", t);
157 TIME (t,mpn_add_n_sub_n(r1p,r2p,r1p,r2p,n));
158 printf (" combined addsub in-place: %.3f\n", t);
159
160 return 0;
161 }
162 #endif
163