Problem Statement

A Mersenne number is a number in the form of 2P-1. If P is prime, the Mersenne number may be a Mersenne prime (if P is not prime, the Mersenne number is also not prime). In the search for Mersenne prime numbers it is advantageous to eliminate exponents by finding a small factor before starting a, potentially lengthy, Lucas-Lehmer test. There are very efficient algorithms for determining if a number divides 2P-1 (or equivalently, if 2P mod (the number) = 1). Some languages already have built-in implementations of this exponent-and-mod operation (called modPow or similar). The following is how to implement this modPow yourself: For example, let's compute 223 mod 47. Convert the exponent 23 to binary, you get 10111. Starting with square = 1, repeatedly square it. Remove the top bit of the exponent, and if it's 1 multiply square by the base of the exponentiation (2), then compute square modulo 47. Use the result of the modulo from the last step as the initial value of square in the next step: Since 223 mod 47 = 1, 47 is a factor of 2P-1. (To see this, subtract 1 from both sides: 223-1 = 0 mod 47.) Since we've shown that 47 is a factor, 223-1 is not prime. Further properties of Mersenne numbers allow us to refine the process even more. Any factor q of 2P-1 must be of the form 2kP+1, k being a positive integer or zero. Furthermore, q must be 1 or 7 mod 8. Finally any potential factor q must be prime. As in other trial division algorithms, the algorithm stops when 2kP+1 > sqrt(N). These primality tests only work on Mersenne numbers where P is prime. For example, M4=15 yields no factors using these techniques, but factors into 3 and 5, neither of which fit 2kP+1.

Using the above method find a factor of 2929-1 (aka M929)

Let's start with the solution:

$ include "seed7_05.s7i";

const func boolean: isPrime (in integer: number) is func
  result
    var boolean: prime is FALSE;
  local
    var integer: upTo is 0;
    var integer: testNum is 3;
  begin
    if number = 2 then
      prime := TRUE;
    elsif odd(number) and number > 2 then
      upTo := sqrt(number);
      while number rem testNum <> 0 and testNum <= upTo do
        testNum +:= 2;
      end while;
      prime := testNum > upTo;
    end if;
  end func;

const func integer: modPow (in var integer: base,
    in var integer: exponent, in integer: modulus) is func
  result
    var integer: power is 1;
  begin
    if exponent < 0 or modulus < 0 then
      raise RANGE_ERROR;
    else
      while exponent > 0 do
        if odd(exponent) then
          power := (power * base) mod modulus;
        end if;
        exponent >>:= 1;
        base := base ** 2 mod modulus;
      end while;
    end if;
  end func;

const func integer: mersenneFactor (in integer: exponent) is func
  result
    var integer: factor is 0;
  local 
    var integer: maxk is 0;
    var integer: k is 1;
    var boolean: searching is TRUE;
  begin
    maxk := 16384 div exponent; # Limit for k to prevent overflow of 32 bit signed integer
    while k <= maxk and searching do
      factor := 2 * exponent * k + 1;
      if (factor mod 8 = 1 or factor mod 8 = 7) and
          isPrime(factor) and modPow(2, exponent, factor) = 1 then
        searching := FALSE;
      end if;
      incr(k);
    end while;
    if searching then
      factor := 0;
    end if;
  end func;

const proc: main is func
  begin
    writeln("Factor of M929: " <& mersenneFactor(929));
  end func;