const func boolean: is_prime (in integer: number) is func
  result
    var boolean: result is FALSE;
  local
    var integer: count is 2;
  begin
    if number = 2 then
      result := TRUE;
    elsif number > 2 then
      while number rem count <> 0 and count * count <= number do
        incr(count);
      end while;
      result := number rem count <> 0;
    end if;
  end func;

const func array boolean: eratosthenes (in integer: n) is func
  result
    var array boolean: sieve is 0 times FALSE;
  local
    var integer: i is 0;
    var integer: j is 0;
  begin 
    sieve := n times TRUE;
    sieve[1] := FALSE;
    for i range 2 to sqrt(n) do
      if sieve[i] then
        for j range i * i to n step i do
          sieve[j] := FALSE;
        end for;
      end if;
    end for;
  end func;

const func array integer: factorise (in var integer: number) is func
  result
    var array integer: result is 0 times 0;
  local
    var integer: checker is 2;
  begin
    while checker * checker <= number do
      if number rem checker = 0 then
        result &:= [](checker);
        number := number div checker;
      else
        incr(checker);
      end if;
    end while;
    if number <> 1 then
      result &:= [](number);
    end if;
  end func;

A Mersenne number is a number that is one less than a power of two (2**p-1). The Lucas-Lehmer test allows to check if a given Mersenne number is prime: For the prime p, the Mersenne number 2**p-1 is prime if and only if S(p-1) rem 2**p-1 = 0 where S(1)=4 and S(n)=S(n-1)**2-2.

const func boolean: mersenne_number_is_prime (in integer: p) is func
  result
    var boolean: result is FALSE;
  local
    var bigInteger: m_p is 0_;
    var bigInteger: s is 4_;
    var integer: i is 0;
  begin
    m_p := 2_ ** p - 1_;
    for i range 2 to pred(p) do
      s := (s ** 2 - 2_) rem m_p;
    end for;
    result := s = 0_;
  end func;

const func integer: gcd (in var integer: a, in var integer: b) is func
  result
    var integer: result is 0;
  local
    var integer: help is 0;
  begin
    while a <> 0 do
      help := b rem a;
      b := a;
      a := help;
    end while;
    result := b;
  end func;

const func integer: lcm (in integer: a, in integer: b) is
  return a * b div gcd(a, b);

const func integer: binaryGcd (in var integer: a, in var integer: b) is func
  result
    var integer: result is 0;
  local
    var integer: shift is 0;
    var integer: diff is 0;
  begin
    if a = 0 then
      result := b;
    elsif b = 0 then
      result := a;
    else
      (* Let shift := log2(K), where K is the greatest power of 2
         dividing both a and b. *)
      while not odd(a) and not odd(b) do
          a >>:= 1;
          b >>:= 1;
        incr(shift);
      end while;

      while not odd(a) do
        a >>:= 1;
      end while;

      (* From here on, a is always odd. *)
      repeat
        while not odd(b) do
          b >>:= 1;
        end while;

        (* Now a and b are both odd, so diff(a, b) is even.
           Let a := min(a, b); b := diff(a, b)/2; *)
        if a < b then
          b -:= a;
        else
          diff := a - b;
          a := b;
          b := diff;
        end if;
        b >>:= 1;
      until b = 0;

      result := a << shift;
    end if;
  end func;

The bigInteger type already has a gcd function defined which can be called with 'gcd(a, b)'.

const func bigInteger: binaryGcd (in var bigInteger: a, in var bigInteger: b) is func
  result
    var bigInteger: result is 0_;
  local
    var integer: lowestSetBitA is 0;
    var integer: lowestSetBitB is 0;
    var integer: shift is 0;
    var bigInteger: diff is 0_;
  begin
    if a = 0_ then
      result := b;
    elsif b = 0_ then
      result := a;
    else
      if a < 0_ then
        a := -a;
      end if;
      if b < 0_ then
        b := -b;
      end if;
      lowestSetBitA := lowestSetBit(a);
      lowestSetBitB := lowestSetBit(b);
      if lowestSetBitA < lowestSetBitB then
        shift := lowestSetBitA;
      else
        shift := lowestSetBitB;
      end if;
      a >>:= lowestSetBitA;
      repeat
        b >>:= lowestSetBit(b);
        if a < b then
          b -:= a;
        else
          diff := a - b;
          a := b;
          b := diff;
        end if;
      until b = 0_;
      result := a << shift;
    end if;
  end func;

This function is part of the "bigint.s7i" library. It returns the modular multiplicative inverse of a modulo b when a and b are coprime (gcd(a, b) <> 1). If a and b are not coprime (gcd(a, b) <> 1) the exception RANGE_ERROR is raised.

const func bigInteger: modInverse (in var bigInteger: a,
    in var bigInteger: b) is func
  result
    var bigInteger: result is 0_;
  local
    var bigInteger: b_bak is 0_;
    var bigInteger: x is 0_;
    var bigInteger: y is 1_;
    var bigInteger: lastx is 1_;
    var bigInteger: lasty is 0_;
    var bigInteger: temp is 0_;
    var bigInteger: quotient is 0_;
  begin
    if b < 0_ then
      raise RANGE_ERROR;
    end if;
    if a < 0_ and b <> 0_ then
      a := a mod b;
    end if;
    b_bak := b;
    while b <> 0_ do
      temp := b;
      quotient := a div b;
      b := a rem b;
      a := temp;

      temp := x;
      x := lastx - quotient * x;
      lastx := temp;

      temp := y;
      y := lasty - quotient * y;
      lasty := temp;
    end while;
    if a = 1_ then
      result := lastx;
      if result < 0_ then
        result +:= b_bak;
      end if;
    else
      raise RANGE_ERROR;
    end if;
  end func;

$ include "seed7_05.s7i";
  include "bigint.s7i";
  include "bigrat.s7i";

# John Machin's formula from 1706:
# PI = 16 * arctan(1 / 5) - 4 * arctan(1 / 239)

# Taylor series for arctan:
# arctan(x) = sum_n_from_0_to_inf((-1) ** n / succ(2 * n) * x ** succ(2 * n))

# Taylor series of John Machin's formula:
# PI = sum_n_from_0_to_inf(16 * (-1) ** n / (succ(2 * n) *   5 ** succ(2 * n)) -
#                           4 * (-1) ** n / (succ(2 * n) * 239 ** succ(2 * n)))

const func bigRational: compute_pi_machin is func
  result
    var bigRational: sum is 0_ / 1_;
  local
    var integer: n is 0;
  begin
    for n range 0 to 713 do
      sum +:= 16_ * (-1_) ** n / (bigInteger conv succ(2 * n) *   5_ ** succ(2 * n)) -
               4_ * (-1_) ** n / (bigInteger conv succ(2 * n) * 239_ ** succ(2 * n));
    end for;
  end func;

 const proc: main is func
   begin
     writeln(compute_pi_machin digits 1000);
   end func;

$ include "seed7_05.s7i";
  include "bigint.s7i";
  include "bigrat.s7i";

# In 1997, David H. Bailey, Peter Borwein and Simon Plouffe published a
# paper (Bailey, 1997) on a new formula for PI as an infinite series:

# PI = sum_n_from_0_to_inf(16 ** (-n) *
#      (4 / (8 * n + 1) - 2 / (8 * n + 4) - 1 / (8 * n + 5) - 1 / (8 * n + 6)))

const func bigRational: compute_pi_bailey_borwein_plouffe is func
  result
    var bigRational: sum is 0_ / 1_;
  local
    var integer: n is 0;
    var bigInteger: k is 0_;
  begin
    for n range 0 to 825 do
      k := bigInteger conv n;
      sum +:= 1_ / 16_ ** n *
          (4_ / (8_ * k + 1_) - 2_ / (8_ * k + 4_) - 1_ / (8_ * k + 5_) - 1_ / (8_ * k + 6_));
    end for;
  end func;

 const proc: main is func
   begin
     writeln(compute_pi_bailey_borwein_plouffe digits 1000);
   end func;

$ include "seed7_05.s7i";

# Newtons formula for PI is:
# PI / 2 = sum_n_from_0_to_inf(n! / (2 * n + 1)!!)

# This can be written as:
# PI / 2 = 1 + 1/3 * (1 + 2/5 * (1 + 3/7 * (1 + 4/9 * (1 + ... ))))

# This algorithm puts 2 * 1000 on the right side and computes everything from inside out.

const integer: SCALE is 10000;
const integer: MAXARR is 3500;
const integer: ARRINIT is 2000;

const proc: main is func
  local
    var integer: i is 0;
    var integer: j is 0;
    var integer: carry is 0;
    var array integer: arr is MAXARR times ARRINIT;
    var integer: sum is 0;
  begin
    for i range MAXARR downto 1 step 14 do
      sum := 0;
      for j range i downto 1 do
        sum := sum*j + SCALE*arr[j];
        arr[j] := sum rem pred(j*2);
        sum := sum div pred(j*2);
      end for;
      write(carry + sum div SCALE lpad0 4);
      carry := sum rem SCALE;
    end for;
    writeln;
  end func;

The bigInteger type already has an integer square root which can be called with 'sqrt(number)'. A hand written version would be:

const func bigInteger: intSqrt (in var bigInteger: number) is func
  result
    var bigInteger: result is 0_;
  local
    var bigInteger: res2 is 0_;
  begin
    if number > 0_ then
      res2 := number;
      repeat
        result := res2;
        res2 := (result + number div result) div 2_;
      until result <= res2;
    end if;
  end func;

The nth root of the number 'a' can be computed with the exponentiation operator: 'a ** (1 / n)'. An alternate function which uses Newton's method is:

const func float: nthRoot (in integer: n, in float: a) is func
  result
    var float: x1 is 0.0;
  local
    var float: x0 is 0.0;
  begin 
    x0 := a;
    x1 := a / flt(n);
    while abs(x1 - x0) >= abs(x0 * 1.0E-9) do
      x0 := x1;
      x1 := (flt(pred(n)) * x0 + a / x0 ** pred(n)) / flt(n);
    end while;
  end func;

The type 'integer' defines an exponentiation operator which can be called with 'base ** exponent'. The Seed7 runtime library contains the function 'intPow' which implements the '**' operator. This function uses the exponentiation by squaring algorithm. The Seed7 version of 'intPow' is:

const func integer: intPow (in var integer: base, in var integer: exponent) is func
  result
    var integer: result is 0;
  begin
    if exponent < 0 then
      raise(NUMERIC_ERROR);
    else
      if odd(exponent) then
        result := base;
      else
        result := 1;
      end if;
      exponent := exponent div 2;
      while exponent <> 0 do
        base *:= base;
        if odd(exponent) then
          result *:= base;
        end if;
        exponent := exponent div 2;
      end while;
    end if;
  end func;

The type 'float' defines an exponentiation operator which can be called with 'base ** exponent'. The Seed7 runtime library contains the function 'fltIPow' which implements the '**' operator. This function uses the exponentiation by squaring algorithm. The Seed7 version of 'fltIPow' is:

const func float: fltIPow (in var float: base, in var integer: exponent) is func
  result
    var float: result is 0.0;
  local
    var boolean: neg_exp is FALSE;
  begin
    if base = 0.0 then
      if exponent < 0 then
        result := Infinity;
      elsif exponent = 0 then
        result := 1.0;
      else
        result := 0.0;
      end if;
    else
      if exponent < 0 then
        exponent := -exponent;
        neg_exp := TRUE;
      end if;
      # In the twos complement representation the most
      # negative number is the only one where both the
      # number and its negation are negative. When the
      # exponent is proven to be to the most negative
      # number fltIPow returns 0.0 .
      if exponent >= 0 then
        if odd(exponent) then
          result := base;
        else
          result := 1.0;
        end if;
        exponent >>:= 1;
        while exponent <> 0 do
          base *:= base;
          if odd(exponent) then
            result *:= base;
          end if;
          exponent >>:= 1;
        end while;
        if neg_exp then
          result := 1.0 / result;
        end if;
      end if;
    end if;
  end func;

const func integer: peasantMult (in var integer: a, in var integer: b) is func
  result
    var integer: result is 0;
  begin
    while a <> 0 do
      if odd(a) then
        result +:= b;
      end if;
      a := a div 2;
      b *:= 2;
    end while;
  end func;

const type: matrix is array array float;

const func matrix: (in matrix: left) + (in matrix: right) is func
  result
    var matrix: result is matrix.value;
  local
    var integer: i is 0;
    var integer: j is 0;
  begin
    if length(left) <> length(right) or
        length(left[1]) <> length(right[1]) then
      raise RANGE_ERROR;
    else
      result := length(left) times length(left[1]) times 0.0;
      for i range 1 to length(left) do
        for j range 1 to length(left[i]) do
          result[i][j] := left[i][j] + right[i][j];
        end for;
      end for;
    end if;
  end func;

const type: matrix is array array float;

const func matrix: (in matrix: left) * (in matrix: right) is func
  result
    var matrix: result is matrix.value;
  local
    var integer: i is 0;
    var integer: j is 0;
    var integer: k is 0;
    var float: accumulator is 0.0;
  begin
    if length(left[1]) <> length(right) then
      raise RANGE_ERROR;
    else
      result := length(left) times length(right[1]) times 0.0;
      for i range 1 to length(left) do
        for j range 1 to length(right) do
          accumulator := 0.0;
          for k range 1 to length(left) do
            accumulator +:= left[i][k] * right[k][j];
          end for;
          result[i][j] := accumulator;
        end for;
      end for;
    end if;
  end func;

This random generator delivers random numbers between 0 and 2**32-1. To get different pseudorandom sequence seed can be initialized with another value.

var bigInteger: seed is 987654321_;

const func bigInteger: rand_32 is func
  result
    var integer: result is 0;
  begin
    seed := (seed * 1073741825_ + 1234567891_) rem 2_ ** 64;
    result := seed div 2_ ** 32;
  end func;

This algorithm just uses the recursive definition of the fibonacci function. It can be used as benchmark to test the performance of recursive calls. For a solution with better scalability use the iterative fibonacci function below.

const func integer: fib (in integer: number) is func
  result
    var integer: result is 1;
  begin
    if number > 2 then
      result := fib(pred(number)) + fib(number - 2);
    elsif number = 0 then
      result := 0;
    end if;
  end func;

const func bigInteger: fib (in integer: number) is func
  result
    var bigInteger: result is 1_;
  local
    var integer: i is 0;
    var bigInteger: a is 0_;
    var bigInteger: c is 0_;
  begin
    for i range 1 to pred(number) do
      c := a;
      a := result;
      result +:= c;
    end for;
  end func;

$ include "seed7_05.s7i";
  include "bigint.s7i";
 
const proc: main is func
  local
    var integer: i is 0;
    var bigInteger: a is 0_;
    var bigInteger: b is 1_;
    var bigInteger: c is 0_;
  begin
    writeln(a);
    for i range 1 to 100 do
      writeln(b);
      c := a;
      a := b;
      b +:= c;
    end for;
  end func;

const func integer: tak (in integer: x, in integer: y, in integer: z) is func
  result
    var integer: result is 0;
  begin
    if y >= x then
      result := z;
    else
      result := tak(tak(pred(x), y, z),
                    tak(pred(y), z, x),
                    tak(pred(z), x, y));
    end if;
  end func;

const func integer: ackermann (in integer: m, in integer: n) is func
  result
    var integer: result is 0;
  begin
    if m = 0 then
      result := succ(n);
    elsif n = 0 then
      result := ackermann(pred(m), 1);
    else
      result := ackermann(pred(m), ackermann(m, pred(n)));
    end if;
  end func;