← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gclcvs/quantal

« back to all changes in this revision

Viewing changes to gmp3/tests/refmpz.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-06-24 15:13:46 UTC
  • Revision ID: james.westby@ubuntu.com-20040624151346-xh0xaaktyyp7aorc
Tags: 2.7.0-26
C_GC_OFFSET is 2 on m68k-linux

Show diffs side-by-side

added added

removed removed

Lines of Context:
gmp3/tests/refmpz.c
 
1
/* Reference mpz functions.
 
2
 
 
3
Copyright 1997, 1999, 2000, 2001 Free Software Foundation, Inc.
 
4
 
 
5
This file is part of the GNU MP Library.
 
6
 
 
7
The GNU MP Library is free software; you can redistribute it and/or modify
 
8
it under the terms of the GNU Lesser General Public License as published by
 
9
the Free Software Foundation; either version 2.1 of the License, or (at your
 
10
option) any later version.
 
11
 
 
12
The GNU MP Library is distributed in the hope that it will be useful, but
 
13
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
14
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
15
License for more details.
 
16
 
 
17
You should have received a copy of the GNU Lesser General Public License
 
18
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 
19
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
20
MA 02111-1307, USA. */
 
21
 
 
22
/* always do assertion checking */
 
23
#define WANT_ASSERT  1
 
24
 
 
25
#include <stdio.h>
 
26
#include <stdlib.h> /* for free */
 
27
#include "gmp.h"
 
28
#include "gmp-impl.h"
 
29
#include "longlong.h"
 
30
#include "tests.h"
 
31
 
 
32
 
 
33
/* Change this to "#define TRACE(x) x" for some traces. */
 
34
#define TRACE(x) 
 
35
 
 
36
 
 
37
unsigned long
 
38
refmpz_hamdist (mpz_srcptr x, mpz_srcptr y)
 
39
{
 
40
  mp_size_t      tsize;
 
41
  mp_ptr         xp, yp;
 
42
  unsigned long  ret;
 
43
 
 
44
  if ((SIZ(x) < 0 && SIZ(y) >= 0)
 
45
      || (SIZ(y) < 0 && SIZ(x) >= 0))
 
46
    return ULONG_MAX;
 
47
 
 
48
  tsize = MAX (ABSIZ(x), ABSIZ(y));
 
49
 
 
50
  xp = refmpn_malloc_limbs (tsize);
 
51
  refmpn_zero (xp, tsize);
 
52
  refmpn_copy (xp, PTR(x), ABSIZ(x));
 
53
  
 
54
  yp = refmpn_malloc_limbs (tsize);
 
55
  refmpn_zero (yp, tsize);
 
56
  refmpn_copy (yp, PTR(y), ABSIZ(y));
 
57
 
 
58
  if (SIZ(x) < 0)
 
59
    refmpn_neg_n (xp, xp, tsize);
 
60
 
 
61
  if (SIZ(x) < 0)
 
62
    refmpn_neg_n (yp, yp, tsize);
 
63
 
 
64
  ret = refmpn_hamdist (xp, yp, tsize);
 
65
 
 
66
  free (xp);
 
67
  free (yp);
 
68
  return ret;
 
69
}
 
70
 
 
71
 
 
72
/* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */
 
73
#define JACOBI_0Z(b)  JACOBI_0LS (PTR(b)[0], SIZ(b))
 
74
 
 
75
/* (a/b) effect due to sign of b: mpz/mpz */
 
76
#define JACOBI_BSGN_ZZ_BIT1(a, b)   JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b))
 
77
 
 
78
/* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd;
 
79
   is (-1/b) if a<0, or +1 if a>=0 */
 
80
#define JACOBI_ASGN_ZZU_BIT1(a, b)  JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0])
 
81
 
 
82
int
 
83
refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig)
 
84
{
 
85
  unsigned long  twos;
 
86
  mpz_t  a, b;
 
87
  int    result_bit1 = 0;
 
88
 
 
89
  if (mpz_sgn (b_orig) == 0)
 
90
    return JACOBI_Z0 (a_orig);  /* (a/0) */
 
91
 
 
92
  if (mpz_sgn (a_orig) == 0)
 
93
    return JACOBI_0Z (b_orig);  /* (0/b) */
 
94
 
 
95
  if (mpz_even_p (a_orig) && mpz_even_p (b_orig))
 
96
    return 0;
 
97
 
 
98
  if (mpz_cmp_ui (b_orig, 1) == 0)
 
99
    return 1;
 
100
 
 
101
  mpz_init_set (a, a_orig);
 
102
  mpz_init_set (b, b_orig);
 
103
 
 
104
  if (mpz_sgn (b) < 0)
 
105
    {
 
106
      result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b);
 
107
      mpz_neg (b, b);
 
108
    }
 
109
  if (mpz_even_p (b))
 
110
    {
 
111
      twos = mpz_scan1 (b, 0L);
 
112
      mpz_tdiv_q_2exp (b, b, twos);
 
113
      result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]);
 
114
    }
 
115
 
 
116
  if (mpz_sgn (a) < 0)
 
117
    {
 
118
      result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]);
 
119
      mpz_neg (a, a);
 
120
    }
 
121
  if (mpz_even_p (a))
 
122
    {
 
123
      twos = mpz_scan1 (a, 0L);
 
124
      mpz_tdiv_q_2exp (a, a, twos);
 
125
      result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
 
126
    }
 
127
 
 
128
  for (;;)
 
129
    {
 
130
      ASSERT (mpz_odd_p (a));
 
131
      ASSERT (mpz_odd_p (b));
 
132
      ASSERT (mpz_sgn (a) > 0);
 
133
      ASSERT (mpz_sgn (b) > 0);
 
134
 
 
135
      TRACE (printf ("top\n");
 
136
             mpz_trace (" a", a);
 
137
             mpz_trace (" b", b));
 
138
 
 
139
      if (mpz_cmp (a, b) < 0)
 
140
        {
 
141
          TRACE (printf ("swap\n"));
 
142
          mpz_swap (a, b);
 
143
          result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]);
 
144
        }
 
145
 
 
146
      if (mpz_cmp_ui (b, 1) == 0)
 
147
        break;
 
148
 
 
149
      mpz_sub (a, a, b);
 
150
      TRACE (printf ("sub\n");
 
151
             mpz_trace (" a", a));
 
152
      if (mpz_sgn (a) == 0)
 
153
        goto zero;
 
154
 
 
155
      twos = mpz_scan1 (a, 0L);
 
156
      mpz_fdiv_q_2exp (a, a, twos);
 
157
      TRACE (printf ("twos %lu\n", twos);
 
158
             mpz_trace (" a", a));
 
159
      result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
 
160
    }
 
161
 
 
162
  mpz_clear (a);
 
163
  mpz_clear (b);
 
164
  return JACOBI_BIT1_TO_PN (result_bit1);
 
165
 
 
166
 zero:
 
167
  mpz_clear (a);
 
168
  mpz_clear (b);
 
169
  return 0;
 
170
}
 
171
 
 
172
/* Same as mpz_kronecker, but ignoring factors of 2 on b */
 
173
int
 
174
refmpz_jacobi (mpz_srcptr a, mpz_srcptr b)
 
175
{
 
176
  mpz_t  b_odd;
 
177
  mpz_init_set (b_odd, b);
 
178
  if (mpz_sgn (b_odd) != 0)
 
179
    mpz_fdiv_q_2exp (b_odd, b_odd, mpz_scan1 (b_odd, 0L));
 
180
  return refmpz_kronecker (a, b_odd);
 
181
}
 
182
 
 
183
int
 
184
refmpz_legendre (mpz_srcptr a, mpz_srcptr b)
 
185
{
 
186
  return refmpz_jacobi (a, b);
 
187
}
 
188
 
 
189
 
 
190
int
 
191
refmpz_kronecker_ui (mpz_srcptr a, unsigned long b)
 
192
{
 
193
  mpz_t  bz;
 
194
  int    ret;
 
195
  mpz_init_set_ui (bz, b);
 
196
  ret = refmpz_kronecker (a, bz);
 
197
  mpz_clear (bz);
 
198
  return ret;
 
199
}
 
200
 
 
201
int
 
202
refmpz_kronecker_si (mpz_srcptr a, long b)
 
203
{
 
204
  mpz_t  bz;
 
205
  int    ret;
 
206
  mpz_init_set_si (bz, b);
 
207
  ret = refmpz_kronecker (a, bz);
 
208
  mpz_clear (bz);
 
209
  return ret;
 
210
}
 
211
 
 
212
int
 
213
refmpz_ui_kronecker (unsigned long a, mpz_srcptr b)
 
214
{
 
215
  mpz_t  az;
 
216
  int    ret;
 
217
  mpz_init_set_ui (az, a);
 
218
  ret = refmpz_kronecker (az, b);
 
219
  mpz_clear (az);
 
220
  return ret;
 
221
}
 
222
 
 
223
int
 
224
refmpz_si_kronecker (long a, mpz_srcptr b)
 
225
{
 
226
  mpz_t  az;
 
227
  int    ret;
 
228
  mpz_init_set_si (az, a);
 
229
  ret = refmpz_kronecker (az, b);
 
230
  mpz_clear (az);
 
231
  return ret;
 
232
}
 
233
 
 
234
 
 
235
void
 
236
refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e)
 
237
{
 
238
  mpz_t          s, t;
 
239
  unsigned long  i;
 
240
 
 
241
  mpz_init_set_ui (t, 1L);
 
242
  mpz_init_set (s, b);
 
243
 
 
244
  if ((e & 1) != 0)
 
245
    mpz_mul (t, t, s);
 
246
 
 
247
  for (i = 2; i <= e; i <<= 1)
 
248
    {
 
249
      mpz_mul (s, s, s);
 
250
      if ((i & e) != 0)
 
251
        mpz_mul (t, t, s);
 
252
    }
 
253
 
 
254
  mpz_set (w, t);
 
255
 
 
256
  mpz_clear (s);
 
257
  mpz_clear (t);
 
258
}