bn_mp_jacobi.c 2.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105
  1. #include <tommath.h>
  2. #ifdef BN_MP_JACOBI_C
  3. /* LibTomMath, multiple-precision integer library -- Tom St Denis
  4. *
  5. * LibTomMath is a library that provides multiple-precision
  6. * integer arithmetic as well as number theoretic functionality.
  7. *
  8. * The library was designed directly after the MPI library by
  9. * Michael Fromberger but has been written from scratch with
  10. * additional optimizations in place.
  11. *
  12. * The library is free for all purposes without any express
  13. * guarantee it works.
  14. *
  15. * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
  16. */
  17. /* computes the jacobi c = (a | n) (or Legendre if n is prime)
  18. * HAC pp. 73 Algorithm 2.149
  19. */
  20. int mp_jacobi (mp_int * a, mp_int * p, int *c)
  21. {
  22. mp_int a1, p1;
  23. int k, s, r, res;
  24. mp_digit residue;
  25. /* if p <= 0 return MP_VAL */
  26. if (mp_cmp_d(p, 0) != MP_GT) {
  27. return MP_VAL;
  28. }
  29. /* step 1. if a == 0, return 0 */
  30. if (mp_iszero (a) == 1) {
  31. *c = 0;
  32. return MP_OKAY;
  33. }
  34. /* step 2. if a == 1, return 1 */
  35. if (mp_cmp_d (a, 1) == MP_EQ) {
  36. *c = 1;
  37. return MP_OKAY;
  38. }
  39. /* default */
  40. s = 0;
  41. /* step 3. write a = a1 * 2**k */
  42. if ((res = mp_init_copy (&a1, a)) != MP_OKAY) {
  43. return res;
  44. }
  45. if ((res = mp_init (&p1)) != MP_OKAY) {
  46. goto LBL_A1;
  47. }
  48. /* divide out larger power of two */
  49. k = mp_cnt_lsb(&a1);
  50. if ((res = mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
  51. goto LBL_P1;
  52. }
  53. /* step 4. if e is even set s=1 */
  54. if ((k & 1) == 0) {
  55. s = 1;
  56. } else {
  57. /* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
  58. residue = p->dp[0] & 7;
  59. if (residue == 1 || residue == 7) {
  60. s = 1;
  61. } else if (residue == 3 || residue == 5) {
  62. s = -1;
  63. }
  64. }
  65. /* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
  66. if ( ((p->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
  67. s = -s;
  68. }
  69. /* if a1 == 1 we're done */
  70. if (mp_cmp_d (&a1, 1) == MP_EQ) {
  71. *c = s;
  72. } else {
  73. /* n1 = n mod a1 */
  74. if ((res = mp_mod (p, &a1, &p1)) != MP_OKAY) {
  75. goto LBL_P1;
  76. }
  77. if ((res = mp_jacobi (&p1, &a1, &r)) != MP_OKAY) {
  78. goto LBL_P1;
  79. }
  80. *c = s * r;
  81. }
  82. /* done */
  83. res = MP_OKAY;
  84. LBL_P1:mp_clear (&p1);
  85. LBL_A1:mp_clear (&a1);
  86. return res;
  87. }
  88. #endif
  89. /* $Source: /cvs/libtom/libtommath/bn_mp_jacobi.c,v $ */
  90. /* $Revision: 1.3 $ */
  91. /* $Date: 2006/03/31 14:18:44 $ */