//
// Euclid-Mullin Sequence
// Using FutureBasic 7.0.35
//
// September 2025, R.W.
//

include "gmp.incl"

_TERMS = 16   // first 16 Euclid–Mullin Sequence

// Return a non-trivial factor of n using Pollard-Rho

void local fn PollardRho( n as mpz_t, out as mpz_t )
  mpz_t x,y,c,d,absDiff
  mpz_inits( x, y, c, d, absDiff, 0 )
  mpz_set_ui( x, 2 ) : mpz_set_ui( y, 2 ) : mpz_set_ui( c, 1 )
  mpz_set_ui( d, 1 )
  while fn mpz_cmp_ui( d, 1 ) == 0
    // x = (x*x + c) mod n
    mpz_mul( x, x, x ) : mpz_add( x, x, c ) : mpz_mod( x, x, n )
    // y = f(f(y))
    mpz_mul( y, y, y ) : mpz_add( y, y, c ) : mpz_mod( y, y, n )
    mpz_mul( y, y, y ) : mpz_add( y, y, c ) : mpz_mod( y, y, n )
    // d = gcd(|x-y|, n)
    mpz_sub( absDiff, x, y ) : mpz_abs( absDiff, absDiff )
    fn mpz_gcd( d, absDiff, n )
  wend
  mpz_set( out, d )
  mpz_clears( x,y,c,d,absDiff,0 )
end fn

void local fn EuclidMullin( terms as long )
  CFStringRef result
  mpz_t total, tmp, factor
  mpz_inits( total, tmp, factor, 0 )
  mpz_set_ui( total, 1 )
  print @
  print @"First ";terms;" numbers of Euclid-Mullin sequence"
  print @
  long k
  for k = 1 to _TERMS
    mpz_add_ui( tmp, total, 1 )        // tmp = total + 1
    fn PollardRho( tmp, factor )       // factor = some factor of tmp
    result = fn mpz_cf2( factor )
    print result
    mpz_mul( total, total, factor )    // total *= factor
  next

  mpz_clears( total, tmp, factor, 0 )
end fn

window 1, @"Euclid-Mullin Sequence"

fn EuclidMullin(_TERMS )

handleEvents