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
7Copyright 1999, 2000, 2001, 2006 Free Software Foundation, Inc.
8
9This file is part of the GNU MP Library.
10
11The GNU MP Library is free software; you can redistribute it and/or modify
12it under the terms of the GNU Lesser General Public License as published by
13the Free Software Foundation; either version 3 of the License, or (at your
14option) any later version.
15
16The GNU MP Library is distributed in the hope that it will be useful, but
17WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
19License for more details.
20
21You should have received a copy of the GNU Lesser General Public License
22along 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.  */
39mp_limb_t
40mpn_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
134long cputime ();
135
136int
137main (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