Merge pull request #1 from skirpichev/continued_fractions-dammina

Continued fractions for dammina

doc/src/modules/ntheory.rst

@@ -106,3 +106,6 @@ Ntheory Functions Reference
 .. autofunction:: legendre_symbol
 
 .. autofunction:: jacobi_symbol
+
+.. automodule:: sympy.ntheory.continued_fraction
+    :members:

sympy/ntheory/__init__.py

@@ -13,3 +13,5 @@
     primitive_root, nthroot_mod, is_nthpow_residue, sqrt_mod_iter
 from .multinomial import binomial_coefficients, binomial_coefficients_list, \
     multinomial_coefficients
+from .continued_fraction import continued_fraction_periodic, \
+    continued_fraction_iterator

sympy/ntheory/continued_fraction.py

@@ -1,192 +1,127 @@
-from sympy.core.compatibility import as_int
-from sympy.core.numbers import Rational
-from sympy.core.numbers import Integer
+from sympy.core.numbers import Integer, Rational
 
-from fractions import gcd
 
+def continued_fraction_periodic(p, q, d=0):
+    r"""
+    Compute the continued fraction expansion of a rational or a
+    quadratic irrational number, i.e. `\frac{p + \sqrt{d}}{q}`, where
+    `p`, `q` and `d \ge 0` are integers.
 
-def continued_fraction(num, den, delta):
-    """
-    Return continued fraction expansion of surd.
-
-    Continued fraction expansion of a rational number is an expression
-    obtained through an iterative process of representing a number as
-    the sum of its integer part and the reciprocal of another number,
-    then writing this other number as the sum of its integer part and
-    another reciprocal, and so on.
+    Returns the continued fraction representation (canonical form) as
+    a list of integers, optionally ending (for quadratic irrationals)
+    with repeating block as the last term of this list.
 
-    Here the three parameters are,
-    * num: the rational part of the number's numerator
-    * delta: the irrational part(discriminator) of the number's numerator
-    * denominator: the denominator of the number
-
-    ex: the golden ratio is (1 + sqrt(5))/2
-    num = 1
-    delta = 5
-    den = 2
+    Parameters
+    ==========
 
-    The denominator of a rational number cannot be zero. So such
-    input will result an error.
+    p : int
+        the rational part of the number's numerator
+    q : int
+        the denominator of the number
+    d : int, optional
+        the irrational part (discriminator) of the number's numerator
 
-    >>> from sympy.ntheory.continued_fraction import continued_fraction
-    >>> continued_fraction(1,0,0)
-    Traceback (most recent call last):
-    ...
-    ValueError: The denominator is zero.
+    Examples
+    ========
 
-    If the discriminator is zero then the number will be a rational number.
+    >>> from sympy.ntheory.continued_fraction import continued_fraction_periodic
+    >>> continued_fraction_periodic(3, 2, 7)
+    [2, [1, 4, 1, 1]]
 
-    >>> from sympy.ntheory.continued_fraction import continued_fraction
-    >>> continued_fraction(4,3,0)
-    [1, 3]
+    Golden ratio has the simplest continued fraction expansion:
 
-    Golden ratio has the simplest continued fraction expansion,
+    >>> continued_fraction_periodic(1, 2, 5)
+    [[1]]
 
-    >>> from sympy.ntheory.continued_fraction import continued_fraction
-    >>> continued_fraction(1,2,5)
-    [[1], [1]]
+    If the discriminator is zero then the number will be a rational number:
 
-    Note: if the length of recurrsive part of the continued part of
-    expansion exceeds 100000 this module will truncate the convergents.
+    >>> continued_fraction_periodic(4, 3, 0)
+    [1, 3]
 
     See Also
     ========
 
-    continued_fraction_rational_number
-    sympy.ntheory.continued_fraction.continued_fraction_rational_number :
-    function which calculate the continued fraction of a rational number.
+    continued_fraction_iterator
 
     References
     ==========
 
-    [1] A. J. van der Poorten, "NOTES ON CONTINUED FRACTIONS AND RECURRENCE
-    SEQUENCES" in Number Theory and Cryptography. New York,
-    USA: Cambridge university press, 2011,
-    ch. 06, pp. 86-96.
-    [2] http://en.wikipedia.org/wiki/Continued_fraction
-    [3] http://www.numbertheory.org/ntw/N4.html#continued_fractions
-    [4] http://www.numbertheory.org/pdfs/CFquadratic.pdf
-    [5] http://maths.mq.edu.au/~alf/www-centre/alfpapers/a117.pdf
+    .. [1] http://en.wikipedia.org/wiki/Periodic_continued_fraction
+    .. [2] K. Rosen. Elementary Number theory and its applications.
+           Addison-Wesley, 3 Sub edition, pages 379-381, January 1992.
+
     """
+    from sympy.core.compatibility import as_int
+    from sympy.functions import sqrt
 
-    list = []
+    p, q, d = list(map(as_int, [p, q, d]))
 
-    # if the denominator is zero the expression cannot be a legal fraction
-    if den == 0:
+    if q == 0:
         raise ValueError("The denominator is zero.")
 
-    # if the discriminator is negative the number is a complex number
-    if delta < 0:
-        raise ValueError("The number is not real, so it does not\n\
-            have continued fraction expansion.")
-
-    # if the discriminator is zero the number is a rational number
-    if delta == 0:
-        list = continued_fraction_rational_number(
-            num, den)
-        return list
-
-    if (delta - (num*num)) % den != 0:
-        delta = delta*den*den
-        num = num*abs(den)
-        den = den*abs(den)
-
-    sqrtDelta = 0 | 1 << (delta.bit_length() + 1)/2
-    sqrtDeltaNext = ((delta/sqrtDelta) + sqrtDelta) >> 1
-
-    while sqrtDelta > sqrtDeltaNext:
-        sqrtDelta = sqrtDeltaNext
-        sqrtDeltaNext = ((delta/sqrtDelta) + sqrtDelta) >> 1
-
-    if sqrtDelta*sqrtDelta == delta:
-        continued_fraction_rational_number(num + sqrtDelta, den)
-        return list
-
-    biP = den
-
-    if biP > 0:
-        biK = sqrtDelta
-    else:
-        biK = sqrtDelta + 1
-
-    biK = biK + num
-
-    if biK > 0:
-        if den > 0:
-            biM = biK/den
-        else:
-            biM = ((den + 1) - biK)/(den*-1)
-    else:
-        if den > 0:
-            biM = ((biK + 1) - den)/den
-        else:
-            biM = (biK*-1)/(den*-1)
-    # appends the integer part of the continued fraction expansion
-    # to the result list
-    list.append(biM)
-    biM = ((biM*den) - num)
-    cont = -1
-    K = -1
-    P = -1
-    L = -1
-    M = -1
-
-    while (cont < 0 or K != P or L != M):
-
-        if (cont < 0 and biP > 0 and biP <= sqrtDelta + biM and
-                biM > 0 and biM <= sqrtDelta):
-            K = P = biP
-            L = M = biM
-            cont = 0
-
-        # both numerator and denominator are positive
-        if cont >= 0:
-            P = (delta - (M*M))/P
-            Z = (sqrtDelta + M)/P
-            M = (Z*P) - M
-            cont += 1
-        else:
-            biP = (delta - (biM*biM))/biP
-            if biP > 0:
-                Z = (sqrtDelta + biM)/biP
-            else:
-                Z = ((sqrtDelta + 1) + biM)/biP
-            biM = (Z*biP) - biM
-
-        # show convergent
-        list.append(Z)
-
-        # too many convergent
-        if cont > 100000:
-            raise NotImplementedError('more than 100,000 convergents')
-
-    if cont >= 1:
-
-        non_recurring_part = list[:len(list) - cont]
-
-        recurring_part = list[len(list) - cont:]
-
-        return [non_recurring_part, recurring_part]
-
-    return list
-
-
-def continued_fraction_rational_number(n, d):
-    """This applies the continued fraction expansion to two numbers
-    numerator/denominator
-
-    >>> from sympy.ntheory.continued_fraction import\
-    continued_fraction_rational_number
-    >>> continued_fraction_rational_number(3, 8)
+    if d < 0:
+        raise ValueError("Delta supposed to be a non-negative "
+                         "integer, got %d" % d)
+    elif d == 0:
+        # the number is a rational number
+        return list(continued_fraction_iterator(Rational(p, q)))
+
+    if (d - p**2)%q:
+        d *= q**2
+        p *= abs(q)
+        q *= abs(q)
+
+    terms = []
+    pq = {}
+    sd = sqrt(d)
+
+    while (p, q) not in pq:
+        pq[(p, q)] = len(terms)
+        terms.append(int((p + sd)/q))
+        p = terms[-1]*q - p
+        q = (d - p**2)/q
+
+    i = pq[(p, q)]
+    return terms[:i] + [terms[i:]]
+
+
+def continued_fraction_iterator(x):
+    """
+    Return continued fraction expansion of x as iterator.
+
+    Examples
+    ========
+
+    >>> from sympy.core import Rational, pi
+    >>> from sympy.ntheory.continued_fraction import continued_fraction_iterator
+
+    >>> list(continued_fraction_iterator(Rational(3, 8)))
     [0, 2, 1, 2]
+
+    >>> for i, v in enumerate(continued_fraction_iterator(pi)):
+    ...    if i > 7:
+    ...        break
+    ...    print(v)
+    3
+    7
+    15
+    1
+    292
+    1
+    1
+    1
+
+    References
+    ==========
+
+    .. [1] http://en.wikipedia.org/wiki/Continued_fraction
+
     """
-    x = Rational(n, d)
-    list = []
+
     while True:
         i = Integer(x)
-        list.append(i)
+        yield i
         x -= i
         if not x:
             break
         x = 1/x
-    return list

sympy/ntheory/tests/test_ntheory.py

@@ -1,7 +1,7 @@
 from collections import defaultdict
 from sympy import Sieve, binomial_coefficients, binomial_coefficients_list, \
     multinomial_coefficients, Mul, S, Pow, sieve, Symbol, summation, Dummy, \
-    factorial as fac, Rational
+    factorial as fac, Rational, pi, GoldenRatio as phi
 from sympy.core.numbers import Integer, igcd
 from sympy.core.compatibility import long
 
@@ -19,7 +19,9 @@
 from sympy.ntheory.primetest import _mr_safe_helper, mr
 from sympy.ntheory.bbp_pi import pi_hex_digits
 from sympy.ntheory.modular import crt, crt1, crt2, solve_congruence
-from sympy.ntheory.continued_fraction import continued_fraction
+from sympy.ntheory.continued_fraction import \
+    (continued_fraction_periodic as cf_p,
+     continued_fraction_iterator as cf_i)
 
 from sympy.polys.domains import ZZ
 
@@ -731,15 +733,31 @@ def test_search():
 
 
 def test_continued_fraction():
-    assert continued_fraction(4,3,0) == [1, 3]
-    assert continued_fraction(0,3,5) == [[0, 1], [2, 1, 12, 1, 2, 2]]
-    assert continued_fraction(1,1,0) == [1]
-    assert continued_fraction(3,4,0) == [0, 1, 3]
-    assert continued_fraction(4,5,0) == [0, 1, 4]
-    assert continued_fraction(5,6,0) == [0, 1, 5]
-    assert continued_fraction(11,13,0) == [0, 1, 5, 2]
-    assert continued_fraction(16,19,0) == [0, 1, 5, 3]
-    assert continued_fraction(27,32,0) == [0, 1, 5, 2, 2]
-    assert continued_fraction(1,2,5) == [[1], [1]]
-    assert continued_fraction(0,1,2) == [[1], [2]]
-    assert continued_fraction(3796,1387,0) == [2, 1, 2, 1, 4]
+    raises(ValueError, lambda: cf_p(1, 0, 0))
+    raises(ValueError, lambda: cf_p(1, 1, -1))
+    assert cf_p(4, 3, 0) == [1, 3]
+    assert cf_p(0, 3, 5) == [0, 1, [2, 1, 12, 1, 2, 2]]
+    assert cf_p(1, 1, 0) == [1]
+    assert cf_p(3, 4, 0) == [0, 1, 3]
+    assert cf_p(4, 5, 0) == [0, 1, 4]
+    assert cf_p(5, 6, 0) == [0, 1, 5]
+    assert cf_p(11, 13, 0) == [0, 1, 5, 2]
+    assert cf_p(16, 19, 0) == [0, 1, 5, 3]
+    assert cf_p(27, 32, 0) == [0, 1, 5, 2, 2]
+    assert cf_p(1, 2, 5) == [[1]]
+    assert cf_p(0, 1, 2) == [1, [2]]
+    assert cf_p(3796, 1387, 0) == [2, 1, 2, 1, 4]
+    assert cf_p(3245, 10000) == [0, 3, 12, 4, 13]
+    assert cf_p(1932, 2568) == [0, 1, 3, 26, 2]
+    assert cf_p(6589, 2569) == [2, 1, 1, 3, 2, 1, 3, 1, 23]
+
+    def take(iterator, n=7):
+        res = []
+        for i, t in enumerate(cf_i(iterator)):
+            if i >= n:
+                break
+            res.append(t)
+        return res
+
+    assert take(phi) == [1, 1, 1, 1, 1, 1, 1]
+    assert take(pi) == [3, 7, 15, 1, 292, 1, 1]

sympy/physics/quantum/shor.py

@@ -15,6 +15,7 @@
 from sympy import Mul, S
 from sympy import log, sqrt
 from sympy.core.numbers import igcd
+from sympy.ntheory import continued_fraction_periodic as continued_fraction
 from sympy.utilities.iterables import variations
 
 from sympy.physics.quantum.gate import Gate
@@ -124,26 +125,6 @@ def ratioize(list, N):
     return list[0] + ratioize(list[1:], N)
 
 
-def continued_fraction(x, y):
-    """This applies the continued fraction expansion to two numbers x/y
-
-    x is the numerator and y is the denominator
-
-    >>> from sympy.physics.quantum.shor import continued_fraction
-    >>> continued_fraction(3, 8)
-    [0, 2, 1, 2]
-    """
-    x = int(x)
-    y = int(y)
-    temp = x//y
-    if temp*y == x:
-        return [temp, ]
-
-    list = continued_fraction(y, x - temp*y)
-    list.insert(0, temp)
-    return list
-
-
 def period_find(a, N):
     """Finds the period of a in modulo N arithmetic