Seed7 - The extensible programming language
Seed7 FAQ Manual Screenshots Examples Libraries Algorithms Download Links
Algorithms Sorting Searching Date & Time String Float Mathematics Message digest Graphics File Puzzles Others
Algorithms
Mathematics
 previous   up   next 

Determine wether a given integer number is prime

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;

Sieve of Eratosthenes

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;

Verify Goldbach's conjecture up to 10000000

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;

Function to calculate the prime factors of a number

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;

Verify if a given Mersenne number is prime

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;

Miller-Rabin primality test

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;

Determine the greatest common divisor of two positive integer numbers

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;

Determine the least common multiple of two positive integer numbers

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

Binary greatest common divisor algorithm for two positive integer numbers

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;

Binary greatest common divisor algorithm for two bigInteger numbers

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;

Return the modular multiplicative inverse of a modulo b

The library "bigint.s7i" defines the bigInteger function modInverse. 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: 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;

Modular exponentiation with the binary method

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;

Compute PI with 1000 decimal digits using John Machin's formula

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;

Compute PI with 1000 decimal digits using the Bailey/Borwein/Plouffe formula

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;

Write PI with 1000 decimal digits using Newtons formula

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

Write the decimal digits of PI with a spigot algorithm

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;

Determine the truncated square root of a big integer number

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;

Function to compute the nth root of a positive float number

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;

Exponentiation function for integers

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;

Exponentiation function for float base and integer exponent

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;

Multiply integers using the peasant multiplication

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;

Binomial coefficient

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;

Gamma function

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;

Matrix addition

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;

Matrix multiplication

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;

Matrix exponentiation with integer exponent

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;

Convert a matrix to row echelon form

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;

Convert a matrix to reduced row echelon form

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;

Dot product

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

Random number generator

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;

Recursive fibonacci function

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;

Iterative fibonacci function

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;

Display 100 numbers of the fibonacci sequence

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

The tak function

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

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;

 previous   up   next