INITIAL COMMIT

LICENSE

@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2017 Louis Abraham, Yassir Akram
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

Makefile

@@ -0,0 +1,26 @@
+define exportdoc
+import sys
+from importlib import import_module
+import pydoc
+pydoc.isdata = lambda _: False
+class MarkdownDoc(pydoc._PlainTextDoc):
+    def getdocloc(self, _):
+        return None
+    def docmodule(self, m):
+        m.__name__ += '\n\n'
+        return '\n'.join(super().docmodule(m).split('\n')[:-4])
+
+renderer = MarkdownDoc()
+for m in sys.argv[1:]:
+	print(renderer.docmodule(import_module(m)),
+				file=open(m + '.txt', 'w'))
+endef
+export exportdoc
+
+doc:
+	@-mkdir doc
+	@path=$$(pwd); \
+	cd doc; \
+	PYTHONPATH=$$path:$$PYTHONPATH python3 -c "$$exportdoc" algnuth algnuth.polynom algnuth.quadratic algnuth.jacobi algnuth.ideals
+
+.PHONY: doc

README.md

@@ -0,0 +1,99 @@
+Algebraic Number Theory package
+===============================
+
+**Louis Abraham** and **Yassir Akram**
+
+Features
+--------
+
+### Jacobi symbol
+
+    >>> from algnuth.jacobi import jacobi
+    >>> jacobi(3763, 20353)
+    -1
+
+### Solovay--Strassen primality test
+
+    >>> from algnuth.jacobi import solovay_strassen
+    >>> p = 12779877140635552275193974526927174906313992988726945426212616053383820179306398832891367199026816638983953765799977121840616466620283861630627224899026453
+    >>> q = 12779877140635552275193974526927174906313992988726945426212616053383820179306398832891367199026816638983953765799977121840616466620283861630627224899027521
+    >>> n = p * q
+    >>> solovay_strassen(p)
+    True
+    >>> solovay_strassen(q)
+    True
+    >>> solovay_strassen(n)
+    False
+
+### Quadratic forms
+
+    >>> from algnuth.quadratic import *
+    >>> display_classes(-44)
+    x^2 + 11y^2
+    2x^2 + 2xy + 6y^2
+    3x^2 - 2xy + 4y^2
+    3x^2 + 2xy + 4y^2
+    >>> display_primitive_forms(-44)
+    x^2 + 11y^2
+    3x^2 - 2xy + 4y^2
+    3x^2 + 2xy + 4y^2
+    >>> display_ambiguous_classes(-44)
+    x^2 + 11y^2
+    2x^2 + 2xy + 6y^2
+    >>> display(*reduced(18, -10, 2))
+    2x^2 + 2xy + 6y^2
+
+### Real polynoms
+
+    >>> from algnuth.polynom import Polynomial
+    >>> P = Polynomial([0] * 10 + [-1, 0, 1])
+    >>> print(P)
+    X^12-X^10
+    >>> P(2)
+    3072
+    >>> P.disc
+    0
+    >>> P.sturm() # Number of distinct real roots
+    3
+    >>> P.r1 # Number of real roots with multiplicity
+    12
+
+### Modular arithmetic
+
+    >>> P = Polynomial([1, 2, 3])
+    >>> Pmodp = P.reduceP(41)
+    >>> print(P ** 3)
+    27⋅X^6+54⋅X^5+63⋅X^4+44⋅X^3+21⋅X^2+6⋅X+1
+
+### Polynomial division
+
+    >>> A = Polynomial([1, 2, 3, 4]).reduceP(7)
+    >>> B = Polynomial([0, 1, 2]).reduceP(7)
+    >>> print(A)
+    4⋅X^3+3⋅X^2+2⋅X+1
+    >>> print(B)
+    2⋅X^2+X
+    >>> print(A % B)
+    5⋅X+1
+    >>> print(A // B)
+    2⋅X+4
+    >>> print((A // B) * B + A % B)
+    4⋅X^3+3⋅X^2+2⋅X+1
+
+### Berlekamp's algorithm
+
+    >>> P = Polynomial([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])
+    >>> Pmodp = P.reduceP(41)
+    >>> print(Polynomial.ppfactors(Pmodp.factor()))
+    12⋅(X+31)⋅X⋅(X^2+40⋅X+24)⋅(X^2+36⋅X+13)⋅(X^6+34⋅X^5+26⋅X^4+13⋅X^3+25⋅X^2+26⋅X+35)
+
+### Unique Factorization of Ideals
+
+    >>> from algnuth.ideals import factorIdeals
+    >>> factorIdeals(Polynomial([4, 0, 0, 1]))
+    X^3+4 mod 2 = X^3
+    (2) = (2, α)^3
+    X^3+4 mod 3 = (X+1)^3
+    (3) = (3, α+1)^3
+    X^3+4 mod 5 = (X+4)⋅(X^2+X+1)
+    (5) = (5, α+4)⋅(5, α^2+α+1)

algnuth/ideals.py

@@ -0,0 +1,49 @@
+"""
+This module provides functions to manipulate
+extension fields given their minimal polynomial.
+"""
+
+from math import pi, factorial
+
+from .polynom import Polynomial
+from .jacobi import sieve
+
+
+def minkowski_bound(P):
+    """
+    Any ideal of the ring of integers
+    of the algebraic number field whose
+    minimal polynomial is P contains
+    an integer N such that
+    1 ≤ N ≤ minkowski_bound(P)
+    """
+    return (4 / pi) ** P.r2 * factorial(P.deg) / P.deg ** P.deg * abs(P.disc) ** .5
+
+
+def idealsContaining(P, p):
+    """
+    Ideals of the extension field of minimal
+    polynomial P containing the prime p
+    """
+    Pmodp = P.reduceP(p)
+    c, Ds = Pmodp.factor()
+    print('%s mod %s = %s' % (P, p, Polynomial.ppfactors((c, Ds))))
+    print("(%s) = " % p +
+          '⋅'.join(("(%s, %s)" % (p, D)
+                    + (v > 1) * ("^%s" % v))
+                   if sum(Ds.values()) > 1 else "(%s)" % p
+                   for D, v in Ds.items()).replace("X", "α"))
+
+
+def factorIdeals(P):
+    """
+    Finds the ideals of the ring of integers
+    of the algebraic number field whose
+    minimal polynomial is P
+    """
+    b = int(minkowski_bound(P))
+    if b == 1:
+        print('Principal!')
+    for p in sieve(b + 1):
+        idealsContaining(P, p)
+    print()

algnuth/jacobi.py

@@ -0,0 +1,90 @@
+"""
+Jacobi symbol, Solovay–Strassen primality test and sieve of Eratosthenes
+"""
+
+from random import randrange
+from math import gcd
+
+
+def expsign(sign, exp):
+    """
+    optimization of sign ** exp
+    """
+    if sign == 1:
+        return 1
+    assert sign == -1
+    return -1 if exp % 2 else 1
+
+
+def jacobi(m, n):
+    """
+    Jacobi's symbol
+    the rule for (-1/n) is not used
+    """
+    assert n % 2
+    if m == 2:
+        if n % 8 in [1, 7]:
+            return 1
+        return -1
+    m %= n
+    q = 0
+    while m & 1 == 0:
+        m >>= 1
+        q += 1
+    if m == 1:
+        return expsign(jacobi(2, n), q)
+    return (expsign(jacobi(2, n), q)
+            * (-1 if (n % 4 == 3) and (m % 4 == 3) else 1)
+            * jacobi(n, m))
+
+
+def solovay_strassen(n, prec=50):
+    """
+    Solovay–Strassen primality test
+    with error probability less than 2^-prec
+    """
+    if n == 1:
+        return False
+    if n % 2 == 0:
+        return n == 2
+    e = (n - 1) // 2
+    for _ in range(prec):
+        x = randrange(1, n)
+        if gcd(x, n) != 1 or pow(x, e, n) != (jacobi(x, n) % n):
+            return False
+    return True
+
+
+def sieve(limit=10**5):
+    """
+    Sieve of Eratosthenes
+    """
+    l = [True] * (limit + 1)
+    l[0] = l[1] = 0
+    next = 2
+    while next < limit:
+        for k in range(2, limit // next + 1):
+            l[k * next] = False
+        next += 1
+        while next < limit and not l[next]:
+            next += 1
+    return [i for i in range(limit) if l[i]]
+
+
+def test_solovay_strassen(limit=10**5):
+    """
+    Runs in ~20s with limit = 10^5
+    """
+    primes = set(sieve(limit))
+    for i in range(limit):
+        assert (i in primes) == solovay_strassen(i)
+
+
+if __name__ == '__main__':
+    test_solovay_strassen(10**3)
+    p = 12779877140635552275193974526927174906313992988726945426212616053383820179306398832891367199026816638983953765799977121840616466620283861630627224899026453
+    q = 12779877140635552275193974526927174906313992988726945426212616053383820179306398832891367199026816638983953765799977121840616466620283861630627224899027521
+    n = p * q
+    assert solovay_strassen(p)
+    assert solovay_strassen(q)
+    assert not solovay_strassen(n)