474 lines
17 KiB
Python
474 lines
17 KiB
Python
|
"""Various algorithms for helping identifying numbers and sequences."""
|
||
|
|
|
||
|
|
|
||
|
|
from sympy.concrete.products import (Product, product)
|
||
|
|
from sympy.core import Function, S
|
||
|
|
from sympy.core.add import Add
|
||
|
|
from sympy.core.numbers import Integer, Rational
|
||
|
|
from sympy.core.symbol import Symbol, symbols
|
||
|
|
from sympy.core.sympify import sympify
|
||
|
|
from sympy.functions.elementary.exponential import exp
|
||
|
|
from sympy.functions.elementary.integers import floor
|
||
|
|
from sympy.integrals.integrals import integrate
|
||
|
|
from sympy.polys.polyfuncs import rational_interpolate as rinterp
|
||
|
|
from sympy.polys.polytools import lcm
|
||
|
|
from sympy.simplify.radsimp import denom
|
||
|
|
from sympy.utilities import public
|
||
|
|
|
||
|
|
|
||
|
|
@public
|
||
|
|
def find_simple_recurrence_vector(l):
|
||
|
|
"""
|
||
|
|
This function is used internally by other functions from the
|
||
|
|
sympy.concrete.guess module. While most users may want to rather use the
|
||
|
|
function find_simple_recurrence when looking for recurrence relations
|
||
|
|
among rational numbers, the current function may still be useful when
|
||
|
|
some post-processing has to be done.
|
||
|
|
|
||
|
|
Explanation
|
||
|
|
===========
|
||
|
|
|
||
|
|
The function returns a vector of length n when a recurrence relation of
|
||
|
|
order n is detected in the sequence of rational numbers v.
|
||
|
|
|
||
|
|
If the returned vector has a length 1, then the returned value is always
|
||
|
|
the list [0], which means that no relation has been found.
|
||
|
|
|
||
|
|
While the functions is intended to be used with rational numbers, it should
|
||
|
|
work for other kinds of real numbers except for some cases involving
|
||
|
|
quadratic numbers; for that reason it should be used with some caution when
|
||
|
|
the argument is not a list of rational numbers.
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy.concrete.guess import find_simple_recurrence_vector
|
||
|
|
>>> from sympy import fibonacci
|
||
|
|
>>> find_simple_recurrence_vector([fibonacci(k) for k in range(12)])
|
||
|
|
[1, -1, -1]
|
||
|
|
|
||
|
|
See Also
|
||
|
|
========
|
||
|
|
|
||
|
|
See the function sympy.concrete.guess.find_simple_recurrence which is more
|
||
|
|
user-friendly.
|
||
|
|
|
||
|
|
"""
|
||
|
|
q1 = [0]
|
||
|
|
q2 = [1]
|
||
|
|
b, z = 0, len(l) >> 1
|
||
|
|
while len(q2) <= z:
|
||
|
|
while l[b]==0:
|
||
|
|
b += 1
|
||
|
|
if b == len(l):
|
||
|
|
c = 1
|
||
|
|
for x in q2:
|
||
|
|
c = lcm(c, denom(x))
|
||
|
|
if q2[0]*c < 0: c = -c
|
||
|
|
for k in range(len(q2)):
|
||
|
|
q2[k] = int(q2[k]*c)
|
||
|
|
return q2
|
||
|
|
a = S.One/l[b]
|
||
|
|
m = [a]
|
||
|
|
for k in range(b+1, len(l)):
|
||
|
|
m.append(-sum(l[j+1]*m[b-j-1] for j in range(b, k))*a)
|
||
|
|
l, m = m, [0] * max(len(q2), b+len(q1))
|
||
|
|
for k, q in enumerate(q2):
|
||
|
|
m[k] = a*q
|
||
|
|
for k, q in enumerate(q1):
|
||
|
|
m[k+b] += q
|
||
|
|
while m[-1]==0: m.pop() # because trailing zeros can occur
|
||
|
|
q1, q2, b = q2, m, 1
|
||
|
|
return [0]
|
||
|
|
|
||
|
|
@public
|
||
|
|
def find_simple_recurrence(v, A=Function('a'), N=Symbol('n')):
|
||
|
|
"""
|
||
|
|
Detects and returns a recurrence relation from a sequence of several integer
|
||
|
|
(or rational) terms. The name of the function in the returned expression is
|
||
|
|
'a' by default; the main variable is 'n' by default. The smallest index in
|
||
|
|
the returned expression is always n (and never n-1, n-2, etc.).
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy.concrete.guess import find_simple_recurrence
|
||
|
|
>>> from sympy import fibonacci
|
||
|
|
>>> find_simple_recurrence([fibonacci(k) for k in range(12)])
|
||
|
|
-a(n) - a(n + 1) + a(n + 2)
|
||
|
|
|
||
|
|
>>> from sympy import Function, Symbol
|
||
|
|
>>> a = [1, 1, 1]
|
||
|
|
>>> for k in range(15): a.append(5*a[-1]-3*a[-2]+8*a[-3])
|
||
|
|
>>> find_simple_recurrence(a, A=Function('f'), N=Symbol('i'))
|
||
|
|
-8*f(i) + 3*f(i + 1) - 5*f(i + 2) + f(i + 3)
|
||
|
|
|
||
|
|
"""
|
||
|
|
p = find_simple_recurrence_vector(v)
|
||
|
|
n = len(p)
|
||
|
|
if n <= 1: return S.Zero
|
||
|
|
|
||
|
|
return Add(*[A(N+n-1-k)*p[k] for k in range(n)])
|
||
|
|
|
||
|
|
|
||
|
|
@public
|
||
|
|
def rationalize(x, maxcoeff=10000):
|
||
|
|
"""
|
||
|
|
Helps identifying a rational number from a float (or mpmath.mpf) value by
|
||
|
|
using a continued fraction. The algorithm stops as soon as a large partial
|
||
|
|
quotient is detected (greater than 10000 by default).
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy.concrete.guess import rationalize
|
||
|
|
>>> from mpmath import cos, pi
|
||
|
|
>>> rationalize(cos(pi/3))
|
||
|
|
1/2
|
||
|
|
|
||
|
|
>>> from mpmath import mpf
|
||
|
|
>>> rationalize(mpf("0.333333333333333"))
|
||
|
|
1/3
|
||
|
|
|
||
|
|
While the function is rather intended to help 'identifying' rational
|
||
|
|
values, it may be used in some cases for approximating real numbers.
|
||
|
|
(Though other functions may be more relevant in that case.)
|
||
|
|
|
||
|
|
>>> rationalize(pi, maxcoeff = 250)
|
||
|
|
355/113
|
||
|
|
|
||
|
|
See Also
|
||
|
|
========
|
||
|
|
|
||
|
|
Several other methods can approximate a real number as a rational, like:
|
||
|
|
|
||
|
|
* fractions.Fraction.from_decimal
|
||
|
|
* fractions.Fraction.from_float
|
||
|
|
* mpmath.identify
|
||
|
|
* mpmath.pslq by using the following syntax: mpmath.pslq([x, 1])
|
||
|
|
* mpmath.findpoly by using the following syntax: mpmath.findpoly(x, 1)
|
||
|
|
* sympy.simplify.nsimplify (which is a more general function)
|
||
|
|
|
||
|
|
The main difference between the current function and all these variants is
|
||
|
|
that control focuses on magnitude of partial quotients here rather than on
|
||
|
|
global precision of the approximation. If the real is "known to be" a
|
||
|
|
rational number, the current function should be able to detect it correctly
|
||
|
|
with the default settings even when denominator is great (unless its
|
||
|
|
expansion contains unusually big partial quotients) which may occur
|
||
|
|
when studying sequences of increasing numbers. If the user cares more
|
||
|
|
on getting simple fractions, other methods may be more convenient.
|
||
|
|
|
||
|
|
"""
|
||
|
|
p0, p1 = 0, 1
|
||
|
|
q0, q1 = 1, 0
|
||
|
|
a = floor(x)
|
||
|
|
while a < maxcoeff or q1==0:
|
||
|
|
p = a*p1 + p0
|
||
|
|
q = a*q1 + q0
|
||
|
|
p0, p1 = p1, p
|
||
|
|
q0, q1 = q1, q
|
||
|
|
if x==a: break
|
||
|
|
x = 1/(x-a)
|
||
|
|
a = floor(x)
|
||
|
|
return sympify(p) / q
|
||
|
|
|
||
|
|
|
||
|
|
@public
|
||
|
|
def guess_generating_function_rational(v, X=Symbol('x')):
|
||
|
|
"""
|
||
|
|
Tries to "guess" a rational generating function for a sequence of rational
|
||
|
|
numbers v.
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy.concrete.guess import guess_generating_function_rational
|
||
|
|
>>> from sympy import fibonacci
|
||
|
|
>>> l = [fibonacci(k) for k in range(5,15)]
|
||
|
|
>>> guess_generating_function_rational(l)
|
||
|
|
(3*x + 5)/(-x**2 - x + 1)
|
||
|
|
|
||
|
|
See Also
|
||
|
|
========
|
||
|
|
|
||
|
|
sympy.series.approximants
|
||
|
|
mpmath.pade
|
||
|
|
|
||
|
|
"""
|
||
|
|
# a) compute the denominator as q
|
||
|
|
q = find_simple_recurrence_vector(v)
|
||
|
|
n = len(q)
|
||
|
|
if n <= 1: return None
|
||
|
|
# b) compute the numerator as p
|
||
|
|
p = [sum(v[i-k]*q[k] for k in range(min(i+1, n)))
|
||
|
|
for i in range(len(v)>>1)]
|
||
|
|
return (sum(p[k]*X**k for k in range(len(p)))
|
||
|
|
/ sum(q[k]*X**k for k in range(n)))
|
||
|
|
|
||
|
|
|
||
|
|
@public
|
||
|
|
def guess_generating_function(v, X=Symbol('x'), types=['all'], maxsqrtn=2):
|
||
|
|
"""
|
||
|
|
Tries to "guess" a generating function for a sequence of rational numbers v.
|
||
|
|
Only a few patterns are implemented yet.
|
||
|
|
|
||
|
|
Explanation
|
||
|
|
===========
|
||
|
|
|
||
|
|
The function returns a dictionary where keys are the name of a given type of
|
||
|
|
generating function. Six types are currently implemented:
|
||
|
|
|
||
|
|
type | formal definition
|
||
|
|
-------+----------------------------------------------------------------
|
||
|
|
ogf | f(x) = Sum( a_k * x^k , k: 0..infinity )
|
||
|
|
egf | f(x) = Sum( a_k * x^k / k! , k: 0..infinity )
|
||
|
|
lgf | f(x) = Sum( (-1)^(k+1) a_k * x^k / k , k: 1..infinity )
|
||
|
|
| (with initial index being hold as 1 rather than 0)
|
||
|
|
hlgf | f(x) = Sum( a_k * x^k / k , k: 1..infinity )
|
||
|
|
| (with initial index being hold as 1 rather than 0)
|
||
|
|
lgdogf | f(x) = derivate( log(Sum( a_k * x^k, k: 0..infinity )), x)
|
||
|
|
lgdegf | f(x) = derivate( log(Sum( a_k * x^k / k!, k: 0..infinity )), x)
|
||
|
|
|
||
|
|
In order to spare time, the user can select only some types of generating
|
||
|
|
functions (default being ['all']). While forgetting to use a list in the
|
||
|
|
case of a single type may seem to work most of the time as in: types='ogf'
|
||
|
|
this (convenient) syntax may lead to unexpected extra results in some cases.
|
||
|
|
|
||
|
|
Discarding a type when calling the function does not mean that the type will
|
||
|
|
not be present in the returned dictionary; it only means that no extra
|
||
|
|
computation will be performed for that type, but the function may still add
|
||
|
|
it in the result when it can be easily converted from another type.
|
||
|
|
|
||
|
|
Two generating functions (lgdogf and lgdegf) are not even computed if the
|
||
|
|
initial term of the sequence is 0; it may be useful in that case to try
|
||
|
|
again after having removed the leading zeros.
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy.concrete.guess import guess_generating_function as ggf
|
||
|
|
>>> ggf([k+1 for k in range(12)], types=['ogf', 'lgf', 'hlgf'])
|
||
|
|
{'hlgf': 1/(1 - x), 'lgf': 1/(x + 1), 'ogf': 1/(x**2 - 2*x + 1)}
|
||
|
|
|
||
|
|
>>> from sympy import sympify
|
||
|
|
>>> l = sympify("[3/2, 11/2, 0, -121/2, -363/2, 121]")
|
||
|
|
>>> ggf(l)
|
||
|
|
{'ogf': (x + 3/2)/(11*x**2 - 3*x + 1)}
|
||
|
|
|
||
|
|
>>> from sympy import fibonacci
|
||
|
|
>>> ggf([fibonacci(k) for k in range(5, 15)], types=['ogf'])
|
||
|
|
{'ogf': (3*x + 5)/(-x**2 - x + 1)}
|
||
|
|
|
||
|
|
>>> from sympy import factorial
|
||
|
|
>>> ggf([factorial(k) for k in range(12)], types=['ogf', 'egf', 'lgf'])
|
||
|
|
{'egf': 1/(1 - x)}
|
||
|
|
|
||
|
|
>>> ggf([k+1 for k in range(12)], types=['egf'])
|
||
|
|
{'egf': (x + 1)*exp(x), 'lgdegf': (x + 2)/(x + 1)}
|
||
|
|
|
||
|
|
N-th root of a rational function can also be detected (below is an example
|
||
|
|
coming from the sequence A108626 from https://oeis.org).
|
||
|
|
The greatest n-th root to be tested is specified as maxsqrtn (default 2).
|
||
|
|
|
||
|
|
>>> ggf([1, 2, 5, 14, 41, 124, 383, 1200, 3799, 12122, 38919])['ogf']
|
||
|
|
sqrt(1/(x**4 + 2*x**2 - 4*x + 1))
|
||
|
|
|
||
|
|
References
|
||
|
|
==========
|
||
|
|
|
||
|
|
.. [1] "Concrete Mathematics", R.L. Graham, D.E. Knuth, O. Patashnik
|
||
|
|
.. [2] https://oeis.org/wiki/Generating_functions
|
||
|
|
|
||
|
|
"""
|
||
|
|
# List of all types of all g.f. known by the algorithm
|
||
|
|
if 'all' in types:
|
||
|
|
types = ('ogf', 'egf', 'lgf', 'hlgf', 'lgdogf', 'lgdegf')
|
||
|
|
|
||
|
|
result = {}
|
||
|
|
|
||
|
|
# Ordinary Generating Function (ogf)
|
||
|
|
if 'ogf' in types:
|
||
|
|
# Perform some convolutions of the sequence with itself
|
||
|
|
t = [1] + [0]*(len(v) - 1)
|
||
|
|
for d in range(max(1, maxsqrtn)):
|
||
|
|
t = [sum(t[n-i]*v[i] for i in range(n+1)) for n in range(len(v))]
|
||
|
|
g = guess_generating_function_rational(t, X=X)
|
||
|
|
if g:
|
||
|
|
result['ogf'] = g**Rational(1, d+1)
|
||
|
|
break
|
||
|
|
|
||
|
|
# Exponential Generating Function (egf)
|
||
|
|
if 'egf' in types:
|
||
|
|
# Transform sequence (division by factorial)
|
||
|
|
w, f = [], S.One
|
||
|
|
for i, k in enumerate(v):
|
||
|
|
f *= i if i else 1
|
||
|
|
w.append(k/f)
|
||
|
|
# Perform some convolutions of the sequence with itself
|
||
|
|
t = [1] + [0]*(len(w) - 1)
|
||
|
|
for d in range(max(1, maxsqrtn)):
|
||
|
|
t = [sum(t[n-i]*w[i] for i in range(n+1)) for n in range(len(w))]
|
||
|
|
g = guess_generating_function_rational(t, X=X)
|
||
|
|
if g:
|
||
|
|
result['egf'] = g**Rational(1, d+1)
|
||
|
|
break
|
||
|
|
|
||
|
|
# Logarithmic Generating Function (lgf)
|
||
|
|
if 'lgf' in types:
|
||
|
|
# Transform sequence (multiplication by (-1)^(n+1) / n)
|
||
|
|
w, f = [], S.NegativeOne
|
||
|
|
for i, k in enumerate(v):
|
||
|
|
f = -f
|
||
|
|
w.append(f*k/Integer(i+1))
|
||
|
|
# Perform some convolutions of the sequence with itself
|
||
|
|
t = [1] + [0]*(len(w) - 1)
|
||
|
|
for d in range(max(1, maxsqrtn)):
|
||
|
|
t = [sum(t[n-i]*w[i] for i in range(n+1)) for n in range(len(w))]
|
||
|
|
g = guess_generating_function_rational(t, X=X)
|
||
|
|
if g:
|
||
|
|
result['lgf'] = g**Rational(1, d+1)
|
||
|
|
break
|
||
|
|
|
||
|
|
# Hyperbolic logarithmic Generating Function (hlgf)
|
||
|
|
if 'hlgf' in types:
|
||
|
|
# Transform sequence (division by n+1)
|
||
|
|
w = []
|
||
|
|
for i, k in enumerate(v):
|
||
|
|
w.append(k/Integer(i+1))
|
||
|
|
# Perform some convolutions of the sequence with itself
|
||
|
|
t = [1] + [0]*(len(w) - 1)
|
||
|
|
for d in range(max(1, maxsqrtn)):
|
||
|
|
t = [sum(t[n-i]*w[i] for i in range(n+1)) for n in range(len(w))]
|
||
|
|
g = guess_generating_function_rational(t, X=X)
|
||
|
|
if g:
|
||
|
|
result['hlgf'] = g**Rational(1, d+1)
|
||
|
|
break
|
||
|
|
|
||
|
|
# Logarithmic derivative of ordinary generating Function (lgdogf)
|
||
|
|
if v[0] != 0 and ('lgdogf' in types
|
||
|
|
or ('ogf' in types and 'ogf' not in result)):
|
||
|
|
# Transform sequence by computing f'(x)/f(x)
|
||
|
|
# because log(f(x)) = integrate( f'(x)/f(x) )
|
||
|
|
a, w = sympify(v[0]), []
|
||
|
|
for n in range(len(v)-1):
|
||
|
|
w.append(
|
||
|
|
(v[n+1]*(n+1) - sum(w[-i-1]*v[i+1] for i in range(n)))/a)
|
||
|
|
# Perform some convolutions of the sequence with itself
|
||
|
|
t = [1] + [0]*(len(w) - 1)
|
||
|
|
for d in range(max(1, maxsqrtn)):
|
||
|
|
t = [sum(t[n-i]*w[i] for i in range(n+1)) for n in range(len(w))]
|
||
|
|
g = guess_generating_function_rational(t, X=X)
|
||
|
|
if g:
|
||
|
|
result['lgdogf'] = g**Rational(1, d+1)
|
||
|
|
if 'ogf' not in result:
|
||
|
|
result['ogf'] = exp(integrate(result['lgdogf'], X))
|
||
|
|
break
|
||
|
|
|
||
|
|
# Logarithmic derivative of exponential generating Function (lgdegf)
|
||
|
|
if v[0] != 0 and ('lgdegf' in types
|
||
|
|
or ('egf' in types and 'egf' not in result)):
|
||
|
|
# Transform sequence / step 1 (division by factorial)
|
||
|
|
z, f = [], S.One
|
||
|
|
for i, k in enumerate(v):
|
||
|
|
f *= i if i else 1
|
||
|
|
z.append(k/f)
|
||
|
|
# Transform sequence / step 2 by computing f'(x)/f(x)
|
||
|
|
# because log(f(x)) = integrate( f'(x)/f(x) )
|
||
|
|
a, w = z[0], []
|
||
|
|
for n in range(len(z)-1):
|
||
|
|
w.append(
|
||
|
|
(z[n+1]*(n+1) - sum(w[-i-1]*z[i+1] for i in range(n)))/a)
|
||
|
|
# Perform some convolutions of the sequence with itself
|
||
|
|
t = [1] + [0]*(len(w) - 1)
|
||
|
|
for d in range(max(1, maxsqrtn)):
|
||
|
|
t = [sum(t[n-i]*w[i] for i in range(n+1)) for n in range(len(w))]
|
||
|
|
g = guess_generating_function_rational(t, X=X)
|
||
|
|
if g:
|
||
|
|
result['lgdegf'] = g**Rational(1, d+1)
|
||
|
|
if 'egf' not in result:
|
||
|
|
result['egf'] = exp(integrate(result['lgdegf'], X))
|
||
|
|
break
|
||
|
|
|
||
|
|
return result
|
||
|
|
|
||
|
|
|
||
|
|
@public
|
||
|
|
def guess(l, all=False, evaluate=True, niter=2, variables=None):
|
||
|
|
"""
|
||
|
|
This function is adapted from the Rate.m package for Mathematica
|
||
|
|
written by Christian Krattenthaler.
|
||
|
|
It tries to guess a formula from a given sequence of rational numbers.
|
||
|
|
|
||
|
|
Explanation
|
||
|
|
===========
|
||
|
|
|
||
|
|
In order to speed up the process, the 'all' variable is set to False by
|
||
|
|
default, stopping the computation as some results are returned during an
|
||
|
|
iteration; the variable can be set to True if more iterations are needed
|
||
|
|
(other formulas may be found; however they may be equivalent to the first
|
||
|
|
ones).
|
||
|
|
|
||
|
|
Another option is the 'evaluate' variable (default is True); setting it
|
||
|
|
to False will leave the involved products unevaluated.
|
||
|
|
|
||
|
|
By default, the number of iterations is set to 2 but a greater value (up
|
||
|
|
to len(l)-1) can be specified with the optional 'niter' variable.
|
||
|
|
More and more convoluted results are found when the order of the
|
||
|
|
iteration gets higher:
|
||
|
|
|
||
|
|
* first iteration returns polynomial or rational functions;
|
||
|
|
* second iteration returns products of rising factorials and their
|
||
|
|
inverses;
|
||
|
|
* third iteration returns products of products of rising factorials
|
||
|
|
and their inverses;
|
||
|
|
* etc.
|
||
|
|
|
||
|
|
The returned formulas contain symbols i0, i1, i2, ... where the main
|
||
|
|
variables is i0 (and auxiliary variables are i1, i2, ...). A list of
|
||
|
|
other symbols can be provided in the 'variables' option; the length of
|
||
|
|
the least should be the value of 'niter' (more is acceptable but only
|
||
|
|
the first symbols will be used); in this case, the main variable will be
|
||
|
|
the first symbol in the list.
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy.concrete.guess import guess
|
||
|
|
>>> guess([1,2,6,24,120], evaluate=False)
|
||
|
|
[Product(i1 + 1, (i1, 1, i0 - 1))]
|
||
|
|
|
||
|
|
>>> from sympy import symbols
|
||
|
|
>>> r = guess([1,2,7,42,429,7436,218348,10850216], niter=4)
|
||
|
|
>>> i0 = symbols("i0")
|
||
|
|
>>> [r[0].subs(i0,n).doit() for n in range(1,10)]
|
||
|
|
[1, 2, 7, 42, 429, 7436, 218348, 10850216, 911835460]
|
||
|
|
"""
|
||
|
|
if any(a==0 for a in l[:-1]):
|
||
|
|
return []
|
||
|
|
N = len(l)
|
||
|
|
niter = min(N-1, niter)
|
||
|
|
myprod = product if evaluate else Product
|
||
|
|
g = []
|
||
|
|
res = []
|
||
|
|
if variables is None:
|
||
|
|
symb = symbols('i:'+str(niter))
|
||
|
|
else:
|
||
|
|
symb = variables
|
||
|
|
for k, s in enumerate(symb):
|
||
|
|
g.append(l)
|
||
|
|
n, r = len(l), []
|
||
|
|
for i in range(n-2-1, -1, -1):
|
||
|
|
ri = rinterp(enumerate(g[k][:-1], start=1), i, X=s)
|
||
|
|
if ((denom(ri).subs({s:n}) != 0)
|
||
|
|
and (ri.subs({s:n}) - g[k][-1] == 0)
|
||
|
|
and ri not in r):
|
||
|
|
r.append(ri)
|
||
|
|
if r:
|
||
|
|
for i in range(k-1, -1, -1):
|
||
|
|
r = [g[i][0]
|
||
|
|
* myprod(v, (symb[i+1], 1, symb[i]-1)) for v in r]
|
||
|
|
if not all: return r
|
||
|
|
res += r
|
||
|
|
l = [Rational(l[i+1], l[i]) for i in range(N-k-1)]
|
||
|
|
return res
|