← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/gmp/mpz/kronsz.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/mpz/kronsz.c
 
1
/* mpz_si_kronecker -- long+mpz Kronecker/Jacobi symbol.
 
2
 
 
3
Copyright 1999, 2000, 2001, 2002 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
#include "gmp.h"
 
23
#include "gmp-impl.h"
 
24
#include "longlong.h"
 
25
 
 
26
 
 
27
/* This implementation depends on BITS_PER_MP_LIMB being even, so that
 
28
   (a/2)^BITS_PER_MP_LIMB = 1 and so there's no need to pay attention to how
 
29
   many low zero limbs are stripped.  */
 
30
#if BITS_PER_MP_LIMB % 2 != 0
 
31
Error, error, unsupported BITS_PER_MP_LIMB
 
32
#endif
 
33
 
 
34
 
 
35
int
 
36
mpz_si_kronecker (long a, mpz_srcptr b)
 
37
{
 
38
  mp_srcptr  b_ptr = PTR(b);
 
39
  mp_limb_t  b_low = b_ptr[0];
 
40
  mp_size_t  b_size = SIZ (b);
 
41
  mp_size_t  b_abs_size = ABS (b_size);
 
42
  mp_limb_t  a_limb, b_rem;
 
43
  unsigned   twos;
 
44
  int        result_bit1;
 
45
 
 
46
#if GMP_NUMB_BITS < BITS_PER_ULONG
 
47
  if (a > GMP_NUMB_MAX || a < -GMP_NUMB_MAX)
 
48
    {
 
49
      mp_limb_t  alimbs[2];
 
50
      mpz_t      az;
 
51
      ALLOC(az) = numberof (alimbs);
 
52
      PTR(az) = alimbs;
 
53
      mpz_set_si (az, a);
 
54
      return mpz_kronecker (az, b);
 
55
    }
 
56
#endif
 
57
 
 
58
  if (b_abs_size == 0)
 
59
    return JACOBI_S0 (a);  /* (a/0) */
 
60
 
 
61
  /* account for the effect of the sign of b, then ignore it */
 
62
  result_bit1 = JACOBI_BSGN_SS_BIT1 (a, b_size);
 
63
 
 
64
  if ((b_low & 1) != 0)
 
65
    {
 
66
      /* b odd */
 
67
 
 
68
      result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low);
 
69
      a_limb = (unsigned long) ABS(a);
 
70
 
 
71
      if ((a_limb & 1) == 0)
 
72
        {
 
73
          /* (0/b)=1 for b=+/-1, 0 otherwise */
 
74
          if (a_limb == 0)
 
75
            return (b_abs_size == 1 && b_low == 1);
 
76
 
 
77
          /* a even, b odd */
 
78
          count_trailing_zeros (twos, a_limb);
 
79
          a_limb >>= twos;
 
80
          /* (a*2^n/b) = (a/b) * twos(n,a) */
 
81
          result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b_low);
 
82
        }
 
83
    }
 
84
  else
 
85
    {
 
86
      /* (even/even)=0, and (0/b)=0 for b!=+/-1 */
 
87
      if ((a & 1) == 0)
 
88
        return 0;
 
89
 
 
90
      /* a odd, b even */
 
91
 
 
92
      /* establish shifted b_low with valid bit1 for use with ASGN and RECIP */
 
93
      MPN_STRIP_LOW_ZEROS_NOT_ZERO (b_ptr, b_abs_size, b_low);
 
94
      if ((b_low & 1) == 0)
 
95
        {
 
96
          if (b_low == GMP_NUMB_HIGHBIT)
 
97
            {
 
98
              if (b_abs_size == 1)
 
99
                {
 
100
                  /* (a/0x80000000) = (a/2)^(BPML-1) */
 
101
                  result_bit1 ^= JACOBI_TWOS_U_BIT1 (GMP_NUMB_BITS - 1, a);
 
102
                  return JACOBI_BIT1_TO_PN (result_bit1);
 
103
                }
 
104
 
 
105
              /* b_abs_size > 1 */
 
106
              b_low = b_ptr[1] << 1;
 
107
            }
 
108
          else
 
109
            {
 
110
              count_trailing_zeros (twos, b_low);
 
111
              b_low >>= twos;
 
112
            }
 
113
        }
 
114
 
 
115
      result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low);
 
116
      a_limb = (unsigned long) ABS(a);
 
117
    }
 
118
 
 
119
  if (a_limb == 1)
 
120
    return JACOBI_BIT1_TO_PN (result_bit1);  /* (1/b)=1 */
 
121
 
 
122
  /* (a/b*2^n) = (b*2^n mod a / a) * recip(a,b) */
 
123
  JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, b_ptr, b_abs_size, a_limb);
 
124
  result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a_limb, b_low);
 
125
  return mpn_jacobi_base (b_rem, a_limb, result_bit1);
 
126
}