<s> """<STR_LIT>""" <EOL> from __future__ import print_function , division <EOL> from collections import defaultdict <EOL> from sympy . core . singleton import S <EOL> from sympy . core . numbers import Rational , I <EOL> from sympy . core . symbol import Symbol , Wild , Dummy <EOL> from sympy . core . relational import Equality <EOL> from sympy . core . add import Add <EOL> from sympy . core . mul import Mul <EOL> from sympy . core import sympify <EOL> from sympy . simplify import simplify , hypersimp , hypersimilar <EOL> from sympy . solvers import solve , solve_undetermined_coeffs <EOL> from sympy . polys import Poly , quo , gcd , lcm , roots , resultant <EOL> from sympy . functions import binomial , factorial , FallingFactorial , RisingFactorial <EOL> from sympy . matrices import Matrix , casoratian <EOL> from sympy . concrete import product <EOL> from sympy . core . compatibility import default_sort_key , range <EOL> from sympy . utilities . iterables import numbered_symbols <EOL> def rsolve_poly ( coeffs , f , n , ** hints ) : <EOL> """<STR_LIT>""" <EOL> f = sympify ( f ) <EOL> if not f . is_polynomial ( n ) : <EOL> return None <EOL> homogeneous = f . is_zero <EOL> r = len ( coeffs ) - <NUM_LIT:1> <EOL> coeffs = [ Poly ( coeff , n ) for coeff in coeffs ] <EOL> polys = [ Poly ( <NUM_LIT:0> , n ) ] * ( r + <NUM_LIT:1> ) <EOL> terms = [ ( S . Zero , S . NegativeInfinity ) ] * ( r + <NUM_LIT:1> ) <EOL> for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> for j in range ( i , r + <NUM_LIT:1> ) : <EOL> polys [ i ] += coeffs [ j ] * binomial ( j , i ) <EOL> if not polys [ i ] . is_zero : <EOL> ( exp , ) , coeff = polys [ i ] . LT ( ) <EOL> terms [ i ] = ( coeff , exp ) <EOL> d = b = terms [ <NUM_LIT:0> ] [ <NUM_LIT:1> ] <EOL> for i in range ( <NUM_LIT:1> , r + <NUM_LIT:1> ) : <EOL> if terms [ i ] [ <NUM_LIT:1> ] > d : <EOL> d = terms [ i ] [ <NUM_LIT:1> ] <EOL> if terms [ i ] [ <NUM_LIT:1> ] - i > b : <EOL> b = terms [ i ] [ <NUM_LIT:1> ] - i <EOL> d , b = int ( d ) , int ( b ) <EOL> x = Dummy ( '<STR_LIT:x>' ) <EOL> degree_poly = S . Zero <EOL> for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> if terms [ i ] [ <NUM_LIT:1> ] - i == b : <EOL> degree_poly += terms [ i ] [ <NUM_LIT:0> ] * FallingFactorial ( x , i ) <EOL> nni_roots = list ( roots ( degree_poly , x , filter = '<STR_LIT>' , <EOL> predicate = lambda r : r >= <NUM_LIT:0> ) . keys ( ) ) <EOL> if nni_roots : <EOL> N = [ max ( nni_roots ) ] <EOL> else : <EOL> N = [ ] <EOL> if homogeneous : <EOL> N += [ - b - <NUM_LIT:1> ] <EOL> else : <EOL> N += [ f . as_poly ( n ) . degree ( ) - b , - b - <NUM_LIT:1> ] <EOL> N = int ( max ( N ) ) <EOL> if N < <NUM_LIT:0> : <EOL> if homogeneous : <EOL> if hints . get ( '<STR_LIT>' , False ) : <EOL> return ( S . Zero , [ ] ) <EOL> else : <EOL> return S . Zero <EOL> else : <EOL> return None <EOL> if N <= r : <EOL> C = [ ] <EOL> y = E = S . Zero <EOL> for i in range ( <NUM_LIT:0> , N + <NUM_LIT:1> ) : <EOL> C . append ( Symbol ( '<STR_LIT:C>' + str ( i ) ) ) <EOL> y += C [ i ] * n ** i <EOL> for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> E += coeffs [ i ] . as_expr ( ) * y . subs ( n , n + i ) <EOL> solutions = solve_undetermined_coeffs ( E - f , C , n ) <EOL> if solutions is not None : <EOL> C = [ c for c in C if ( c not in solutions ) ] <EOL> result = y . subs ( solutions ) <EOL> else : <EOL> return None <EOL> else : <EOL> A = r <EOL> U = N + A + b + <NUM_LIT:1> <EOL> nni_roots = list ( roots ( polys [ r ] , filter = '<STR_LIT>' , <EOL> predicate = lambda r : r >= <NUM_LIT:0> ) . keys ( ) ) <EOL> if nni_roots != [ ] : <EOL> a = max ( nni_roots ) + <NUM_LIT:1> <EOL> else : <EOL> a = S . Zero <EOL> def _zero_vector ( k ) : <EOL> return [ S . Zero ] * k <EOL> def _one_vector ( k ) : <EOL> return [ S . One ] * k <EOL> def _delta ( p , k ) : <EOL> B = S . One <EOL> D = p . subs ( n , a + k ) <EOL> for i in range ( <NUM_LIT:1> , k + <NUM_LIT:1> ) : <EOL> B *= - Rational ( k - i + <NUM_LIT:1> , i ) <EOL> D += B * p . subs ( n , a + k - i ) <EOL> return D <EOL> alpha = { } <EOL> for i in range ( - A , d + <NUM_LIT:1> ) : <EOL> I = _one_vector ( d + <NUM_LIT:1> ) <EOL> for k in range ( <NUM_LIT:1> , d + <NUM_LIT:1> ) : <EOL> I [ k ] = I [ k - <NUM_LIT:1> ] * ( x + i - k + <NUM_LIT:1> ) / k <EOL> alpha [ i ] = S . Zero <EOL> for j in range ( <NUM_LIT:0> , A + <NUM_LIT:1> ) : <EOL> for k in range ( <NUM_LIT:0> , d + <NUM_LIT:1> ) : <EOL> B = binomial ( k , i + j ) <EOL> D = _delta ( polys [ j ] . as_expr ( ) , k ) <EOL> alpha [ i ] += I [ k ] * B * D <EOL> V = Matrix ( U , A , lambda i , j : int ( i == j ) ) <EOL> if homogeneous : <EOL> for i in range ( A , U ) : <EOL> v = _zero_vector ( A ) <EOL> for k in range ( <NUM_LIT:1> , A + b + <NUM_LIT:1> ) : <EOL> if i - k < <NUM_LIT:0> : <EOL> break <EOL> B = alpha [ k - A ] . subs ( x , i - k ) <EOL> for j in range ( <NUM_LIT:0> , A ) : <EOL> v [ j ] += B * V [ i - k , j ] <EOL> denom = alpha [ - A ] . subs ( x , i ) <EOL> for j in range ( <NUM_LIT:0> , A ) : <EOL> V [ i , j ] = - v [ j ] / denom <EOL> else : <EOL> G = _zero_vector ( U ) <EOL> for i in range ( A , U ) : <EOL> v = _zero_vector ( A ) <EOL> g = S . Zero <EOL> for k in range ( <NUM_LIT:1> , A + b + <NUM_LIT:1> ) : <EOL> if i - k < <NUM_LIT:0> : <EOL> break <EOL> B = alpha [ k - A ] . subs ( x , i - k ) <EOL> for j in range ( <NUM_LIT:0> , A ) : <EOL> v [ j ] += B * V [ i - k , j ] <EOL> g += B * G [ i - k ] <EOL> denom = alpha [ - A ] . subs ( x , i ) <EOL> for j in range ( <NUM_LIT:0> , A ) : <EOL> V [ i , j ] = - v [ j ] / denom <EOL> G [ i ] = ( _delta ( f , i - A ) - g ) / denom <EOL> P , Q = _one_vector ( U ) , _zero_vector ( A ) <EOL> for i in range ( <NUM_LIT:1> , U ) : <EOL> P [ i ] = ( P [ i - <NUM_LIT:1> ] * ( n - a - i + <NUM_LIT:1> ) / i ) . expand ( ) <EOL> for i in range ( <NUM_LIT:0> , A ) : <EOL> Q [ i ] = Add ( * [ ( v * p ) . expand ( ) for v , p in zip ( V [ : , i ] , P ) ] ) <EOL> if not homogeneous : <EOL> h = Add ( * [ ( g * p ) . expand ( ) for g , p in zip ( G , P ) ] ) <EOL> C = [ Symbol ( '<STR_LIT:C>' + str ( i ) ) for i in range ( <NUM_LIT:0> , A ) ] <EOL> g = lambda i : Add ( * [ c * _delta ( q , i ) for c , q in zip ( C , Q ) ] ) <EOL> if homogeneous : <EOL> E = [ g ( i ) for i in range ( N + <NUM_LIT:1> , U ) ] <EOL> else : <EOL> E = [ g ( i ) + _delta ( h , i ) for i in range ( N + <NUM_LIT:1> , U ) ] <EOL> if E != [ ] : <EOL> solutions = solve ( E , * C ) <EOL> if not solutions : <EOL> if homogeneous : <EOL> if hints . get ( '<STR_LIT>' , False ) : <EOL> return ( S . Zero , [ ] ) <EOL> else : <EOL> return S . Zero <EOL> else : <EOL> return None <EOL> else : <EOL> solutions = { } <EOL> if homogeneous : <EOL> result = S . Zero <EOL> else : <EOL> result = h <EOL> for c , q in list ( zip ( C , Q ) ) : <EOL> if c in solutions : <EOL> s = solutions [ c ] * q <EOL> C . remove ( c ) <EOL> else : <EOL> s = c * q <EOL> result += s . expand ( ) <EOL> if hints . get ( '<STR_LIT>' , False ) : <EOL> return ( result , C ) <EOL> else : <EOL> return result <EOL> def rsolve_ratio ( coeffs , f , n , ** hints ) : <EOL> """<STR_LIT>""" <EOL> f = sympify ( f ) <EOL> if not f . is_polynomial ( n ) : <EOL> return None <EOL> coeffs = list ( map ( sympify , coeffs ) ) <EOL> r = len ( coeffs ) - <NUM_LIT:1> <EOL> A , B = coeffs [ r ] , coeffs [ <NUM_LIT:0> ] <EOL> A = A . subs ( n , n - r ) . expand ( ) <EOL> h = Dummy ( '<STR_LIT:h>' ) <EOL> res = resultant ( A , B . subs ( n , n + h ) , n ) <EOL> if not res . is_polynomial ( h ) : <EOL> p , q = res . as_numer_denom ( ) <EOL> res = quo ( p , q , h ) <EOL> nni_roots = list ( roots ( res , h , filter = '<STR_LIT>' , <EOL> predicate = lambda r : r >= <NUM_LIT:0> ) . keys ( ) ) <EOL> if not nni_roots : <EOL> return rsolve_poly ( coeffs , f , n , ** hints ) <EOL> else : <EOL> C , numers = S . One , [ S . Zero ] * ( r + <NUM_LIT:1> ) <EOL> for i in range ( int ( max ( nni_roots ) ) , - <NUM_LIT:1> , - <NUM_LIT:1> ) : <EOL> d = gcd ( A , B . subs ( n , n + i ) , n ) <EOL> A = quo ( A , d , n ) <EOL> B = quo ( B , d . subs ( n , n - i ) , n ) <EOL> C *= Mul ( * [ d . subs ( n , n - j ) for j in range ( <NUM_LIT:0> , i + <NUM_LIT:1> ) ] ) <EOL> denoms = [ C . subs ( n , n + i ) for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) ] <EOL> for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> g = gcd ( coeffs [ i ] , denoms [ i ] , n ) <EOL> numers [ i ] = quo ( coeffs [ i ] , g , n ) <EOL> denoms [ i ] = quo ( denoms [ i ] , g , n ) <EOL> for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> numers [ i ] *= Mul ( * ( denoms [ : i ] + denoms [ i + <NUM_LIT:1> : ] ) ) <EOL> result = rsolve_poly ( numers , f * Mul ( * denoms ) , n , ** hints ) <EOL> if result is not None : <EOL> if hints . get ( '<STR_LIT>' , False ) : <EOL> return ( simplify ( result [ <NUM_LIT:0> ] / C ) , result [ <NUM_LIT:1> ] ) <EOL> else : <EOL> return simplify ( result / C ) <EOL> else : <EOL> return None <EOL> def rsolve_hyper ( coeffs , f , n , ** hints ) : <EOL> """<STR_LIT>""" <EOL> coeffs = list ( map ( sympify , coeffs ) ) <EOL> f = sympify ( f ) <EOL> r , kernel , symbols = len ( coeffs ) - <NUM_LIT:1> , [ ] , set ( ) <EOL> if not f . is_zero : <EOL> if f . is_Add : <EOL> similar = { } <EOL> for g in f . expand ( ) . args : <EOL> if not g . is_hypergeometric ( n ) : <EOL> return None <EOL> for h in similar . keys ( ) : <EOL> if hypersimilar ( g , h , n ) : <EOL> similar [ h ] += g <EOL> break <EOL> else : <EOL> similar [ g ] = S . Zero <EOL> inhomogeneous = [ ] <EOL> for g , h in similar . items ( ) : <EOL> inhomogeneous . append ( g + h ) <EOL> elif f . is_hypergeometric ( n ) : <EOL> inhomogeneous = [ f ] <EOL> else : <EOL> return None <EOL> for i , g in enumerate ( inhomogeneous ) : <EOL> coeff , polys = S . One , coeffs [ : ] <EOL> denoms = [ S . One ] * ( r + <NUM_LIT:1> ) <EOL> s = hypersimp ( g , n ) <EOL> for j in range ( <NUM_LIT:1> , r + <NUM_LIT:1> ) : <EOL> coeff *= s . subs ( n , n + j - <NUM_LIT:1> ) <EOL> p , q = coeff . as_numer_denom ( ) <EOL> polys [ j ] *= p <EOL> denoms [ j ] = q <EOL> for j in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> polys [ j ] *= Mul ( * ( denoms [ : j ] + denoms [ j + <NUM_LIT:1> : ] ) ) <EOL> R = rsolve_poly ( polys , Mul ( * denoms ) , n ) <EOL> if not ( R is None or R is S . Zero ) : <EOL> inhomogeneous [ i ] *= R <EOL> else : <EOL> return None <EOL> result = Add ( * inhomogeneous ) <EOL> else : <EOL> result = S . Zero <EOL> Z = Dummy ( '<STR_LIT>' ) <EOL> p , q = coeffs [ <NUM_LIT:0> ] , coeffs [ r ] . subs ( n , n - r + <NUM_LIT:1> ) <EOL> p_factors = [ z for z in roots ( p , n ) . keys ( ) ] <EOL> q_factors = [ z for z in roots ( q , n ) . keys ( ) ] <EOL> factors = [ ( S . One , S . One ) ] <EOL> for p in p_factors : <EOL> for q in q_factors : <EOL> if p . is_integer and q . is_integer and p <= q : <EOL> continue <EOL> else : <EOL> factors += [ ( n - p , n - q ) ] <EOL> p = [ ( n - p , S . One ) for p in p_factors ] <EOL> q = [ ( S . One , n - q ) for q in q_factors ] <EOL> factors = p + factors + q <EOL> for A , B in factors : <EOL> polys , degrees = [ ] , [ ] <EOL> D = A * B . subs ( n , n + r - <NUM_LIT:1> ) <EOL> for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> a = Mul ( * [ A . subs ( n , n + j ) for j in range ( <NUM_LIT:0> , i ) ] ) <EOL> b = Mul ( * [ B . subs ( n , n + j ) for j in range ( i , r ) ] ) <EOL> poly = quo ( coeffs [ i ] * a * b , D , n ) <EOL> polys . append ( poly . as_poly ( n ) ) <EOL> if not poly . is_zero : <EOL> degrees . append ( polys [ i ] . degree ( ) ) <EOL> if degrees : <EOL> d , poly = max ( degrees ) , S . Zero <EOL> else : <EOL> return None <EOL> for i in range ( <NUM_LIT:0> , r + <NUM_LIT:1> ) : <EOL> coeff = polys [ i ] . nth ( d ) <EOL> if coeff is not S . Zero : <EOL> poly += coeff * Z ** i <EOL> for z in roots ( poly , Z ) . keys ( ) : <EOL> if z . is_zero : <EOL> continue <EOL> ( C , s ) = rsolve_poly ( [ polys [ i ] * z ** i for i in range ( r + <NUM_LIT:1> ) ] , <NUM_LIT:0> , n , symbols = True ) <EOL> if C is not None and C is not S . Zero : <EOL> symbols |= set ( s ) <EOL> ratio = z * A * C . subs ( n , n + <NUM_LIT:1> ) / B / C <EOL> ratio = simplify ( ratio ) <EOL> n0 = <NUM_LIT:0> <EOL> for n_root in roots ( ratio . as_numer_denom ( ) [ <NUM_LIT:1> ] , n ) . keys ( ) : <EOL> if n_root . has ( I ) : <EOL> return None <EOL> elif ( n0 < ( n_root + <NUM_LIT:1> ) ) == True : <EOL> n0 = n_root + <NUM_LIT:1> <EOL> K = product ( ratio , ( n , n0 , n - <NUM_LIT:1> ) ) <EOL> if K . has ( factorial , FallingFactorial , RisingFactorial ) : <EOL> K = simplify ( K ) <EOL> if casoratian ( kernel + [ K ] , n , zero = False ) != <NUM_LIT:0> : <EOL> kernel . append ( K ) <EOL> kernel . sort ( key = default_sort_key ) <EOL> sk = list ( zip ( numbered_symbols ( '<STR_LIT:C>' ) , kernel ) ) <EOL> if sk : <EOL> for C , ker in sk : <EOL> result += C * ker <EOL> else : <EOL> return None <EOL> if hints . get ( '<STR_LIT>' , False ) : <EOL> symbols |= { s for s , k in sk } <EOL> return ( result , list ( symbols ) ) <EOL> else : <EOL> return result <EOL> def rsolve ( f , y , init = None ) : <EOL> """<STR_LIT>""" <EOL> if isinstance ( f , Equality ) : <EOL> f = f . lhs - f . rhs <EOL> n = y . args [ <NUM_LIT:0> ] <EOL> k = Wild ( '<STR_LIT:k>' , exclude = ( n , ) ) <EOL> f = f . expand ( ) . collect ( y . func ( Wild ( '<STR_LIT:m>' , integer = True ) ) ) <EOL> h_part = defaultdict ( lambda : S . Zero ) <EOL> i_part = S . Zero <EOL> for g in Add . make_args ( f ) : <EOL> coeff = S . One <EOL> kspec = None <EOL> for h in Mul . make_args ( g ) : <EOL> if h . is_Function : <EOL> if h . func == y . func : <EOL> result = h . args [ <NUM_LIT:0> ] . match ( n + k ) <EOL> if result is not None : <EOL> kspec = int ( result [ k ] ) <EOL> else : <EOL> raise ValueError ( <EOL> "<STR_LIT>" % ( y . func , n , h ) ) <EOL> else : <EOL> raise ValueError ( <EOL> "<STR_LIT>" % ( y . func , h . func ) ) <EOL> else : <EOL> coeff *= h <EOL> if kspec is not None : <EOL> h_part [ kspec ] += coeff <EOL> else : <EOL> i_part += coeff <EOL> for k , coeff in h_part . items ( ) : <EOL> h_part [ k ] = simplify ( coeff ) <EOL> common = S . One <EOL> for coeff in h_part . values ( ) : <EOL> if coeff . is_rational_function ( n ) : <EOL> if not coeff . is_polynomial ( n ) : <EOL> common = lcm ( common , coeff . as_numer_denom ( ) [ <NUM_LIT:1> ] , n ) <EOL> else : <EOL> raise ValueError ( <EOL> "<STR_LIT>" % coeff ) <EOL> i_numer , i_denom = i_part . as_numer_denom ( ) <EOL> if i_denom . is_polynomial ( n ) : <EOL> common = lcm ( common , i_denom , n ) <EOL> if common is not S . One : <EOL> for k , coeff in h_part . items ( ) : <EOL> numer , denom = coeff . as_numer_denom ( ) <EOL> h_part [ k ] = numer * quo ( common , denom , n ) <EOL> i_part = i_numer * quo ( common , i_denom , n ) <EOL> K_min = min ( h_part . keys ( ) ) <EOL> if K_min < <NUM_LIT:0> : <EOL> K = abs ( K_min ) <EOL> H_part = defaultdict ( lambda : S . Zero ) <EOL> i_part = i_part . subs ( n , n + K ) . expand ( ) <EOL> common = common . subs ( n , n + K ) . expand ( ) <EOL> for k , coeff in h_part . items ( ) : <EOL> H_part [ k + K ] = coeff . subs ( n , n + K ) . expand ( ) <EOL> else : <EOL> H_part = h_part <EOL> K_max = max ( H_part . keys ( ) ) <EOL> coeffs = [ H_part [ i ] for i in range ( K_max + <NUM_LIT:1> ) ] <EOL> result = rsolve_hyper ( coeffs , - i_part , n , symbols = True ) <EOL> if result is None : <EOL> return None <EOL> solution , symbols = result <EOL> if init == { } or init == [ ] : <EOL> init = None <EOL> if symbols and init is not None : <EOL> if type ( init ) is list : <EOL> init = { i : init [ i ] for i in range ( len ( init ) ) } <EOL> equations = [ ] <EOL> for k , v in init . items ( ) : <EOL> try : <EOL> i = int ( k ) <EOL> except TypeError : <EOL> if k . is_Function and k . func == y . func : <EOL> i = int ( k . args [ <NUM_LIT:0> ] ) <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" % k ) <EOL> try : <EOL> eq = solution . limit ( n , i ) - v <EOL> except NotImplementedError : <EOL> eq = solution . subs ( n , i ) - v <EOL> equations . append ( eq ) <EOL> result = solve ( equations , * symbols ) <EOL> if not result : <EOL> return None <EOL> else : <EOL> solution = solution . subs ( result ) <EOL> return solution </s>
