mp_kronecker.c

   1 #include "tommath_private.h"
   2 #ifdef MP_KRONECKER_C
   3
   4 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   5 /* SPDX-License-Identifier: Unlicense */
   6
   7 /*
   8    Kronecker symbol (a|p)
   9    Straightforward implementation of algorithm 1.4.10 in
  10    Henri Cohen: "A Course in Computational Algebraic Number Theory"
  11
  12    @book{cohen2013course,
  13      title={A course in computational algebraic number theory},
  14      author={Cohen, Henri},
  15      volume={138},
  16      year={2013},
  17      publisher={Springer Science \& Business Media}
  18     }
  19  */
  20 mp_err mp_kronecker(const mp_int *a, const mp_int *p, int *c)
  21 {
  22    mp_int a1, p1, r;
  23    mp_err err;
  24    int v, k;
  25
  26    static const char table[] = {0, 1, 0, -1, 0, -1, 0, 1};
  27
  28    if (mp_iszero(p)) {
  29       if ((a->used == 1) && (a->dp[0] == 1u)) {
  30          *c = 1;
  31       } else {
  32          *c = 0;
  33       }
  34       return MP_OKAY;
  35    }
  36
  37    if (mp_iseven(a) && mp_iseven(p)) {
  38       *c = 0;
  39       return MP_OKAY;
  40    }
  41
  42    if ((err = mp_init_copy(&a1, a)) != MP_OKAY) {
  43       return err;
  44    }
  45    if ((err = mp_init_copy(&p1, p)) != MP_OKAY) {
  46       goto LBL_KRON_0;
  47    }
  48
  49    v = mp_cnt_lsb(&p1);
  50    if ((err = mp_div_2d(&p1, v, &p1, NULL)) != MP_OKAY) {
  51       goto LBL_KRON_1;
  52    }
  53
  54    if ((v & 1) == 0) {
  55       k = 1;
  56    } else {
  57       k = table[a->dp[0] & 7u];
  58    }
  59
  60    if (mp_isneg(&p1)) {
  61       p1.sign = MP_ZPOS;
  62       if (mp_isneg(&a1)) {
  63          k = -k;
  64       }
  65    }
  66
  67    if ((err = mp_init(&r)) != MP_OKAY) {
  68       goto LBL_KRON_1;
  69    }
  70
  71    for (;;) {
  72       if (mp_iszero(&a1)) {
  73          if (mp_cmp_d(&p1, 1uL) == MP_EQ) {
  74             *c = k;
  75             goto LBL_KRON;
  76          } else {
  77             *c = 0;
  78             goto LBL_KRON;
  79          }
  80       }
  81
  82       v = mp_cnt_lsb(&a1);
  83       if ((err = mp_div_2d(&a1, v, &a1, NULL)) != MP_OKAY) {
  84          goto LBL_KRON;
  85       }
  86
  87       if ((v & 1) == 1) {
  88          k = k * table[p1.dp[0] & 7u];
  89       }
  90
  91       if (mp_isneg(&a1)) {
  92          /*
  93           * Compute k = (-1)^((a1)*(p1-1)/4) * k
  94           * a1.dp[0] + 1 cannot overflow because the MSB
  95           * of the type mp_digit is not set by definition
  96           */
  97          if (((a1.dp[0] + 1u) & p1.dp[0] & 2u) != 0u) {
  98             k = -k;
  99          }
 100       } else {
 101          /* compute k = (-1)^((a1-1)*(p1-1)/4) * k */
 102          if ((a1.dp[0] & p1.dp[0] & 2u) != 0u) {
 103             k = -k;
 104          }
 105       }
 106
 107       if ((err = mp_copy(&a1, &r)) != MP_OKAY) {
 108          goto LBL_KRON;
 109       }
 110       r.sign = MP_ZPOS;
 111       if ((err = mp_mod(&p1, &r, &a1)) != MP_OKAY) {
 112          goto LBL_KRON;
 113       }
 114       if ((err = mp_copy(&r, &p1)) != MP_OKAY) {
 115          goto LBL_KRON;
 116       }
 117    }
 118
 119 LBL_KRON:
 120    mp_clear(&r);
 121 LBL_KRON_1:
 122    mp_clear(&p1);
 123 LBL_KRON_0:
 124    mp_clear(&a1);
 125
 126    return err;
 127 }
 128
 129 #endif