source/libs/gmp/gmp-6.1.0/mpn/generic/jacobi_2.c

   1 /* jacobi_2.c
   2
   3    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
   4    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   5    GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
   6
   7 Copyright 1996, 1998, 2000-2004, 2008, 2010 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 #include "longlong.h"
  38
  39 #ifndef JACOBI_2_METHOD
  40 #define JACOBI_2_METHOD 2
  41 #endif
  42
  43 /* Computes (a / b) where b is odd, and a and b are otherwise arbitrary
  44    two-limb numbers. */
  45 #if JACOBI_2_METHOD == 1
  46 int
  47 mpn_jacobi_2 (mp_srcptr ap, mp_srcptr bp, unsigned bit)
  48 {
  49   mp_limb_t ah, al, bh, bl;
  50   int c;
  51
  52   al = ap[0];
  53   ah = ap[1];
  54   bl = bp[0];
  55   bh = bp[1];
  56
  57   ASSERT (bl & 1);
  58
  59   bl = ((bh << (GMP_NUMB_BITS - 1)) & GMP_NUMB_MASK) | (bl >> 1);
  60   bh >>= 1;
  61
  62   if ( (bh | bl) == 0)
  63     return 1 - 2*(bit & 1);
  64
  65   if ( (ah | al) == 0)
  66     return 0;
  67
  68   if (al == 0)
  69     {
  70       al = ah;
  71       ah = 0;
  72       bit ^= GMP_NUMB_BITS & (bl ^ (bl >> 1));
  73     }
  74   count_trailing_zeros (c, al);
  75   bit ^= c & (bl ^ (bl >> 1));
  76
  77   c++;
  78   if (UNLIKELY (c == GMP_NUMB_BITS))
  79     {
  80       al = ah;
  81       ah = 0;
  82     }
  83   else
  84     {
  85       al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
  86       ah >>= c;
  87     }
  88   while ( (ah | bh) > 0)
  89     {
  90       mp_limb_t th, tl;
  91       mp_limb_t bgta;
  92
  93       sub_ddmmss (th, tl, ah, al, bh, bl);
  94       if ( (tl | th) == 0)
  95         return 0;
  96
  97       bgta = LIMB_HIGHBIT_TO_MASK (th);
  98
  99       /* If b > a, invoke reciprocity */
 100       bit ^= (bgta & al & bl);
 101
 102       /* b <-- min (a, b) */
 103       add_ssaaaa (bh, bl, bh, bl, th & bgta, tl & bgta);
 104
 105       if ( (bh | bl) == 0)
 106         return 1 - 2*(bit & 1);
 107
 108       /* a <-- |a - b| */
 109       al = (bgta ^ tl) - bgta;
 110       ah = (bgta ^ th);
 111
 112       if (UNLIKELY (al == 0))
 113         {
 114           /* If b > a, al == 0 implies that we have a carry to
 115              propagate. */
 116           al = ah - bgta;
 117           ah = 0;
 118           bit ^= GMP_NUMB_BITS & (bl ^ (bl >> 1));
 119         }
 120       count_trailing_zeros (c, al);
 121       c++;
 122       bit ^= c & (bl ^ (bl >> 1));
 123
 124       if (UNLIKELY (c == GMP_NUMB_BITS))
 125         {
 126           al = ah;
 127           ah = 0;
 128         }
 129       else
 130         {
 131           al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
 132           ah >>= c;
 133         }
 134     }
 135
 136   ASSERT (bl > 0);
 137
 138   while ( (al | bl) & GMP_LIMB_HIGHBIT)
 139     {
 140       /* Need an extra comparison to get the mask. */
 141       mp_limb_t t = al - bl;
 142       mp_limb_t bgta = - (bl > al);
 143
 144       if (t == 0)
 145         return 0;
 146
 147       /* If b > a, invoke reciprocity */
 148       bit ^= (bgta & al & bl);
 149
 150       /* b <-- min (a, b) */
 151       bl += (bgta & t);
 152
 153       /* a <-- |a - b| */
 154       al = (t ^ bgta) - bgta;
 155
 156       /* Number of trailing zeros is the same no matter if we look at
 157        * t or a, but using t gives more parallelism. */
 158       count_trailing_zeros (c, t);
 159       c ++;
 160       /* (2/b) = -1 if b = 3 or 5 mod 8 */
 161       bit ^= c & (bl ^ (bl >> 1));
 162
 163       if (UNLIKELY (c == GMP_NUMB_BITS))
 164         return 1 - 2*(bit & 1);
 165
 166       al >>= c;
 167     }
 168
 169   /* Here we have a little impedance mismatch. Better to inline it? */
 170   return mpn_jacobi_base (2*al+1, 2*bl+1, bit << 1);
 171 }
 172 #elif JACOBI_2_METHOD == 2
 173 int
 174 mpn_jacobi_2 (mp_srcptr ap, mp_srcptr bp, unsigned bit)
 175 {
 176   mp_limb_t ah, al, bh, bl;
 177   int c;
 178
 179   al = ap[0];
 180   ah = ap[1];
 181   bl = bp[0];
 182   bh = bp[1];
 183
 184   ASSERT (bl & 1);
 185
 186   /* Use bit 1. */
 187   bit <<= 1;
 188
 189   if (bh == 0 && bl == 1)
 190     /* (a|1) = 1 */
 191     return 1 - (bit & 2);
 192
 193   if (al == 0)
 194     {
 195       if (ah == 0)
 196         /* (0|b) = 0, b > 1 */
 197         return 0;
 198
 199       count_trailing_zeros (c, ah);
 200       bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1));
 201
 202       al = bl;
 203       bl = ah >> c;
 204
 205       if (bl == 1)
 206         /* (1|b) = 1 */
 207         return 1 - (bit & 2);
 208
 209       ah = bh;
 210
 211       bit ^= al & bl;
 212
 213       goto b_reduced;
 214     }
 215   if ( (al & 1) == 0)
 216     {
 217       count_trailing_zeros (c, al);
 218
 219       al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
 220       ah >>= c;
 221       bit ^= (c << 1) & (bl ^ (bl >> 1));
 222     }
 223   if (ah == 0)
 224     {
 225       if (bh > 0)
 226         {
 227           bit ^= al & bl;
 228           MP_LIMB_T_SWAP (al, bl);
 229           ah = bh;
 230           goto b_reduced;
 231         }
 232       goto ab_reduced;
 233     }
 234
 235   while (bh > 0)
 236     {
 237       /* Compute (a|b) */
 238       while (ah > bh)
 239         {
 240           sub_ddmmss (ah, al, ah, al, bh, bl);
 241           if (al == 0)
 242             {
 243               count_trailing_zeros (c, ah);
 244               bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1));
 245
 246               al = bl;
 247               bl = ah >> c;
 248               ah = bh;
 249
 250               bit ^= al & bl;
 251               goto b_reduced;
 252             }
 253           count_trailing_zeros (c, al);
 254           bit ^= (c << 1) & (bl ^ (bl >> 1));
 255           al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
 256           ah >>= c;
 257         }
 258       if (ah == bh)
 259         goto cancel_hi;
 260
 261       if (ah == 0)
 262         {
 263           bit ^= al & bl;
 264           MP_LIMB_T_SWAP (al, bl);
 265           ah = bh;
 266           break;
 267         }
 268
 269       bit ^= al & bl;
 270
 271       /* Compute (b|a) */
 272       while (bh > ah)
 273         {
 274           sub_ddmmss (bh, bl, bh, bl, ah, al);
 275           if (bl == 0)
 276             {
 277               count_trailing_zeros (c, bh);
 278               bit ^= ((GMP_NUMB_BITS + c) << 1) & (al ^ (al >> 1));
 279
 280               bl = bh >> c;
 281               bit ^= al & bl;
 282               goto b_reduced;
 283             }
 284           count_trailing_zeros (c, bl);
 285           bit ^= (c << 1) & (al ^ (al >> 1));
 286           bl = ((bh << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (bl >> c);
 287           bh >>= c;
 288         }
 289       bit ^= al & bl;
 290
 291       /* Compute (a|b) */
 292       if (ah == bh)
 293         {
 294         cancel_hi:
 295           if (al < bl)
 296             {
 297               MP_LIMB_T_SWAP (al, bl);
 298               bit ^= al & bl;
 299             }
 300           al -= bl;
 301           if (al == 0)
 302             return 0;
 303
 304           count_trailing_zeros (c, al);
 305           bit ^= (c << 1) & (bl ^ (bl >> 1));
 306           al >>= c;
 307
 308           if (al == 1)
 309             return 1 - (bit & 2);
 310
 311           MP_LIMB_T_SWAP (al, bl);
 312           bit ^= al & bl;
 313           break;
 314         }
 315     }
 316
 317  b_reduced:
 318   /* Compute (a|b), with b a single limb. */
 319   ASSERT (bl & 1);
 320
 321   if (bl == 1)
 322     /* (a|1) = 1 */
 323     return 1 - (bit & 2);
 324
 325   while (ah > 0)
 326     {
 327       ah -= (al < bl);
 328       al -= bl;
 329       if (al == 0)
 330         {
 331           if (ah == 0)
 332             return 0;
 333           count_trailing_zeros (c, ah);
 334           bit ^= ((GMP_NUMB_BITS + c) << 1) & (bl ^ (bl >> 1));
 335           al = ah >> c;
 336           goto ab_reduced;
 337         }
 338       count_trailing_zeros (c, al);
 339
 340       al = ((ah << (GMP_NUMB_BITS - c)) & GMP_NUMB_MASK) | (al >> c);
 341       ah >>= c;
 342       bit ^= (c << 1) & (bl ^ (bl >> 1));
 343     }
 344  ab_reduced:
 345   ASSERT (bl & 1);
 346   ASSERT (bl > 1);
 347
 348   return mpn_jacobi_base (al, bl, bit);
 349 }
 350 #else
 351 #error Unsupported value for JACOBI_2_METHOD
 352 #endif