#! /usr/bin/env python3
#
def is_prime1 ( n ):

#*****************************************************************************80
#
## is_prime1() reports whether an integer is prime.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 January 2023
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer n: the value to be tested.
#
#  Output:
#
#    logical value: true if n is a prime.
#
#    integer d: the number of times the mod function was used.
#
  if ( n != int ( n ) ):
    raise Exception ( 'is_prime1(): input n is not an integer!' )

  d = 0

  if ( n <= 0 ):
    value = False
    return value, d

  if ( n == 1 ):
    value = False
    return value, d

  value = True
  for i in range ( 2, n ):
    d = d + 1
    if ( ( n % i ) == 0 ):
      value = False

  return value, d

def is_prime2 ( n ):

#*****************************************************************************80
#
## is_prime2() reports whether an integer is prime.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 January 2023
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer n: the value to be tested.
#
#  Output:
#
#    logical value: true if n is a prime.
#
#    integer d: the number of times the mod function was used.
#
  d = 0

  if ( n != int ( n ) ):
    raise Exception ( 'is_prime2(): input n is not an integer!' )

  if ( n <= 0 ):
    value = False
    return value, d

  if ( n == 1 ):
    value = False
    return value, d

  value = True
  for i in range ( 2, n ):
    d = d + 1
    if ( ( n % i ) == 0 ):
      value = False
      break

  return value, d

def is_prime3 ( n ):

#*****************************************************************************80
#
## is_prime3() reports whether an integer is prime.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 January 2023
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer n: the value to be tested.
#
#  Output:
#
#    logical value: true if n is a prime.
#
#    integer d: the number of times the mod function was used.
#
  from math import sqrt

  d = 0

  if ( n != int ( n ) ):
    raise Exception ( 'is_prime3(): input n is not an integer!' )

  if ( n <= 0 ):
    value = False
    return value, d

  if ( n == 1 ):
    value = False
    return value, d

  n_sqrt = int ( sqrt ( n ) )

  for i in range ( 2, n_sqrt + 1 ):
    d = d + 1
    if ( ( n % i ) == 0 ):
      value = False
      return value, d

  value = True

  return value, d

def is_prime4 ( n ):

#*****************************************************************************80
#
## is_prime4() reports whether an integer is prime.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 January 2023
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer n: the value to be tested.
#
#  Output:
#
#    logical value: true if n is a prime.
#
#    integer d: the number of times the mod function was used.
#
  from math import sqrt

  d = 0

  if ( n != int ( n ) ):
    raise Exception ( 'is_prime4(): input n is not an integer!' )

  if ( n <= 1 ):
    value = False
    return value, d

  if ( ( n == 2 ) or ( n == 3 ) ):
    value = True
    return value, d

  d = d + 1
  if ( ( n % 2 ) == 0 ):
    value = False
    return value, d

  d = d + 1
  if ( ( n % 3 ) == 0 ):
    value = False
    return value, d

  sqrt_n = int ( sqrt ( n ) )

  for i in range ( 5, sqrt_n + 1, 6 ):
    d = d + 1
    if ( ( n % i ) == 0 ):
      value = False
      return value, d
    d = d + 1
    if ( ( n % ( i + 2 ) ) == 0 ):
      value = False
      return value, d

  value = True

  return value, d
  
def problem0 ( ):

#*****************************************************************************80
#
## problem0(): Which numbers are prime?
#
  from sympy import isprime

  print ( '' )
  print ( 'problem0():' )
  print ( '  Which numbers are prime?' )
  print ( '' )
  
  n_list = [ 31, 331, 3331, 33331, 333331, 3333331, 33333331, 333333331 ]
  
  for n in n_list:
    print ( '  n = ', n, ' is prime? = ', isprime ( n ) )
    
  return

def problem1 ( ):

#*****************************************************************************80
#
## problem1(): Is 1,000,009 prime.
#
  from sympy import isprime
  
  print ( '' )
  print ( 'problem1():' )
  print ( '  Is 1,000,009 prime?' )
  
  n = 1000009
  print ( '' )
  print ( '  n = ', n, ' is prime? = ', isprime ( n ) )
  
  return
  
def problem2 ( ):

#*****************************************************************************80
#
## problem2(): Euler's formula for primes.
#
  from sympy import isprime
  
  print ( '' )
  print ( 'problem2():' )
  print ( '  Check Euler\'s formula for primes.' )
  
  n = 0
  fn = n**2 + n + 41
  while ( isprime ( fn ) ):
    print ( '  n = ', n, '  f(n) =', fn, ' is prime.' )
    n = n + 1 
    fn = n**2 + n + 41

  print ( '' )
  print ( '  Euler\'s formula fails for n = ', n )
  print ( '  fn = ', fn, ' is not prime.' )

  return
  
def problem3 ( ):

#*****************************************************************************80
#
## problem3(): count calls to mod for is_prime tests.
#  
  print ( '' )
  print ( 'problem3():' )
  print ( '  Count calls to mod function for prime tests' )
  print ( '' )
  
  n = 27644437

  value, d = is_prime1 ( n )
  print ( '  is_prime1(', n, ') = ', value, '  calls = ', d )
  value, d = is_prime2 ( n )
  print ( '  is_prime2(', n, ') = ', value, '  calls = ', d )
  value, d = is_prime3 ( n )
  print ( '  is_prime3(', n, ') = ', value, '  calls = ', d )
  value, d = is_prime4 ( n )
  print ( '  is_prime4(', n, ') = ', value, '  calls = ', d )
  
  return
  
def problem4 ( ):

#*****************************************************************************80
#
## problem4(): prime factorization of n.
#
  from math import sqrt
  from sympy import isprime
  
  print ( '' )
  print ( 'problem4():' )
  print ( '  prime factorization of n' )
 
  for n in [ 1120, 2023, 314159265 ]:
  
    print ( '' )
    print ( '  Factors of ', n )
    
    factors = []
    n_hi = int ( sqrt ( n ) ) + 1
  
    for d in range ( 2, n_hi ):
      if ( isprime ( d ) ):
        while ( n % d == 0 ):
          factors.append ( d )
          n = n // d

    print ( factors )

  return

def problem5 ( ):

#*****************************************************************************80
#
## problem5(): prime_pi()
# 
  print ( '' )
  print ( 'problem5():' )
  print ( '  prime_pi(n): number of primes up to n.' )
  print ( '' )
  for n in [ 64, 256, 1024 ]:
    value = prime_pi ( n )
    print ( '  prime_pi(', n, ') = ', value )
    
  return
  
def prime_pi ( n ):

#*****************************************************************************80
#
## prime_pi() counts primes up to n.
#
  from sympy import isprime
  
  value = 0
  for d in range ( 2, n + 1 ):
    if ( isprime ( d ) ):
      value = value + 1

  return value

def problem6 ( ):

#*****************************************************************************80
#
## problem6(): Estimating prime_pi()
#
  from math import log
  
  print ( '' )
  print ( 'problem6():' )
  print ( '  Compare prime_pi(n) and n/logn.' )
 
  for k in range ( 8, 13 ):
    n = 2**k
    value = prime_pi ( n )
    gauss = n / log ( n )
    print ( '  n = ', n, '  prime_pi(n) = ', value, '  gauss = ', gauss )
    
  return
  
def problem7 ( ):

#*****************************************************************************80
#
## problem7(): sums distinct divisors
#  
  print ( '' )
  print ( 'problem7():' )
  print ( '  sigma(n): sum the divisors of n.' )
  print ( '' )
  for n in [ 617, 816, 1000, 1024 ]:
    value = sigma ( n )
    print ( '  sigma(', n, ') = ', value )
    
  return
  
def sigma ( n ):

#*****************************************************************************80
#
## sigma() counts divisors.
# 
  value = 0
  for d in range ( 1, n + 1 ):
    if ( n % d == 0 ):
      value = value + d
      
  return value
  
def problem8 ( ):

#*****************************************************************************80
#
## problem8(): count divisors.
#  
  print ( '' )
  print ( 'problem8():' )
  print ( '  tau(n): count the divisors of n.' )
  print ( '' )
  for n in [ 512, 610, 832, 960 ]:
    value = tau ( n )
    print ( '  tau(', n, ') = ', value )
  
  return

def tau ( n ):

#*****************************************************************************80
#
## tau() counts divisors.
# 
  value = 0
  for d in range ( 1, n + 1 ):
    if ( n % d == 0 ):
      value = value + 1
      
  return value
  
def problem9 ( ):

#*****************************************************************************80
#
## problem9(): Greatest prime factor.
#
  print ( '' )
  print ( 'problem9():' )
  print ( '  gpf(n): greatest prime factor of n.' )
  print ( '' )
  for n in [ 25, 698, 751, 364, 526 ]:
    value = gpf ( n )
    print ( '  gpf (', n, ') = ', value )
  
def gpf ( n ):

#*****************************************************************************80
#
## gpf() greatest prime factor of n.
# 
  from sympy import isprime
  
  value = 1
  for d in range ( 1, n + 1 ):
    if ( isprime ( d ) ):
      if ( n % d == 0 ):
        value = d
      
  return value

def problem10 ( ):

#*****************************************************************************80
#
## problem10(): multiply big values
# 
  print ( '' )
  print ( 'problem10():' )
  print ( '  Multiply big numbers' )
  
  p1 =    3_490_529_510_847_650_949_147_849_619_903_898_133_417_764_638_493_387_843_990_820_577
            
  p2 =   32_769_132_993_266_709_549_961_988_190_834_461_413_177_642_967_992_942_539_798_288_533
            
  key = 114_381_625_757_888_867_669_235_779_976_146_612_010_218_296_721_242_362_562_561_842_935_706_935_245_733_897_830_597_123_563_958_705_058_989_075_147_599_290_026_879_543_541
               
  diff = key - p1 * p2
  print ( '  diff = ', diff )
  
if ( __name__ == "__main__" ):
  problem0 ( )
  problem1 ( )
  problem2 ( )
  problem3 ( )
  problem4 ( )
  problem5 ( )
  problem6 ( )
  problem7 ( )
  problem8 ( )
  problem9 ( )
  problem10 ( )