← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/mpn/generic/jacbase.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2006-05-17 02:46:26 UTC
  • Revision ID: james.westby@ubuntu.com-20060517024626-lljr08ftv9g9vefl
Tags: upstream-0.9h-20060510
Import upstream version 0.9h-20060510

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpn/generic/jacbase.c
 
1
/* mpn_jacobi_base -- limb/limb Jacobi symbol with restricted arguments.
 
2
 
 
3
   THIS INTERFACE IS PRELIMINARY AND MIGHT DISAPPEAR OR BE SUBJECT TO
 
4
   INCOMPATIBLE CHANGES IN A FUTURE RELEASE OF GMP.
 
5
 
 
6
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
 
7
 
 
8
This file is part of the GNU MP Library.
 
9
 
 
10
The GNU MP Library is free software; you can redistribute it and/or modify
 
11
it under the terms of the GNU Lesser General Public License as published by
 
12
the Free Software Foundation; either version 2.1 of the License, or (at your
 
13
option) any later version.
 
14
 
 
15
The GNU MP Library is distributed in the hope that it will be useful, but
 
16
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
17
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
18
License for more details.
 
19
 
 
20
You should have received a copy of the GNU Lesser General Public License
 
21
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 
22
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
23
MA 02111-1307, USA. */
 
24
 
 
25
#include "gmp.h"
 
26
#include "gmp-impl.h"
 
27
#include "longlong.h"
 
28
 
 
29
 
 
30
/* Use the simple loop by default.  The generic count_trailing_zeros is not
 
31
   very fast, and the extra trickery of method 3 has proven to be less use
 
32
   than might have been though.  */
 
33
#ifndef JACOBI_BASE_METHOD
 
34
#define JACOBI_BASE_METHOD  2
 
35
#endif
 
36
 
 
37
 
 
38
/* Use count_trailing_zeros.  */
 
39
#if JACOBI_BASE_METHOD == 1
 
40
#define PROCESS_TWOS_ANY                                \
 
41
  {                                                     \
 
42
    mp_limb_t  twos;                                    \
 
43
    count_trailing_zeros (twos, a);                     \
 
44
    result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b);        \
 
45
    a >>= twos;                                         \
 
46
  }
 
47
#define PROCESS_TWOS_EVEN  PROCESS_TWOS_ANY
 
48
#endif
 
49
 
 
50
/* Use a simple loop.  A disadvantage of this is that there's a branch on a
 
51
   50/50 chance of a 0 or 1 low bit.  */
 
52
#if JACOBI_BASE_METHOD == 2
 
53
#define PROCESS_TWOS_EVEN               \
 
54
  {                                     \
 
55
    int  two;                           \
 
56
    two = JACOBI_TWO_U_BIT1 (b);        \
 
57
    do                                  \
 
58
      {                                 \
 
59
        a >>= 1;                        \
 
60
        result_bit1 ^= two;             \
 
61
        ASSERT (a != 0);                \
 
62
      }                                 \
 
63
    while ((a & 1) == 0);               \
 
64
  }
 
65
#define PROCESS_TWOS_ANY        \
 
66
  if ((a & 1) == 0)             \
 
67
    PROCESS_TWOS_EVEN;
 
68
#endif
 
69
 
 
70
/* Process one bit arithmetically, then a simple loop.  This cuts the loop
 
71
   condition down to a 25/75 chance, which should branch predict better.
 
72
   The CPU will need a reasonable variable left shift.  */
 
73
#if JACOBI_BASE_METHOD == 3
 
74
#define PROCESS_TWOS_EVEN               \
 
75
  {                                     \
 
76
    int  two, mask, shift;              \
 
77
                                        \
 
78
    two = JACOBI_TWO_U_BIT1 (b);        \
 
79
    mask = (~a & 2);                    \
 
80
    a >>= 1;                            \
 
81
                                        \
 
82
    shift = (~a & 1);                   \
 
83
    a >>= shift;                        \
 
84
    result_bit1 ^= two ^ (two & mask);  \
 
85
                                        \
 
86
    while ((a & 1) == 0)                \
 
87
      {                                 \
 
88
        a >>= 1;                        \
 
89
        result_bit1 ^= two;             \
 
90
        ASSERT (a != 0);                \
 
91
      }                                 \
 
92
  }
 
93
#define PROCESS_TWOS_ANY                \
 
94
  {                                     \
 
95
    int  two, mask, shift;              \
 
96
                                        \
 
97
    two = JACOBI_TWO_U_BIT1 (b);        \
 
98
    shift = (~a & 1);                   \
 
99
    a >>= shift;                        \
 
100
                                        \
 
101
    mask = shift << 1;                  \
 
102
    result_bit1 ^= (two & mask);        \
 
103
                                        \
 
104
    while ((a & 1) == 0)                \
 
105
      {                                 \
 
106
        a >>= 1;                        \
 
107
        result_bit1 ^= two;             \
 
108
        ASSERT (a != 0);                \
 
109
      }                                 \
 
110
  }
 
111
#endif
 
112
 
 
113
 
 
114
/* Calculate the value of the Jacobi symbol (a/b) of two mp_limb_t's, but
 
115
   with a restricted range of inputs accepted, namely b>1, b odd, and a<=b.
 
116
 
 
117
   The initial result_bit1 is taken as a parameter for the convenience of
 
118
   mpz_kronecker_ui() et al.  The sign changes both here and in those
 
119
   routines accumulate nicely in bit 1, see the JACOBI macros.
 
120
 
 
121
   The return value here is the normal +1, 0, or -1.  Note that +1 and -1
 
122
   have bit 1 in the "BIT1" sense, which could be useful if the caller is
 
123
   accumulating it into some extended calculation.
 
124
 
 
125
   Duplicating the loop body to avoid the MP_LIMB_T_SWAP(a,b) would be
 
126
   possible, but a couple of tests suggest it's not a significant speedup,
 
127
   and may even be a slowdown, so what's here is good enough for now.
 
128
 
 
129
   Future: The code doesn't demand a<=b actually, so maybe this could be
 
130
   relaxed.  All the places this is used currently call with a<=b though.  */
 
131
 
 
132
int
 
133
mpn_jacobi_base (mp_limb_t a, mp_limb_t b, int result_bit1)
 
134
{
 
135
  ASSERT (b & 1);  /* b odd */
 
136
  ASSERT (b != 1);
 
137
  ASSERT (a <= b);
 
138
 
 
139
  if (a == 0)
 
140
    return 0;
 
141
 
 
142
  PROCESS_TWOS_ANY;
 
143
  if (a == 1)
 
144
    goto done;
 
145
 
 
146
  for (;;)
 
147
    {
 
148
      result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a, b);
 
149
      MP_LIMB_T_SWAP (a, b);
 
150
 
 
151
      do
 
152
        {
 
153
          /* working on (a/b), a,b odd, a>=b */
 
154
          ASSERT (a & 1);
 
155
          ASSERT (b & 1);
 
156
          ASSERT (a >= b);
 
157
 
 
158
          if ((a -= b) == 0)
 
159
            return 0;
 
160
 
 
161
          PROCESS_TWOS_EVEN;
 
162
          if (a == 1)
 
163
            goto done;
 
164
        }
 
165
      while (a >= b);
 
166
    }
 
167
 
 
168
 done:
 
169
  return JACOBI_BIT1_TO_PN (result_bit1);
 
170
}