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