# -*- coding: UTF-8 -*-

# mathGUIde  Version 2.1"
# Copyright © 2004-2011 Hartmut Ring"
# http://www.math.uni-siegen.de/ring/mathGUIde/
# This file is for internal use only by the mathGUIde GUI.
#
# mathGUIde is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3 of the License, or (at your option) any later
# version.
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with
# this program; if not, see http://www.gnu.org/licenses/.""")
#............................................................


import operator, random, math, re, inspect, functools

class _MathGUIdeError(Exception):
    def __init__(self, value):
        self.value = value
    def __str__(self):
        return repr(self.value)

_pow = pow # import everything from math except pow:
from math import *
pow = _pow

# make original Python functions floor and ceil invisible:
_floor, _ceil = floor, ceil

def _cmp(a,b):
    return int(a > b) - int(a < b)

#............................................................
#begin

# Change floor and ceil so that they map from R-->Z
# instead of R-->R as given by Python.

def floor(x):
    """ <b>floor(x)</b><br/>
        <i>highest integer n &lt;= x</i><br/>
    ----de----
        <b>floor(x)</b><br/>
        <i>größte ganze Zahl n &lt;= x (Gauß-Klammer)</i><br/>
    """
    return int(_floor(x))

def ceil(x):
    """ <b>floor(x)</b><br/>
        <i>smallest integer n >= x</i><br/>
    ----de----
        <b>ceil(x)</b><br/>
        <i>kleinste ganze Zahl n >= x</i><br/>
    """
    return int(_ceil(x))

#............................................................

def copy(x):
    """ <b>copy(x)</b><br/>
        <i>independent copy ("clone") of x</i>
    ----de----
        <b>copy(x)</b><br/>
        <i>unabhängige Kopie</i>
    """
    try:
        c = x.copy()
    except:
        c = x
    return c

#===============================================================

def fromTo(a, b, step=1):
    """ <b>fromTo(a, b, step=1)</b><br/>
        <i>Iterator over equidistant integers from <var>a</var> to <var>b</var></i><br/>
        If no <var>step</var> is given, <var>step</var> is set to 1.<br/>
        <b>Examples</b>:<pre>
    list(fromTo(3, 6))               # --->  [3,4,5,6]
    list(fromTo(3, 6, 2))            # --->  [3,5]
    list(fromTo(6, 3))               # --->  []
    list(fromTo(6, 3, -2))           # --->  [6,4]
    [2*i for i in fromTo(6, 3, -2)]  # --->  [12,8]
    sum(fromTo(1, 10))               # --->  55</pre>
    ----de----
        <b>fromTo(a, b, step=1)</b><br/>
        <i>Iterator über die ganzen Zahlen von a bis einschl. b im Abstand step</i><br/>
        Wenn <var>step</var> nicht übergeben wird, wird <var>step</var> auf 1 gesetzt.<br/>
        <b>Beispiele</b>:<pre>
    list(fromTo(3, 6))               # --->  [3,4,5,6]
    list(fromTo(3, 6, 2))            # --->  [3,5]
    list(fromTo(6, 3))               # --->  []
    list(fromTo(6, 3, -2))           # --->  [6,4]
    [2*i for i in fromTo(6, 3, -2)]  # --->  [12,8]
    sum(fromTo(1, 10))               # --->  55</pre>
    """
    if step > 0:
        return range(a, b+1, step)
    else:
        return range(a, b-1, step)

#===============================================================
#  Elementary Functions
#---------------------------------------------------------------

def isInteger(n):
    """ <b>isInteger(n)</b><br/>
        returns <code>True</code> if and only if n is of type <code>int</code>.
    ----de----
        <b>isInteger(n)</b><br/>
        gibt <code>True</code> zurück genau dann, wenn n vom Typ <code>int</code> ist.
    """
    return type(n) == int or type(n) == Rational and n.q == 1

def isNatural(n):
    return isInteger(n) and n > 0

#............................................................

def sum(obj, start=None, end=None, step=1):
    """ <b>sum(obj, start=None, end=None, step=1)</b><br/>
        <i>The sum of the elements of iterable obj or of all obj(i) for i in fromTo(start,end,step)</i><br/>
        The element type must define the operator + (e.g. numbers or strings).<br/>
        <b>Examples</b>:<br/>
        <code>sum([i^2 for i n fromTo(1,10)])</code><br/>
        <code>sum(fibonacci, 0, 10)</code><br/>
        <code>sum(lambda i:i^2, 2, 10, 2)</code>
    ----de----
        <b>sum(obj, start=None, end=None, step=1)</b><br/>
        <i>Summe der Elemente des iterierbaren Objekts a bzw. aller obj(i) für i in fromTo(start,end,step)</i><br/>
        Die Elemente müssen mit dem Operator + verknüpfbar sein
        (z. B. Zahlen oder Strings).<br/>
        <b>Beispiele</b>:<br/>
        <code>sum([i^2 for i n fromTo(1,10)])</code><br/>
        <code>sum(fibonacci, 0, 10)</code><br/>
        <code>sum(lambda i:i^2, 2, 10, 2)</code>
    """
    s = 0 # for an empty iteration the element type ist assumed to be numeric.
    if str(type(obj)) != "<class 'function'>":
        # obj is sequence or iterator or generator
        for i in obj:
            if s == 0: # adding to 0 would not work for types with neutral element != 0 (e.g. for a list of strings).
                s = i
            else:
                s += i
    else: # obj is a function
        for i in fromTo(start, end, step):
            if s == None:
                s = obj(i)
            else:
                s += obj(i)
    return s

#............................................................

def product(obj, start=None, end=None, step=1):
    """ <b>product(obj, start=None, end=None, step=1)</b><br/>
        <i>The product of the elements of iterable obj or of all obj(i) for i in fromTo(start,end,step)</i><br/>
        The element type must define the operator *.<br/>
        <b>Examples</b>:<br/>
        <code>product([i^2 for i n fromTo(1,10)])</code><br/>
        <code>product(fibonacci, 0, 10)</code><br/>
        <code>product(lambda i:i^2, 2, 10, 2)</code>
    ----de----
        <b>product(obj, start=None, end=None, step=1)</b><br/>
        <i>Product der Elemente des iterierbaren Objekts a bzw. aller obj(i) für i in fromTo(start,end,step)</i><br/>
        Die Elemente müssen mit dem Operator * verknüpfbar sein<br/>
        <b>Beispiele</b>:<br/>
        <code>product([i^2 for i n fromTo(1,10)])</code><br/>
        <code>product(fibonacci, 0, 10)</code><br/>
        <code>product(lambda i:i^2, 2, 10, 2)</code>
    """
    s = 1 # for an empty iteration the element type ist assumed to be numeric.
    if str(type(obj)) != "<class 'function'>":
        # obj is sequence or iterator or generator
        for i in obj:
            if s == 1:  # multiplying to 1 would not work for types with neutral element != 1
                s = i
            else:
                s *= i
    else: # obj is a function
        for i in fromTo(start, end, step):
            if s == None:
                s = obj(i)
            else:
                s *= obj(i)
    return s

#............................................................

def closure(elements, operation):
    """ <b>closure(elements, operation)</b><br/>
        <i>The closure of the elements under the binary operation</i><br/>
        operation must be given as a function with two parameters.
    ----de----
        <b>closure(elements, operation)</b><br/>
        <i>Der Abschluss der Menge elements unter der binären Operation operation</i><br/>
        operation muss eine function mit zwei Parametern sein.
    """
    closureList = list(elements)
    stringList = [str(x) for x in closureList]
    more = True
    while more:
        more = False
        for x in closureList:
            for y in closureList:
                e = operation(x,y)
                s = str(e)
                if s not in stringList:
                    closureList.append(e)
                    stringList.append(s)
                    more = True
    closureList.sort(key=str)
    return closureList

#............................................................

def join(l, filler=""):
    """ <b>join(l, filler="")</b><br/>
        <i>Joins the elements of <code>l</code> to a string</i><br/>
        Between each two adjacent elements <code>filler</code> is inserted.<br/>
        This is a replacement for the string method <code>join()</code>.<br/>
        In addition to the string method the list elements are converted to strings.<br/>
        <b>Example</b>:<br/>
        <code>join(fromTo(1,5), "-")</code>
        returns <code>"1-2-3-4-5"</code>.
    ----de----
        <b>join(l, filler="")</b><br/>
        <i>verbindet die Elemente von <code>l</code> zu einem String</i><br/>
        Zwischen je zwei Elemente wird <code>filler</code> eingefügt.<br/>
        Ersatz für die schwer lesbare String-Methode join()
        bzw. für die Funktion string.join().<br/>
        Im Gegensatz zu diesen Funktionen werden zusätzlich
        die Listenelemente zu Strings konvertiert.<br/>
        <b>Beispiel</b>:<br/>
        <code>join(fromTo(1,5), "-")</code>
        liefert <code>"1-2-3-4-5"</code>.
    """
    return filler.join([str(x) for x in l])


#===============================================================
#  Algorithms of Elementary Number Theory
#---------------------------------------------------------------

#---------------------------------------------------------------
#  Large random numbers

def rand(n):
    """ <b>rand(n)</b><br/>
        Integer random number from the range 0 ... <var>n</var>-1<br/>
        <var>n</var> must be a natural number (int).
    ----de----
        <b>rand(n)</b><br/>
        Zufallszahl im Bereich 0 ... <var>n</var>-1<br/>
        <var>n</var> kann eine beliebige natürliche Zahl sein (int).
    """
    return random.randrange(n)

#---------------------------------------------------------------
#  Conversion between positional notation systems

class NumeralSystem: # conversions between numeral systems # Umrechnung zwischen Stellenwertsystemen
    """ Class <code>NumeralSystem</code>: Numeral Systems<br/>
        This is a static class (only static methods available).
    ----de----
        Klasse <code>NumeralSystem</code>: Stellenwertsysteme<br/>
        Statische Klasse: nur Klassenmethoden.
    """
    #............................................................
    _staticMethods = ("toBase", "binary", "fromBase")

    @staticmethod
    def toBase(n, b, ascending=False):
        """ <b>toBase(n, b)</b><br/>
            <i><var>n</var> in the base <var>b</var> numeral system</i><br/>
            The result is returned as a list of the digits.<br/>
            <b>See also</b> <a href="#binary"><code>binary</code></a>,
                            <a href="#fromBase"><code>fromBase</code></a>.
        ----de----
            <b>toBase(n, b)</b><br/>
            <i><var>n</var> im Stellenwertsystem zur Basis <var>b</var></i><br/>
            Das Ergebnis wird als Liste der Stellen ausgegeben.<br/>
            <b>Siehe auch</b> <a href="#binary"><code>binary</code></a>,
                              <a href="#fromBase"><code>fromBase</code></a>.
        """
        v = []
        while n > 0:
            if ascending:
                v.append(n % b)    # insert n % b at the end
            else:
                v.insert(0, n % b) # insert n % b at the beginning
            n //= b
        return v

    @staticmethod
    def binary(n, ascending=False):
        """ <b>binary(n)</b><br/>
            <i>Binary number representation of <var>n</var> as list</i><br/>
            Abbreviation for <code>toBase(n,2)</code><br/>
            <b>Example</b>:<br/>
            <code>NumeralSystem.binary(11)</code> returns <code>[1, 0, 1, 1]</code>.<br/>
            <b>See also</b> <a href="#toBase"><code>toBase</code></a>,
                            <a href="#fromBase"><code>fromBase</code></a>.
        ----de----
            <b>binary(n)</b><br/>
            <i>Binärdarstellung von <var>n</var> als Liste</i><br/>
            Kurzform für <code>toBase(n,2)</code><br/>
            <b>Beispiel</b>:<br/>
            <code>NumeralSystem.binary(11)</code> gibt <code>[1, 0, 1, 1]</code> zurück.<br/>
            <b>Siehe auch</b> <a href="#toBase"><code>toBase</code></a>,
                              <a href="#fromBase"><code>fromBase</code></a>.
        """
        return NumeralSystem.toBase(n, 2, ascending)

    @staticmethod
    def fromBase(v, b, ascending=False):
        """ <b>fromBase(v, b)</b><br/>
            <i>Converts <var>v</var> from the base <var>b</var> numeral system to a number</i><br/>
            <var>v</var> is expected as list of digits.<br/>
            <b>See also</b> <a href="#toBase"><code>toBase</code></a>,
                            <a href="#binary"><code>binary</code></a>.
        ----de----
            <b>fromBase(v, b)</b><br/>
            <i><var>v</var> aus dem Stellenwertsystem zur Basis <var>b</var> als Zahl</i><br/>
            <var>v</var> wird als Liste der Stellen erwartet.<br/>
            <b>Siehe auch</b> <a href="#toBase"><code>toBase</code></a>,
                              <a href="#binary"><code>binary</code></a>.
        """
        n = 0
        if ascending:
            for i in fromTo(len(v)-1, 0, -1):
                n = b * n + i
        else:
            for i in v:
                n = b * n + i
        return n

#---------------------------------------------------------------
#  Integer functions

def binomial(n, k):
    """ <b>binomial(n, k)</b><br/>
        <i>Binomial coefficient</i><br/>
        n must be a nonnegative integer and k in the range 0 .. n.
    ----de----
        <b>binomial(n, k)</b><br/>
        <i>Binomialkoeffizient (n über k)</i><br/>
        n muss nichtnegative ganze Zahl sein, k im Bereich 0 .. n.
    """
    if k > n//2:
        k = n - k
    b = 1
    for i in fromTo(1, k):
        b = b * (n+1-i) // i
    return b

#............................................................

def factorial(n):
    """ <b>factorial(n)</b><br/>
        <i>n! = 1*2*3* ... *n (factorial of n)</i>
    ----de----
        <b>factorial(n)</b><br/>
        <i>n! = 1*2*3* ... *n (n Fakultät)</i>
    """
    f = 1
    for i in fromTo(2, n):
        f *= i
    return f

#............................................................

def fibonacci(n):
    """ <b>fibonacci(n)</b><br/>
        <i><var>n</var>-th element of the Fibonacci series</i>
        fibonacci(0) = 0<br/>
        fibonacci(1) = 1<br/>
        fibonacci(n) = fibonacci(n-1) + fibonacci(n-2) für n >= 2
    ----de----
        <b>fibonacci(n)</b><br/>
        <i><var>n</var>-tes Element der Fibonacci-Folge</i>
        fibonacci(0) = 0<br/>
        fibonacci(1) = 1<br/>
        fibonacci(n) = fibonacci(n-1) + fibonacci(n-2) für n >= 2
    """
    a, b = 0, 1
    while n > 0:
        a, b = b, a+b
        n -= 1
    return a

#---------------------------------------------------------------
#  Elementary number theoretic functions
#  Funktionen der elementaren Zahlentheorie

def gcd(a,b):
    """ <b>gcd(a,b)</b><br/>
        <i>greatest common divisor of a and b</i><br/>
        Calculated using the Euclidean algorithm<br/>
        <b>See also</b> <a href="#lcm"><code>lcm</code></a>
    ----de----
        <b>gcd(a,b)</b><br/>
        <i>Größter gemeinsamer Teiler von a und b</i><br/>
        Berechnung mit dem euklidischen Algorithmus<br/>
        <b>Siehe auch</b> <a href="#lcm"><code>lcm</code></a>
    """
    while b != 0:
        a, b = b, a % b
    return a

#............................................................

def _doc_gcdBinary(a,b):
    """ Binary Euclidean algorithm
        For optimization replace:
        x //= 2     by  x >>= 1
        x *= 2      by  x <<= 1
        x % 2 == 0  by  x &amp; 1 == 0
    ----de----
        Binärer Euklidischer Algorithmus
        Zur Optimierung sollte ersetzt werden:
        x //= 2     durch  x >>= 1
        x *= 2      durch  x <<= 1
        x % 2 == 0  durch  x &amp; 1 == 0
    """
    assert a >= b
    g = 1
    while a % 2 == 0 and b % 2 == 0:
        # größte gemeinsame Zweierpotenz:
        a //= 2
        b //= 2
        g *= 2
    while a != 0:
        while a % 2 == 0:
            a //= 2
        while b % 2 == 0:
            b //= 2
        t = abs(a-b) // 2
        if a >= b:
            a = t
        else:
            b = t
    return g * b

#............................................................

def _doc_gcdRecursive(a,b):
    """ recursive version of gcd
    ----de----
        rekursive Variante von gcd
    """
    if b == 0:
        return a
    else:
        return _doc_gcdRecursive(b, a % b)

#............................................................

def gcdExt(a,b):
    """ <b>gcdExt(a,b)</b><br/>
        <i>Extended Euclidean algorithm.</i><br/>
        Returns (d, x, y) with d == gcd(a,b) and x*a + y*b == d
    ----de----
        <b>gcdExt(a,b)</b><br/>
        <i>Erweiterter Euklidischer Algorithmus.</i><br/>
        Gibt (d, x, y) zurück mit d == gcd(a,b) und  x*a + y*b == d
    """
    if b == 0:                 # Basis for Math. Induction over b
        return a, 1, 0         # gcd(a,b) == 1*a + 0*b
    else:
        q, r = divmod(a, b)    # r < b
        d, z, x = gcdExt(b, r) # d == z*b + x*r (Induction assumption)
                               #   == z*b + x*(a - b*q)
                               #   == x*a + (z - x*q)*b
        y = z - x * q          # d == x*a + y*b
        return d, x, y

#............................................................

def lcm(a, b):
    """ <b>lcm(a,b)</b><br/>
        <i>Least common multiple of <var>a</var> and <var>b</var></i><br/>
        <b>See also</b> <a href="#gcd"><code>gcd</code></a>
    ----de----
        <b>lcm(a,b)</b><br/>
        <i>Kleinstes gemeinsames Vielfaches von <var>a</var> und <var>b</var></i><br/>
        <b>Siehe auch</b> <a href="#gcd"><code>gcd</code></a>
    """
    assert a > 0 and b > 0
    return (a * b) // gcd(a, b)

#............................................................

def _chinese2(mod1, mod2):
    m, n = mod1.m, mod2.m
    a, b = mod1.n, mod2.n
    d, x, y = gcdExt(m,n)
    assert d == 1    # m*x + n*y == 1
    x1 = x % n       # m*x1 % n == 1
    y1 = y % m       # n*y1 % m == 1
    c = n * y1       # c % m == 1
                     # c % n == 0
    d = m * x1       # d % n == 1
                     # d % m == 0
                     # (a*c + b*d) % m == a
                     # (a*c + b*d) % n == b
    return Mod(a*c + b*d, m*n)

def chinese(mod1, *mod2):
    """ <b>chinese(mod1, mod2,...)</b><br/>
        <i>Mod object from mod1, mod2, ... according to the Chinese remainder theorem</i><br/>
        Requires that mod1, mod2, ... are Mod objects with pairwise coprime moduli.<br/>
        <code>chinese(Mod(a1, m1), Mod(a2, m2), ...)</code><br/>
        returns Mod(x, m1*m2*...) with 0 &lt;= <code>x</code> &lt; m1*m2*... and<br/>
        <code>x % m1 == a1</code><br/>
        <code>x % m2 == a2</code><br/>, ...<br/>
        <b>Example</b>:<pre>
    chinese(Mod(1,5), Mod(2,7))  # ---> Mod(16,35)</pre>
        (16 mod 5 = 1, 16 mod 7 = 2).<br/>
        <b>See also</b> <a href="Mod.html"><code>Mod</code></a>
    ----de----
        <b>chinese(mod1, mod2,...)</b><br/>
        <i>Mod-Objekt zu mod1 mod2, ... nach Chinesischem Restsatz</i><br/>
        Voraussetzung: mod1, mod2, ... sind Mod-Objekte mit paarweise
        teilerfremden Moduln.<br/>
        <code>chinese(Mod(a1, m1), Mod(a2, m2, ...)</code><br/>
        berechnet Mod(x, m1*m2*...) mit 0 &lt;= <code>x</code> &lt; m1*m2*... und<br/>
        <code>x % m1 == a1</code><br/>
        <code>x % m2 == a2</code><br/>, ...<br/>
        <b>Beispiel</b>:<pre>
    chinese(Mod(1,5), Mod(2,7))  # ---> Mod(16,35)</pre>
        (16 mod 5 = 1, 16 mod 7 = 2).<br/>
        <b>Siehe auch</b> <a href="Mod.html"><code>Mod</code></a>
    """
    r = mod1.copy()
    for m in mod2:
        r = _chinese2(r, m)
    return r

#............................................................

def eulerPhi(n):
    """ <b>eulerPhi(n)</b><br/>
        <i>Euler's totient function</i><br/>
        Number of positive integers less than or equal to <var>n</var> that are coprime to <var>n</var>
    ----de----
        <b>eulerPhi(n)</b><br/>
        <i>Eulersche Phi-Funktion</i><br/>
        Anzahl der zu <var>n</var> teilerfremden Zahlen im Bereich von 1 bis <var>n</var>
    """
    #--- directly following the definition:
    # return sum([int(gcd(i,n) == 1) for i in Range(n)])

    #--- more efficient:
    for p in factors(n):
        n = n*(p-1)//p
    return n

#............................................................

def expMod(a, x, m):
    """ <b>expMod(a, x, m)</b><br/>
        <i>Calculates (a**x) % m efficiently.</i><br/>
        Equivalent to the Python standard function <code>pow</code>.
    ----de----
        <b>expMod(a, x, m)</b><br/>
        <i>Berechnet (a**x) % m effizient.</i><br/>
        Synonym für die Python-Standardfunktion <code>pow</code>.
    """
    return pow(a, x, m)  # Python standard function pow

#............................................................
# The following two implementations are for documentation only:
# Die folgenden beiden Implementierungen dienen nur der Dokumentation:

def _doc_expMod_recursive(a, x, n):
    """ recursive implementation of expMod
        (Stack overflow possible!)
    ----de----
        rekursive Implementierung von expMod
        (Stack-Überlauf möglich!)
    """
    if x == 1:
        return a % n
    elif x % 2 == 1:
        return (a * _doc_expMod_recursive(a,x-1,n)) % n
    else:
        t = _doc_expMod_recursive(a, x//2, n)
        return (t * t) % n

def _doc_expMod_iter(a,x,n):
    """ iterative implementation of expMod
    ----de----
        iterative Implementierung von expMod
    """
    b, c = 1, a      # TODO: verification
    while x > 1:
        if x % 2 == 1:
            b = (b*c) % n
        c = (c*c) % n
        x = x // 2
    return (b*c) % n

#............................................................

def _divisor(n):
    """ Pollard-Rho Algorithm (see Cormen, p. 845)
    """
    i, k = 1, 2
    x = rand(n)
    y = x
    while i < n:
        i += 1
        x = (x*x-1) % n
        d = gcd(y-x, n)
        if d != 1 and d != n:
            return d
        if i == k:
            y = x
            k = 2 * k

#............................................................

def isPrime(n,s=20):
    """ <b>isPrime(n)</b><br/>
        <i>Returns bool (is <var>n</var> a prime number?)</i><br/>
        Implementation using Miller Rabin test<br/>
        (see. Cormen, Leiserson, Rivest: Algorithms, p. 841)<br/>
        <b>See also</b> <a href="#prime"><code>prime</code></a>,
                        <a href="#nextPrime"><code>nextPrime</code></a>.
    ----de----
        <b>isPrime(n)</b><br/>
        <i>Rückgabewert bool (ist <var>n</var> Primzahl?)</i><br/>
        Implementierung mit Miller-Rabin-Test<br/>
        (vgl. Cormen, Leiserson, Rivest: Algorithms, S. 841)<br/>
        <b>Siehe auch</b> <a href="#prime"><code>prime</code></a>,
                          <a href="#nextPrime"><code>nextPrime</code></a>.
    """
    def _witness(a,n):
        """ witness for n being a product (see. Cormen p. 840)
            Local helper function
        ----de----
            Zeuge dafür, dass n zusammengesetzt ist (vgl. Cormen S. 840)
            Lokale Hilfsfunktion
        """
        b = NumeralSystem.binary(n-1)
        d = 1
        for i in b:
            x = d
            d = (d*d) % n
            if d == 1 and x != 1 and x != n-1:
                return True
            if i == 1:
                d = (d*a) % n
        return d != 1

    if n < 4:
        return n > 1
    for i in fromTo(1, s):
        a = 1 + rand(n-1)
        if _witness(a,n):
            return False
    return True

#............................................................
# nextPrime

def nextPrime(n):
    """ <b>nextPrime(n)</b><br/>
        <i>next prime larger than or equal to <var>n</var></i><br/>
        <b>See also</b> <a href="#isPrime"><code>isPrime</code></a>,
                        <a href="#prime"><code>prime</code></a>.
    ----de----
        <b>nextPrime(n)</b><br/>
        <i>nächste Primzahl ab (einschließlich) <var>n</var></i><br/>
        <b>Siehe auch</b> <a href="#isPrime"><code>isPrime</code></a>,
                          <a href="#prime"><code>prime</code></a>.
    """
    if n < 3:
        return 2
    if n % 2 == 0:
        n += 1
    while not isPrime(n):
        n += 2
    return n

#............................................................

# pre-calculated primes:
#        prime(0), prime(1000), prime(2000), ..., prime(99000),
_P1000 = [     0,    7919,   17389,   27449,   37813,   48611,   59359,   70657,   81799,   93179,
          104729,  116447,  128189,  139901,  151703,  163841,  176081,  187963,  200183,  212369,
          224737,  237203,  249439,  262139,  274529,  287117,  300023,  312583,  324949,  337541,
          350377,  363269,  376127,  389171,  401987,  414977,  427991,  440723,  453889,  467213,
          479909,  493127,  506131,  519227,  532333,  545747,  559081,  572311,  585493,  598687,
          611953,  625187,  638977,  652429,  665659,  679277,  692543,  706019,  719639,  732923,
          746773,  760267,  773603,  787207,  800573,  814279,  827719,  841459,  855359,  868771,
          882377,  896009,  909683,  923591,  937379,  951161,  965113,  978947,  993107, 1006721,
         1020379, 1034221, 1048129, 1062511, 1076143, 1090373, 1103923, 1117579, 1131617, 1145689,
         1159523, 1173301, 1187003, 1200949, 1215133, 1229269, 1243709, 1257517, 1271293, 1285517]

def prime(n):
    """ <b>prime(n)</b><br/>
        <i><var>n</var>-th prime number</i><br/>
        <b>See also</b> <a href="#isPrime"><code>isPrime</code></a>,
                        <a href="#nextPrime"><code>nextPrime</code></a>.
    ----de----
        <b>prime(n)</b><br/>
        <i><var>n</var>-te Primzahl</i><br/>
        <b>Siehe auch</b> <a href="#isPrime"><code>isPrime</code></a>,
                          <a href="#nextPrime"><code>nextPrime</code></a>.
    """
    if n < 100000:
        t, e = divmod(n, 1000)
        p = _P1000[t]
        for i in range(e):
            p = nextPrime(p+1)
        return p
    else:
        raise _MathGUIdeError("Argument is to large (must be &lt; 100000)!")

#............................................................

def factor(n):
    """ <b>factor(n)</b><br/>
        <i>Prime factorization</i><br/>
        Returns a list of pairs [factor, multiplicity]<br/>
        <b>Example</b>:<br/>
        <code>factor(24)</code><br/>
        returns <code>[[2,3],[3,1]]</code>.<br/>
        <b>See also</b> <a href="#factors"><code>factors</code></a>,
                        <a href="#allFactors"><code>allFactors</code></a>.
    ----de----
        <b>factor(n)</b><br/>
        <i>Primfaktorenzerlegung</i><br/>
        Rückgabewert: Liste mit Paaren [Primfaktor, Exponent]<br/>
        <b>Beispiel</b>:<br/>
        <code>factor(24)</code><br/>
        gibt <code>[[2,3],[3,1]]</code> zurück.<br/>
        <b>Siehe auch</b> <a href="#factors"><code>factors</code></a>,
                          <a href="#allFactors"><code>allFactors</code></a>.
    """
    p = 2
    i = 0
    r = []
    while p*p <= n:
        while n % p == 0:
            i += 1
            n //= p
        if i > 0:
            r.append([p,i])
        p = nextPrime(p+1)
        i = 0
    if n > 1:
        r.append([n,1])
    return r

#............................................................

def factors(n):
    """ <b>factors(n)</b><br/>
        <i>List of prime factors of <var>n</var> (without multiplicity)</i><br/>
        <b>Example</b>:<br/>
        <code>factors(24)</code>
        returns <code>[2, 3]</code>.<br/>
        <b>See also</b> <a href="#factor"><code>factor</code></a>,
                        <a href="#allFactors"><code>allFactors</code></a>.
    ----de----
        <b>factors(n)</b><br/>
        <i>Liste der Primfaktoren von <var>n</var> (jeder nur einmal)</i><br/>
        <b>Beispiel</b>:<br/>
        <code>factors(24)</code>
        gibt <code>[2, 3]</code> zurück.<br/>
        <b>Siehe auch</b> <a href="#factor"><code>factor</code></a>,
                            <a href="#allFactors"><code>allFactors</code></a>.
    """
    return [f[0] for f in factor(n)]

#............................................................

def allFactors(n):
    """ <b>allFactors(n)</b><br/>
        <i>List of all prime factors of <var>n</var></i><br/>
        <b>Example</b>: <code>allFactors(24)</code>
        returns <code>[2, 2, 2, 3]</code>.<br/>
        <b>See also</b> <a href="#factor"><code>factor</code></a>,
                        <a href="#factors"><code>factors</code></a>.
    ----de----
        <b>allFactors(n)</b><br/>
        <i>Liste aller Primfaktoren von <var>n</var></i><br/>
        <b>Beispiel</b>: <code>allFactors(24)</code>
        gibt <code>[2, 2, 2, 3]</code> zurück.<br/>
        <b>Siehe auch</b> <a href="#factor"><code>factor</code></a>,
                          <a href="#factors"><code>factors</code></a>.
    """
    p = 2
    r = []
    while p*p <= n:
        while n % p == 0:
            r.append(p)
            n //= p
        p = nextPrime(p+1)
    if n > 1:
        r.append(n)
    return r

#............................................................

def divisors(n):
    """ <b>divisors(n)</b><br/>
        <i>List of all divisors of <var>n</var></i>
    ----de----
        <b>divisors(n)</b><br/>
        <i>Liste der Teiler von <var>n</var></i>
    """
    # TODO: more efficient implementation!
    d = []
    for i in fromTo(1, n//2):
        if n % i == 0:
            d.append(i)
    d.append(n)
    return d


#===============================================================
#  Rational numbers and continued fractions
#  Rationale Zahlen, Kettenbrüche
#---------------------------------------------------------------

class Rational: # Rational numbers # Rationale Zahlen
    """ Class <code>Rational</code>: Rational numbers<br/>
        <b>Rational(p, q)</b> Example: <code>Rational(3,17)</code><br/>
        <b>Rational(n)</b> Example: <code>Rational(3)</code><br/>
        <b>Rational(s)</b> Example: <code>Rational("3/17")</code><br/>
        For <code>Rational</code> objects the arithmetical operators
        +, -, *, / are defined.
    ----de----
        Klasse <code>Rational</code>: Rationale Zahl<br/>
        <b>Rational(p, q)</b> Beispiel: <code>Rational(3,17)</code><br/>
        <b>Rational(n)</b> Beispiel: <code>Rational(3)</code><br/>
        <b>Rational(s)</b> Beispiel: <code>Rational("3/17")</code><br/>
        <code>Rational</code>-Objekte können mit den
        arithmetischen Operatoren +, -, *, / verknüpft werden.
    """
    def operators(self):
        """ <b>Rational.operators()</b><br/>
            <i>For documentation only</i><br/>
            The following operators are defined in the class <b>Rational</b>:<br/>
            <table>
             <tr><th>Op.</th><th>Function</th><th>Examples</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Subtraction</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>Multiplication</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>/</th><td>Division</td><td><code>a / b; a /= b</code></td></tr>
             <tr><th>-</th><td>Unary minus</td><td><code>-a</code></td></tr>
             <tr><th>&lt; > &lt;= >=</th><td>Comparision operators</td><td></td></tr>
            </table>
            <font color="#000080" size="-2">This method is for documentation only. It has no effect.</font>
        ----de----
            <b>Rational.operators()</b><br/>
            <i>Nur zur Dokumentation</i><br/>
            Die Klasse <b>Rational</b> erlaubt folgende Operatoren:<br/>
            <table>
             <tr><th>Op.</th><th>Funktion</th><th>Beispiele</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Subtraktion</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>Multiplikation</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>/</th><td>Division</td><td><code>a / b; a /= b</code></td></tr>
             <tr><th>-</th><td>unäres Minus</td><td><code>-a</code></td></tr>
             <tr><th>&lt; > &lt;= >=</th><td>Vergleichsoperatoren</td><td></td></tr>
            </table>
            <font color="#000080" size="-2">Diese Methode dient nur zur Dokumentation. Sie hat keine Wirkung.</font>
        """
        pass

    def __init__(self, p, q=1):
        """ Constructor
        ----de----
            Konstruktor
        """
        if type(p) == float:
            r = ContFrac(p,12).toRational()
            self.p, self.q = r.p, r.q
            return
        if type(p) == Rational:
            self.p, self.q = p.p, p.q
            return
        if type(p) == str:   #  e.g. "3/17"
            l = p.split("/")
            p = int(l[0])
            if len(l) > 1:
                q = int(l[1])
        if q < 0:
            p,q = -p, -q  # denominator always positive
        d = gcd(abs(p), q)
        if d > 1:         # reduce fraction if possible
            p //= d
            q //= d
        self.p = p
        self.q = q

    def copy(self):
        """ <b>r.copy()</b><br/>
            <i>independent copy ("clone") of r</i>
        ----de----
            <b>r.copy()</b><br/>
            <i>unabhängige Kopie von r</i>
        """
        return Rational(self.p, self.q)

    def __abs__(self):
        """ absolute value
        ----de----
            Absolutbetrag
        """
        return Rational(abs(self.p), self.q)

    def abs(self):
        """ <b>r.abs()</b><br/>
            <i>absolute value</i><br/>
            Also as global function: <b>abs(r)</b><br/>
        ----de----
            <b>r.abs()</b><br/>
            <i>Absolutbetrag</i><br/>
            Auch globale Funktion: <b>abs(r)</b><br/>
        """
        return Rational(abs(self.p), self.q)

    def __float__(self):
        """ Conversion to floating point number
        ----de----
            Umwandlung in Gleitkommazahl
        """
        return float(self.p) / float(self.q)

    def toFloat(self):
        """ <b>r.toFloat()</b><br/>
            <i>Conversion to floating point number</i><br/>
            Also as global function: <b>float(r)</b><br/>
        ----de----
            <b>r.toFloat()</b><br/>
            <i>Konvertierung in Gleitkommazahl</i><br/>
            Auch globale Funktion: <b>float(r)</b><br/>
        """
        return float(self.p) / float(self.q)

    def __int__(self):
        """ Round to the nearest integer
        ----de----
            Auf-/Abrundung zur nächsten Ganzzahl
        """
        assert self.q > 0
        return (self.p + self.q/2) // self.q

    def isInteger(self):
        return self.q == 1

    def toInt(self):
        """ <b>r.toInt()</b><br/>
            <i>Round to the nearest integer</i><br/>
            Also as global function: <b>int(r)</b><br/>
        ----de----
            <b>r.toInt()</b><br/>
            <i>Rundung zur nächsten ganzen Zahl</i><br/>
            Auch globale Funktion: <b>int(r)</b><br/>
        """

    def __long__(self):
        """ Round to the nearest integer
        ----de----
            Auf-/Abrundung zur nächsten Ganzzahl
        """
        assert self.q > 0
        return (self.p + self.q/2) // self.q

    def _unify(self, other):
        return (self, Rational(other))

    def __repr__(self):
        """ Object representation as string
        """
        if self.q == 1:
            return str(self.p)
        else:
            return "{0}/{1}".format(self.p, self.q)

    def toStr(self):
        """ <b>r.toStr()</b><br/>
            <i>Representation as string</i><br/>
            Also as global function: <b>str(r)</b><br/>
        ----de----
            <b>r.toStr()</b><br/>
            <i>Darstellung als String</i><br/>
            Auch globale Funktion: <b>str(r)</b><br/>
        """

    def _cmp(a, b):
        if type(b) == float:
            return _cmp(float(a), b)
        if not isinstance(b, Rational):
            a,b = a._unify(b)
        return _cmp(a.p*b.q, b.p*a.q)

    def __lt__(a,b): return a._cmp(b) < 0
    def __le__(a,b): return a._cmp(b) <= 0
    def __eq__(a,b): return a._cmp(b) == 0
    def __ge__(a,b): return a._cmp(b) >= 0
    def __gt__(a,b): return a._cmp(b) > 0
    def __ne__(a,b): return a._cmp(b) != 0

    def __neg__(self):
        """ Unary operator -
            -a --> a.__neg__()
        """
        return Rational(-self.p, self.q)

    def __add__(a, b):
        """ <b>a.__add__(b)</b><br/>
            <i>Addition</i><br/>
            Also operator notation: <b>a + b</b><br/>
            Instead of <code>a = a + b</code> you can write <code>a += b</code>.
        """
        if type(b) == float:
            return float(a) + b
        if not isinstance(b, Rational):
            a, b = a._unify(b)
        return Rational(a.p*b.q + b.p*a.q, a.q*b.q)

    def __radd__(a, b):
        """ Addition, see __add__
        """
        return a + b

    def __sub__(a, b):
        """ Subtraction:
            a - b  --> a.__sub__(b),  if a is of type Rational
                       b.__rsub__(a), otherwise
        """
        if type(b) == float:
            return float(a) - b
        if not isinstance(b, Rational):
            a,b = a._unify(b)
        return Rational(a.p*b.q - b.p*a.q, a.q*b.q)

    def __rsub__(b, a):
        """ Subtraction, see __sub__
        """
        if type(a) == float:
            return a - float(b)
        if not isinstance(a, Rational):
            b,a = b._unify(a)
        return a - b

    def __mul__(a, b):
        """ Multiplication:
            a * b  --> a.__mul__(b),  if a is of type Rational
                       b.__rmul__(a), otherwise
        """
        if type(b) == float:
            return float(a) * b
        if type(b) in (Matrix, Poly):
            return b * a
        if not isinstance(b, Rational):
            a,b = a._unify(b)
        return Rational(a.p*b.p, a.q*b.q)

    def __rmul__(b, a):
        """ Multiplication, see __mul__
        """
        return b * a

    def __div__(a, b):
        """ exact Division
            a / b  --> a.__div__(b),  if a is of type Rational
                       b.__rdiv__(a), otherwise
        """
        if type(b) == float:
            return float(a) / b
        if not isinstance(b, Rational):
            a,b = a._unify(b)
        return Rational(a.p*b.q, a.q*b.p)

    def __rdiv__(b, a):
        """ exact Division, see __div__
        """
        if type(a) == float:
            return a / float(b)
        if not isinstance(a, Rational):
            b,a = b._unify(a)
        return Rational(a.p*b.q, a.q*b.p)

    def __truediv__(a, b):   return a.__div__(b)
    def __rtruediv__(a, b):  return a.__rdiv__(b)
    def __floordiv__(a, b):  return a.__div__(b)
    def __rfloordiv__(a, b): return a.__rdiv__(b)

def toNumber(s):
    """ <b>toNumber(s)</b><br/>
        <i>converts the string s into a rational or floating point number</i><br/>
        <b>Examples</b>:<br/>
        <code>toNumber("3.1")</code><br/>
        <code>toNumber("7/3")</code>
        ----de----
        <b>toNumber(s)</b><br/>
        <i>verwandelt den String s in eine rationale oder Gleitkommazahl</i><br/>
        <b>Beispiele</b>:<br/>
        <code>toNumber("3.1")</code><br/>
        <code>toNumber("7/3")</code>
    """
    if "." in s or "e" in s or "E" in s:
        return float(s)
    else:
        return Rational(s)

#---------------------------------------------------------------

class ContFrac: # Continued fractions # Kettenbrüche
    """ <b>ContFrac(number, steps=5)</b><br/>
        <i>Continued fraction from <b>number</b> with maximal length <b>steps</b></i><br/>
        <b>Example</b>:<br/>
        <code>ContFrac(pi, 3)</code> returns
        <code>ContFrac(3, 7, 15)</code>.<br/>
        This represent the number 3 + 1/(7 + 1/15).
        ----de----
        <b>ContFrac(number, steps=5)</b><br/>
        <i>Kettenbruch aus <b>number</b> mit Maximallänge <b>steps</b></i><br/>
        <b>Beispiel</b>:<br/>
        <code>ContFrac(pi, 3)</code> ergibt
        <code>ContFrac(3, 7, 15)</code>.<br/>
        Das steht für die Zahl 3 + 1/(7 + 1/15).
    """
    def __init__(self, number, steps=5):
        """ Constructor
        ----de----
            Konstruktor
        """
        self._l = []
        for i in range(steps):
            i = floor(number)
            self._l.append(i)
            if abs(number-i) < 0.00000001:
                break
            number = 1 / (number-i)

    def __repr__(self):
        """ Object representation as string
        """
        return "ContFrac({0})".format(join(self._l, ","))

    def toRational(self, n=0):
        """ <b>toRational(n=0)</b><br/>
            <i>Rational number defined by the first n list elements</i><br/>
            Converts the continued fraction int a rational number.<br/>
            <b>Example</b>:<br/>
            <code>ContFrac(pi).toRational(2)</code> returns 22/7.
        ----de----
            <b>toRational(n=0)</b><br/>
            <i>Rationale Zahl aus den ersten n Gliedern</i><br/>
            Wandelt den Kettenbruch in eine rationale Zahl um.<br/>
            <b>Beispiel</b>:<br/>
            <code>ContFrac(pi).toRational(2)</code> ergibt 22/7.
        """
        def contFracListToRat(cf):
            if len(cf) == 0: return Rational(0)
            if len(cf) == 1: return Rational(cf[0])
            return cf[0] + Rational(1) / contFracListToRat(cf[1:])

        if n == 0 or n > len(self._l):
            n = len(self._l)
        return contFracListToRat(self._l[0:n])

#===============================================================
#  Lineare Algebra
#  Klassen Vector und Matrix
#---------------------------------------------------------------

class Vector (list): # Vectors # Vektoren
    """ Class <code>Vector</code> (with arithmetic operators)<br/>
        <b>Vector(v)</b> (v Liste or Tuple)<br/>
        The Indices of the elements of Vector are counted from 0 (as in Python lists etc.).
        <b>Example:</b><br/>
        <code>Vector([1,2,3])</code>
        ----de----
        Klasse <code>Vector</code> (mit arithmetischen Operatoren)<br/>
        <b>Vector(v)</b> (v Liste oder Tupel)<br/>
        Die Indizes der Elemente von Vector werden (wie in Python-Listen etc.) ab 0 gezählt.
        <b>Beispiel:</b><br/>
        <code>Vector([1,2,3])</code>
    """
    #............................................................
    _staticMethods = ("fromFunction", "fromString")

    @staticmethod
    def fromFunction(n, fn, offset=0):
        """ <b>Vector.fromFunction(n, fn, offset=0)</b><br/>
            <i>Vector, the elements of which are to be calculated using the function fn</i><br/>
            Returns: <code>Vector([fn(i+offset) for i in range(n)])</code><br/>
            <b>Example:</b><br/>
            <code>Vector.fromFunction(5, sqrt, 1)</code> returns Vector([1.0, 1.414, 1.732, 2.0, 2.236])
        ----de----
            <b>Vector.fromFunction(n, fn, offset=0)</b><br/>
            <i>Vector, dessen Elemente mit der Funktion fn berechnet werden</i><br/>
            Rückgabewert: <code>Vector([fn(i+offset) for i in range(n)])</code><br/>
            <b>Beispiel:</b><br/>
            <code>Vector.fromFunction(5, sqrt, 1)</code> ergibt Vector([1.0, 1.414, 1.732, 2.0, 2.236])
        """
        return Vector([fn(i+offset) for i in range(n)])

    @staticmethod
    def fromString(s):
        """ <b>Vector.fromString(s)</b><br/>
            <i>Vector from String</i><br/>
            Elements (may also bbe Rationals) must be given as comma separated string.<br/>
            <b>Example:</b><br/>
            <code>Vector.fromString("1, 2/3, 2")</code>
        ----de----
            <b>Vector.fromString(s)</b><br/>
            <i>Vector aus String</i><br/>
            Elemente (auch rationale Zahlen) müssen im String durch Kommas getrennt werden.<br/>
            <b>Beispiel:</b><br/>
            <code>Vector.fromString("1, 2/3, 2")</code>
        """
        return Vector([toNumber(x) for x in s.split(",")])
    #............................................................

    def operators(self):
        """ <b>Vector.operators()</b><br/>
            <i>For documentation only</i><br/>
            The following operators are defined in the class <b>Vector</b>:<br/>
            <table>
             <tr><th>Op.</th><th>Function</th><th>Examples</th></tr>
             <tr><th>+</th><td>elementwise Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>elementwise Subtraction</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>if both operators are Vectors: inner product,
               otherwise: scalar multiplication (elementwise)</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>-</th><td>Unary minus</td><td><code>-a</code></td></tr>
             <tr><th>[ ]</th><td>Index operator</td><td><code>A[i]</code> (Element)</td></tr>
            </table>
            <font color="#000080" size="-2">This method is for documentation only. It has no effect.</font>
        ----de----
            <b>Vector.operators()</b><br/>
            <i>Nur zur Dokumentation</i><br/>
            Die Klasse <b>Vector</b> erlaubt folgende Operatoren:<br/>
            <table>
             <tr><th>Op.</th><th>Funktion</th><th>Beispiele</th></tr>
             <tr><th>+</th><td>Elementweise Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Elementweise Subtraktion</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>wenn beide Operanden Vektoren sind: Skalarprodukt,
                sonst Skalarmultiplikation (elementweise)</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>-</th><td>unäres Minus</td><td><code>-a</code></td></tr>
             <tr><th>[ ]</th><td>Indexoperator</td><td><code>A[i]</code> (Element)</td></tr>
            </table>
            <font color="#000080" size="-2">Diese Methode dient nur zur Dokumentation. Sie hat keine Wirkung.</font>
        """
        pass

    def __init__(v, n):
        """ Constructor
        ----de----
            Konstruktor
        """
        if isinstance(n, list):
            list.__init__(v, n)
        elif isinstance(n, tuple):
            list.__init__(v, list(n))

    def __repr__(A):
        """ Object representation as string
        """
        return "Vector({0})".format(list.__repr__(A))

    def __invert__(v):
        """ Operator ~ (transposition)
        """
        return v.transp()

    def transp(v):
        """ <b>v.transp()</b><br/>
            <i>berechnet die transponierte Matrix (Spaltenvektor)</i><br/>
            kürzer: ~v
        ----de----
            <b>v.transp()</b><br/>
            <i>calculates the transposed Matrix (column vektor)</i><br/>
            shorter: ~v
        """
        return Matrix([[a] for a in v])

    def join(x, y):
        """ <b>v.join(w)</b><br/>
            <i>Concatenation with Vector w</i>
        ----de----
            <b>v.join(w)</b><br/>
            <i>Verkettung mit dem Vektor w</i>
        """
        return Vector(list(x)+list(y))

    def __neg__(v):
        """ unary oparator - (negative Vektor)
        """
        return Vector([-x for x in v])

    def __add__(x, y):
        """ operators + and += (Elementwise addition)
        """
        #   a + b  --> a.__add__ (b), if a is of type Vector
        #              b.__radd__(a), otherwise
        n = len(x)
        assert n == len(y)
        return Vector([x[i]+y[i] for i in range(n)])

    def __sub__(x, y):
        """ operators - and -= (elementwise subtraction)
        """
        return x + (-y)

    def __mul__(x, y):
        """ Operators * and *=
            Multiplikation
            if both operands are vektors sind: inner product,
            otherwise scalar multiplication (elementwise)
        """
        if isinstance(y, Vector): # Skalarprodukt
            n = len(x)
            assert n == len(y)
            return sum([x[i]*y[i] for i in range(n)])
        else:
            return Vector([y*x[i] for i in range(len(x))])

    def __rmul__(x, c):
        return x * c

    def __imul__(x, c):
        for i in range(len(x)):
            x[i] *= c
        return x

    def norm(x):
        """ <b>v.norm()</b><br/>
            <i>Norm des Vektors v</i>
        ----de----
            <b>v.norm()</b><br/>
            <i>Norm of the Vector v</i>
        """
        return math.sqrt(x*x)



class Matrix (Vector): # Matrices # Matrizen
    """ class <code>Matrix</code><br/>
        <i>see also: Menu: Insert -- Matrix</i><br/>
        <b>Matrix(v)</b>: v list of lists with uniform length.<br/>
        <b>Example</b>:<br/>
        <code>Matrix([[11,12,13],[21,22,23]])</code><br/>
        The class <code>Matrix</code> ist derived from <code>Vector</code>:<br/>
        A <code>Matrix</code> object ist a <code>Vector</code>
        of <code>Vector</code> objects with uniform length.<br/>
        The elements can be addressed using double indices (e.g.: <code>A[i,k]</code>).
        Indices are counted from 0 (as in Python lists etc.).<br/>
        see also. class methods with <code>Matrix.</code>
        ----de----
        Klasse <code>Matrix</code><br/>
        <i>vgl. Menü: Einfügen -- Matrix</i><br/>
        <b>Matrix(v)</b>: v Liste von gleichlangen Listen,<br/>
        <b>Beispiel</b>:<br/>
        <code>Matrix([[11,12,13],[21,22,23]])</code><br/>
        Die Klasse <code>Matrix</code> ist abgeleitet von <code>Vector</code>:<br/>
        Ein <code>Matrix</code>-Objekt ist ein <code>Vector</code>
        von gleichlangen <code>Vector</code>-Objekten.<br/>
        Die Elemente können mit Doppelindizes (z.B.: <code>A[i,k]</code>) angesprochen werden.
        Indizes werden (wie in Python-Listen etc.) ab 0 gezählt.<br/>
        vgl. Klassenmethoden mit <code>Matrix.</code>
    """
    #............................................................
    _staticMethods = ("null", "identity", "fromFunction", "fromString", "random")

    @staticmethod
    def null(m, n=0):
        """ <b>Matrix.null(m, n=0)</b><br/>
            <i>m*n-Null matrix</i><br/>
            Square Matrix, if n is omitted.<br/>
        ----de----
            <b>Matrix.null(m, n=0)</b><br/>
            <i>m*n-Nullmatrix</i><br/>
            Quadratische Matrix, falls n weggelassen wird.<br/>
        """
        if n == 0:
            n = m
        return Matrix([ [0 for k in range(n)]
                        for i in range(m)])

    @staticmethod
    def identity(n):
        """ <b>Matrix.identity(n)</b><br/>
            <i>n*n Identity Matrix</i><br/>
        ----de----
            <b>Matrix.identity(n)</b><br/>
            <i>n*n-Einheitsmatrix</i><br/>
        """
        return Matrix([ [Rational(int(i==k)) for k in range(n)]
                        for i in range(n)])

    @staticmethod
    def fromFunction(m, n, fn, offset=0):
        """ <b>Matrix.fromFunction(m,n, fn, offset=0)</b><br/>
            <i><code>m*n</code>-Matrix defined by function <code>fn</code></i><br/>
            <code>fn</code> must be a binary function.<br/>
            The Matrix element <code>A[i,k]</code> is defined as <code>fn(i+offset,k+offset)</code>.<br/>
            <b>Example</b>:<br/>
            <code>Matrix.fromFunction(2,4, pow, 1)</code><br/>
            returns the Matrix<pre>
     / 1 1 1  1 \
     |          |
     \ 2 4 8 16 /</pre>
        ----de----
            <b>Matrix.fromFunction(m,n, fn, offset=0)</b><br/>
            <i>Durch Funktion <code>fn</code> berechnete <code>m*n</code>-Matrix</i><br/>
            <code>fn</code> muss eine zweistellige Funktion sein.<br/>
            Das allgemeine Element <code>A[i,k]</code> wird mit <code>fn(i+offset,k+offset)</code> berechnet.<br/>
            <b>Beispiel</b>:<br/>
            <code>Matrix.fromFunction(2,4, pow, 1)</code><br/>
            liefert die Matrix<pre>
     / 1 1 1  1 \
     |          |
     \ 2 4 8 16 /</pre>
        """
        return Matrix([ [fn(i+offset,k+offset) for k in range(n)]
                        for i in range(m)])

    @staticmethod
    def fromString(s):
        """ <b>Matrix.fromString(s)</b><br/>
            <i>Matrix given by a string</i><br/>
            The string s must contains the elementa row by row.<br/>
            The rows are divided by semicolon, the elemente within a row by comma.<br/>
            Elemente may also be rational (notated with slash).<br/>
            <b>Example</b>:<br/>
            <code>Matrix.fromString("1, 2, 3.14; 4/5, 5, 6")</code>
            returns the matrix<pre>
     / 1    2  3.14 \
     |              |
     \ 4/5  5   6   /</pre>
        ----de----
            <b>Matrix.fromString(s)</b><br/>
            <i>Matrix aus String</i><br/>
            Im String s werden die Elemente der Matrix zeilenweise angegeben.<br/>
            Die Zeilen werden durch Semikolon getrennt, die Elemente innerhalb der Zeilen durch Komma.<br/>
            Elemente können auch rational sein (mit Schrägstrich).<br/>
            <b>Beispiel</b>:<br/>
            <code>Matrix.fromString("1, 2, 3.14; 4/5, 5, 6")</code>
            liefert die Matrix<pre>
     / 1    2  3.14 \
     |              |
     \ 4/5  5   6   /</pre>
        """
        return Matrix([ [toNumber(a) for a in row.split(",")]
                        for row in s.split(";")])

    @staticmethod
    def random(m,n=0, r=100):
        """ <b>Matrix.random(m,n=0, r=100)</b><br/>
            <i>m*n matrix with random elements in 0..r-1</i>
        ----de----
            <b>Matrix.random(m,n=0, r=100)</b><br/>
            <i>m*n-Matrix mit Zufallswerten 0..r-1</i>
        """
        if n == 0:
            n = m
        return Matrix([ [Rational(rand(r)) for k in range(n)]
                        for i in range(m)])
    #............................................................

    def operators(self):
        """ <b>Matrix.operators()</b><br/>
            <i>For documentation only</i><br/>
            The following operators are defined in the class <b>Matrix</b>:<br/>
            <table>
             <tr><th>Op.</th><th>Function</th><th width="30%">Examples</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Subtraction</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>Matrizenmultiplikation if both operands are Matrices, otherwise: Scalar multiplication</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>-</th><td>Unary minus</td><td><code>-a</code></td></tr>
             <tr><th>~</th><td>Transposed Matrix</td><td><code>~A</code></td></tr>
             <tr><th>|</th><td>Concatenation</td><td><code>A | B; A |= B</code></td></tr>
             <tr><th>[ ]</th><td>Index operator</td><td><code>A[i,k]</code><br/><code>A[i]</code> (Row Vector)</td></tr>
            </table>
            <font color="#000080" size="-2">This method is for documentation only. It has no effect.</font>
        ----de----
            <b>Matrix.operators()</b><br/>
            <i>Nur zur Dokumentation</i><br/>
            Die Klasse <b>Matrix</b> erlaubt folgende Operatoren:<br/>
            <table>
             <tr><th>Op.</th><th>Funktion</th><th width="30%">Beispiele</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>A + B; A += B</code></td></tr>
             <tr><th>-</th><td>Subtraktion</td><td><code>A - B; A -= B</code></td></tr>
             <tr><th>*</th><td>Matrizenmultiplikation (wenn beide Operanden Matrizen sind) bzw. Skalarmultiplikation (sonst)</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>-</th><td>unäres Minus</td><td><code>-A</code></td></tr>
             <tr><th>~</th><td>Transponierte Matrix</td><td><code>~A</code></td></tr>
             <tr><th>|</th><td>Verkettung</td><td><code>A | B; A |= B</code></td></tr>
             <tr><th>[ ]</th><td>Indexoperator</td><td><code>A[i,k]</code><br/><code>A[i]</code> (Zeilenvektor)</td></tr>
            </table>
            <font color="#000080" size="-2">Diese Methode dient nur zur Dokumentation. Sie hat keine Wirkung.</font>
        """
        pass

    def __init__(A, v):
        """ Constructor
        ----de----
            Konstruktor
        """
        Vector.__init__(A, [Vector(v[i]) for i in range(len(v))])
        A._makeRational()

    def _makeRational(A):
        for i in A.rowRange():
            for k in A.colRange():
                if isInteger(A[i,k]):
                    A[i,k] = Rational(A[i,k])

    def float(A):
        """ <b>A.float()</b><br/>
            <i>returns the Matrix A with all elements converted to <code>float</code></i>
        ----de----
            <b>A.float()</b><br/>
            <i>gibt die Matrix A zurück, in der alle Elemente nach <code>float</code> konvertiert sind</i>
        """
        return Matrix([ [float(A[i,k]) for k in A.colRange()]
                        for i in A.rowRange()])

    def copy(A):
        """ <b>A.copy()</b><br/>
            <i>independent copy of the matrix ("clone") A</i>
        ----de----
            <b>A.copy()</b><br/>
            <i>unabhängige Kopie der Matrix A</i>
        """
        return Matrix([ [copy(A[i,k]) for k in A.colRange()]
                        for i in A.rowRange()])

    def isSquare(A):
        """ <b>A.isSquare()</b><br/>
            <i>Truth value (A is square matrix)</i>
        ----de----
            <b>A.isSquare()</b><br/>
            <i>Wahrheitswert (A quadratisch)</i>
        """
        return A.height() == A.width()

    def str(A, i,k):
        """ <b>A.str(i,k)</b><br/>
            <i>string representation of the element [i,k]</i>
        ----de----
            <b>A.str(i,k)</b><br/>
            <i>String-Darstellung des Elements [i,k]</i>
        """
        if type(A[i,k]) == float:
            return "%g" % A[i,k]
        else:
            return str(A[i,k])

    def _pp(A, name=""):
        n = A.height()
        s = """{{html}}<table><tr valign="middle">
                 <td width="30"><font color="#3060C0">{0}</font></td>
                 <td><table cellspacing="0" cellpadding="-1">""".format(name)
        for i in range(n):
            left, right = "21", "22"
            if n > 1:
                if   i == 0:   left, right = "11", "12"
                elif i == n-1: left, right = "31", "32"
            else:
                left, right = "01", "02"
            s += """<tr><td align="center" valign="middle"><img src=":/img/matrix/{0}.png"/></td>""".format(left)
            spc = {False: "&nbsp;", True: ""}
            for k in range(A.width()):
                s += """<td align="center" valign="middle"><font color="#3060C0">{1}{0}{2}</font></td>""".format(A[i,k], spc[k==0], spc[k==A.width()-1])
            s += """<td align="center" valign="middle"><img src=":/img/matrix/{0}.png"/></td></tr>""".format(right)
        s += """</table></td></table>{/html}"""
        return s

    def __repr__(A):
        """ Object representation as string
        """
        return A._pp()

    def __getitem__(A, i):
        """ Definition of the Index operator []
            Counts from 0
            A[i]   --> i-th row vector
            A[i,k] --> element in i-th row, k-th column
        """
        if isinstance(i, tuple):
            return A[i[0]][i[1]]
        else:
            return Vector.__getitem__(A, i)

    def __setitem__(A, i, x):
        """ Definition of the Index operator [] for assignments
            Counts from 0
            A[i]   --> i-th row vector
            A[i,k] --> element in i-th row, k-th column
        """
        if isinstance(i, tuple):
            A[i[0]][i[1]] = x
        else:
            Vector.__setitem__(A, i, x)

    def height(A):
        """ <b>A.height()</b><br/>
            <i>Number of rows of the Matrix A</i>
        ----de----
            <b>A.height()</b><br/>
            <i>Anzahl der Zeilen der Matrix A</i>
        """
        return len(A)

    def width(A):
        """ <b>A.height()</b><br/>
            <i>Number of columns of the Matrix A</i>
        ----de----
            <b>A.width()</b><br/>
            <i>Anzahl der Spalten der Matrix A</i>
        """
        return len(A[0])

    def rowRange(A):
        """ <b>A.rowRange()</b><br/>
            <i>List of all row indices</i>
        ----de----
            <b>A.rowRange()</b><br/>
            <i>Liste aller Zeilenindizes</i>
        """
        return range(A.height())

    def colRange(A):
        """ <b>A.rowRange()</b><br/>
            <i>List of all column indices</i>
        ----de----
            <b>A.rowRange()</b><br/>
            <i>Liste aller Spaltenindizes</i>
        """
        return range(A.width())

    def __neg__(A):
        """ unary minus operator
        """
        return Matrix([ [-A[i,k] for k in A.colRange()]
                        for i in A.rowRange()])

    def joinRight(A, B):
        """ <b>A.joinRight(B)</b><br/>
            <i>Horizontal concatenation of A with B</i><br/>
            <b>Requires</b>: Both Matrices must have the same row count.
        ----de----
            <b>A.concat(B)</b><br/>
            <i>Horizontale Verkettung von A mit B</i><br/>
            <b>Voraussetzung</b>: Die beiden Matrizen müssen gleich viele Zeilen haben.
        """
        assert A.height() == B.height()
        return Matrix([ list(A[i])+list(B[i]) for i in A.rowRange()])

    def joinBottom(A, B):
        """ <b>A.joinBottom(B)</b><br/>
            <i>Vertical concatenation of A with B</i><br/>
            <b>Requires</b>: Both Matrices must have the same column count.
        ----de----
            <b>A.concat(B)</b><br/>
            <i>Vertikale Verkettung von A mit B</i><br/>
            <b>Voraussetzung</b>: Die beiden Matrizen müssen gleich viele Spalten haben.
        """
        assert A.width() == B.width()
        def el(i,k):
            if i < A.height(): return A[i,k]
            else:              return B[i-A.height(), k]
        return Matrix.fromFunction(A.height()+B.height(), A.width(), el)

    def __or__(A, B):
        """ <b>A.__or__(B)</b><br/>
            <i>Horizontal concatenation with B</i><br/>
            A | B is the same as A.joinRight(B)<br/>
            Instead of A = A | B you may write A |= B.
            <b>Requires</b>: Both matrices must have the same row count.
        """
        assert A.height() == B.height()
        return Matrix([ list(A[i])+list(B[i]) for i in A.rowRange()])

    def __add__(A, B):
        """ <b>A.__add__(B)</b><br/>
            <i>Matrix addition operator</i><br/>
            Instead of A = A + B you may write A += B.
        """
        #   a + b  --> a.__add__ (b), wenn a vom Typ Matrix ist
        #              b.__radd__(a), sonst
        assert A.height() == B.height() and A.width() == B.width()
        return Matrix([ [A[i,k]+B[i,k] for k in A.colRange()]
                        for i in A.rowRange()
                      ])

    def __sub__(A, B):
        """ <b>A.__sub__(B)</b><br/>
            <i>Matrix subtraction operator</i><br/>
            Instead of A = A - B you may write A -= B.
        """
        return A + (-B)

    def __mul__(A, B):
        """ <b>A.__sub__(B)</b><br/>
            <i>Scalar or matrix multiplication</i><br/>
            Instead of A = A * B you may write A *= B.
        """
        # If both operands are matrices, matrix multiplication is applied,
        # otherwise scalar multiplikation
        if isinstance(B, Matrix): # matrix multiplication
            assert A.width() == B.height()
            return Matrix(
                [ [ sum([A[i,j] * B[j,k] for j in A.colRange()])
                    for k in B.colRange()]
                  for i in A.rowRange()])
        else: # scalar multiplikation
            return Matrix([ [B*A[i,k] for k in A.colRange()]
                            for i in A.rowRange()])

    def __rmul__(A, B):
        assert not isinstance(B, Matrix)
        return A * B

    def transp(A):
        """ <b>A.transp()</b><br/>
            <i>returns the transposed Matrix</i><br/>
            short form: ~A
        ----de----
            <b>A.transp()</b><br/>
            <i>berechnet die Transponierte der Matrix A</i><br/>
            kürzer: ~A
        """
        return Matrix([ [A[i,k] for i in A.rowRange()]
                        for k in A.colRange()])

    def __invert__(A):
        """ Operator ~ (transposition)
        """
        return A.transp()

    def submatrix(A, i, k, m, n):
        """ <b>A.submatrix(i, k, m, n)</b><br/>
            <i>m*n Submatrix (rows  i..i+m-1, columns k..k+n-1)</i><br/>
        ----de----
            <b>A.submatrix(i, k, m, n)</b><br/>
            <i>m*n-Untermatrix (Zeilen i..i+m-1, Spalten k..k+n-1)</i><br/>
        """
        assert i+m-1 < A.height() and k+n-1 < A.width()
        return Matrix([ [A[i1,k1] for k1 in range(k,k+n)]
                        for i1 in range(i,i+m)])

    def complement(A, i, k):
        """ <b>A.complement(i, k)</b><br/>
            <i>Matrix without row i and column k</i><br/>
            (Algebraic complement)
        ----de----
            <b>A.complement(i, k)</b><br/>
            <i>Matrix ohne die Zeile i und die Spalte k</i><br/>
            (Algebraisches Komplement)
        """
        return Matrix([ [A[i1,k1] for k1 in A.colRange() if k1 != k]
                        for i1 in A.rowRange() if i1 != i])

    def minor(A, i, k):
        """ <b>A.minor(i, k)</b><br/>
            <i>Minor for row i, column k</i><br/>
            (Determinant of the algebraic complement)
        ----de----
            <b>A.minor(i, k)</b><br/>
            <i>Minor zur Zeile i und der Spalte k</i><br/>
            (Determinante des Algebraischen Komplements)
        """
        return A.complement(i,k).det()

    def cofactor(A, i, k):
        """ <b>A.cofactor(i, k)</b><br/>
            <i>Cofactor for row i, columns k</i><br/>
            (signed minor)
        ----de----
            <b>A.cofactor(i, k)</b><br/>
            <i>Kofaktor zur Zeile i und der Spalte k</i><br/>
            (Minor mit Vorzeichen)
        """
        return (-1)**(i+k) * A.minor(i,k)

    def adjoint(A):
        """ <b>A.adjoint()</b><br/>
            <i>Adjoint matrix</i><br/>
        ----de----
            <b>A.adjoint()</b><br/>
            <i>Adjungierte Matrix</i><br/>
        """
        assert A.isSquare()
        return Matrix([ [A.cofactor(k,i) for k in A.rowRange()]
                        for i in A.colRange()])
    def det(A):
        """ <b>A.det()</b><br/>
            <i>Determinant of the matrix A</i><br/>
        ----de----
            <b>A.det()</b><br/>
            <i>Determinante der Matrix A</i><br/>
        """
        if not A.isSquare():
            raise _MathGUIdeError("A is not a square matrix")
        n = A.height()
        if n == 1:
            return A[0,0]
        else:
            return sum([A[0,k] * A.cofactor(0,k)
                        for k in A.colRange()])

    def _gaussElim0(A, jordan=True):
        """ unsichere Vorstufe zu gaussElim
        """
        m, n = A.height(), A.width()
        k = 0
        for i0 in A.rowRange():
            k = i0
            # oberste Zeile für führende Eins passend durchmultiplizieren
            A[i0] *= 1/A[i0,k]

            # Passende Vielfache zu anderen Zeilen addieren, so dass
            # unterhalb der führenden Eins Nullen entstehen
            for i in range(i0+1, m):
                A[i] -= A[i0] * A[i,k]
            if jordan:
                # Gauß-Jordan: dto. auch oberhalb der führenden Eins
                for i in range(i0):
                    A[i] -= A[i0] * A[i,k]

    def gaussElim(self, jordan=True):
        """ <b>A.gaussElim(jordan=True)</b><br/>
            <i>Row echelon form of the matrix A</i><br/>
            If <var>jordan</var> is <code>True</code>, the reduced row echelon form
            is returned (Gaussian elimination), otherwise the normal row echelon form.
        ----de----
            <b>A.gaussElim(jordan=True)</b><br/>
            <i>Zeilenstufenform der Matrix A.</i><br/>
            Wenn <var>jordan</var> <code>True</code> ist, wird die reduzierte Zeilenstufenform
            (Gauß-Jordan-Elimination) zurückgegeben, sonst die einfache Zeilenstufenform (Gauß-Elimination).
        """
        A = self.copy()
        m, n = A.height(), A.width()
        k = -1
        for i0 in A.rowRange():
            # Bestimme die am weitesten links stehende Spalte k, die
            # (ab Zeile i0) von Null verschiedene Elemente enthält:
            k += 1
            aMax, iMax = 0, i0
            while aMax == 0 and k < n:
                for i in range(i0,m):
                    if abs(A[i,k]) > aMax:
                        aMax, iMax = abs(A[i,k]), i
                if aMax == 0:
                    k += 1
            if k < n:
                # Oberste Zeile mit der vertauschen, die das (dem Betrag nach)
                # größte Element in Spalte i0 enthält
                A[iMax], A[i0] = A[i0], A[iMax]

                # oberste Zeile für führende Eins passend durchmultiplizieren
                A[i0] *= 1/A[i0,k]

                # Passende Vielfache zu anderen Zeilen addieren, so dass
                # unterhalb der führenden Eins Nullen entstehen
                for i in range(i0+1, m):
                    A[i] -= A[i0] * A[i,k]
                if jordan:
                    # Gauß-Jordan: dto. auch oberhalb der führenden Eins
                    for i in range(i0):
                        A[i] -= A[i0] * A[i,k]
        return A

    def rank(A):
        """ <b>A.rank()</b><br/>
            <i>Rank of the matrix A</i><br/>
        ----de----
            <b>A.rank()</b><br/>
            <i>Rang der Matrix A</i><br/>
        """
        B = A.gaussElim()
        r = 0
        while r < B.height() and B[r].norm() != 0:
            r += 1
        return r

    def inverse(A):
        """ <b>A.inverse()</b><br/>
            <i>Inverse og the matrix A</i><br/>
        ----de----
            <b>A.inverse()</b><br/>
            <i>Inverse Matrix A</i><br/>
        """
        if not A.isSquare():
            raise _MathGUIdeError("A is not a square matrix")
        n = A.height()
        I = Matrix.identity(n)
        B = A.joinRight(I).gaussElim()
        if B.submatrix(0,0,n,n) != I:
            raise _MathGUIdeError("Matrix is not invertible")
        return B.submatrix(0,n,n,n)

    def gramSchmidt(self, normalize=False):
        """ <b>A.gramSchmidt(normalize=False)</b><br/>
            <i>Orthogonalizing using the Gram–Schmidt process</i><br/>
            <code>A</code> must be the matrix of the (linearly independent) vectors
            to be orthogonalized (given as rows).<br/>
            If <var>normalize</var> is True, the vectors are normalized.<br/>
            return value: Matrix consisting of the orthogonalen vectors.<br/>
            <b>Example</b>:<br/><pre>
     A = Matrix([[1,1],[2,0]])
     A.gramSchmidt()</pre>
            returns the matrix<pre>
     / 1   1 \
     |       |
     \ 1  -1 /</pre>
        ----de----
            <b>A.gramSchmidt(normalize=False)</b><br/>
            <i>Gram-Schmidt'sches Orthogonalisierungsverfahren.</i><br/>
            <code>A</code> muss Matrix aus den (linear unabhängigen)
            zu orthogonalisierenden Vektoren (als Zeilen) sein.<br/>
            Wenn <var>normalize</var> True ist, werden die Vektoren normalisiert.<br/>
            Rückgabewert: Matrix mit den orthogonalen Vektoren.<br/>
            <b>Beispiel</b><pre>
     A = Matrix([[1,1],[2,0]])
     A.gramSchmidt()</pre>
            liefert die Matrix<pre>
     / 1   1 \
     |       |
     \ 1  -1 /</pre>
        """
        A = self.copy()
        n = A.height()
        for i in range(1,n):
            for k in range(i):
                A[i] -= A[k] * ((A[i]*A[k])/(A[k]*A[k]))
        if normalize:
            for v in A:
                v *= 1/v.norm()
        return A

    def _luRecursive(A):
        """ Vorstufe: funktioniert nicht für jede reguläre Matrix!
        """
        assert A.isSquare()
        n = A.height()
        if n == 1: # Rekursionsabbruch: 1·1-Matrix
            return Matrix.identity(1), A.copy()

        #     /a11 :  w  \   / 1      :  0 \   /a11:     w       \
        #     |..........|   |.............|   |................ |
        # A = |    :     | = |  1     :    | * |   :      1      |
        #     | ~v :  A1 |   | --- ~v :  I |   | 0 : A1- --- ~v*w|
        #     \    :     /   \ a11    :    /   \   :     a11     /
        v = A.submatrix(1,0, n-1,1)
        w = A.submatrix(0,1, 1,n-1)
        A1 = A.submatrix(1,1, n-1,n-1)
        B = A1 - (1/A[0,0]) * v * w
        L1, U1 = B._luRecursive() # Rekursiver Aufruf für (n-1)·(n-1)-Matrix
        L_top = Matrix.identity(1).joinRight(Matrix.null(1,n-1))
        L_btm = ((1/A[0,0]) * v).joinRight(L1)
        L = L_top.joinBottom(L_btm)
        U_top = A.submatrix(0,0, 1,1).joinRight(w)
        U_btm = Matrix.null(n-1,1).joinRight(U1)
        U = U_top.joinBottom(U_btm)
        return L,U

    def _luIter(self):
        A = self.copy()
        assert A.isSquare()
        n = A.height()
        for k in range(n-1):
            assert A[k,k] != 0 # sonst scheitert der Algorithmus!
            for i in range(k+1,n):
                A[i,k] /= A[k,k]
                for j in range(k+1,n):
                    A[i,j] -= A[i,k] * A[k,j]
        # A in linke und rechte Dreiecksmatrix zerlegen
        L = Matrix.identity(n)
        U = Matrix.identity(n)
        for i in range(n):
            for k in range(n):
                if k < i:  L[i,k] = A[i,k]
                else:      U[i,k] = A[i,k]
        return L, U

    def lup(self):
        """ <b>A.lup()</b><br/>
            <i>LUP decomposition of the matrix</i><br/>
            Returns a triplet (L, U, P) with P*A = L*R, where:<br/>
            P is a Permutation matrix,<br/>
            L is a Lower triangular matrix where all diagonal elements are one.<br/>
            R is an upper triangular matrix.<br/>
            <bRequires</b>: Matrix ist invertible.
        ----de----
        <b>.lup()</b><br/>
            <i>LRP-Zerlegung (LUP decomposition) der Matrix</i><br/>
            Rückgabewert: Tripel (L, U, P) mit P*A = L*R, wobei:<br/>
            P Permutationsmatrix,<br/>
            L linke (lower) Dreiecksmatrix mit Einsen auf Diagonale,<br/>
            U rechte (upper) Dreiecksmatrix ist.<br/>
            <b>Voraussetzung</b>: Matrix ist invertierbar.
        """
        A = self.copy()
        assert A.isSquare()
        n = A.height()
        P = Matrix.identity(n)
        for k in range(n-1):
            # Suche das dem Betrag nach größte Element
            # der k-ten Spalte ab der k-ten Zeile und seinen Zeilenindex
            absMax, kMax = max([(abs(A[i,k]), -i) for i in range(k,n)])
            kMax = -kMax # -i im Paar --> bei mehreren gleichen erstes
            assert absMax > 0

            # Vertausche die Zeile des größten Elements
            # mit der k-ten Zeile (Permutation und Matrix)
            P[k], P[kMax] = P[kMax], P[k]
            A[k], A[kMax] = A[kMax], A[k]

            for i in range(k+1,n):
                A[i,k] /= A[k,k]
                for j in range(k+1,n):
                    A[i,j] -= A[i,k] * A[k,j]

        # A in linke und rechte Dreiecksmatrix zerlegen
        L = Matrix.identity(n)
        U = Matrix.identity(n)
        for i in range(n):
            for k in range(n):
                if k < i:  L[i,k] = A[i,k]
                else:      U[i,k] = A[i,k]
        return L, U, P

    def solve(A, b):
        """ <b>A.solve(b)</b><br/>
            <i>Solution of the linear equation system Ax = b</i><br/>
            b may be Vektor or column Matrix.
        ----de----
            <b>A.solve(b)</b><br/>
            <i>Lösung x des linearen Gleichungssystems Ax = b</i><br/>
            b kann Vektor oder Spaltenmatrix sein.
        """
        if not isinstance(b, Matrix):
            assert isinstance(b, Vector)
            b = ~b
        return A.inverse() * b

    def leastSquares(A, b):
        """ <b>A.leastSquares(b)</b><br/>
            <i>Solution of the normal system</i><br/>
            Returns: <code>(~A*A).solve(~A*b)</code>
        ----de----
            <b>A.leastSquares(b)</b><br/>
            <i>Lösung der Normalgleichung</i><br/>
            Rückgabewert: <code>(~A*A).solve(~A*b)</code>
        """
        return (~A*A).solve(~A*b)


def fit(data, values):
    """ <b>fit(data, values)</b><br/>
        <i>Interpolation using the Gaussian method of least squares</i><br/>
        <code>fit(data, values)</code>
        calculates (using the method of least squares)
        a linear combination of the functions mapping x
        to the terms in <code>values</code>.<br/>
        (corresponds to the Mathematica function <code>Fit</code>).<br/>
        <b>Example</b>:<pre>
 fit([(-1,2), (1,1), (2,1), (3,0), (5,3)],
     ["1", "x", "x**2"])</pre>
        calculates <code>[c1, c2, c3]</code> such that the function<br/>
        <code>f(x) = c1*1  + c2*x + c3 * x**2</code><br/>
        approximates the conditions f(-1)=2, f(1)=1, f(2)=1, f(3)=0, f(5)=3 .<br/>
        Result: "(6/5 * 1) + (-53/70 * x) + (3/14 * x**2)"
        ----de----
        <b>fit(data, values)</b><br/>
        <i>Interpolation mit der Gaußschen Methode der kleinsten Quadrate</i><br/>
        <code>fit(data, values)</code>
        berechnet nach der Methode der kleinsten Quadrate
        eine Linearkombination der Funktionen, die x auf die
        Terme in <code>values</code> abbilden.<br/>
        (entspricht der Mathematica-Funktion <code>Fit</code>).<br/>
        <b>Beispiel</b>:<pre>
 fit([(-1,2), (1,1), (2,1), (3,0), (5,3)],
     ["1", "x", "x**2"])</pre>
        berechnet <code>[c1, c2, c3]</code> so dass die Funktion<br/>
        <code>f(x) = c1*1  + c2*x + c3 * x**2</code><br/>
        die Bedingungen f(-1)=2, f(1)=1, f(2)=1, f(3)=0, f(5)=3 approximiert.<br/>
        Ergebnis: "(6/5 * 1) + (-53/70 * x) + (3/14 * x**2)"
    """
                              # from the example:
    m = len(data)             # m = 5
    n = len(values)           # n = 3
    functions = [eval("lambda x:{}".format(val)) for val in values]
    A = Matrix.fromFunction(m,n, lambda i,k: functions[k](data[i][0]))
    b = ~Vector.fromFunction(m, lambda i: data[i][1])
                              #         / 1  -1   1 \         /  2  \
                              #        |  1   1   1  |        |  1  |
                              # A =    |  1   2   4  |    b = |  1  |
                              #        |  1   3   9  |        |  0  |
                              #         \ 1   5  25 /         \  3  /
    y = A.leastSquares(b)     # Solution of the normal system (column matrix):
    s = ""                    # y = [[6/5], [-53/70], [3/14]]
    for i in y.rowRange():    # Create linear combination:
        s += "({0} * {1})".format(str(y[i,0]), values[i])
        if i < y.height()-1:
            s += " + "
    return s                  # "(6/5 * 1) + (-53/70 * x) + (3/14 * x**2)"

#---------------------------------------------------------------
#  RSA-Chiffrierung

class Rsa: # RSA Encryption # RSA-Verschlüsselung
    """ Class <code>Rsa</code>: RSA Encryption<br/>
        This is a static class (only static methods available).
    ----de----
        Klasse <code>Rsa</code>: RSA-Verschlüsselung<br/>
        Statische Klasse: nur Klassenmethoden.
    """

    _staticMethods = ("createKeys", "encrypt")

    @staticmethod
    def createKeys(bits=768):
        """ <b>createKeys(bits=768)</b><br/>
            <i>returns an RSA key pair</i><br/>
            The pair is returned as tuple <code>((e,n),(d,n))</code>.<br/>
            <b>Example</b>:<br/>
            <code>public, private = Rsa.createKeys(24)<br/>
                  m = 12345 # message<br/>
                  c = Rsa.encrypt(m, public) # encrypted message<br/>
                  assert Rsa.encrypt(c, private) == m</code><br/>
            <b>See also</b> <a href="#Rsa.encrypt"><code>Rsa.encrypt</code></a>.
        ----de----
            <b>createKeys(bits=768)</b><br/>
            <i>liefert ein RSA-Schlüsselpaar</i><br/>
            Das Paar wird als Tupel <code>((e,n),(d,n))</code> zurückgegeben.<br/>
            <b>Beispiel</b>:<br/>
            <code>public, private = Rsa.createKeys(24)<br/>
                  m = 12345 # Nachricht<br/>
                  c = Rsa.encrypt(m, public) # verschlüsselte Nachricht<br/>
                  assert Rsa.encrypt(c, private) == m</code><br/>
            <b>Siehe auch</b> <a href="#Rsa.encrypt"><code>Rsa.encrypt</code></a>.
        """
        m = 2 ** (bits//2 - 1)
        # wähle zwei Primzahlen p, q
        p = nextPrime(m + rand(m))
        q = nextPrime(m + rand(m))
        n = p * q
        phi = (p-1)*(q-1)       # phi(n) (Eulersche Phi-Funktion)
        d = phi                 # wähle d aus 1..phi(n)
        while gcd(d,phi) != 1:  # mit d, phi(n) teilerfremd
            d = rand(phi)
        t,x,y = gcdExt(d,phi)
        e = x % phi
        publicKey = (e,n)
        privateKey = (d,n)
        return (publicKey, privateKey)

    @staticmethod
    def encrypt(m, key):
        """ <b>encrypt(m, key)</b><br/>
            <i>RSA encryption</i><br/>
            encrypts m using the mit dem key <code>key = (e,n)</code><br/>
            <b>See also</b> <a href="#Rsa.createKeys"><code>Rsa.createKeys</code></a>.
        ----de----
            <b>encrypt(m, key)</b><br/>
            <i>RSA-Verschlüsselung</i><br/>
            Verschlüsselt m mit dem Schlüssel key = (e,n)<br/>
            <b>Siehe auch</b> <a href="#Rsa.createKeys"><code>Rsa.createKeys</code></a>.
        """
        return pow(m, key[0], key[1])


#===============================================================
#  Permutationen, Permutationsgruppen
#---------------------------------------------------------------

def permute(v, fn, prefix=[]):
    """ <b>permute(v, fn)</b><br/>
        <i>applies procedure <code>fn</code> to all permutations of liste v.</i><br/>
        If v is sorted, the permutations are also generated in sorted order.
        ----de----
        <b>permute(v, fn)</b><br/>
        <i>wendet Prozedur <code>fn</code> auf alle Permutationen von Liste v an.</i><br/>
        Wenn v sortiert ist, werden auch die Permutationen in
        sortierter Reihenfolge erzeugt.
    """
    if len(v) == 0: # nichts zu permutieren!
        fn(prefix)
    else:
        for i in v:
            p = prefix + [i]  # i an Praefix anhängen
            w = list(v)       # Duplikat von w, ...
            w.remove(i)       # ... ohne i
            permute(w, fn, p)

def permutations(v, prefix=[]):
    """ <b>permutations(v)</b><br/>
        <i>List of all permutations of v</i><br/>
        If v is sorted, the permutations are also generated in sorted order.
        ----de----
        <b>permutations(v)</b><br/>
        <i>Liste aller Permutationen von v</i><br/>
        Wenn v sortiert ist, werden auch die Permutationen in
        sortierter Reihenfolge erzeugt.
    """
    if len(v) == 0: # nichts zu permutieren!
        yield prefix
    else:
        for i in v:
            p = prefix + [i]  # i an Praefix anhängen
            w = list(v)       # Duplikat von w, ...
            w.remove(i)       # ... ohne i
            for i in permutations(w, p):
                yield i


class Perm: # Permutations (from dictionary or cycle representation) # Permutationen (aus Dictionary oder Zyklendarstellung)
    """ Class <code>Perm</code><br/>
        <i>Permutations (from dictionary or cycle representation)</i>
        <b>Examples</b> (three times the same permutation:<br/>
        <code>Perm({1:4, 2:3, 3:2, 4:1})</code><br/>
        <code>Perm("(1,4)(2,3)")</code><br/>
        <code>Perm([[1,4],[2,3]])</code>
        ----de----
        Klasse <code>Perm</code><br/>
        <i>Permutationen (aus Dictionary oder Zyklendarstellung)</i>
        <b>Beispiele</b> (dreimal die gleiche Permutation:<br/>
        <code>Perm({1:4, 2:3, 3:2, 4:1})</code><br/>
        <code>Perm("(1,4)(2,3)")</code><br/>
        <code>Perm([[1,4],[2,3]])</code>
    """
    _staticMethods = []

    def __init__(self, p):
        """ Constructor
        ----de----
            Konstruktor
        """
        if type(p) == dict:
            # Zuordnungstafel
            # Beispiel: Perm({1:4, 2:3, 3:2, 4:1})
            self._d = p
        elif type(p) == str:
            # Zyklendarstellung als String
            # Der String darf nur Zahlen enthalten
            # Beispiel: Perm("(1,4)(2,3)")
            v = []
            p = p.replace(" ", "") # Leerzeichen entfernen
            assert p[0] == "(" and p[-1] == ")"
            p = p[1:-1] # Klammern vorn und hinten entfernen
            for cycle in p.split(")("):
                v.append([eval(x) for x in cycle.split(",")])
            self.fromVector(v)
        else:
            # Zyklendarstellung
            # Beispiel: Perm([[1,4],[2,3]])
            self.fromVector(p)
        self.reduce()

    def fromVector(self, v):
        """ <b>p.fromVector(v)</b><br/>
            <i>Generation from vector v with argument/value pairs.</i><br/>
            Example: <code>p.fromVector([[1,4],[2,3]])</code>
        ----de----
            <b>p.fromVector(v)</b><br/>
            <i>Erzeugung aus Vektor v mit Argument-Wert-Paaren.</i><br/>
            Beispiel: <code>p.fromVector([[1,4],[2,3]])</code>
        """
        assert type(v) == list
        self._d = {}
        for cycle in v:
            assert type(cycle) == list
            for i in range(len(cycle)):
                if i < len(cycle) - 1:
                    self._d[cycle[i]] = cycle[i+1]
                else:
                    self._d[cycle[i]] = cycle[0]

    def __repr__(self):
        """ Object representation as string (as cycles)
        """
        if len(self._d) == 0:
            return "identity"
        else:
            return self.cycles()

    def __getitem__(self, i):
        return self._d.get(i, i)

    def __mul__(p1, p2):
        """ product of two Permutations (left first, then right!)
        """
        p = {}
        m = p1.members()
        for i in m:
            p[i] = p2[p1[i]]
        for i in p2.members():
            if i not in m:
                p[i] = p2[i]
        return Perm(p)

    def __eq__(p1,p2): return p1.cycles() == p2.cycles()
    def __ne__(p1,p2): return p1.cycles() != p2.cycles()

    def members(self):
        """ <b>p.members()</b><br/>
            <i>List of the permuted elements</i>
        ----de----
            <b>p.members()</b><br/>
            <i>Liste der permutierten Elemente</i>
        """
        m = list(self._d.keys())
        m.sort()
        return m

    def inverse(self):
        """ <b>p.inverse()</b><br/>
            <i>inverse permutation</i>
        ----de----
            <b>p.inverse()</b><br/>
            <i>Umkehrpermutation</i>
        """
        d = {}
        for i in self._d:
            d[self._d[i]] = i
        return Perm(d)

    def cycles(self):
        """ <b>p.cycles()</b><br/>
            <i>Cycle representation</i>
        ----de----
            <b>p.cycles()</b><br/>
            <i>Zyklendarstellung</i>
        """
        s = ""
        m = self.members()
        while len(m) > 0:
            first = m[0]
            m.remove(first)
            i = first
            s += "({0}".format(str(first))
            while self._d[i] != first:
                i = self._d[i]
                m.remove(i)
                s += ",{0}".format(str(i))
            s += ")"
        return s

    def check(self):
        """ <b>p.check()</b><br/>
            <i>check for consistency</i>
        ----de----
            <b>p.check()</b><br/>
            <i>Konsistenzprüfung</i>
        """
        m = self.members()
        v = self._d.values()
        v.sort()
        assert(k == v)

    def reduce(self):
        """ <b>p.reduce()</b><br/>
            <i>removes fixed points</i>
        ----de----
            <b>p.reduce()</b><br/>
            <i>entfernt Fixpunkte</i>
        """
        id = []
        for i in self._d:
            if self._d[i] == i:
                id.append(i)
        for i in id:
            del self._d[i]

    def inversions(self):
        """ <b>p.inversions()</b><br/>
            <i>Inversion number (count of "swapped pairs")</i>
        ----de----
            <b>p.inversions()</b><br/>
            <i>Inversionszahl (Anzahl der "vertauschten Paare")</i>
        """
        m = self.members()
        s = 0
        n = len(m)
        for i in range(n):
            for j in range(i+1, n):
                if self[m[i]] > self[m[j]]:
                    s += 1
        return s

    def sign(self):
        """ <b>p.sign()</b><br/>
            <i>signum</i><br/>
            1 for even, -1 for odd permutations
        ----de----
            <b>p.sign()</b><br/>
            <i>Signum</i><br/>
            1 bzw. -1 für gerade bzw. ungerade Permutation
        """
        return 1 - 2*(self.inversions() % 2)

    def odd(self):
        """ <b>p.odd()</b><br/>
            <i>odd permutation?</i> (bool)
        ----de----
            <b>p.odd()</b><br/>
            <i>ungerade Permutation?</i> (bool)
        """
        return self.inversions() % 2 == 1


#===============================================================
# Modulo-Arithmetik und Rechnen in Polynomringen
#---------------------------------------------------------------

class Mod: # Elements of the residue class rings Z[m] # Elemente der Restklassenringe Z[m]
    """ Class <code>Mod</code><br/>
        <i>Elements of the residue class rings Z[m]</i><br/>
        <code>Mod(n, m)</code> returns element n des of the residue class ring Z[m]<br/>
        Arithmetic operators may be applied to the objects.<br/>
        Some operations are defined only,
        if Z[m] is a field (&lt;==> m is prime number)<br/>
        <b>See also</b> <a href="mathguide.html#mod"><code>Mod</code></a>
        ----de----
        Klasse <code>Mod</code><br/>
        <i>Elemente der Restklassenringe Z[m]</i><br/>
        <code>Mod(n, m)</code> liefert Element n des Restklassenrings Z[m]<br/>
        Die Objekte können mit arithmetischen Operatoren verknüpft werden.<br/>
        Einige Operationen sind nur definiert,
        wenn Z[m] ein Körper ist (&lt;==> m ist Primzahl)<br/>
        <b>Siehe auch</b> <a href="mathguide.html#mod"><code>Mod</code></a>
    """
    _staticMethods = []

    def operators(self):
        """ <b>Mod.operators()</b><br/>
            <i>For documentation only</i><br/>
            The following operators are defined in the class <b>Mod</b>:<br/>
            <table>
             <tr><th>Op.</th><th>Function</th><th>Examples</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Subtraction</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>Multiplication</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>/</th><td>Division</td><td><code>a / b; a /= b</code></td></tr>
             <tr><th>%</th><td>Modulo (in fields only , i.e. modulus is prime)</td><td><code>a &amp; b; a &amp;= b</code></td></tr>
             <tr><th>-</th><td>Unary minus</td><td><code>-a</code></td></tr>
            </table>
            <font color="#000080" size="-2">This method is for documentation only. It has no effect.</font>
        ----de----
            <b>Mod.operators()</b><br/>
            <i>Nur zur Dokumentation</i><br/>
            Die Klasse <b>Mod</b> erlaubt folgende Operatoren:<br/>
            <table>
             <tr><th>Op.</th><th>Funktion</th><th>Beispiele</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Subtraktion</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>Multiplikation</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>/</th><td>Division (nur in Körpern, d.h. Primzahlmodul)</td><td><code>a / b; a /= b</code></td></tr>
             <tr><th>%</th><td>Modulo (nur in Körpern, d.h. Primzahlmodul)</td><td><code>a &amp; b; a &amp;= b</code></td></tr>
             <tr><th>-</th><td>unäres Minus</td><td>-a</td></tr>
            </table>
            <font color="#000080" size="-2">Diese Methode dient nur zur Dokumentation. Sie hat keine Wirkung.</font>
        """
        pass

    def __init__(self, n, m):
        """ Constructor
        ----de----
            Konstruktor
        """
        #assert isInteger(n) and isNatural(m)
        self.n = n % m
        self.m = m

    def copy(self):
        """ <b>A.copy()</b><br/>
            <i>independent copy ("clone") of A</i>
        ----de----
            <b>A.copy()</b><br/>
            <i>unabhängige Kopie von A</i>
        """
        return Mod(self.n, self.m)

    def __int__(self):
        return self.n

    def __repr__(self):
        """ Object representation as string
        """
        return "Mod({0}, {1})".format(str(self.n), str(self.m))

    def __cmp__(A, B):
        if not isinstance(B, Mod):
            B = Mod(B, A.m)
        return _cmp(A.n, B.n)

    def inverse(self):
        """ <b>p.inverse()</b><br/>
            <i>Multiplicative inverse</i><br/>
            defined in fields only. Therefore modulus must be prime number!
        ----de----
            <b>p.inverse()</b><br/>
            <i>Multiplikatives Inverses</i><br/>
            Nur in Körpern definiert. Deshalb muss der Modul Primzahl sein!
        """
        def gcdExt(a,b):
            if b == 0:
                return a, 1, 0
            else:
                q, r = divmod(a, b)
                d, z, x = gcdExt(b, r)
                return d, x, z - x * q

        d, inv, y = gcdExt(self.n, self.m)
        assert d == 1
        # inv * n + y * m == 1 ==> (inv * n) % m == 1
        assert (inv * self.n) % self.m == 1
        return Mod(inv, self.m)

    def __eq__(a,b):
        """ Test for equality (operator ==)
            with other Mod object or 0
        """
        if type(b) == int and b == 0:
            return a.n == 0
        return (a.n == b.n and a.m  == b.m)

    def __abs__(self):
        """ Absolute value
        """
        return self

    def __ne__(a,b):
        """ Test for unequality (Operator !=)
            with other Mod object or 0
        """
        return not a.__eq__(b)

    def __neg__(self):
        """ Unary Operator -
            -a --> a.__neg__()
        """
        return Mod((self.m - self.n) % self.m, self.m)

    def __add__(a, b):
        """ <b>a.__add__(b)</b><br/>
            <i>Addition</i><br/>
            Also operator notation: <b>a + b</b><br/>
            Instead of <code>a = a + b</code> you can write <code>a += b</code>.
        """
        #   a + b  --> a.__add__(b), wenn a vom Typ Mod ist
        #              b.__radd__(a), sonst
        if isinstance(b, Mod):
            d = gcd(a.m, b.m)
            return Mod(a.n + b.n, d)
        else:
            return Mod(a.n + b, a.m)

    def __radd__(b, a):
        """ Addition, see __add__ """
        return b.__add__(a)

    def __sub__(a, b):
        """ Subtraction:
            a - b  --> a.__sub__(b), if a is of type Mod
                       b.__rsub__(a), otherwise
        """
        if isinstance(b, Mod):
            d = gcd(a.m, b.m)
            return Mod(a.n - b.n, d)
        else:
            return Mod(a.n - b, a.m)

    def __rsub__(b, a):
        """ Subtraction, see __sub__
        """
        return Mod(a - b.m, b.m)

    def __mul__(a, b):
        """ Multiplication:
            a * b  --> a.__mul__(b), if a is of type Mod
                       b.__rmul__(a), otherwise
        """
        if isinstance(b, Mod):
            d = gcd(a.m, b.m)
            return Mod(a.n * b.n, d)
        else:
            return Mod(a.n * b, a.m)

    def __rmul__(b, a):
        """ Multiplikation, see __mul__
        """
        return b.__mul__(a)

    def __truediv__(a, b):
        """ exact division:
            a / b  --> a.__div__(b), if a is of type Mod
                       b.__rdiv__(a), otherwise
        """
        if not isinstance(b, Mod):
            b = Mod(b, a.m)
        assert a.m == b.m
        return a * b.inverse()

    def __rtruediv__(a, b):
        """ exact division, see __div__
        """
        assert not isinstance(b, Mod)
        b = Mod(b, a.m)
        return b.__div__(a)

    def __div__(a, b): return a.__truediv__(b)
    def __rdiv__(a, b): return a.__rtruediv__(b)

#===============================================================
# Poly
#---------------------------------------------------------------

class Poly: # Polynomials over arbitrary rings # Polynome über beliebigen Ringen
    """ Class <code>Poly</code><br/>
        <i>Polynomials</i><br/>
        <code>Poly(list)</code> returns a polynomial with the coefficients from <code>list</code> (ascending).<br/>
        The coefficients must comply with the ring axioms,<br/>
        for division and modulo operator (//, %, divmod) also the field axioms,<br/>
        and implement the equality with 0.<br/>
        Short form for polynomial with Mod objektes:<br/>
        <b>Example</b><br/>
        <code>Poly([Mod(1,2),Mod(0,2),Mod(1,2)])</code><br/>
        Short form:<br/>
        <code>Poly([1,0,1,1],2)</code>
        ----de----
        Klasse <code>Poly</code><br/>
        <i>Polynome</i><br/>
        <code>Poly(list)</code> liefert Polynom mit den Koeffizienten aus <code>list</code> (aufsteigend).<br/>
        Die Koeffizienten müssen die Ringaxiome erfüllen,<br/>
        für Division und Modulo-Operator (//, %, divmod) auch Körperaxiome,<br/>
        und die Abfrage auf Gleichheit mit 0 implementieren.<br/>
        Kurzform für Polynom mit Mod-Objekten:<br/>
        <b>Beispiel</b><br/>
        <code>Poly([Mod(1,2),Mod(0,2),Mod(1,2)])</code><br/>
        Dafür auch Kurzschreibweise:<br/>
        <code>Poly([1,0,1,1],2)</code>
    """
    #............................................................
    _staticMethods = ("gcd", "gcdExt", "irreducibles")

    @staticmethod
    def gcd(a,b):
        """ <b>gcd(a,b)</b><br/>
            <i>greatest common divisor of a and b</i><br/>
            Calculated using the Euclidean algorithm
        ----de----
            <b>gcd(a,b)</b><br/>
            <i>Größter gemeinsamer Teiler von a und b</i><br/>
            Berechnung mit dem Euklidischen Algorithmus
        """
        p = a.copy()
        q = b.copy()
        while not q.isZero():
            p, q = q, p % q
        return p.normalized()

    @staticmethod
    def gcdExt(a,b):
        """ <b>gcdExt(a,b)</b><br/>
            <i>Extended Euclidean algorithm.</i><br/>
            Returns (d, x, y) with d == gcd(a,b) and x*a + y*b == d
        ----de----
            <b>gcdExt(a,b)</b><br/>
            <i>Erweiterter Euklidischer Algorithmus.</i><br/>
            Gibt (d, x, y) zurück mit d == gcd(a,b) und  x*a + y*b == d
        """
        assert type(a[0]) == Mod
        m = a[0].m
        if b.isZero():
            return a.copy(), Poly([1],m), Poly([0],m)
        else:
            q, r = divmod(a, b)
            d, z, x = Poly.gcdExt(b, r)
            y = z - x * q
            return d, x, y

    @staticmethod
    def irreducibles(maxDegree, mod=2):
        """ <b>Poly.irreducibles(maxDegree, mod=2)</b><br/>
            <i>Irreducible polynomials</i><br/>
            Returns: List of irreducible polynomials
            over Z[mod] up to maxDegree
        ----de----
            <b>Poly.irreducibles(maxDegree, mod=2)</b><br/>
            <i>Irreduzible Polynome</i><br/>
            Rückgabewert: Liste der irreduziblen Polynome
            über Z[mod] bis zum Grad maxDegree
        """
        primeList = [Poly([0,1],mod)]
        for i in fromTo(mod+1, mod**(maxDegree+1)-1):
            f = Poly(NumeralSystem.toBase(i, mod, True), mod)
            divisible = False
            for p in primeList:
                if f % p == 0:
                    divisible = True
                    break
            if not divisible:
                primeList.append(f)
        return primeList
    #............................................................

    def operators(self):
        """ <b>Poly.operators()</b><br/>
            <i>For documentation only</i><br/>
            The following operators are defined in the class <b>Poly</b>:<br/>
            <table>
             <tr><th>Op.</th><th>Function</th><th>Examples</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Subtraction</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>Multiplication</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>/</th><td>Polynomial division (with remainder)</td><td><code>a / b; a /= b</code></td></tr>
             <tr><th>divmod</th><td>Division and remainder simultaneously</td><td><code>d,m = divmod(a,b)</code></td></tr>
             <tr><th>-</th><td>Unary minus</td><td><code>-a</code></td></tr>
            </table>
            <font color="#000080" size="-2">This method is for documentation only. It has no effect.</font>
        ----de----
            <b>Poly.operators()</b><br/>
            <i>Nur zur Dokumentation</i><br/>
            Die Klasse <b>Poly</b> erlaubt folgende Operatoren:<br/>
            <table>
             <tr><th>Op.</th><th>Funktion</th><th>Beispiele</th></tr>
             <tr><th>+</th><td>Addition</td><td><code>a + b; a += b</code></td></tr>
             <tr><th>-</th><td>Subtraktion</td><td><code>a - b; a -= b</code></td></tr>
             <tr><th>*</th><td>Multiplikation</td><td><code>a * b; a *= b</code></td></tr>
             <tr><th>//</th><td>Polynomdivision (mit Rest)</td><td><code>a // b; a //= b</code></td></tr>
             <tr><th>divmod</th><td>Division und Modulo simultan</td><td><code>d,m = divmod(a,b)</code></td></tr>
             <tr><th>-</th><td>unäres Minus</td><td>-a</td></tr>
            </table>
            <font color="#000080" size="-2">Diese Methode dient nur zur Dokumentation. Sie hat keine Wirkung.</font>
        """
        pass
    def __init__(self, coeff=[], mod=0):
        """ Constructor
        ----de----
            Konstruktor
        """
        if mod != 0:
            self.__p = [Mod(i,mod) for i in coeff]
            if len(coeff) == 0: # zero polynomial
                self.__p = [Mod(0, mod)]
        else:
            # make all coefficients rational numbers
            self.__p = [Rational(c) for c in coeff]
            if len(coeff) == 0: # zero polynomial
                self.__p = [Rational(0)]

        # drop leading zero coefficients
        n = len(self.__p)
        while n > 0 and self[n-1] == 0:
            n -= 1
        self.__p = self.__p[0:n]

    def copy(self):
        """ <b>p.copy()</b><br/>
            <i>independent copy ("clone") of p</i>
        ----de----
            <b>p.copy()</b><br/>
            <i>unabhängige Kopie von p</i>
        """
        c = Poly()
        c.__p = [copy(x) for x in self.__p]
        return c

    def isZero(self):
        return (len(self.__p) == 0
            or len(self.__p) == 1 and self.__p[0] == 0)

    def isOne(self):
        if len(self.__p) != 1:
            return False
        if type(self.__p[0]) == Mod:
            return self.__p[0].n == 1
        return self.__p[0] == 1

    def degree(self):
        """ <b>p.degree()</b><br/>
            <i>Grad des Polynoms</i>
        ----de----
            <b>p.degree()</b><br/>
            <i>degree of the polynomial</i>
        """
        if self.isZero():
            return -1
        else:
            return len(self.__p) - 1

    def value(self, x):
        """ <b>p.value(x)</b><br/>
            <i>Value of the polynomial for argument x</i>
        ----de----
            <b>p.value(x)</b><br/>
            <i>Wert des Polynoms für Argument x</i>
        """
        val = 0
        n = self.degree()
        for i in fromTo(0, self.degree()):
            val = (x * val) + self.__p[n-i]
        return val

    def __normalize(self):
        """ reduce to true degree
        """
        n = len(self.__p)
        while n > 0 and self[n-1] == 0:
            n -= 1
        self.__p = self.__p[0:n]

    def toDigits(f):
        """ <b>p.toDigits()</b><br/>
            <i>String of coefficients</i><br/>
            Short form, <b>nur für singled digit</b> Mod coefficients:<br/>
            Concatenation of the coefficient digits (without mod).<br/>
            Highest coeffizient left!<br/>
            <b>Example</b><br/>
            <code>p=Poly([1,0,1,1],2)</code><br/>
            <code>p.toDigits()</code><br/>
            Result: <code>"1101"</code>
        ----de----
            <b>p.toDigits()</b><br/>
            <i>String aus Koeffizienten</i><br/>
            Kurzform, <b>nur für einstellige</b> Mod-Koeffizienten:<br/>
            Verkettung der Koeffizienten-Ziffern (ohne mod-Angabe).<br/>
            Höchster Koeffizient links!<br/>
            <b>Beispiel</b><br/>
            <code>p=Poly([1,0,1,1],2)</code><br/>
            <code>p.toDigits()</code><br/>
            Ergebnis: <code>"1101"</code>
        """
        assert isinstance(f[0], Mod) and f[0].n < 10
        l = [str(x.n) for x in f.__p]
        l.reverse()
        return join(l)

    def __repr__(f):
        """ Object representation as string
        """
        # prüfen, ob mind. ein Koeffizien vom Typ Mod ist:
        m = 0
        for x in f.__p:
            if isinstance(x, Mod):
                m = x.m
                break
        if m > 0:
            # Voraussetzung: Mod-Objekte kompatibel!
            return "Poly({0}, {1})".format(str([int(x) for x in f.__p]), m)
        else:
            return "Poly([{0}])".format(join(f.__p, ", "))

    def expand(f, x="x"):
        """ <b>p.expand(x="x")</b><br/>
            <i>Polynomial written for the variable x</i><br/>
        ----de----
            <b>p.expand(x="x")</b><br/>
            <i>ausgeschriebenes Polynom für die Variable x</i><br/>
        """
        def pot(i):
            if i == 0:
                return 1
            elif i == 1:
                return x
            else:
                return "{0}^{1}".format(x, i)
        l = ["({0} * {1})".format(str(f[i]), pot(i)) for i in fromTo(0, f.degree())]
        return join(l, " + ")

    def pretty(f, width=0):
        """ <b>p.pretty(width=0)</b><br/>
            <i>Human readable written representation</i><br/>
        ----de----
            <b>p.pretty(width=0)</b><br/>
            <i>Lesbare Darstellung des Polynoms</i><br/>
        """
        if f.isZero():
            pp = "0"
        else:
            x = "x"
            def pot(i):
                if   i == 0:  return ""
                elif i == 1:  return x
                elif i == 2:  return "{0}²".format(x)
                elif i == 3:  return "{0}³".format(x)
                else:         return "{0}^{1}".format(x, i)
            def sign(i):
                if f[i] < 0: return "- "
                elif i > 0: return "+ "
                else: return ""
            pp = ""
            for i in fromTo(f.degree(), 0, -1):
                # sign
                sign = ""
                if type(f[i]) != Mod and f[i] < 0: sign = " - "
                elif i < f.degree(): sign = " + "

                # abs value
                a = str(abs(f[i]))
                if type(f[i]) == Mod:
                    a = f[i].n
                if i > 0 and type(f[i]) == Mod and f[i].m == 2:
                    a = ""

                if f[i] != 0 and pot != "":
                    pp += "{0}{1}{2}".format(sign, a, pot(i))
        if len(pp) < width:
            width -= len(pp)
            pp = "{0}{1}{2}".format((width//2)*" ", pp, (width-width//2)*" ")
        return pp

    def __getitem__(f, i):
        """ Definition of the index operator []
            f[i] returns the i-th coefficient of the polynomial.
        """
        if i <= f.degree():
            return f.__p[i]
        else:
            return 0 * f.__p[i]

    #--- unäre Operatoren: ---

    def __neg__(f):
        """ unärer Operator -
        """
        r = Poly()
        r.__p = [-copy(x) for x in f.coeff]
        return r

    #--- binäre Operatoren: ---

    def __eq__(f,g):
        """ Test for equality (operator ==)
            with other polynomial or ring element
        """
        if not isinstance(g, Poly):
            g = Poly([g])
        return f.__p == g.__p

    def __ne__(f,g):
        """ Test for unequality (operator !=)
            with other polynomial or ring element
        """
        return not f.__eq__(g)

    def __add__(f, g):
        """ <b>a.__add__(b)</b><br/>
            <i>Addition</i><br/>
            Also operator notation: <b>a + b</b><br/>
            Instead of <code>a = a + b</code> you can write <code>a += b</code>.
        """
        #   a + b  --> a.__add__ (b), wenn a vom Typ Poly ist
        #              b.__radd__(a), sonst
        if not isinstance(g, Poly):
            g = Poly([g])
        fDeg, gDeg = f.degree(), g.degree()
        minDeg = min(fDeg, gDeg)
        maxDeg = max(fDeg, gDeg)
        s = Poly()
        s.__p = (maxDeg + 1) * [0] # list with (maxDeg + 1) zeroes
        for i in fromTo(0, maxDeg):
            if i > fDeg:
                s.__p[i] = copy(g.__p[i])
            elif i > gDeg:
                s.__p[i] = copy(f.__p[i])
            else:
                s.__p[i] = f.__p[i] + g.__p[i]
        s.__normalize()
        return s

    def __radd__(f, g):
        """ Addition, see __add__
        """
        return f + g

    def __sub__(f, g):
        """ Subtraction (operators - and -=):
            a - b  --> a.__sub__(b), if a is of type Poly
                       b.__rsub__(a), otherwise
        """
        if not isinstance(g, Poly):
            g = Poly([g])
        fDeg, gDeg = f.degree(), g.degree()
        minDeg = min(fDeg, gDeg)
        maxDeg = max(fDeg, gDeg)
        s = Poly()
        s.__p = (maxDeg + 1) * [0] # list with (maxDeg + 1) zeroes
        for i in fromTo(0, maxDeg):
            if i > fDeg:
                s.__p[i] = -g.__p[i]
            elif i > gDeg:
                s.__p[i] = copy(f.__p[i])
            else:
                s.__p[i] = f.__p[i] - g.__p[i]
        s.__normalize()
        return s

    def __rsub__(f, g):
        """ Subtraktion, see __sub__
        """
        return Poly([g]) - f

    def __mul__(f, g):
        """ Multiplication (operators * and *=):
            a * b  --> a.__mul__(b), if a is of type Poly
                       b.__rmul__(a), otherwise
        """
        if not isinstance(g, Poly):
            g = Poly([g])
        if f.isZero():
            return f.copy()
        if g.isZero():
            return g.copy()
        n, m = f.degree(), g.degree()

        p = Poly()
        p.__p = (n + m + 1) * [0] # list with (n + m + 1) zeroes
        for i in fromTo(0, n):
            for j in fromTo(0, m):
                p.__p[i+j] += f.__p[i] * g.__p[j]
        p.__normalize()
        return p

    def __rmul__(f, g):
        """ Multiplikation, see __mul__
        """
        return f * g

    def __rshift__(f, n):
        """ f >> n -- Multiplication of Polynomial by x**n
        """
        r = Poly()
        r.__p = n * [0] + [copy(x) for x in f.__p]
        return r

    def __divmod__(f, g):
        """ Polynomial division with remainder. (global function divmod)
            divmod(f, g) --> f.__divmod__(g), if f is of type Poly
                             g.__rdivmod__(f), otherwise

            Returns the pair (q,r)
            with:  f == g * q + r
            and:   r.degree() &lt; g.degree() or r == 0
        """
        if not isinstance(g, Poly):
            g = Poly([g])
        r = f.copy()   # Kopie: Zählerpolynom wird zu Restpolynom
        q = Poly()     # Divisionsergebnis: beginnen mit Nullpolynom
        while r.degree() >= g.degree() and not r.isZero():
            assert f == g * q + r
            qCoeff = r[r.degree()] / g[g.degree()]  # Quotient der hoechsten Koeff.
            qPoly = Poly()
            qPoly.__p = [qCoeff]
            qPoly >>= (r.degree() - g.degree())
            q += qPoly
            assert f == g * (q-qPoly) + r
            deg0 = r.degree()
            r -= qPoly * g # grad(r) muss jetzt kleiner werden
            assert r.degree() < deg0
            assert f == g * (q-qPoly) + (r + qPoly*g)
            assert f == g * q + r
        return q, r

    def __rdivmod__(f, g):
        """ Polynomial division with remainder. See __divmod__
        """
        assert not isinstance(g, Poly)
        g = Poly([g])
        return g.__divmod__(f)

    def __div__(f, g):
        """ Polynomial division with remainder (operators // und //=):
            f // g  --> f.__floordiv__ (g), if f is of type Poly
                        g.__rfloordiv__(f), otherwise
        """
        return divmod(f, g)[0]

    def __rdiv__(f, g):
        """ Polynomial division with remainder. See __floordiv__
        """
        assert not isinstance(g, Poly)
        g = Poly([g])
        return divmod(g, f)[0]

    def __truediv__(a, b):   return a.__div__(b)
    def __rtruediv__(a, b):  return a.__rdiv__(b)
    def __floordiv__(a, b):  return a.__div__(b)
    def __rfloordiv__(a, b): return a.__rdiv__(b)

    #-- Modulo-Operator ---
    # a & b  --> __mod__ (a, b), wenn a ein Poly-Objekt ist
    #            __rmod__(b, a), sonst

    def __mod__(f, d):
        """ Modulo operator for polynomials (operators % und %=):
            a % b  --> a.__mod__ (b), if f is of type Poly
                       b.__rmod__(a), otherwise
        """
        return divmod(f, d)[1]

    def __rmod__(f, g):
        """ Modulo operator for polynomials. See __mod__
        """
        assert not isinstance(g, Poly)
        g = Poly([g])
        return divmod(g, f)[1]

    def __pow__(a, x):
        assert isNatural(x)
        if x == 1:
            return a
        elif x % 2 == 1:
            return a * a**(x-1)
        else:
            t = a ** (x//2)
            return t * t

    def normalized(self):
        """ Divides the polynomial by the highest coefficient
        ----de----
            Dividiert das Polynom durch den höchsten Koeffizienten
        """
        if self.degree() >= 0:
            r = self.copy()
            r.__p = [c/r.__p[-1] for c in r.__p]
            return r
        else:
            return self.copy()

    def inverse(self, mod):
        """ <b>p.inverse(m)</b><br/>
            <i>Multiplicative inverse modulo polynomial m</i><br/>
            The polynomial mod must be relative prime to p!
        ----de----
            <b>p.inverse(m)</b><br/>
            <i>Multiplikatives Inverses modulo Polynom m</i><br/>
            Das Polynom mod muss teilerfremd zu p sein!
        """
        assert type(self.__p[0]) == Mod
        m = self.__p[0].m
        d, inv, y = Poly.gcdExt(self, mod)
        if not d.isOne():
            raise _MathGUIdeError("Polynomial is not invertible")
        # inv * self + y * mod == d == 1
        # ==> (inv * self) % mod == 1
        return inv % mod


#===============================================================
# Set
#---------------------------------------------------------------

class Set (list): # Polymorphic sets # Polymorphe Mengen
    """ Class <code>Set</code><br/>
        <i>Polymorphic sets</i><br/>
        <b>Example</b>: <code>A = Set([1,2,3])</code><br/>
        The class Set allows the following operations:<br/>
        <ul>
          <li>Intersection (operator &amp;)</li>
          <li>Union (operator +)</li>
          <li>Set difference (operator -)</li>
          <li>Symmetric difference (operator ^)</li>
          <li>Power set (method powerset)</li>
          <li>Cartesian product (Operator *)</li>
          <li>Element test (Operator <code>in</code>)</li>
          <li>Cardinality (Python function <code>len</code>)</li>
        </ul>
    ----de----
        Klasse <code>Set</code><br/>
        <i>Polymorphe Mengen</i><br/>
        <b>Beispiel</b>: <code>A = Set([1,2,3])</code><br/>
        Die Klasse Set erlaubt folgende Operationen:<br/>
        <ul>
          <li>Durchschnitt (Operator &amp;)</li>
          <li>Vereinigung (Operator +)</li>
          <li>Mengendifferenz (Operator -)</li>
          <li>Symmetrische Differenz (Operator ^)</li>
          <li>Potenzmenge (powerset)</li>
          <li>Kartesisches Produkt (Operator *)</li>
          <li>Elementtest (Operator <code>in</code>)</li>
          <li>Kardinalzahl (Python-Funktion <code>len</code>)</li>
        </ul>
    """
    _staticMethods = []

    def operators(self):
        """ <b>Set.operators()</b><br/>
            <i>For documentation only</i><br/>
            The following operators are defined in the class <b>Set</b>:<br/>
            <table>
             <tr><th>Op.</th><th>Function</th><th>Examples</th></tr>
             <tr><th>&amp;</th><td>Intersection</td><td><code>A &amp; B; A &amp;= B</code></td></tr>
             <tr><th>+</th><td>Union</td><td><code>A + B; A += B</code></td></tr>
             <tr><th>-</th><td>Set difference</td><td><code>A - B; A -= B</code></td></tr>
             <tr><th>^</th><td>Symmetric difference</td><td><code>A ^ B; A ^= B</code></td></tr>
             <tr><th>*</th><td>Cartesian product</td><td><code>A * B; A *= B</code></td></tr>
             <tr><th>in</th><td>is element of</td><td><code>x in A</code></td></tr>
            </table>
            <font color="#000080" size="-2">This method is for documentation only. It has no effect.</font>
        ----de----
            <b>Set.operators()</b><br/>
            <i>Nur zur Dokumentation</i><br/>
            Die Klasse <b>Set</b> erlaubt folgende Operatoren:<br/>
            <table>
             <tr><th>Op.</th><th>Funktion</th><th>Beispiele</th></tr>
             <tr><th>&amp;</th><td>Durchschnitt</td><td><code>A &amp; B; A &amp;= B</code></td></tr>
             <tr><th>+</th><td>Vereinigung</td><td><code>A + B; A += B</code></td></tr>
             <tr><th>-</th><td>Mengendifferenz</td><td><code>A - B; A -= B</code></td></tr>
             <tr><th>^</th><td>Symmetrische Differenz</td><td><code>A ^ B; A ^= B</code></td></tr>
             <tr><th>*</th><td>Kartesisches Produkt</td><td><code>A * B; A *= B</code></td></tr>
             <tr><th>in</th><td>ist Element von</td><td><code>x in A</code></td></tr>
            </table>
            <font color="#000080" size="-2">Diese Methode dient nur zur Dokumentation. Sie hat keine Wirkung.</font>
        """
        pass

    def __init__(self, el):
        """ Constructor
        ----de----
            Konstruktor
        """
        list.__init__(self, list(el))
        self.__normalize()

    def __normalize(self):
        self.sort()
        for i in range(len(self)-1, 0, -1):
            if self[i] == self[i-1]:
                del self[i]

    def __repr__(self):
        """ Object representation as string (in curly braces)
        """
        r = "{"
        for e in self:
            r += str(e) + ","
        return (r+"}").replace(",}","}")

    # In der Basisklasse list sind die "rich comparison"
    # Methoden definiert und haben Vorrang vor __cmp___
    # Deshalb werden sie hier wieder auf __cmp__ geleitet:
    def __lt__(A, B):
        return A.__cmp__(B) < 0

    def __le__(A, B):
        return A.__cmp__(B) <= 0

    def __eq__(A, B):
        return A.__cmp__(B) == 0

    def __ne__(A, B):
        return A.__cmp__(B) != 0

    def __gt__(A, B):
        return A.__cmp__(B) > 0

    def __ge__(A, B):
        return A.__cmp__(B) >= 0

    def __cmp__(A, B):
        """ For ordering sets
            Elements, wihch are sets themself come first and are sorted by length (cardinality)

            simultaneous overloading of  &lt; &lt;=, >, >=, == and !=

            If A is a Set and b not, A.__cmp__(b)  is called for comparing A with b (or vice versa)!

        """
        if not isinstance(B, Set):
            return -1           # erst Mengen, dann Nichtmengen
        elif len(A) != len(B):  # zwei Mengen: nach Größe ordnen
            return _cmp(len(A), len(B))
        elif len(A) == 0:
            return 0            # zwei leere Mengen
        else:                   # zwei gleich große nichtleere Mengen
            if A[0] != B[0]:
                return _cmp(A[0], B[0])
            else:
                return _cmp(Set(A[1:]), Set(B[1:]))

    def insert(self, e):
        """ <b>A.insert(e)</b><br/>
            <i>Inserts an additional element ti the set</i><br/>
            A.insert(e) is equivalent to A += Set(e)
        ----de----
            <b>A.insert(e)</b><br/>
            <i>Fügt ein zusätzliches Element in die Menge ein</i><br/>
            A.insert(e) ist äquivalent zu A += Set(e)
        """
        self.append(e)
        self.__normalize()

    def issubset(A, B):
        """ <b>A.issubset(B)</b><br/>
            <i>Is A subset of (oder equal to) B?</i><br/>
            Return value: bool
        ----de----
            <b>A.issubset(B)</b><br/>
            <i>Ist A Teilmenge (oder gleich) B?</i><br/>
            Rückgabewert: bool
        """
        for e in A:
            if e not in B:
                return False
        return True

    def issuperset(A, B):
        """ <b>A.issuperset(B)</b><br/>
            <i>Is A superset of (oder equal to) B?</i><br/>
            Return value: bool
        ----de----
            <b>A.issuperset(B)</b><br/>
            <i>Ist A Obermenge (oder gleich) B</i><br/>
            Rückgabewert: bool
        """
        return B.issubset(A)

    def __add__(A, B):
        """ Union (operators + and +=):
            a + b  --> __add__ (a, b), if a is of type Set
                       __radd__(b, a), otherwise
        """
        return Set(list(A) + list(B))

    def union(A,B):
        """ <b>A.union(B)</b><br/>
            <i>Union set</i><br/>
            Also operator notation: <b>A + B</b><br/>
            Instead of <code>A = A + B</code> you can write <code>A += B</code>.
        ----de----
            <b>A.union(B)</b><br/>
            <i>Vereinigungsmenge</i><br/>
            Auch Operatorschreibweise: <b>A + B</b><br/>
            Statt A = A + B kann auch A += B geschrieben werden.
        """
        return A + B

    def __sub__(A, B):
        """ Set difference (operators - and -=)
        """
        result = list(A)
        for e in B:
            if e in result:
                result.remove(e)
        return Set(result)

    def difference(A,B):
        """ <b>A.difference(B)</b><br/>
            <i>Set difference</i><br/>
            Also operator notation: <b>A - B</b><br/>
            Instead of <code>A = A - B</code> you can write <code>A -= B</code>.
        ----de----
            <b>A.difference(B)</b><br/>
            <i>Mengendifferenz</i><br/>
            Auch Operatorschreibweise: <b>A - B</b><br/>
            Statt A = A - B kann auch A -= B geschrieben werden.
        """
        return A - B

    def __and__(A, B):
        """ Intersection (operators &amp; and &amp;=)
        """
        result = list(A)
        for e in A:
            if e not in B:
                result.remove(e)
        return Set(result)

    def intersection(A,B):
        """ <b>A.intersection(B)</b><br/>
            <i>Intersection</i><br/>
            Also operator notation: <b>A - B</b><br/>
            Instead of <code>A = A &amp; B</code> you can write <code>A &amp;= B</code>.
        ----de----
            <b>A.intersection(B)</b><br/>
            <i>Schnittmenge</i><br/>
            Auch Operatorschreibweise: <b>A &amp; B</b><br/>
            Statt A = A &amp; B kann auch A &amp;= B geschrieben werden.
        """
        return A & B

    def __xor__(A, B):
        """ Symmetric difference (operators ^ and ^=)
        """
        return (A + B) - (A & B)

    def __pow__(A, B):
        """ Symmetric difference (operators ^ and ^=)
        """
        return (A + B) - (A & B)

    def symmDiff(A,B):
        """ <b>A.symmDiff(B)</b><br/>
            <i>Symmetric difference</i><br/>
            Also operator notation: <b>A ^ B</b><br/>
            Instead of <code>A = A ^ B</code> you can write <code>A ^= B</code>.
        ----de----
            <b>A.symmDiff(B)</b><br/>
            <i>Symmetrische Differenz</i><br/>
            Auch Operatorschreibweise: <b>A ^ B</b><br/>
            Statt A = A ^ B kann auch A ^= B geschrieben werden.
        """
        return A^B

    def __mul__(A, B):
        """ Cartesian product (operator *)
        """
        return Set([(x,y) for x in A for y in B])

    def product(A,B):
        """ <b>A.product(B)</b><br/>
            <i>Cartesian product</i><br/>
            Also operator notation: <b>A * B</b><br/>
            Instead of <code>A = A * B</code> you can write <code>A *= B</code>.
        ----de----
            <b>A.product(B)</b><br/>
            <i>Kartesisches Produkt</i><br/>
            Auch Operatorschreibweise: <b>A * B</b><br/>
            Statt A = A * B kann auch A *= B geschrieben werden.
        """
        return A*B

    def powerset(self):
        """ <b>A.powerset()</b><br/>
            <i>Powerset (set of alle subsets)</i><br/>
        ----de----
            <b>A.powerset()</b><br/>
            <i>Potenzmenge (Menge aller Teilmengen)</i><br/>
        """
        if len(self) == 0:
            return Set([Set([])])
        else:
            first = self[0]
            rest = self[1:]
            powRest = Set(list(rest)).powerset()
            pl = [Set([first])+x for x in powRest]
            return powRest + Set(pl)

    def card(self):
        """ <b>A.card()</b><br/>
            <i>Cardinality</i><br/>
            Same as len(A)
        ----de----
            <b>A.card()</b><br/>
            <i>Kardinalzahl</i><br/>
            Synonym zu len(A)
        """
        return len(self)

###"Repräsentative Mengen"
##
##               # 1 2 12
##A1 = Set(1,12) # + - +
##A2 = Set(2,12) # - + +
##
##                      # 1 2 3 12 13 23 123
##B1 = Set(1,12,13,123) # + - - +  +  -   +
##B2 = Set(2,12,23,123) # - + - +  -  +   +
##B3 = Set(3,13,23,123) # - - + -  +  +   +
##
##C1 = Set(1,12,13,14, 123, 124, 134, 1234)
##C2 = Set(2,12,23,24, 123, 124, 234, 1234)
##C3 = Set(3,13,23,34, 123, 134, 234, 1234)
##C4 = Set(4,14,24,34, 124, 134, 234, 1234)
##
##D1 = Set(1, 12,13,14,15, 123,124,125,134,135,145, 1234,1235,1245,1345, 12345)
##D2 = Set(2, 12,23,24,25, 123,124,125,234,235,245, 1234,1235,1245,2345, 12345)
##D3 = Set(3, 13,23,34,35, 123,134,135,234,234,345, 1234,1235,1345,2345, 12345)
##D4 = Set(4, 14,24,34,45, 124,134,145,234,245,345, 1234,1245,1345,2345, 12345)
##D5 = Set(5, 15,25,35,45, 125,135,145,235,245,345, 1235,1245,1345,2345, 12345)


#===============================================================
#  Propositional Logic
#---------------------------------------------------------------

class Logic: # propositional logic # Aussagenlogik
    """ Class <code>Logic</code>: Propositional Logic<br/>
        This is a static class (only static methods available)<br/>
        Formulae of propositional logic, passed to the static methods,
        must be composed as follows:<br/>
        <b>Variables</b> must be single letters.<br/>
        <b>Junctors</b>: "<code>not</code>", "<code>and</code>", "<code>or</code>", "<code>xor</code>", "<code>-></code>", "<code>&lt;-></code>".
    ----de----
        Klasse <code>Logic</code>: Aussagenlogik<br/>
        Statische Klasse: nur Klassenmethoden<br/>
        Aussagenlogische Formeln, die den Klassenmethoden übergeben werden,
        müssen wie folgt aufgebaut sein:<br/>
        <b>Variablen</b> müssen einzelne Buchstaben sein.<br/>
        <b>Junktoren</b>: "<code>not</code>", "<code>and</code>", "<code>or</code>", "<code>xor</code>", "<code>-></code>", "<code>&lt;-></code>".
    """
    _staticMethods = ("truthTables", "satisfies", "valid", "satisfiable", "implies", "equivalent")

    @staticmethod
    def truthTables(formula):
        """ <b>truthTables(formula)</b><br/>
            <i>All possible truth tables for variables in formula</i><br/>
            The result is returned as a generator object (for a single iteration).<br/>
            <b>Example</b><br/>
            <code>for t in Logic.truthTables("A and B"): print(t)</code><br/>
            returns:<br/>
            <code>{'A': False, 'B': False}<br/>
            {'A': False, 'B': True}<br/>
            {'A': True, 'B': False}<br/>
            {'A': True, 'B': True}</code>
        ----de----
            <b>truthTables(formula)</b><br/>
            <i>Alle möglichen Wahrheitswerttabellen für im Ausdruck
            <b>formula</b> vorkommenden Variablen</i><br/>
            Das Ergebnis wird als Generatorobjekt (zur einmaligen Iteration) ausgegeben.<br/>
            <b>Beispiel</b><br/>
            <code>for t in Logic.truthTables("A and B"): print(t)</code><br/>
            gibt folgendes aus:<br/>
            <code>{'A': False, 'B': False}<br/>
            {'A': False, 'B': True}<br/>
            {'A': True, 'B': False}<br/>
            {'A': True, 'B': True}</code>
        """
        table = {}
        for op in ("not", "and", "or", "xor"):
            formula = formula.replace(op, " ")
        for c in formula:
            if c.isalpha(): # Variable!
                table[c] = False
        n = len(table)
        vars = list(table.keys())
        while True:
            yield table
            k = n-1
            table[vars[k]] = not table[vars[k]]
            while not table[vars[k]]: # "Übertrag"
                k -= 1
                if k < 0:
                    return
                table[vars[k]] = not table[vars[k]]

    @staticmethod
    def satisfies(f, truthTable):
        """ <b>satisfies(f, truthTable)</b><br/>
            <i>Is the formula <b>f</b> true under the <b>truthTable</b>?</i><br/>
            <b>Example</b><br/>
            <code>Logic.satisfies("A or B", {"A":True, "B":False})</code>
            returns <code>True<code>.
        ----de----
            <b>satisfies(f, truthTable)</b><br/>
            <i>Wird die Formel <b>f</b> unter der Belegung <b>truthTable</b> wahr?</i><br/>
            <b>Beispiel</b><br/>
            <code>Logic.satisfies("A or B", {"A":True, "B":False})</code>
            gibt <code>True</code> zurück.
        """
        f = f.replace("not", "@0").replace("and", "@1").replace("or", "@2").replace("xor", "@3")
        for var in truthTable:
            f = f.replace(var, str(truthTable[var]))
        f = f.replace("@0", "not").replace("@1", "and").replace("@2", "or").replace("@3", "xor")
        return eval(f)

    @staticmethod
    def valid(f):
        """ <b>satisfies(f, truthTable)</b><br/>
            <i>Is the formula <b>f</b> valid</b>?</i><br/>
            <b>Example</b><br/>
            <code>Logic.valid("A or not A")</code>
            returns <code>True<code>.
        ----de----
            <b>satisfies(f, truthTable)</b><br/>
            <i>Ist die Formel <b>f</b> allgemeingültig?</i><br/>
            <b>Beispiel</b><br/>
            <code>Logic.valid("A or not A")</code>
            gibt <code>True</code> zurück.
        """
        f = f.replace("->", "<=").replace("<->", "==").replace("xor", "!=").replace("<<=", "==")
        for t in Logic.truthTables(f):
            if not Logic.satisfies(f, t):
                return False
        return True

    @staticmethod
    def satisfiable(f):
        """ <b>satisfiable(f)</b><br/>
            <i>Is the formula <b>f</b> satisfiable</b>?</i><br/>
            <b>Example</b><br/>
            <code>Logic.satisfiable("A and not A")</code>
            returns <code>False<code>.
        ----de----
            <b>satisfiable(f)</b><br/>
            <i>Ist die Formel <b>f</b> erfüllbar?</i><br/>
            <b>Beispiel</b><br/>
            <code>Logic.satisfiable("A and not A")</code>
            gibt <code>False</code> zurück.
        """
        return not Logic.valid("not({0})".format(f))

    @staticmethod
    def implies(f1, f2):
        """ <b>implies(f1, f2)</b><br/>
            <i>implies the formula <b>f1</b> the formula <b>f2</b>?</b>?</i><br/>
            <b>Example</b><br/>
            <code>Logic.implies("A and B", "A or B")</code>
            returns <code>True<code>.
        ----de----
            <b>implies(f1, f2)</b><br/>
            <i>Folgt aus der Formel <b>f1</b> die Formel <b>f2</b>?</i><br/>
            <b>Beispiel</b><br/>
            <code>Logic.implies("A and B", "A or B")</code>
            gibt <code>True</code> zurück.
        """
        return Logic.valid("not({0}) or ({1})".format(f1, f2))

    @staticmethod
    def equivalent(f1, f2):
        """ <b>equivalent(f1, f2)</b><br/>
            <i>Are the formulae <b>f1</b> and <b>f2</b> equivalent?</b>?</i><br/>
            <b>Example</b><br/>
            <code>Logic.equivalent("not (A and B)", "not A or not B"))</code>
            returns <code>True<code>.
        ----de----
            <b>equivalent(f1, f2)</b><br/>
            <i>Sind die Formeln <b>f1</b> und <b>f2</b> äquivalent?</i><br/>
            <b>Beispiel</b><br/>
            <code>Logic.equivalent("not (A and B)", "not A or not B"))</code>
            gibt <code>True</code> zurück.
        """
        return Logic.valid("({0}) == ({1})".format(f1, f2))


#===============================================================
#  Graphs and Trees
#---------------------------------------------------------------

class Edge: # Edges of a Graph # Kanten eines Graphen
    """ Class <code>Edge</code>: Edge of a <a href="Graph.html">Graph</a><br/>
    ----de----
        Klasse <code>Edge</code>: Kante eines Graphen (Klasse <a href="Graph.html">Graph</a>)<br/>
    """
    _staticMethods = []

    def __init__(self, source, target, weight=0, highlight=False):
        """ Constructor
        ----de----
            Konstruktor
        """
        self.source = source
        self.target = target
        self.weight = weight
        self.highlight = highlight

    def __repr__(self):
        """ Object representation as string
        """
        if self.weight == 0:
            return "{{{0},{1}}}".format(self.source, self.target)
        else:
            return "({{{0},{1}}}, w={2})".format(self.source, self.target, self.weight)

    @staticmethod
    def _fromShortRepr(s):
        p = s.split(",")
        n0, n1 = p[0:2]
        w = 0
        m = len(p) >= 4
        if len(p) >= 3:
            w = p[2]
        return Edge(int(n0), int(n1), int(w), m)

    def isIncidentNode(self, node):
        return node in (self.source, self.target)

    def adjacentNode(self, node):
        if node == self.source:
            return self.target
        elif node == self.target:
            return self.source
        else:
            return None

    def _shortRepr(self):
        """ short representation as string
        """
        if self.highlight:
            return "{0},{1},{2},1".format(self.source, self.target, self.weight)
        elif self.weight != 0:
            return "{0},{1},{2}".format(self.source, self.target, self.weight)
        else:
            return "{0},{1}".format(self.source, self.target)


class Graph: # Graphs # Graphen
    """ Class <code>Graph</code>: Graph (in the sense of graph theory)<br/>
        A Graph consists of nodes and edges.
        An edge (class <a href="Edge.html">Edge</a>) connects two nodes.<br/>
        Edges optionally may have a weight.<br/>
        <p class="remark">Although a pure Graph does not carry geometrical data, mathGUIde allows to attach
        positions to nodes corresponding to a Graph representation in the Graph window.</p>
    ----de----
        Klasse <code>Graph</code>: Graph (im Sinne der Graphentheorie)<br/>
        Ein Graph besteht aus Knoten und Kanten.
        Eine Kante (Klasse <a href="Edge.html">Edge</a>) verbindet zwei Knoten.<br/>
        Kanten können auch mit einem Gewicht versehen werden.<br/>
        <p class="remark">Obwohl zur Graphendefinition keine geometrischen Informationen gehören,
        kann mathGUIde die Positionen von im Graphenfenster angeordneten Knoten speichern.</p>
    """
    #............................................................
    _staticMethods = ("fromShortRepr", "example1", "complete", "completeBipartite", "petersen")

    @staticmethod
    def fromShortRepr(s):
        """ <b>Graph.fromShortRepr(s)</b><br/>
            <i>Makes the Graph from the short representation s</i>
        ----de----
            <b>Graph.fromShortRepr(s)</b><br/>
            <i>Erzeugt den Graphen aus der Kurzdarstellung s</i>
        """
        g = Graph()
        parts = s.split("|")
        if len(parts[0]) > 0:
            # nodes
            for n in parts[0].split(";"):
                p = n.split(",")
                g.nodes.append(int(p[0]))
                if len(p) >= 3:
                    g.locations[int(p[0])] = float(p[1]), float(p[2])
        if len(parts[1]) > 0:
            # edges
            for e in parts[1].split(";"):
                g.edges.append(Edge._fromShortRepr(e))
        return g

    @staticmethod
    def example1():
        """ <b>Graph.example1()</b>
        ----de----
            <b>Graph.example1()</b>
        """
        return Graph.fromShortRepr("1,80,20;2,240,20;3,0,100;4,160,100;5,320,100;6,0,180;7,160,180;8,320,180;9,80,260;10,240,260|1,2,36;1,4,15;2,4,69;2,5,1;4,5,53;3,4,92;1,3,99;3,6,37;6,7,21;7,8,88;4,6,27;4,7,32;5,7,10;5,8,66;6,9,3;7,9,18;7,10,49;8,10,43;9,10,74")

    @staticmethod
    def complete(n, r=130, margin=20):
        """ <b>Graph.complete(n, r=130, margin=20)</b><br/>
            <i>Complete Graph of order n, displayed on a circle with radius r.</i>
        ----de----
            <b>Graph.complete(n, r=130, margin=20)</b><br/>
            <i>Vollständiger Graph der Ordnung n angeordnet auf einem Kreis mit Radius r.</i>
        """
        c = margin + r # center x and y
        nodes = ";".join(["{0},{1},{2}".format(i+1, c+r*sin(i*2*pi/n), c-r*cos(i*2*pi/n)) for i in range(n)])
        edges = ";".join(["{0},{1},{2}".format(i+1,j+1,0) for i in range(n) for j in range(n) if i < j])
        return Graph.fromShortRepr("{0}|{1}".format(nodes, edges))

    @staticmethod
    def completeBipartite(n1, n2):
        """ <b>Graph.completeBipartite(n1,n2)</b><br/>
            <i>Complete bipartite Graph with partition sizes n1 and n2.</i>
        ----de----
            <b>Graph.completeBipartite(n1,n2)</b><br/>
            <i>Vollständiger bipartiter Graph mit den Partitionsgrößen n1 und n2.</i>
        """
        dx, dy, margin = 60, 160, 20
        margin1 = margin2 = margin
        if n1 < n2:
            margin1 += (n2-n1) * dx / 2
        elif n2 < n1:
            margin2 += (n1-n2) * dx / 2
        nodes1 = ";".join(["{0},{1},{2}".format(i+1, margin1+i*dx, margin) for i in range(n1)])
        nodes2 = ";".join(["{0},{1},{2}".format(n1+i+1, margin2+i*dx, margin+dy) for i in range(n2)])
        edges = ";".join(["{0},{1},{2}".format(i+1,n1+j+1,0) for i in range(n1) for j in range(n2)])
        return Graph.fromShortRepr("{0};{1}|{2}".format(nodes1, nodes2, edges))
        
    @staticmethod
    def petersen(r=150, margin=20):
        """ <b>Graph.petersen()</b><br/>
            <i>The Peterson Graph</i><br/>
            It serves as an example and counterexample for many problems in graph theory.
        ----de----
            <b>Graph.completeBipartite()</b><br/>
            <i>Der Peterson-Graph</i><br/>
            Er dient als Beispiel oder Gegenbeispiel für viele Probleme der Graphentheorie.
        """
        c = margin + r # center x and y
        r0 = r/2
        nodes1 = ";".join(["{0},{1},{2}".format(i+1, c+r*sin(i*2*pi/5),  c-r*cos(i*2*pi/5))  for i in range(5)])
        nodes2 = ";".join(["{0},{1},{2}".format(i+6, c+r0*sin(i*2*pi/5), c-r0*cos(i*2*pi/5)) for i in range(5)])
        edges = "1,2;1,5;1,6;2,3;2,7;3,4;3,8;4,5;4,9;5,10;6,8;6,9;7,9;7,10;8,10"
        return Graph.fromShortRepr("{0};{1}|{2}".format(nodes1, nodes2, edges))
    
    #............................................................

    def __init__(self):
        """ Constructor
        ----de----
            Konstruktor
        """
        self.nodes = []
        self.edges = []
        self.locations = {}

    def __repr__(self):
        """ Object representation as string
        """
        return "Graph(Nodes = {{{0}}}, Edges = {{{1}}})".format(
                join([str(n) for n in self.nodes], ","),
                join([str(e) for e in self.edges], ","))

    def shortRepr(self):
        """ <b>G.shortRepr()</b><br/>
            <i>short representation of Graph G as string</i><br/>
        ----de----
            <b>.shortRepr()</b><br/>
            <i>Kurzdarstellung des Graphen G (als String)</i>
        """
        # example: "1;2;3;4|1,2,36;1,3,99;1,4,15;2,4,69"
        return "{0}|{1}".format(join([self._shortNodeRepr(n) for n in self.nodes], ";"),
                                join([e._shortRepr() for e in self.edges], ";"))

    def _shortNodeRepr(self, node):
        """ short representation of a Node as string
        """
        if node in self.locations:
            return "{0},{1},{2}".format(node, self.locations[node][0], self.locations[node][1])
        else:
            return "{0}".format(node)

    def order(self):
        """ <b>G.order()</b><br/>
            <i>Order (number of vertices) of the Graph G.</i>
        ----de----
            <b>G.order()</b><br/>
            <i>Ordnung (Anzahl der Knoten) des Graphen G.</i>
        """
        return len(self.nodes)

    def complement(self):
        """ <b>G.complement()</b><br/>
            <i>Returns the graph with the same vertices as G, such that two vertices are adjacent if and only if they are not adjacent in G.</i>
        ----de----
            <b>G.complement()</b><br/>
            <i>Gibt den Graphen mit den gleichen Knoten wie G zurück, in dem sich zwischen zwei Knoten genau dann eine Kante befindet, wenn G dort keine Kante hat.</i>
        """
        g = Graph()
        g.nodes = list(self.nodes) # independent copy
        g.edges = []
        g.locations = dict(self.locations)
        for n1 in g.nodes:
            for n2 in g.nodes:
                if n1 < n2:
                    if not self.adjacent(n1, n2):
                        g.addEdges(Edge(n1, n2))
        return g

    def connectedNodes(self, node):
        """ <b>G.componentsCount()</b><br/>
            <i>Returns a list of all nodes of G which are connected with node (including node).</i>
        ----de----
            <b>G.componentsCount()</b><br/>
            <i>Gibt eine Liste aller mit node durch einen Kantenzug verbundenen Knoten in G zurück (einschließlich node).</i>
        """
        checked = 0
        nodes = [node]
        while checked < len(nodes):
            n = nodes[checked]
            for m in self.adjacentNodes(n):
                if (m not in nodes):
                    nodes.append(m)
            checked += 1
        nodes.sort()
        return nodes

    def componentsCount(self):
        """ <b>G.componentsCount()</b><br/>
            <i>The number connectivity components of the Graph G</i>
        ----de----
            <b>G.componentsCount()</b><br/>
            <i>Anzahl der Zusammenhangskomponenten des Graphen G</i>
        """
        count = 0
        uncheckedNodes = list(self.nodes) # independent copy
        while len(uncheckedNodes) > 0:
            count += 1
            n = uncheckedNodes[0]
            for m in self.connectedNodes(n):
                uncheckedNodes.remove(m)
        return count

    def isConnected(self):
        """ <b>G.isConnected()</b><br/>
            <i>Returns a truth value indicating wether the Graph G is connected</i>
        ----de----
            <b>G.isConnected()</b><br/>
            <i>Gibt als Wahrheitswert zurück, ob der Graph G zusammenhängend ist</i>
        """
        if len(self.nodes) == 0: # an empty graph is connected!
            return True
        else:
            return (self.componentsCount() == 1)

    def isTree(self):
        """ <b>G.isTree()</b><br/>
            <i>Returns a truth value indicating wether the Graph G is a tree</i>
        ----de----
            <b>G.isTree()</b><br/>
            <i>Gibt als Wahrheitswert zurück, ob der Graph G ein Baum ist</i>
        """
        return self.isConnected() and len(self.edges) == len(self.nodes) - 1

    def adjacentNodes(self, node):
        """ <b>G.adjacentNodes(node)</b><br/>
            <i>Returns a sorted list of all nodes of G adjacent to <b>node</b>.</i>
        ----de----
            <b>G.adjacentNodes(node)</b><br/>
            <i>Gibt eine sortierte Liste aller zu <b>node</b> adjazenten Knoten von G zurück.</i>
        """
        adjNodes = []
        for e in self.edges:
            if node == e.source:
                adjNodes.append(e.target)
            elif node == e.target:
                adjNodes.append(e.source)
        adjNodes.sort()
        return adjNodes

    def adjacent(self, node1, node2):
        """ <b>G.adjacentNodes(node1, node2)</b><br/>
            <i>Determines wether node1 and node2 are adjacent (bool return value) in G.</i>
        ----de----
            <b>G.adjacentNodes(node1, node2)</b><br/>
            <i>Gibt an, ob node1 und node2 in G adjazent sind (bool).</i>
        """
        return node2 in self.adjacentNodes(node1)

    def incidentEdges(self, node):
        incEdges = []
        for e in self.edges:
            if e.isIncidentNode(node):
                incEdges.append(e)
        return incEdges

    def degree(self, node):
        """ <b>G.degree(node)</b><br/>
            <i>Returns the degree of node (count of adjacent nodes)</i>
        ----de----
            <b>G.degree(node)</b><br/>
            <i>Gibt den Grad des Knotens node (Anzahl der adjazenten Knoten) in G zurück</i>
        """
        return len(self.incidentEdges(node))

    def edge(self, node1, node2):
        """ <b>G.edge(node1, node2)</b><br/>
            <i>The connecting edge between the two nodes</i><br/>
            If node1 and node2 are adjacent, the connecting edge is returned, otherwise <code>None</code>.
        ----de----
            <b>G.edge(node1, node2)</b><br/>
            <i>Die Verbindungskante zwischen den zwei Knoten</i><br/>
            Wenn node1 und node2 adjazent sind, wird die the Verbindungskante zurückgegeben, sonst <code>None</code>.
        """
        for e in self.edges:
            if (e.source == node1 and e.target == node2) or (e.source == node2 and e.target == node1):
                return e
        return None

    def display(self, n=0):
        """ <b>G.display(n=0)</b><br/>
            <i>Displays the Graph G in the n-th Graph tab.</i><br/>
            if n == 0, the Graph is displayed in the active Graph window.
        ----de----
            <b>G.display(n=0)</b><br/>
            <i>Zeigt den Graphen G im n-ten Graphen-Fenster an.</i>
            Ist n = 0, wird der Graph im aktiven Graphen-Fenster angezeigt.
        """
        return "@G{0}:{1}".format(n, self.shortRepr())

    def addNodes(self, *nodes):
        self.nodes += nodes

    def addEdges(self, *edges):
        for e in edges:
            if e.source not in self.nodes:
                self.nodes.append(e.source)
            if e.target not in self.nodes:
                self.nodes.append(e.target)
        self.edges += edges

    def copy(self):
        """ <b>G.copy()</b><br/>
            <i>independent copy ("clone") of G</i>
        ----de----
            <b>G.copy()</b><br/>
            <i>unabhängige Kopie von G</i>
        """
        return Graph.fromShortRepr(self.shortRepr())

    def reduceToHighlightedEdges(self):
        """ <b>G.reduceToHighlightedEdges()</b><br/>
            <i>Removes all not highlighted edges from the Graph G.</i><br/>
            The highlighting of the remaining edges is turned off.
        ----de----
            <b>G.reduceToHighlightedEdges()</b><br/>
            <i>Entfernt alle nicht hervorgehobenen Kanten aus dem Graphen G.</i><br/>
            Die Hervorhebung der verbleibenden Kanten wird zurückgenommen.
        """
        highlightedEdges = [e for e in self.edges if e.highlight]
        otherEdges = [e for e in self.edges if not e.highlight]
        for e in highlightedEdges:
            e.highlight = False
        for e in otherEdges:
            self.edges.remove(e)

    def cancelHighlight(self):
        """ <b>G.cancelHighlight()</b><br/>
            <i>Resets all highlighted edges to normal state.</i>
        ----de----
            <b>G.cancelHighlight()</b><br/>
            <i>Setzt alle hervorgehobenen Kanten in den Normalzustand zurück.</i>
        """
        for e in self.edges:
            e.highlight = False

    def highlightShortestPath(self, startNode, goalNode):
        """ <b>G.highlightShortestPath(startNode, goalNode)</b><br/>
            <i>Highlightes a shortest path (if existent) from startNode to goalNode using the Dijkstra algorithm.</i><br/>
            Returns the length of the path (or inf if there is no connection).<br/>
            Use the method <a href="Graph.display"><code>display</code></a> in order
            to make the highlighting visible.
        ----de----
            <b>G.highlightShortestPath(startNode, goalNode)</b><br/>
            <i>Hebt einen kürzesten Pfad von startNode nach goalNode hervor (mit dem Dijkstra-Algorithmus).</i><br/>
            Die Länge des Pfades wird zurückgegeben (oder inf, falls keine Verbindung vorhanden).<br/>
            Mit der Methode <a href="Graph.display"><code>display</code></a> können Sie
            die Hervorhebung sichtbar machen.
        """
        self.cancelHighlight()
        predecessorNode = {}
        unfinishedNodes = list(self.nodes) # independent copy

        upperBound = {} # dictionary with upper bounds for path length from startNode to node
        infinite = float("inf")
        for n in self.nodes:
            upperBound[n] = infinite
        upperBound[startNode] = 0

        activeNode = startNode
        while True:
            unfinishedNodes.remove(activeNode)
            fail = True
            for e in self.incidentEdges(activeNode):
                n = e.adjacentNode(activeNode)
                if n in unfinishedNodes:
                    fail = False
                    # add 1 to weights in order to count path lengths
                    # for unweighted Graphs (default weight = 0):
                    dist = upperBound[activeNode] + e.weight + 1
                    if dist < upperBound[n]: # a shorter path to n discovered
                        upperBound[n] = dist
                        predecessorNode[n] = activeNode
            if fail:
                break
            best, newNode = min([(upperBound[n], n) for n in unfinishedNodes])
            activeNode = newNode

        if goalNode in unfinishedNodes:
            return infinite
        # go back from goalNode to startNode:
        n = goalNode
        length = 0
        while n != startNode:
            p = predecessorNode[n]
            e = self.edge(p, n)
            e.highlight = True
            length += e.weight
            n = p
        return length

    def highlightMinimalSpanTree(g, trace=False):
        """ <b>G.minimalSpanTree(trace=False)</b><br/>
            <i>Highlightes a minimal spanning tree using the Kruskal algorithm.</i><br/>
            If <code>trace</code> is <code>True</code>, all steps are printed.
        ----de----
            <b>G.minimalSpanTree(trace=False)</b><br/>
            <i>Hebt einen minimalen aufspannenden Baum hervor (mit dem Kruskal-Algorithmus).</i><br/>
            Wenn <code>trace</code> den Wert <code>True</code> hat, werden alle Schritte angezeigt.
        """

        nodeColor = {}
        marked = {}
        for n in g.nodes:
            nodeColor[n] = 0
        def weight(edge):
            return edge.weight
        g.edges.sort(key=weight)

        curColor = 1
        spanTreeSize = 0

        for edge in g.edges:
            marked[edge] = False

        for edge in g.edges:
            choose = True
            n0, n1 = edge.source, edge.target
            if nodeColor[n0] == nodeColor[n1]:
            # both nodes have tha same color
                if nodeColor[n0] == 0:
                    nodeColor[n0] = curColor   # both nodes uncolored (no edge):
                    nodeColor[n1] = curColor   # use current color,
                    curColor += 1              # new acurrent color
                else:
                    choose = False  # coonecting would result in a circle ==> skip edge
            else:
            # both nodes have different colors
                if min(nodeColor[n0], nodeColor[n1]) == 0:
                # one node uncolored: use color of other node
                    if nodeColor[n0] == 0:
                        nodeColor[n0] = nodeColor[n1]
                    else:
                        nodeColor[n1] = nodeColor[n0]
                else:
                # both nodes colored: join two components (change all nodes with lower color to higher color)
                    lowColor  = min(nodeColor[n0], nodeColor[n1])
                    highColor = max(nodeColor[n0], nodeColor[n1])
                    for n in nodeColor:
                        if nodeColor[n] == lowColor:
                            nodeColor[n] = highColor
            marked[edge] = choose
            if trace:
                if choose:
                    spanTreeSize += 1
                    print("{0} Choose".format(spanTreeSize), end=" ")
                else:
                    print("Skip", end=" ")
                print("Edge {0}-{1} (Weight {2})".format(edge.source, edge.target, edge.weight))
            if spanTreeSize == len(nodeColor) - 1:
                break
        for edge in marked:
            edge.highlight = marked[edge]

    def isomorphic(G,H):
        """ <b>G.isomorphic(H)</b><br/>
            <i>Determines wether G and H are isomorphic Graphs.</i><br/>
        ----de----
            <b>G.isomorphic(H)</b><br/>
            <i>Entscheidet, ob G und H isomorphe Graphen sind.</i><br/>
        """
        # Very inefficient implementation!
        # As a first optimization we could permute only nodes of equal degree
        g, h = G.nodes, H.nodes
        n = len(g)
        if n != len(h):
            return False
        if n > 8:
            raise _MathGUIdeError("max. 8 nodes!")
        for p in permutations(g):
            match = True
            for i in fromTo(0, n-2):
                for k in fromTo(i+1, n-1):
                    if G.adjacent(p[i], p[k]) != H.adjacent(h[i], h[k]):
                        match = False
                        break
                if not match:
                    break
            if match:
                return True
        return False

#===============================================================


def printTable(height, width, cellFunction, **options):
    """ <b>printTable(height, width, cellFunction, rowHeadFn=None, colHeadFn=None, corner="")</b><br/>
        <i>displays a table of <var>height</var> rows and <var>width</var> columns, whose cells are calculated by <var>cellFunction</var>.</i><br/>
        <var>cellFunction</var> must be a function with two parameters (row, column).<br/>
        The following optional named parameters are supported:<br/>
        <var>rowHeadFn</var>: if given, row header cells are added and calculated using this function<br/>
        <var>colHeadFn</var>: if given, column header cells are added and calculated using this function<br/>
        <var>headAlign</var>: "left" (default), "right" or "center"<br/>
        <var>cellAlign</var>: "left", "right" (default) or "center"<br/>
        <var>title</var>: if given, a title line is added, using this text<br/>
        <var>corner</var>: if given (and rowHeadFn and colHeadFn are given), this text is displayed in the top left table corner.<br/>
        This function is used by the menu command  Insert -- Table.<br/>
        <b>Example</b>:<pre>
 printTable(7, 7, lambda i,k: i*k % 7,
    rowHeadFn = lambda i: i,
    colHeadFn = lambda k: k,
    title = "Multiplication mod 7",
    corner = "*")</pre>
        prints a multiplication table modulo 7.
    ----de----
        <b>printTable(height, width, cellFunction, rowHeadFn=None, colHeadFn=None, corner="")</b><br/>
        <i>gibt eine Tabelle mit <var>height</var> Zeilen und <var>width</var> Spalten aus, deren Zellen mit <var>cellFunction</var> berechnet werden.</i><br/>
        <var>cellFunction</var> muss eine Funktion mit zwei Parametern (Zeile, Spalte) sein.<br/>
        Folgende optionalen benannten Parameter werden unterstützt:<br/>
        <var>rowHeadFn</var>: Falls vorhanden, werden Zeilenköpfe angezeigt und mit dieser Funktion berechnet.<br/>
        <var>colHeadFn</var>: Falls vorhanden, werden Spaltenköpfe angezeigt und mit dieser Funktion berechnet.<br/>
        <var>headAlign</var>: "left" (Standard), "right" oder "center"<br/>
        <var>cellAlign</var>: "left", "right" (Standard) oder "center"<br/>
        <var>title</var>: Falls vorhanden, wird eine Titelzeile mit diesem Text über der Tabelle angezeigt<br/>
        <var>corner</var>: Falls vorhanden (und auch rowHeadFn und colHeadFn), wird dieser Text
        in der oberen linken Tabellenecke angezeigt.<br/>
        Diese Funktion wird vom Menübefehl Einfügen -- Tabelle verwendet.<br/>
        <b>Beispiel</b>:<pre>
 printTable(7, 7, lambda i,k: i*k % 7,
    rowHeadFn = lambda i: i,
    colHeadFn = lambda k: k,
    title = "Multiplikation mod 7",
    corner = "*")</pre>
        Gibt eine Multiplikationstabelle modulo 7 aus.
    """
    rowHeadFn = options.get("rowHeadFn", None)
    colHeadFn = options.get("colHeadFn", None)
    headAlign = options.get("headAlign", "left")
    cellAlign = options.get("cellAlign", "right")
    title = options.get("title", "")
    titleAlign = options.get("titleAlign", "center")
    corner = options.get("corner", "")

    tdHead  = '<td bgcolor="#C6D9F1" align="{0}">'.format(headAlign)
    tdValue = '<td bgcolor="#F2F7FC" align="{0}">'.format(cellAlign)
    print('{html}<table cellspacing="6">') # don't convert to HTML again !
    if len(title) > 0:
        colspan = width
        if rowHeadFn:
            colspan = width + 1
        print('<tr><td colspan="{0}" align="{1}" bgcolor="#376092"><font color="white">{2}</font></td></tr>'.format(colspan, titleAlign, title))
    if colHeadFn:
        print("<tr>", end="") # heading
        if rowHeadFn:
            print("<td>{0}</td>".format(corner), end="") # top left corner
        for k in range(width):
            print(tdHead, colHeadFn(k), "</td>", end="")
        print("</tr>")
    for i in range(height):
        print("<tr>", end="") # heading
        if rowHeadFn:
            print(tdHead, rowHeadFn(i), "</td>", end="")
        for k in range(width):
            print(tdValue, cellFunction(i,k), "</td>", end="")
        print("</tr>")
    print("</table>{/html}")


def printValueTable(values, variable="n", start=0, end=10, increment=1, vertical=True):
    """ <b>printValueTable(values, variable="n", start=0, end=10, increment=1, vertical=True)</b><br/>
        <i>Displays a value table.</i><br/>
        This function is used by the menu command  Insert -- Value Table.<br/>
        <b>Example</b>:<pre>
printValueTable([(lambda n:n^2, "n²"), (lambda n:2^n, "2<sup>n</sup>")],
                 "n", 0, 10, True)</pre>
        displays the values of n² and 2<sup>n</sup> for n from 0 to 10.
    ----de----
        <b>printValueTable(values, variable="n", start=0, end=10, increment=1, vertical=True)</b><br/>
        <i>Gibt eine Wertetafel aus.</i><br/>
        Diese Funktion wird vom Menübefehl Einfügen -- Wertetafel verwendet.<br/>
        <b>Beispiel</b>:<pre>
printValueTable([(lambda n:n^2, "n²"), (lambda n:2^n, "2<sup>n</sup>")],
                 "n", 0, 10, True)</pre>
        Listet die Werte n² und 2<sup>n</sup> für n von 0 bis 10 auf.
    """

    titles = [variable] + [v[1] for v in values]
    functions = [lambda n:n] + [v[0] for v in values]

    tdHead  = '<td bgcolor="#C6D9F1" align="left">'
    tdValue = '<td bgcolor="#F2F7FC" align="right">'

    def args():
        n = start
        if increment > 0:
            while n <= end:
                yield n
                n += increment
        else:
            while n >= end:
                yield n
                n += increment

    print('{html}<table cellspacing="6">') # don't convert to HTML again !
    if vertical:
        print("<tr>", end="") # heading
        for t in titles:
            print(tdHead, t, "</td>", end="")
        print("</tr>")

        for n in args():
            print("<tr>", end="") # heading
            for f in functions:
                print(tdValue, f(n), "</td>", end="")
            print("</tr>")
    else:
        for i in range(len(functions)):
            print("<tr>", end="")
            print(tdHead, titles[i], "</td>", end="")
            for n in args():
                print(tdValue, functions[i](n), "</td>", end="")
    print("</table>{/html}")

def plot(terms, x0, x1, **options):
    termList = [terms]
    if type(terms) in (tuple, list):
        termList = terms
    nSteps = options.get("nSteps", 0)
    dots = options.get("dots", False)
    if nSteps == 0 and isInteger(x0) and isInteger(x1) and x0 < x1:
        bars = options.get("bars", True)
        lines = options.get("lines", False)
        nSteps = x1 - x0
    else:
        nSteps = options.get("nSteps", 200)
        bars = options.get("bars", False)
        lines = options.get("lines", True)
    optStr = ""
    if dots:  optStr += "d"
    if bars:  optStr += "b"
    if not lines: optStr += "n"
    if nSteps == x1 - x0:
        d = 1
    else:
        d = (x1-x0)/float(nSteps)
    val = []
    for i in range(nSteps+1):
        x = x0 + i*d
        try:
            s = str(float(eval(termList[0])))
        except:
            s = "0"
        for j in range(1, len(termList)):
            try:
                s += "&" + str(float(eval(termList[j])))
            except:
                s += "&0"
        val.append(s)
    mark = ["{0},{1}".format(p[0], p[1]) for p in options.get("mark", [])]
    return "@P{0}|{1}|{2}|{3}|{4}".format(optStr, str(x0), str(x1), ";".join(val), ";".join(mark))

def plotFunctions(functions, x0, x1, **options):
    if type(functions) not in (tuple, list):
        functions = [functions]
    nSteps = options.get("nSteps", 0)
    dots = options.get("dots", False)
    if nSteps == 0 and isInteger(x0) and isInteger(x1) and x0 < x1:
        bars = options.get("bars", True)
        lines = options.get("lines", False)
        nSteps = x1 - x0
    else:
        nSteps = options.get("nSteps", 200)
        bars = options.get("bars", False)
        lines = options.get("lines", True)
    optStr = ""
    if dots:  optStr += "d"
    if bars:  optStr += "b"
    if not lines: optStr += "n"
    if nSteps == x1 - x0:
        d = 1
    else:
        d = (x1-x0)/float(nSteps)
    val = []
    for i in range(nSteps+1):
        x = x0 + i*d
        try:
            s = str(float(functions[0](x)))
        except:
            s = "0"
        for j in range(1, len(functions)):
            try:
                s += "&" + str(float(functions[j](x)))
            except:
                s += "&0"
        val.append(s)
    mark = ["{0},{1}".format(p[0], p[1]) for p in options.get("mark", [])]
    return "@P{0}|{1}|{2}|{3}|{4}".format(optStr, str(x0), str(x1), ";".join(val), ";".join(mark))


#===============================================================
#  Introspection
#  (internal functions for the mathGUIde GUI)
#---------------------------------------------------------------

def _x_getFunc(glob, loc):
    all = [x.strip() for x in glob.keys() if not "_" in x]
    classes = [x+"C" for x in all if str(glob[x]).startswith("<class") ]
    functions = [x+"F" for x in all if str(glob[x]).startswith("<function") ]
    all = classes + functions
    all.sort(key=str.lower)
    return "L" + ",".join(all)

def _x_getMeth(c, glob, loc):
    s_static = [x.strip() for x in dir(eval(c, glob, loc)) if x[0] == "_"]
    all = [x.strip() for x in dir(eval(c, glob, loc)) if x[0] != "_"]
    typeDescr = str(type(eval(c, glob, loc)))
    elements = []
    if typeDescr == "<class 'type'>":
        # class: display static methods
        try:
            elements = [m+"F" for m in eval(c+"._staticMethods")]
        except:
            pass
    elif typeDescr[:6] == "<class":
        # object: display members
        for x in all:
            typeStr = str(type(eval(c+"."+x, glob, loc)))
            if "function" not in typeStr:
                if "'method'" in typeStr:
                    elements.append(x+"F")
                else:
                    elements.append(x+"E")
    else:
        # global functions
        elements = [x+"()" for x in all if "function" in str(type(eval(c+"."+x, glob, loc)))]
    elements.sort(key=str.lower)
    return "L" + ",".join(elements)

def _x_getFnPar(fn, glob, loc):
    fn = eval(fn, glob, loc)
    descr = str(fn)
    s = inspect.getdoc(fn)
    return "S" + s

def _x_getMethPar(obj, fn, glob, loc):
    fn = eval(obj + "." + fn, glob, loc)
    s = inspect.getdoc(fn)
    return "S" + s

def _x_parse(s):
    op = s[-1]
    s = s[:-1]
    if op == "(":                               # \w --> word characters [a-zA-Z0-9_]
        m = re.search(r"\w+\s*\.\s*\w+\s*$", s) # \s --> whitespace
        if m:
            return [x.strip() for x in s[m.start():].split(".")] + ["("]
    m = re.search(r"\w+\s*$", s)
    if m:
        return ["", s[m.start():].strip(), op]

def _contextInfo(s, glob, loc):
    """ Erwartet den Inhalt der Cursorzeile links vom Cursor.
        Rückgabewert:
          - erstes Zeichen "S": danach Info-String für Statuszeile
          - erstes Zeichen "L": danach durch Komma getrennte Methoden
    """
    if s == "@@@": # globale Klassen und Funktionen:
        try:
            r = _x_getFunc(glob, loc)
        except:
            r = "***"
        return r

    p = _x_parse(s)
    try:
        if p[-1] == ".":
            r = _x_getMeth(p[-2], glob, loc)
        elif p[-1] == "(":
            if p[-3] != "":
                r = _x_getMethPar(p[-3], p[-2], glob, loc)
            else:
                r = _x_getFnPar(p[-2], glob, loc)
    except:
        r = "***"
    return r

def _x_classMethInfo(s, glob, loc):
    c, m, p = _x_parse(s)
    assert p == "("
    if c == "": # global
        return "global.{0}".format(m)
    else:
        t = str(type(eval(c, glob, loc)))
        if not t == "<class 'type'>":
            c = t[18:-2] # e.g. "<class 'mathguide.Matrix'>" --> "Matrix"
        return "{0}.{1}".format(c, m)

#end