/* Copyright 2005-2006, Voltage Security, all rights reserved.
 */

#include "vibecrypto.h"
#include "environment.h"
#include "base.h"
#include "libctx.h"
#include "mpint.h"
#include "prime.h"
#include "errorctx.h"

/* Find the Lucas Test D value (called D in X9.31). It is a number from
 * the sequence
 *   5, -7, 9, -11, 13, -15, 17, ...
 * It is the smallest from the set of numbers for which the Jacobi
 * is -1.
 * <p>The function will check to see if Jacobi (5) is -1. If not, it
 * will try -7 (actually primeCandidate - 7). If -7 fails, it will try
 * 9. And so on. However, the function will not try every possible
 * number. At some point it will give up. The likelihood of giving up
 * in the real world is extremely small, but the code is there to
 * prevent a "near-infinite" loop. See the function to find out where
 * it gives up. If the function gives up, it will return 0 and set
 * valueFound to 0 (no/false).
 * <p>What will almost certainly happen every time in the real world,
 * is that the function will find D (the vast majority of the time it
 * will be 5 or -7) and set valueFound to 1 (yes/true).
 * <p>The function also calls for some temp MpInts, just so we don't
 * create and destroy MpInts again and again.
 *
 * @param primeCandidate
 * @param DValue The created mpInt that this function will set with the
 * found D value.
 * @param temp1
 * @param temp2
 * @param valueFound The address where the routine will deposit an
 * unsigned int indicating whether it found a D or not.
 * @return an int, 0 if the function completed successfully or a
 * non-zero error code.
 */
static int VOLT_CALLING_CONV FindLucasDValue VOLT_PROTO_LIST ((
   VoltMpInt *primeCandidate,
   VoltMpInt *DValue,
   VoltMpInt *temp1,
   VoltMpInt *temp2,
   unsigned int *valueFound
));

/* Find the Jacobi (a/n).
 * <p>This function follows the algorithm outlined in X9.42.
 * <p>The function will set the int at jacobi to 0, +1, or -1.
 *
 * @param nVal
 * @param aVal The number for which the Jacobi is requested.
 * @param jacobi The address where the routine will deposit an int, the
 * result, either 0, +1, or -1.
 * @return an int, 0 if the function completed successfully or a
 * non-zero error code.
 */
static int VOLT_CALLING_CONV FindJacobi VOLT_PROTO_LIST ((
   VoltMpInt *aVal,
   VoltMpInt *nVal,
   VoltMpInt *temp1,
   VoltMpInt *temp2,
   int *jacobi
));

int VoltLucasTest (
   VoltMpInt *primeCandidate,
   unsigned int *isPrime
   )
{
  int status, getZero;
  unsigned int index, bitLen, valueFound, getBit;
  VoltMpIntCtx *mpCtx = (VoltMpIntCtx *)(primeCandidate->mpCtx);
  VoltMpInt *NPlus1 = (VoltMpInt *)0;
  VoltMpInt *DValue = (VoltMpInt *)0;
  VoltMpInt *UVal = (VoltMpInt *)0;
  VoltMpInt *VVal = (VoltMpInt *)0;
  VoltMpInt *newU = (VoltMpInt *)0;
  VoltMpInt *newV = (VoltMpInt *)0;
  VOLT_DECLARE_ERROR_TYPE (errorType)
  VOLT_DECLARE_FNCT_LINE (fnctLine)

  /* Init the isPrime return to not prime, change it if it passes the
   * tests.
   */
  *isPrime = 0;

  do
  {
    VOLT_SET_ERROR_TYPE (errorType, VT_ERROR_TYPE_PRIMARY)

    /* In X9.31, the candidate is called N.
     * We need N + 1, we'll want the bit length of this value, and
     * later on, the actual work in each iteration of the loop will
     * depend on each of the bits of N + 1.
     */
    VOLT_SET_ERROR_TYPE (errorType, 0)
    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->CreateMpInt ((Pointer)mpCtx, &NPlus1);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->IntToMpInt (0, 1, NPlus1);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->Add (primeCandidate, NPlus1, NPlus1);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->GetBitLength (NPlus1, &bitLen);
    if (status != 0)
      break;

    /* Find D. It is from the sequence
     *   5, -7, 9, -11, 13, -15, 17, ...
     * It is the smallest from the set of numbers for which the Jacobi
     * is -1. We'll call FindLucasDValue to do that. It needs a created
     * but empty mpInt.
     * See the comments on the FindLucasDValue function for more
     * info.
     * The FindD function requires a couple temp variables. Use U and V.
     */
    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->CreateMpInt ((Pointer)mpCtx, &DValue);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->CreateMpInt ((Pointer)mpCtx, &UVal);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->CreateMpInt ((Pointer)mpCtx, &VVal);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = FindLucasDValue (
      primeCandidate, DValue, UVal, VVal, &valueFound);
    if (status != 0)
      break;

    /* If we reach this point and valueFound is 0 (status is 0), we
     * need to give up on this candidate. isPrime was initialized to
     * false.
     */
    if (valueFound == 0)
      break;

    /* For each bit in N + 1, compute U and V.
     */
    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->CreateMpInt ((Pointer)mpCtx, &newU);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->CreateMpInt ((Pointer)mpCtx, &newV);
    if (status != 0)
      break;

    /* Step 3, initialize U to 1 and V to P (which is fixed at 1).
     */
    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->IntToMpInt (0, 1, UVal);
    if (status != 0)
      break;

    VOLT_SET_FNCT_LINE (fnctLine)
    status = mpCtx->IntToMpInt (0, 1, VVal);
    if (status != 0)
      break;

    for (index = bitLen - 1; index > 0; --index)
    {
      /* Step 4: First, newU = UV mod N
       *                newV = (V^2 + DU^2) / 2  mod N
       * Do newV first so we can use newU as a temp variable.
       */
      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->Square (UVal, newU);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->ModReduce (newU, primeCandidate, newV);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->Multiply (DValue, newV, newU);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->Square (VVal, newV);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->Add (newV, newU, newU);
      if (status != 0)
        break;

      /* If newU is now even, divide by 2 by shifting right by 1 bit.
       * If newU is odd, add modulus (primeCandidate) and then divide
       * by 2.
       */
      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->GetBit (newU, 0, &getBit);
      if (status != 0)
        break;

      if (getBit == 1)
      {
        VOLT_SET_FNCT_LINE (fnctLine)
        status = mpCtx->Add (primeCandidate, newU, newU);
        if (status != 0)
          break;
      }

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->ShiftRightBits (newU, 1);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->ModReduce (newU, primeCandidate, newV);
      if (status != 0)
        break;

      /* Now compute the newU.
       */
      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->Multiply (UVal, VVal, newU);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->ModReduce (newU, primeCandidate, UVal);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->MpIntToMpInt (newV, VVal);
      if (status != 0)
        break;

      /* Get the current bit from N + 1 (index is "off" by one so we
       * can have a loop limit of 0).
       */
      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->GetBit (NPlus1, index - 1, &getBit);
      if (status != 0)
        break;

      /* If this bit is 0, we're done for this loop.
       */
      if (getBit == 0)
        continue;

      /* If this bit is set, we need to continue
       * Step 4: newU = (U + V) / 2 mod N
       *         newV = (V + DU) / 2 mod N
       * Once again, do newV first so we can use newU as a temp
       * variable.
       */
      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->Multiply (DValue, UVal, newU);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->Add (VVal, newU, newU);
      if (status != 0)
        break;

      /* If newU is now even, divide by 2 by shifting right by 1 bit.
       * If newU is odd, add modulus (primeCandidate) and then divide
       * by 2.
       */
      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->GetBit (newU, 0, &getBit);
      if (status != 0)
        break;

      if (getBit == 1)
      {
        VOLT_SET_FNCT_LINE (fnctLine)
        status = mpCtx->Add (primeCandidate, newU, newU);
        if (status != 0)
          break;
      }

      VOLT_SET_FNCT_LINE (fnctLine)
      status = mpCtx->ShiftRightBits (newU, 1);
      if (status != 0)
        break;

      VOLT_SET_FNCT_LINE (fnctLine)