Hex Artifact Content
Not logged in

Artifact 824d104d2a91001812f5d4c600112f1fcb37330a:


0000: 2f 2a 20 6d 70 7a 5f 6b 72 6f 6e 65 63 6b 65 72  /* mpz_kronecker
0010: 5f 75 69 20 2d 2d 20 6d 70 7a 2b 75 6c 6f 6e 67  _ui -- mpz+ulong
0020: 20 4b 72 6f 6e 65 63 6b 65 72 2f 4a 61 63 6f 62   Kronecker/Jacob
0030: 69 20 73 79 6d 62 6f 6c 2e 0a 0a 43 6f 70 79 72  i symbol...Copyr
0040: 69 67 68 74 20 31 39 39 39 2c 20 32 30 30 30 2c  ight 1999, 2000,
0050: 20 32 30 30 31 2c 20 32 30 30 32 20 46 72 65 65   2001, 2002 Free
0060: 20 53 6f 66 74 77 61 72 65 20 46 6f 75 6e 64 61   Software Founda
0070: 74 69 6f 6e 2c 20 49 6e 63 2e 0a 0a 54 68 69 73  tion, Inc...This
0080: 20 66 69 6c 65 20 69 73 20 70 61 72 74 20 6f 66   file is part of
0090: 20 74 68 65 20 47 4e 55 20 4d 50 20 4c 69 62 72   the GNU MP Libr
00a0: 61 72 79 2e 0a 0a 54 68 65 20 47 4e 55 20 4d 50  ary...The GNU MP
00b0: 20 4c 69 62 72 61 72 79 20 69 73 20 66 72 65 65   Library is free
00c0: 20 73 6f 66 74 77 61 72 65 3b 20 79 6f 75 20 63   software; you c
00d0: 61 6e 20 72 65 64 69 73 74 72 69 62 75 74 65 20  an redistribute 
00e0: 69 74 20 61 6e 64 2f 6f 72 20 6d 6f 64 69 66 79  it and/or modify
00f0: 0a 69 74 20 75 6e 64 65 72 20 74 68 65 20 74 65  .it under the te
0100: 72 6d 73 20 6f 66 20 74 68 65 20 47 4e 55 20 4c  rms of the GNU L
0110: 65 73 73 65 72 20 47 65 6e 65 72 61 6c 20 50 75  esser General Pu
0120: 62 6c 69 63 20 4c 69 63 65 6e 73 65 20 61 73 20  blic License as 
0130: 70 75 62 6c 69 73 68 65 64 20 62 79 0a 74 68 65  published by.the
0140: 20 46 72 65 65 20 53 6f 66 74 77 61 72 65 20 46   Free Software F
0150: 6f 75 6e 64 61 74 69 6f 6e 3b 20 65 69 74 68 65  oundation; eithe
0160: 72 20 76 65 72 73 69 6f 6e 20 32 2e 31 20 6f 66  r version 2.1 of
0170: 20 74 68 65 20 4c 69 63 65 6e 73 65 2c 20 6f 72   the License, or
0180: 20 28 61 74 20 79 6f 75 72 0a 6f 70 74 69 6f 6e   (at your.option
0190: 29 20 61 6e 79 20 6c 61 74 65 72 20 76 65 72 73  ) any later vers
01a0: 69 6f 6e 2e 0a 0a 54 68 65 20 47 4e 55 20 4d 50  ion...The GNU MP
01b0: 20 4c 69 62 72 61 72 79 20 69 73 20 64 69 73 74   Library is dist
01c0: 72 69 62 75 74 65 64 20 69 6e 20 74 68 65 20 68  ributed in the h
01d0: 6f 70 65 20 74 68 61 74 20 69 74 20 77 69 6c 6c  ope that it will
01e0: 20 62 65 20 75 73 65 66 75 6c 2c 20 62 75 74 0a   be useful, but.
01f0: 57 49 54 48 4f 55 54 20 41 4e 59 20 57 41 52 52  WITHOUT ANY WARR
0200: 41 4e 54 59 3b 20 77 69 74 68 6f 75 74 20 65 76  ANTY; without ev
0210: 65 6e 20 74 68 65 20 69 6d 70 6c 69 65 64 20 77  en the implied w
0220: 61 72 72 61 6e 74 79 20 6f 66 20 4d 45 52 43 48  arranty of MERCH
0230: 41 4e 54 41 42 49 4c 49 54 59 0a 6f 72 20 46 49  ANTABILITY.or FI
0240: 54 4e 45 53 53 20 46 4f 52 20 41 20 50 41 52 54  TNESS FOR A PART
0250: 49 43 55 4c 41 52 20 50 55 52 50 4f 53 45 2e 20  ICULAR PURPOSE. 
0260: 20 53 65 65 20 74 68 65 20 47 4e 55 20 4c 65 73   See the GNU Les
0270: 73 65 72 20 47 65 6e 65 72 61 6c 20 50 75 62 6c  ser General Publ
0280: 69 63 0a 4c 69 63 65 6e 73 65 20 66 6f 72 20 6d  ic.License for m
0290: 6f 72 65 20 64 65 74 61 69 6c 73 2e 0a 0a 59 6f  ore details...Yo
02a0: 75 20 73 68 6f 75 6c 64 20 68 61 76 65 20 72 65  u should have re
02b0: 63 65 69 76 65 64 20 61 20 63 6f 70 79 20 6f 66  ceived a copy of
02c0: 20 74 68 65 20 47 4e 55 20 4c 65 73 73 65 72 20   the GNU Lesser 
02d0: 47 65 6e 65 72 61 6c 20 50 75 62 6c 69 63 20 4c  General Public L
02e0: 69 63 65 6e 73 65 0a 61 6c 6f 6e 67 20 77 69 74  icense.along wit
02f0: 68 20 74 68 65 20 47 4e 55 20 4d 50 20 4c 69 62  h the GNU MP Lib
0300: 72 61 72 79 3b 20 73 65 65 20 74 68 65 20 66 69  rary; see the fi
0310: 6c 65 20 43 4f 50 59 49 4e 47 2e 4c 49 42 2e 20  le COPYING.LIB. 
0320: 20 49 66 20 6e 6f 74 2c 20 77 72 69 74 65 20 74   If not, write t
0330: 6f 0a 74 68 65 20 46 72 65 65 20 53 6f 66 74 77  o.the Free Softw
0340: 61 72 65 20 46 6f 75 6e 64 61 74 69 6f 6e 2c 20  are Foundation, 
0350: 49 6e 63 2e 2c 20 35 39 20 54 65 6d 70 6c 65 20  Inc., 59 Temple 
0360: 50 6c 61 63 65 20 2d 20 53 75 69 74 65 20 33 33  Place - Suite 33
0370: 30 2c 20 42 6f 73 74 6f 6e 2c 0a 4d 41 20 30 32  0, Boston,.MA 02
0380: 31 31 31 2d 31 33 30 37 2c 20 55 53 41 2e 20 2a  111-1307, USA. *
0390: 2f 0a 0a 23 69 6e 63 6c 75 64 65 20 22 67 6d 70  /..#include "gmp
03a0: 2e 68 22 0a 23 69 6e 63 6c 75 64 65 20 22 67 6d  .h".#include "gm
03b0: 70 2d 69 6d 70 6c 2e 68 22 0a 23 69 6e 63 6c 75  p-impl.h".#inclu
03c0: 64 65 20 22 6c 6f 6e 67 6c 6f 6e 67 2e 68 22 0a  de "longlong.h".
03d0: 0a 0a 2f 2a 20 54 68 69 73 20 69 6d 70 6c 65 6d  ../* This implem
03e0: 65 6e 74 61 74 69 6f 6e 20 64 65 70 65 6e 64 73  entation depends
03f0: 20 6f 6e 20 42 49 54 53 5f 50 45 52 5f 4d 50 5f   on BITS_PER_MP_
0400: 4c 49 4d 42 20 62 65 69 6e 67 20 65 76 65 6e 2c  LIMB being even,
0410: 20 73 6f 20 74 68 61 74 0a 20 20 20 28 61 2f 32   so that.   (a/2
0420: 29 5e 42 49 54 53 5f 50 45 52 5f 4d 50 5f 4c 49  )^BITS_PER_MP_LI
0430: 4d 42 20 3d 20 31 20 61 6e 64 20 73 6f 20 74 68  MB = 1 and so th
0440: 65 72 65 27 73 20 6e 6f 20 6e 65 65 64 20 74 6f  ere's no need to
0450: 20 70 61 79 20 61 74 74 65 6e 74 69 6f 6e 20 74   pay attention t
0460: 6f 20 68 6f 77 0a 20 20 20 6d 61 6e 79 20 6c 6f  o how.   many lo
0470: 77 20 7a 65 72 6f 20 6c 69 6d 62 73 20 61 72 65  w zero limbs are
0480: 20 73 74 72 69 70 70 65 64 2e 20 20 2a 2f 0a 23   stripped.  */.#
0490: 69 66 20 42 49 54 53 5f 50 45 52 5f 4d 50 5f 4c  if BITS_PER_MP_L
04a0: 49 4d 42 20 25 20 32 20 21 3d 20 30 0a 45 72 72  IMB % 2 != 0.Err
04b0: 6f 72 2c 20 65 72 72 6f 72 2c 20 75 6e 73 75 70  or, error, unsup
04c0: 70 6f 72 74 65 64 20 42 49 54 53 5f 50 45 52 5f  ported BITS_PER_
04d0: 4d 50 5f 4c 49 4d 42 0a 23 65 6e 64 69 66 0a 0a  MP_LIMB.#endif..
04e0: 0a 69 6e 74 0a 6d 70 7a 5f 6b 72 6f 6e 65 63 6b  .int.mpz_kroneck
04f0: 65 72 5f 75 69 20 28 6d 70 7a 5f 73 72 63 70 74  er_ui (mpz_srcpt
0500: 72 20 61 2c 20 75 6e 73 69 67 6e 65 64 20 6c 6f  r a, unsigned lo
0510: 6e 67 20 62 29 0a 7b 0a 20 20 6d 70 5f 73 72 63  ng b).{.  mp_src
0520: 70 74 72 20 20 61 5f 70 74 72 20 3d 20 50 54 52  ptr  a_ptr = PTR
0530: 28 61 29 3b 0a 20 20 6d 70 5f 73 69 7a 65 5f 74  (a);.  mp_size_t
0540: 20 20 61 5f 73 69 7a 65 3b 0a 20 20 6d 70 5f 6c    a_size;.  mp_l
0550: 69 6d 62 5f 74 20 20 61 5f 72 65 6d 3b 0a 20 20  imb_t  a_rem;.  
0560: 69 6e 74 20 20 20 20 20 20 20 20 72 65 73 75 6c  int        resul
0570: 74 5f 62 69 74 31 3b 0a 0a 20 20 61 5f 73 69 7a  t_bit1;..  a_siz
0580: 65 20 3d 20 53 49 5a 28 61 29 3b 0a 20 20 69 66  e = SIZ(a);.  if
0590: 20 28 61 5f 73 69 7a 65 20 3d 3d 20 30 29 0a 20   (a_size == 0). 
05a0: 20 20 20 72 65 74 75 72 6e 20 4a 41 43 4f 42 49     return JACOBI
05b0: 5f 30 55 20 28 62 29 3b 0a 0a 20 20 69 66 20 28  _0U (b);..  if (
05c0: 62 20 3e 20 47 4d 50 5f 4e 55 4d 42 5f 4d 41 58  b > GMP_NUMB_MAX
05d0: 29 0a 20 20 20 20 7b 0a 20 20 20 20 20 20 6d 70  ).    {.      mp
05e0: 5f 6c 69 6d 62 5f 74 20 20 62 6c 69 6d 62 73 5b  _limb_t  blimbs[
05f0: 32 5d 3b 0a 20 20 20 20 20 20 6d 70 7a 5f 74 20  2];.      mpz_t 
0600: 20 20 20 20 20 62 7a 3b 0a 20 20 20 20 20 20 41       bz;.      A
0610: 4c 4c 4f 43 28 62 7a 29 20 3d 20 6e 75 6d 62 65  LLOC(bz) = numbe
0620: 72 6f 66 20 28 62 6c 69 6d 62 73 29 3b 0a 20 20  rof (blimbs);.  
0630: 20 20 20 20 50 54 52 28 62 7a 29 20 3d 20 62 6c      PTR(bz) = bl
0640: 69 6d 62 73 3b 0a 20 20 20 20 20 20 6d 70 7a 5f  imbs;.      mpz_
0650: 73 65 74 5f 75 69 20 28 62 7a 2c 20 62 29 3b 0a  set_ui (bz, b);.
0660: 20 20 20 20 20 20 72 65 74 75 72 6e 20 6d 70 7a        return mpz
0670: 5f 6b 72 6f 6e 65 63 6b 65 72 20 28 61 2c 20 62  _kronecker (a, b
0680: 7a 29 3b 0a 20 20 20 20 7d 0a 0a 20 20 69 66 20  z);.    }..  if 
0690: 28 28 62 20 26 20 31 29 20 21 3d 20 30 29 0a 20  ((b & 1) != 0). 
06a0: 20 20 20 7b 0a 20 20 20 20 20 20 72 65 73 75 6c     {.      resul
06b0: 74 5f 62 69 74 31 20 3d 20 4a 41 43 4f 42 49 5f  t_bit1 = JACOBI_
06c0: 41 53 47 4e 5f 53 55 5f 42 49 54 31 20 28 61 5f  ASGN_SU_BIT1 (a_
06d0: 73 69 7a 65 2c 20 62 29 3b 0a 20 20 20 20 7d 0a  size, b);.    }.
06e0: 20 20 65 6c 73 65 0a 20 20 20 20 7b 0a 20 20 20    else.    {.   
06f0: 20 20 20 6d 70 5f 6c 69 6d 62 5f 74 20 20 61 5f     mp_limb_t  a_
0700: 6c 6f 77 20 3d 20 61 5f 70 74 72 5b 30 5d 3b 0a  low = a_ptr[0];.
0710: 20 20 20 20 20 20 69 6e 74 20 20 20 20 20 20 20        int       
0720: 20 74 77 6f 73 3b 0a 0a 20 20 20 20 20 20 69 66   twos;..      if
0730: 20 28 62 20 3d 3d 20 30 29 0a 20 20 20 20 20 20   (b == 0).      
0740: 20 20 72 65 74 75 72 6e 20 4a 41 43 4f 42 49 5f    return JACOBI_
0750: 4c 53 30 20 28 61 5f 6c 6f 77 2c 20 61 5f 73 69  LS0 (a_low, a_si
0760: 7a 65 29 3b 20 20 20 2f 2a 20 28 61 2f 30 29 20  ze);   /* (a/0) 
0770: 2a 2f 0a 0a 20 20 20 20 20 20 69 66 20 28 21 20  */..      if (! 
0780: 28 61 5f 6c 6f 77 20 26 20 31 29 29 0a 20 20 20  (a_low & 1)).   
0790: 20 20 20 20 20 72 65 74 75 72 6e 20 30 3b 20 20       return 0;  
07a0: 2f 2a 20 28 65 76 65 6e 2f 65 76 65 6e 29 3d 30  /* (even/even)=0
07b0: 20 2a 2f 0a 0a 20 20 20 20 20 20 2f 2a 20 28 61   */..      /* (a
07c0: 2f 32 29 3d 28 32 2f 61 29 20 66 6f 72 20 61 20  /2)=(2/a) for a 
07d0: 6f 64 64 20 2a 2f 0a 20 20 20 20 20 20 63 6f 75  odd */.      cou
07e0: 6e 74 5f 74 72 61 69 6c 69 6e 67 5f 7a 65 72 6f  nt_trailing_zero
07f0: 73 20 28 74 77 6f 73 2c 20 62 29 3b 0a 20 20 20  s (twos, b);.   
0800: 20 20 20 62 20 3e 3e 3d 20 74 77 6f 73 3b 0a 20     b >>= twos;. 
0810: 20 20 20 20 20 72 65 73 75 6c 74 5f 62 69 74 31       result_bit1
0820: 20 3d 20 28 4a 41 43 4f 42 49 5f 54 57 4f 53 5f   = (JACOBI_TWOS_
0830: 55 5f 42 49 54 31 20 28 74 77 6f 73 2c 20 61 5f  U_BIT1 (twos, a_
0840: 6c 6f 77 29 0a 20 20 20 20 20 20 20 20 20 20 20  low).           
0850: 20 20 20 20 20 20 20 20 20 20 5e 20 4a 41 43 4f            ^ JACO
0860: 42 49 5f 41 53 47 4e 5f 53 55 5f 42 49 54 31 20  BI_ASGN_SU_BIT1 
0870: 28 61 5f 73 69 7a 65 2c 20 62 29 29 3b 0a 20 20  (a_size, b));.  
0880: 20 20 7d 0a 0a 20 20 69 66 20 28 62 20 3d 3d 20    }..  if (b == 
0890: 31 29 0a 20 20 20 20 72 65 74 75 72 6e 20 4a 41  1).    return JA
08a0: 43 4f 42 49 5f 42 49 54 31 5f 54 4f 5f 50 4e 20  COBI_BIT1_TO_PN 
08b0: 28 72 65 73 75 6c 74 5f 62 69 74 31 29 3b 20 20  (result_bit1);  
08c0: 2f 2a 20 28 61 2f 31 29 3d 31 20 66 6f 72 20 61  /* (a/1)=1 for a
08d0: 6e 79 20 61 20 2a 2f 0a 0a 20 20 61 5f 73 69 7a  ny a */..  a_siz
08e0: 65 20 3d 20 41 42 53 28 61 5f 73 69 7a 65 29 3b  e = ABS(a_size);
08f0: 0a 0a 20 20 2f 2a 20 28 61 2f 62 29 20 3d 20 28  ..  /* (a/b) = (
0900: 61 20 6d 6f 64 20 62 20 2f 20 62 29 20 2a 2f 0a  a mod b / b) */.
0910: 20 20 4a 41 43 4f 42 49 5f 4d 4f 44 5f 4f 52 5f    JACOBI_MOD_OR_
0920: 4d 4f 44 45 58 41 43 54 5f 31 5f 4f 44 44 20 28  MODEXACT_1_ODD (
0930: 72 65 73 75 6c 74 5f 62 69 74 31 2c 20 61 5f 72  result_bit1, a_r
0940: 65 6d 2c 20 61 5f 70 74 72 2c 20 61 5f 73 69 7a  em, a_ptr, a_siz
0950: 65 2c 20 62 29 3b 0a 20 20 72 65 74 75 72 6e 20  e, b);.  return 
0960: 6d 70 6e 5f 6a 61 63 6f 62 69 5f 62 61 73 65 20  mpn_jacobi_base 
0970: 28 61 5f 72 65 6d 2c 20 28 6d 70 5f 6c 69 6d 62  (a_rem, (mp_limb
0980: 5f 74 29 20 62 2c 20 72 65 73 75 6c 74 5f 62 69  _t) b, result_bi
0990: 74 31 29 3b 0a 7d 0a                             t1);.}.