Line data    Source code

       1             : #include "tommath_private.h"
       2             : #ifdef BN_MP_PRIME_FROBENIUS_UNDERWOOD_C
       3             : 
       4             : /* LibTomMath, multiple-precision integer library -- Tom St Denis */
       5             : /* SPDX-License-Identifier: Unlicense */
       6             : 
       7             : /*
       8             :  *  See file bn_mp_prime_is_prime.c or the documentation in doc/bn.tex for the details
       9             :  */
      10             : #ifndef LTM_USE_ONLY_MR
      11             : 
      12             : #ifdef MP_8BIT
      13             : /*
      14             :  * floor of positive solution of
      15             :  * (2^16)-1 = (a+4)*(2*a+5)
      16             :  * TODO: Both values are smaller than N^(1/4), would have to use a bigint
      17             :  *       for a instead but any a biger than about 120 are already so rare that
      18             :  *       it is possible to ignore them and still get enough pseudoprimes.
      19             :  *       But it is still a restriction of the set of available pseudoprimes
      20             :  *       which makes this implementation less secure if used stand-alone.
      21             :  */
      22             : #define LTM_FROBENIUS_UNDERWOOD_A 177
      23             : #else
      24             : #define LTM_FROBENIUS_UNDERWOOD_A 32764
      25             : #endif
      26           0 : mp_err mp_prime_frobenius_underwood(const mp_int *N, mp_bool *result)
      27             : {
      28             :    mp_int T1z, T2z, Np1z, sz, tz;
      29             : 
      30             :    int a, ap2, length, i, j;
      31             :    mp_err err;
      32             : 
      33           0 :    *result = MP_NO;
      34             : 
      35           0 :    if ((err = mp_init_multi(&T1z, &T2z, &Np1z, &sz, &tz, NULL)) != MP_OKAY) {
      36           0 :       return err;
      37             :    }
      38             : 
      39           0 :    for (a = 0; a < LTM_FROBENIUS_UNDERWOOD_A; a++) {
      40             :       /* TODO: That's ugly! No, really, it is! */
      41           0 :       if ((a==2) || (a==4) || (a==7) || (a==8) || (a==10) ||
      42           0 :           (a==14) || (a==18) || (a==23) || (a==26) || (a==28)) {
      43           0 :          continue;
      44             :       }
      45             :       /* (32764^2 - 4) < 2^31, no bigint for >MP_8BIT needed) */
      46           0 :       mp_set_u32(&T1z, (uint32_t)a);
      47             : 
      48           0 :       if ((err = mp_sqr(&T1z, &T1z)) != MP_OKAY)                  goto LBL_FU_ERR;
      49             : 
      50           0 :       if ((err = mp_sub_d(&T1z, 4uL, &T1z)) != MP_OKAY)           goto LBL_FU_ERR;
      51             : 
      52           0 :       if ((err = mp_kronecker(&T1z, N, &j)) != MP_OKAY)           goto LBL_FU_ERR;
      53             : 
      54           0 :       if (j == -1) {
      55           0 :          break;
      56             :       }
      57             : 
      58           0 :       if (j == 0) {
      59             :          /* composite */
      60           0 :          goto LBL_FU_ERR;
      61             :       }
      62             :    }
      63             :    /* Tell it a composite and set return value accordingly */
      64           0 :    if (a >= LTM_FROBENIUS_UNDERWOOD_A) {
      65           0 :       err = MP_ITER;
      66           0 :       goto LBL_FU_ERR;
      67             :    }
      68             :    /* Composite if N and (a+4)*(2*a+5) are not coprime */
      69           0 :    mp_set_u32(&T1z, (uint32_t)((a+4)*((2*a)+5)));
      70             : 
      71           0 :    if ((err = mp_gcd(N, &T1z, &T1z)) != MP_OKAY)                  goto LBL_FU_ERR;
      72             : 
      73           0 :    if (!((T1z.used == 1) && (T1z.dp[0] == 1u)))                   goto LBL_FU_ERR;
      74             : 
      75           0 :    ap2 = a + 2;
      76           0 :    if ((err = mp_add_d(N, 1uL, &Np1z)) != MP_OKAY)                goto LBL_FU_ERR;
      77             : 
      78           0 :    mp_set(&sz, 1uL);
      79           0 :    mp_set(&tz, 2uL);
      80           0 :    length = mp_count_bits(&Np1z);
      81             : 
      82           0 :    for (i = length - 2; i >= 0; i--) {
      83             :       /*
      84             :        * temp = (sz*(a*sz+2*tz))%N;
      85             :        * tz   = ((tz-sz)*(tz+sz))%N;
      86             :        * sz   = temp;
      87             :        */
      88           0 :       if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY)                 goto LBL_FU_ERR;
      89             : 
      90             :       /* a = 0 at about 50% of the cases (non-square and odd input) */
      91           0 :       if (a != 0) {
      92           0 :          if ((err = mp_mul_d(&sz, (mp_digit)a, &T1z)) != MP_OKAY) goto LBL_FU_ERR;
      93           0 :          if ((err = mp_add(&T1z, &T2z, &T2z)) != MP_OKAY)         goto LBL_FU_ERR;
      94             :       }
      95             : 
      96           0 :       if ((err = mp_mul(&T2z, &sz, &T1z)) != MP_OKAY)             goto LBL_FU_ERR;
      97           0 :       if ((err = mp_sub(&tz, &sz, &T2z)) != MP_OKAY)              goto LBL_FU_ERR;
      98           0 :       if ((err = mp_add(&sz, &tz, &sz)) != MP_OKAY)               goto LBL_FU_ERR;
      99           0 :       if ((err = mp_mul(&sz, &T2z, &tz)) != MP_OKAY)              goto LBL_FU_ERR;
     100           0 :       if ((err = mp_mod(&tz, N, &tz)) != MP_OKAY)                 goto LBL_FU_ERR;
     101           0 :       if ((err = mp_mod(&T1z, N, &sz)) != MP_OKAY)                goto LBL_FU_ERR;
     102           0 :       if (s_mp_get_bit(&Np1z, (unsigned int)i) == MP_YES) {
     103             :          /*
     104             :           *  temp = (a+2) * sz + tz
     105             :           *  tz   = 2 * tz - sz
     106             :           *  sz   = temp
     107             :           */
     108           0 :          if (a == 0) {
     109           0 :             if ((err = mp_mul_2(&sz, &T1z)) != MP_OKAY)           goto LBL_FU_ERR;
     110             :          } else {
     111           0 :             if ((err = mp_mul_d(&sz, (mp_digit)ap2, &T1z)) != MP_OKAY) goto LBL_FU_ERR;
     112             :          }
     113           0 :          if ((err = mp_add(&T1z, &tz, &T1z)) != MP_OKAY)          goto LBL_FU_ERR;
     114           0 :          if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY)              goto LBL_FU_ERR;
     115           0 :          if ((err = mp_sub(&T2z, &sz, &tz)) != MP_OKAY)           goto LBL_FU_ERR;
     116           0 :          mp_exch(&sz, &T1z);
     117             :       }
     118             :    }
     119             : 
     120           0 :    mp_set_u32(&T1z, (uint32_t)((2 * a) + 5));
     121           0 :    if ((err = mp_mod(&T1z, N, &T1z)) != MP_OKAY)                  goto LBL_FU_ERR;
     122           0 :    if (MP_IS_ZERO(&sz) && (mp_cmp(&tz, &T1z) == MP_EQ)) {
     123           0 :       *result = MP_YES;
     124             :    }
     125             : 
     126           0 : LBL_FU_ERR:
     127           0 :    mp_clear_multi(&tz, &sz, &Np1z, &T2z, &T1z, NULL);
     128           0 :    return err;
     129             : }
     130             : 
     131             : #endif
     132             : #endif