# include <cstdlib>
# include <ctime>
# include <iostream>

using namespace std;

# include "prime_fermat.hpp"

//****************************************************************************80

int gcd ( int a, int b )

//****************************************************************************80
//
//  Purpose:
//
//    gcd() calculates the greatest common divisor of two integers.
//
//  Discussion:
//
//    This function uses recursion.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    20 September 2024
//
//  Author:
//
//    Original C++ code on Geeks for Geeks.
//    This version by John Burkardt.
//
//  Input:
//
//    int a, b: two numbers.
//
//  Output:
//
//    int gcd: the greatest common divisor of a and b.
//
{
  if ( a < b )
  {
    return gcd ( b, a );
  }
  else if ( a % b == 0 )
  {
    return b;
  }
  else
  {
    return gcd ( b, a % b );  
  }
}
//****************************************************************************80

bool is_prime ( unsigned int n, int k )

//****************************************************************************80
//
//  Purpose:
//
//    is_prime() applies Fermat's test for primality.
//
//  Discussion:
//
//    The test is only probabilistic.
//    If the input n is a prime, is_prime(n) will be true.
//    If n is not prime, then is_prime(n) is probably false.
//    Reliability can be increased by increasing the value of k, which
//    repeats the test with different bases.
//    However, there are some values of n which are not prime, but for
//    which is_prime(n) will always return true.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    20 September 2024
//
//  Author:
//
//    Original C++ code on Geeks for Geeks.
//    This version by John Burkardt.
//
//  Reference:
//
//    Thomas Cormen, Charles Leiserson, Ronald Rivest, Clifford Stein,
//    Introduction to Algorithms,
//    Section 31.8: Primality testing,
//    MIT Press; McGraw-Hill, 2001.
//
//  Input:
//
//    unsigned int n: the value to be tested.
//
//    int k: the number of times the test is to be carried out.
//
//  Output:
//
//    bool is_prime: true if n passed the test.
//
{
  int a;
  int i;

  if ( n <= 1 || n == 4 )
  {
    return false;
  }
  if ( n <= 3 )
  {
    return true;
  }
//
//  Do the test k times.
//
  for ( i = 0; i < k; k++ )
  {
//
//  Pick a random number in [2..n-2]        
//
    a = 2 + rand() % ( n - 4 );  
// 
//  Checking that a and n are co-prime.
//
    if ( gcd ( n, a ) != 1 )
    {
      return false;
    }
//
//  Check Fermat's little theorem.
//
    if ( power ( a, n - 1, n ) != 1 )
    {
      return false;
    }
  }
 
  return true;
}
//****************************************************************************80

int power ( int a, unsigned int n, int p )

//****************************************************************************80
//
//  Purpose:
//
//    power() computes a^n mod p.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    20 September 2024
//
//  Author:
//
//    Original C++ code on Geeks for Geeks.
//    This version by John Burkardt.
//
//  Input:
//
//    int a: the base.
//
//    unsigned int n: the power.
//
//    int p: the calculation is to be done mod p.
//
//  Output:
//
//    int power: the value of a^n mod p.
//
{
  int res = 1;
//
//  Reduce a if p <= a.
//
  a = a % p;
 
  while ( 0 < n )
  {
//
//  If n is odd, multiply result by a.
//
    if ( n & 1 )
    {
      res = ( res * a ) % p;
    } 
//
//  n must be even now.  Divide it by 2.
//
    n = n >> 1;
    a = ( a * a ) % p;
  }
  return res;
}