mpn/generic/toom3_sqr.c

   1 /* mpn_toom3_sqr -- Square {ap,an}.
   2
   3    Contributed to the GNU project by Torbjorn Granlund.
   4    Additional improvements by Marco Bodrato.
   5
   6    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
   7    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   8    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
   9
  10 Copyright 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
  11
  12 This file is part of the GNU MP Library.
  13
  14 The GNU MP Library is free software; you can redistribute it and/or modify
  15 it under the terms of the GNU Lesser General Public License as published by
  16 the Free Software Foundation; either version 3 of the License, or (at your
  17 option) any later version.
  18
  19 The GNU MP Library is distributed in the hope that it will be useful, but
  20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  21 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  22 License for more details.
  23
  24 You should have received a copy of the GNU Lesser General Public License
  25 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  26
  27
  28 #include "gmp.h"
  29 #include "gmp-impl.h"
  30
  31 /* Evaluate in: -1, 0, +1, +2, +inf
  32
  33   <-s--><--n--><--n-->
  34    ____ ______ ______
  35   |_a2_|___a1_|___a0_|
  36
  37   v0  =  a0         ^2 #   A(0)^2
  38   v1  = (a0+ a1+ a2)^2 #   A(1)^2    ah  <= 2
  39   vm1 = (a0- a1+ a2)^2 #  A(-1)^2   |ah| <= 1
  40   v2  = (a0+2a1+4a2)^2 #   A(2)^2    ah  <= 6
  41   vinf=          a2 ^2 # A(inf)^2
  42 */
  43
  44 #if TUNE_PROGRAM_BUILD
  45 #define MAYBE_sqr_basecase 1
  46 #define MAYBE_sqr_toom3   1
  47 #else
  48 #define MAYBE_sqr_basecase                                              \
  49   (SQR_TOOM3_THRESHOLD < 3 * SQR_TOOM2_THRESHOLD)
  50 #define MAYBE_sqr_toom3                                                 \
  51   (SQR_TOOM4_THRESHOLD >= 3 * SQR_TOOM3_THRESHOLD)
  52 #endif
  53
  54 #define TOOM3_SQR_REC(p, a, n, ws)                                      \
  55   do {                                                                  \
  56     if (MAYBE_sqr_basecase                                              \
  57         && BELOW_THRESHOLD (n, SQR_TOOM2_THRESHOLD))                    \
  58       mpn_sqr_basecase (p, a, n);                                       \
  59     else if (! MAYBE_sqr_toom3                                          \
  60              || BELOW_THRESHOLD (n, SQR_TOOM3_THRESHOLD))               \
  61       mpn_toom2_sqr (p, a, n, ws);                                      \
  62     else                                                                \
  63       mpn_toom3_sqr (p, a, n, ws);                                      \
  64   } while (0)
  65
  66 void
  67 mpn_toom3_sqr (mp_ptr pp,
  68                mp_srcptr ap, mp_size_t an,
  69                mp_ptr scratch)
  70 {
  71   mp_size_t n, s;
  72   mp_limb_t cy, vinf0;
  73   mp_ptr gp;
  74   mp_ptr as1, asm1, as2;
  75
  76 #define a0  ap
  77 #define a1  (ap + n)
  78 #define a2  (ap + 2*n)
  79
  80   n = (an + 2) / (size_t) 3;
  81
  82   s = an - 2 * n;
  83
  84   ASSERT (0 < s && s <= n);
  85
  86   as1 = scratch + 4 * n + 4;
  87   asm1 = scratch + 2 * n + 2;
  88   as2 = pp + n + 1;
  89
  90   gp = scratch;
  91
  92   /* Compute as1 and asm1.  */
  93   cy = mpn_add (gp, a0, n, a2, s);
  94 #if HAVE_NATIVE_mpn_add_n_sub_n
  95   if (cy == 0 && mpn_cmp (gp, a1, n) < 0)
  96     {
  97       cy = mpn_add_n_sub_n (as1, asm1, a1, gp, n);
  98       as1[n] = 0;
  99       asm1[n] = 0;
 100     }
 101   else
 102     {
 103       cy2 = mpn_add_n_sub_n (as1, asm1, gp, a1, n);
 104       as1[n] = cy + (cy2 >> 1);
 105       asm1[n] = cy - (cy & 1);
 106     }
 107 #else
 108   as1[n] = cy + mpn_add_n (as1, gp, a1, n);
 109   if (cy == 0 && mpn_cmp (gp, a1, n) < 0)
 110     {
 111       mpn_sub_n (asm1, a1, gp, n);
 112       asm1[n] = 0;
 113     }
 114   else
 115     {
 116       cy -= mpn_sub_n (asm1, gp, a1, n);
 117       asm1[n] = cy;
 118     }
 119 #endif
 120
 121   /* Compute as2.  */
 122 #if HAVE_NATIVE_mpn_rsblsh1_n
 123   cy = mpn_add_n (as2, a2, as1, s);
 124   if (s != n)
 125     cy = mpn_add_1 (as2 + s, as1 + s, n - s, cy);
 126   cy += as1[n];
 127   cy = 2 * cy + mpn_rsblsh1_n (as2, a0, as2, n);
 128 #else
 129 #if HAVE_NATIVE_mpn_addlsh1_n
 130   cy  = mpn_addlsh1_n (as2, a1, a2, s);
 131   if (s != n)
 132     cy = mpn_add_1 (as2 + s, a1 + s, n - s, cy);
 133   cy = 2 * cy + mpn_addlsh1_n (as2, a0, as2, n);
 134 #else
 135   cy = mpn_add_n (as2, a2, as1, s);
 136   if (s != n)
 137     cy = mpn_add_1 (as2 + s, as1 + s, n - s, cy);
 138   cy += as1[n];
 139   cy = 2 * cy + mpn_lshift (as2, as2, n, 1);
 140   cy -= mpn_sub_n (as2, as2, a0, n);
 141 #endif
 142 #endif
 143   as2[n] = cy;
 144
 145   ASSERT (as1[n] <= 2);
 146   ASSERT (asm1[n] <= 1);
 147
 148 #define v0    pp                                /* 2n */
 149 #define v1    (pp + 2 * n)                      /* 2n+1 */
 150 #define vinf  (pp + 4 * n)                      /* s+s */
 151 #define vm1   scratch                           /* 2n+1 */
 152 #define v2    (scratch + 2 * n + 1)             /* 2n+2 */
 153 #define scratch_out  (scratch + 5 * n + 5)
 154
 155   /* vm1, 2n+1 limbs */
 156 #ifdef SMALLER_RECURSION
 157   TOOM3_SQR_REC (vm1, asm1, n, scratch_out);
 158   cy = 0;
 159   if (asm1[n] != 0)
 160     cy = asm1[n] + mpn_add_n (vm1 + n, vm1 + n, asm1, n);
 161   if (asm1[n] != 0)
 162     cy += mpn_add_n (vm1 + n, vm1 + n, asm1, n);
 163   vm1[2 * n] = cy;
 164 #else
 165   TOOM3_SQR_REC (vm1, asm1, n + 1, scratch_out);
 166 #endif
 167
 168   TOOM3_SQR_REC (v2, as2, n + 1, scratch_out);  /* v2, 2n+1 limbs */
 169
 170   TOOM3_SQR_REC (vinf, a2, s, scratch_out);     /* vinf, s+s limbs */
 171
 172   vinf0 = vinf[0];                              /* v1 overlaps with this */
 173
 174 #ifdef SMALLER_RECURSION
 175   /* v1, 2n+1 limbs */
 176   TOOM3_SQR_REC (v1, as1, n, scratch_out);
 177   if (as1[n] == 1)
 178     {
 179       cy = as1[n] + mpn_add_n (v1 + n, v1 + n, as1, n);
 180     }
 181   else if (as1[n] != 0)
 182     {
 183 #if HAVE_NATIVE_mpn_addlsh1_n
 184       cy = 2 * as1[n] + mpn_addlsh1_n (v1 + n, v1 + n, as1, n);
 185 #else
 186       cy = 2 * as1[n] + mpn_addmul_1 (v1 + n, as1, n, CNST_LIMB(2));
 187 #endif
 188     }
 189   else
 190     cy = 0;
 191   if (as1[n] == 1)
 192     {
 193       cy += mpn_add_n (v1 + n, v1 + n, as1, n);
 194     }
 195   else if (as1[n] != 0)
 196     {
 197 #if HAVE_NATIVE_mpn_addlsh1_n
 198       cy += mpn_addlsh1_n (v1 + n, v1 + n, as1, n);
 199 #else
 200       cy += mpn_addmul_1 (v1 + n, as1, n, CNST_LIMB(2));
 201 #endif
 202     }
 203   v1[2 * n] = cy;
 204 #else
 205   cy = vinf[1];
 206   TOOM3_SQR_REC (v1, as1, n + 1, scratch_out);
 207   vinf[1] = cy;
 208 #endif
 209
 210   TOOM3_SQR_REC (v0, ap, n, scratch_out);       /* v0, 2n limbs */
 211
 212   mpn_toom_interpolate_5pts (pp, v2, vm1, n, s + s, 0, vinf0);
 213 }