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;

The program below computes the number of primes between 1 and 10000000:

$ include "seed7_05.s7i";

const func set of integer: eratosthenes (in integer: n) is func
  result
    var set of integer: sieve is EMPTY_SET;
  local
    var integer: i is 0;
    var integer: j is 0;
  begin
    sieve := {2 .. n};
    for i range 2 to sqrt(n) do
      if i in sieve then
        for j range i ** 2 to n step i do
          excl(sieve, j);
        end for;
      end if;
    end for;
  end func;

const proc: main is func
  begin
    writeln(card(eratosthenes(10000000)));
  end func;

Goldbach's conjecture states that all even integers greater than 2 can be expressed as the sum of two primes. Goldbach's conjecture is one of the oldest unsolved problems. The program below checks all even numbers between 4 and 10000000:

$ include "seed7_05.s7i";

const func set of integer: eratosthenes (in integer: n) is func
  result
    var set of integer: sieve is EMPTY_SET;
  local
    var integer: i is 0;
    var integer: j is 0;
  begin
    sieve := {2 .. n};
    for i range 2 to sqrt(n) do
      if i in sieve then
        for j range i ** 2 to n step i do
          excl(sieve, j);
        end for;
      end if;
    end for;
  end func;

const integer: limit is 10000000;

const func boolean: isGoldbachNumber (in integer: number) is func
  result
    var boolean: isGoldbachNumber is FALSE;
  local
    const set of integer: primes is eratosthenes(limit);
    var integer: testNum is 0;
  begin
    for testNum range primes until isGoldbachNumber do
      isGoldbachNumber := (number - testNum) in primes;
    end for;
  end func;

const proc: main is func
  local
    var integer: number is 0;
  begin
    for number range 4 to limit step 2 do
      if not isGoldbachNumber(number) then
        writeln(number <& " is not a Goldbach number");
      end if;
      if number rem 16384 = 0 then
        write(".");
        flush(OUT);
      end if;
    end for;
    writeln;
  end func;

const func array integer: factorise (in var integer: number) is func
  result
    var array integer: factors is 0 times 0;
  local
    var integer: checker is 2;
  begin
    while checker * checker <= number do
      if number rem checker = 0 then
        factors &:= [](checker);
        number := number div checker;
      else
        incr(checker);
      end if;
    end while;
    if number <> 1 then
      factors &:= [](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 odd 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: lucasLehmerTest (in integer: p) is func
  result
    var boolean: prime is TRUE;
  local
    var bigInteger: m_p is 0_;
    var bigInteger: s is 4_;
    var integer: i is 0;
  begin
    if p <> 2 then
      m_p := 2_ ** p - 1_;
      for i range 2 to pred(p) do
        s := (s ** 2 - 2_) rem m_p;
      end for;
      prime := s = 0_;
    end if;
  end func;

The Miller-Rabin primality test is a probabilistic algorithm, which determines whether a given number is prime or not. The number of tests is specified with the parameter k.

const func boolean: millerRabin (in bigInteger: n, in integer: k) is func
  result
    var boolean: probablyPrime is TRUE;
  local
    var bigInteger: d is 0_;
    var integer: r is 0;
    var integer: s is 0;
    var bigInteger: a is 0_;
    var bigInteger: x is 0_;
    var integer: tests is 0;
  begin
    if n < 2_ or (n > 2_ and not odd(n)) then
      probablyPrime := FALSE;
    elsif n > 3_ then
      d := pred(n);
      s := lowestSetBit(d);
      d >>:= s;
      while tests < k and probablyPrime do
        a := rand(2_, pred(n));
        x := modPow(a, d, n);
        if x <> 1_ and x <> pred(n) then
          r := 1;
          while r < s and x <> 1_ and x <> pred(n) do
            x := modPow(x, 2_, n);
            incr(r);
          end while;
          probablyPrime := x = pred(n);
        end if;
        incr(tests);
      end while;
    end if;
  end func;

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

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

const func integer: binaryGcd (in var integer: a, in var integer: b) is func
  result
    var integer: gcd is 0;
  local
    var integer: shift is 0;
    var integer: diff is 0;
  begin
    if a = 0 then
      gcd := b;
    elsif b = 0 then
      gcd := 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;

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

The library "bigint.s7i" defines the bigInteger function gcd. The function 'binaryGcd' below uses the binary greatest common divisor algorithm:

const func bigInteger: binaryGcd (in var bigInteger: a, in var bigInteger: b) is func
  result
    var bigInteger: gcd is 0_;
  local
    var integer: lowestSetBitA is 0;
    var integer: shift is 0;
    var bigInteger: diff is 0_;
  begin
    if a = 0_ then
      gcd := b;
    elsif b = 0_ then
      gcd := a;
    else
      if a < 0_ then
        a := -a;
      end if;
      if b < 0_ then
        b := -b;
      end if;
      lowestSetBitA := lowestSetBit(a);
      shift := lowestSetBit(b);
      if lowestSetBitA < shift then
        shift := lowestSetBitA;
      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_;
      gcd := a << shift;
    end if;
  end func;

The library "bigint.s7i" defines the bigInteger function modInverse. It returns the modular multiplicative inverse of a modulo b if 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: inverse 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
      inverse := lastx;
      if inverse < 0_ then
        inverse +:= b_bak;
      end if;
    else
      raise RANGE_ERROR;
    end if;
  end func;

The function modPow is part of the "bigint.s7i" library.

const func bigInteger: modPow (in var bigInteger: base,
    in var bigInteger: exponent, in bigInteger: modulus) is func
  result
    var bigInteger: 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;

The function 'compute_pi_machin' approximates PI as bigRational number. To compute PI with more digits the upper bound of the for-loop needs to be adjusted, together with the right operand of the digits operator. E.g.: For 2000 digits the upper bound of the for-loop must be 1429 instead of 713.

$ 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;

The function 'compute_pi_bailey_borwein_plouffe' approximates PI as bigRational number. To compute PI with more digits the upper bound of the for-loop needs to be adjusted, together with the right operand of the digits operator. E.g.: For 2000 digits the upper bound of the for-loop must be 1655 instead of 825.

$ 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: k8 is 0_;
  begin
    for n range 0 to 825 do
      k8 := bigInteger conv (8 * n);
      sum +:= 1_ / 16_ ** n *
          (4_ / (k8 + 1_) - 2_ / (k8 + 4_) - 1_ / (k8 + 5_) - 1_ / (k8 + 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 program below uses the unbounded spigot algorithm published by Jeremy Gibbons. The digits of PI are calculated and written in succession. To stop the program press ctrl-C followed by * and ENTER.

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

const proc: main is func
  local
    var bigInteger: q is 1_;
    var bigInteger: r is 0_;
    var bigInteger: t is 1_;
    var bigInteger: k is 1_;
    var bigInteger: n is 3_;
    var bigInteger: l is 3_;
    var bigInteger: nn is 0_;
    var bigInteger: nr is 0_;
    var boolean: first is TRUE;
  begin
    while TRUE do
      if 4_ * q + r - t < n * t then
        write(n);
        if first then
          write(".");
          first := FALSE;
        end if;
        nr := 10_ * (r - n * t);
        n := 10_ * (3_ * q + r) div t - 10_ * n;
        q *:= 10_;
        r := nr;
        flush(OUT);
      else
        nr := (2_ * q + r) * l;
        nn := (q * (7_ * k + 2_) + r * l) div (t * l);
        q *:= k;
        t *:= l;
        l +:= 2_;
        incr(k);
        n := nn;
        r := nr;
      end if;
    end while;
  end func;

The library "bigint.s7i" defines the bigInteger function sqrt.

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

The nth root of the number 'a' can be computed with the exponentiation operator: 'a ** (1.0 / flt(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: power is 0;
  begin
    if exponent < 0 then
      raise(NUMERIC_ERROR);
    else
      if odd(exponent) then
        power := base;
      else
        power := 1;
      end if;
      exponent := exponent div 2;
      while exponent <> 0 do
        base *:= base;
        if odd(exponent) then
          power *:= 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: power is 1.0;
  local
    var integer: stop is 0;
  begin
    if base = 0.0 then
      if exponent < 0 then
        power := Infinity;
      elsif exponent > 0 then
        power := 0.0;
      end if;
    else
      if exponent < 0 then
        stop := -1;
      end if;
      if odd(exponent) then
        power := base;
      end if;
      exponent >>:= 1;
      while exponent <> stop do
        base *:= base;
        if odd(exponent) then
          power *:= base;
        end if;
        exponent >>:= 1;
      end while;
      if stop = -1 then
        power := 1.0 / power;
      end if;
    end if;
  end func;

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

The type integer defines the infix operator ! to compute the binomial coefficient. Below is the definition of the bigInteger infix operator ! from the library "bigint.s7i".

const func bigInteger: (in bigInteger: n) ! (in var bigInteger: k) is func
  result
    var bigInteger: binom is 0_;
  local
    var bigInteger: numerator is 0_;
    var bigInteger: denominator is 0_;
  begin
    if n >= 0_ and k > n >> 1 then
      k := n - k;
    end if;
    if k < 0_ then
      binom := 0_;
    elsif k = 0_ then
      binom := 1_;
    else
      binom := n;
      numerator := pred(n);
      denominator := 2_;
      while denominator <= k do
        binom *:= numerator;
        binom := binom div denominator;
        decr(numerator);
        incr(denominator);
      end while;
    end if;
  end func;

const func float: gamma (in float: X) is func
  result
    var float: gamma is 0.0;
  local
    const array float: A is [] (
         1.00000000000000000000,  0.57721566490153286061,
        -0.65587807152025388108, -0.04200263503409523553,
         0.16653861138229148950, -0.04219773455554433675,
        -0.00962197152787697356,  0.00721894324666309954,
        -0.00116516759185906511, -0.00021524167411495097,
         0.00012805028238811619, -0.00002013485478078824,
        -0.00000125049348214267,  0.00000113302723198170,
        -0.00000020563384169776,  0.00000000611609510448,
         0.00000000500200764447, -0.00000000118127457049,
         0.00000000010434267117,  0.00000000000778226344,
        -0.00000000000369680562,  0.00000000000051003703,
        -0.00000000000002058326, -0.00000000000000534812,
         0.00000000000000122678, -0.00000000000000011813,
         0.00000000000000000119,  0.00000000000000000141,
        -0.00000000000000000023,  0.00000000000000000002);
    var float: Y is 0.0;
    var float: Sum is A[maxIdx(A)];
    var integer: N is 0;
  begin
    Y := X - 1.0;
    for N range pred(maxIdx(A)) downto minIdx(A) do
      Sum := Sum * Y + A[N];
    end for;
    gamma := 1.0 / Sum;
  end func;

const type: matrix is array array float;

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

const type: matrix is array array float;

const func matrix: (in matrix: factor1) * (in matrix: factor2) is func
  result
    var matrix: product 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(factor1[1]) <> length(factor2) then
      raise RANGE_ERROR;
    else
      product := length(factor1) times length(factor2[1]) times 0.0;
      for i range 1 to length(factor1) do
        for j range 1 to length(factor2) do
          accumulator := 0.0;
          for k range 1 to length(factor1) do
            accumulator +:= factor1[i][k] * factor2[k][j];
          end for;
          product[i][j] := accumulator;
        end for;
      end for;
    end if;
  end func;

Based on the matrix multiplication and the type matrix a matrix exponentiation with an integer exponent can be defined:

const func matrix: (in var matrix: base) ** (in var integer: exponent) is func
  result
    var matrix: power is matrix.value;
  local
    var integer: row is 0;
    var integer: column is 0;
  begin
    if exponent < 0 then
      raise NUMERIC_ERROR;
    else
      if odd(exponent) then
        power := base;
      else
        # Create identity matrix
        power := length(base) times length(base) times 0.0;
        for row range 1 to length(base) do
          for column range 1 to length(base) do
            if row = column then
              power[row][column] := 1.0;
            end if;
          end for;
        end for;
      end if;
      exponent := exponent div 2;
      while exponent > 0 do
        base := base * base;
        if odd(exponent) then
          power := power * base;
        end if;
        exponent := exponent div 2;
      end while;
    end if;
  end func;

A matrix is in row echelon form if

• All rows with at least one nonzero element are above any rows with all zeroes.
• The first nonzero number from the left, of a nonzero row is always strictly to the right of the first nonzero number from the left, in the row above it.

The function below converts a matrix like

2.0		1.0		-1.0		8.0
-3.0		-1.0		2.0		-11.0
-2.0		1.0		2.0		-3.0

into its row echelon form

2.0		1.0		-1.0		8.0
0.0		0.5		0.5		1.0
0.0		0.0		-1.0		1.0

const type: matrix is array array float;

const proc: toRowEchelonForm (inout matrix: mat) is func
  local
    var integer: numRows is 0;
    var integer: numColumns is 0;
    var integer: row is 0;
    var integer: column is 0;
    var integer: pivot is 0;
    var float: factor is 0.0;
  begin
    numRows := length(mat);
    numColumns := length(mat[1]);
    for row range numRows downto 1 do
      column := 1;
      while column <= numColumns and mat[row][column] = 0.0 do
        incr(column);
      end while;
      if column > numColumns then
        # Empty rows are moved to the bottom
        mat := mat[.. pred(row)] & mat[succ(row) ..] & [] (mat[row]);
        decr(numRows);
      end if;
    end for;
    for pivot range 1 to numRows do
      if mat[pivot][pivot] = 0.0 then
        # Find a row were the pivot column is not zero
        row := 1;
        while row <= numRows and mat[row][pivot] = 0.0 do
          incr(row);
        end while;
        # Add row were the pivot column is not zero
        for column range 1 to numColumns do
          mat[pivot][column] +:= mat[row][column];
        end for;
      end if;
      for row range succ(pivot) to numRows do
        if mat[row][pivot] <> 0.0 then
          # Make sure that the pivot column contains zero in all rows below the pivot row
          factor := -mat[row][pivot] / mat[pivot][pivot];
          for column range pivot to numColumns do
            mat[row][column] +:= mat[pivot][column] * factor;
          end for;
        end if;
      end for;
    end for;
  end func;

A matrix is in reduced row echelon form if

• All rows with at least one nonzero element are above any rows with all zeroes.
• The first nonzero number from the left, of a nonzero row is always strictly to the right of the first nonzero number from the left, in the row above it.
• Every leading coefficient is 1 and is the only nonzero entry in its column

The function below converts a matrix like

2.0		1.0		-1.0		8.0
-3.0		-1.0		2.0		-11.0
-2.0		1.0		2.0		-3.0

into its reduced row echelon form

1.0		0.0		0.0		2.0
0.0		1.0		0.0		3.0
0.0		0.0		1.0		-1.0

const type: matrix is array array float;

const proc: toReducedRowEchelonForm (inout matrix: mat) is func
  local
    var integer: numRows is 0;
    var integer: numColumns is 0;
    var integer: row is 0;
    var integer: column is 0;
    var integer: pivot is 0;
    var float: factor is 0.0;
  begin
    numRows := length(mat);
    numColumns := length(mat[1]);
    for row range numRows downto 1 do
      column := 1;
      while column <= numColumns and mat[row][column] = 0.0 do
        incr(column);
      end while;
      if column > numColumns then
        # Empty rows are moved to the bottom
        mat := mat[.. pred(row)] & mat[succ(row) ..] & [] (mat[row]);
        decr(numRows);
      end if;
    end for;
    for pivot range 1 to numRows do
      if mat[pivot][pivot] = 0.0 then
        # Find a row were the pivot column is not zero
        row := 1;
        while row <= numRows and mat[row][pivot] = 0.0 do
          incr(row);
        end while;
        # Add row were the pivot column is not zero
        for column range 1 to numColumns do
          mat[pivot][column] +:= mat[row][column];
        end for;
      end if;
      if mat[pivot][pivot] <> 1.0 then
        # Make sure that the pivot element is 1.0
        factor := 1.0 / mat[pivot][pivot];
        for column range pivot to numColumns do
          mat[pivot][column] := mat[pivot][column] * factor;
        end for;
      end if;
      for row range 1 to numRows do
        if row <> pivot and mat[row][pivot] <> 0.0 then
          # Make sure that in all other rows the pivot column contains zero
          factor := -mat[row][pivot];
          for column range pivot to numColumns do
            mat[row][column] +:= mat[pivot][column] * factor;
          end for;
        end if;
      end for;
    end for;
  end func;

$ include "seed7_05.s7i";

$ syntax expr: .().dot.() is  -> 6;  # priority of dot operator

const func integer: (in array integer: a) dot (in array integer: b) is func
  result
    var integer: sum is 0;
  local
    var integer: index is 0;
  begin
    if length(a) <> length(b) then
      raise RANGE_ERROR;
    else
      for index range 1 to length(a) do
        sum +:= a[index] * b[index];
      end for;
    end if;
  end func;

const proc: main is func
  begin
    writeln([](1, 3, -5) dot [](4, -2, -1));
  end func;

The function 'rand', which creates a random value between a lower and an upper bound, is predefined for many types (boolean, integer, bigInteger, float, char). The random generator below delivers random numbers between 0 and 2**32-1. To get a different pseudorandom sequence the variable 'seed' can be initialized with another value.

var bigInteger: seed is 987654321_;

const func bigInteger: rand_32 is func
  result
    var integer: randomNumber is 0;
  begin
    seed := (seed * 1103515245_ + 12345_ rem 2_ ** 64;
    randomNumber := 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: fib is 1;
  begin
    if number > 2 then
      fib := fib(pred(number)) + fib(number - 2);
    elsif number = 0 then
      fib := 0;
    end if;
  end func;

const func bigInteger: fib (in integer: number) is func
  result
    var bigInteger: fib 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 := fib;
      fib +:= 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: tak is 0;
  begin
    if y >= x then
      tak := z;
    else
      tak := tak(tak(pred(x), y, z),
                 tak(pred(y), z, x),
                 tak(pred(z), x, y));
    end if;
  end func;

The ackermann function grows rapidly, even for small inputs:

ackermann(4, 1) is 65533

ackermann(4, 2) has 19729 decimal digits.

Below is a simple implementation of the ackermann function, that follows the original specification. The recursion depth of ackermann grows so fast that computing ackermann(4, 1), with the implementation below, might trigger a stack overflow. In practice the bigInteger ackermann function implementation further below should be preferred.

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

This implementation of the ackermann function uses bigInteger instead of integer. Additionally there are optimizations for special cases, to reduce the recursion depth. Because of these improvements the ackermann implementation below can compute ackermann(4, 2) without problems.

const func bigInteger: ackermann (in bigInteger: m, in bigInteger: n) is func
  result
    var bigInteger: ackermann is 0_;
  begin
    case m of
      when {0_}: ackermann := succ(n);
      when {1_}: ackermann := n + 2_;
      when {2_}: ackermann := 3_ + 2_ * n;
      when {3_}: ackermann := 5_ + 8_ * pred(2_ ** ord(n));
      otherwise:
        if n = 0_ then
          ackermann := ackermann(pred(m), 1_);
        else
          ackermann := ackermann(pred(m), ackermann(m, pred(n)));
        end if;
    end case;
  end func;