From d2815161b48ec427c4fc2e2cf64c7f38c357d9d3 Mon Sep 17 00:00:00 2001 From: mattpap Date: Sat, 28 Jul 2007 22:21:00 +0000 Subject: [PATCH] Made Bernoulli numbers and polynomials work with new core and moved to integer_sequences.py. Refactored and bug fixed log(). Bug fixed floor() and ceiling(). Changed string representation of rf() and ff(). Added Euler gamma, golden ratio and catalan symbolic numbers and refactored others. --- sympy/core/defined_functions.py | 60 ++++--- sympy/core/integer_sequences.py | 176 +++++++++++++++------ sympy/core/numbers.py | 81 ++++++++-- sympy/modules/__init__.py | 2 +- .../specfun/tests/test_orthogonal_polynomials.py | 2 +- 5 files changed, 241 insertions(+), 80 deletions(-) rewrite sympy/core/integer_sequences.py (71%) diff --git a/sympy/core/defined_functions.py b/sympy/core/defined_functions.py index 99808df..0419d26 100644 --- a/sympy/core/defined_functions.py +++ b/sympy/core/defined_functions.py @@ -196,7 +196,7 @@ class Log(DefinedFunction): """ Log() -> log """ is_comparable = True - nofargs = 1 + #nofargs = 1 def fdiff(self, argindex=1): if argindex == 1: @@ -209,8 +209,6 @@ class Log(DefinedFunction): return Exp() def _eval_apply(self, arg, base=None): - I = Basic.ImaginaryUnit() - if base is not None: base = Basic.sympify(base) @@ -221,19 +219,19 @@ class Log(DefinedFunction): if isinstance(arg, Basic.Number): if isinstance(arg, Basic.Zero): - return Basic.NegativeInfinity() + return S.NegativeInfinity elif isinstance(arg, Basic.One): - return Basic.Zero() + return S.Zero elif isinstance(arg, Basic.Infinity): - return Basic.Infinity() + return S.Infinity elif isinstance(arg, Basic.NegativeInfinity): - return Basic.Infinity() + return S.Infinity elif isinstance(arg, Basic.NaN): - return Basic.NaN() + return S.NaN elif arg.is_negative: - return Basic.Pi() * I + self(-arg) + return S.Pi * S.ImaginaryUnit + self(-arg) elif isinstance(arg, Basic.Exp1): - return Basic.One() + return S.One elif isinstance(arg, ApplyExp) and arg.args[0].is_real: return arg.args[0] elif isinstance(arg, Basic.Pow): @@ -243,25 +241,21 @@ class Log(DefinedFunction): elif isinstance(arg, Basic.Mul) and arg.is_real: return Basic.Add(*[self(a) for a in arg]) elif not isinstance(arg, Basic.Add): - x = Basic.Wild('x') - - coeff = arg.match(x*I) + coeff = arg.as_coefficient(S.ImaginaryUnit) if coeff is not None: - coeff = coeff[x] - if isinstance(coeff, Basic.Infinity): - return Basic.Infinity() + return S.Infinity elif isinstance(coeff, Basic.NegativeInfinity): - return Basic.Infinity() + return S.Infinity elif isinstance(coeff, Basic.Rational): if coeff.is_nonnegative: - return Basic.Pi() * I * Basic.Half() + self(coeff) + return S.Pi * S.ImaginaryUnit * S.Half + self(coeff) else: - return -Basic.Pi() * I * Basic.Half() + self(-coeff) + return -S.Pi * S.ImaginaryUnit * S.Half + self(-coeff) def as_base_exp(self): - return Exp(),Basic.Integer(-1) + return Exp(), Basic.Integer(-1) def _eval_apply_evalf(self, arg): arg = arg.evalf() @@ -719,7 +713,7 @@ class Floor(DefinedFunction): elif isinstance(arg, Basic.Rational): return Basic.Integer(arg.p // arg.q) elif isinstance(arg, Basic.Real): - return arg.floor() + return Basic.Integer(int(arg.floor())) elif isinstance(arg, Basic.NumberSymbol): return arg.approximation_interval(Basic.Integer)[0] elif arg.has(Basic.ImaginaryUnit()): @@ -812,7 +806,7 @@ class Ceiling(DefinedFunction): elif isinstance(arg, Basic.Rational): return Basic.Integer(arg.p // arg.q + 1) elif isinstance(arg, Basic.Real): - return arg.ceiling() + return Basic.Integer(int(arg.ceiling())) elif isinstance(arg, Basic.NumberSymbol): return arg.approximation_interval(Basic.Integer)[1] elif arg.has(Basic.ImaginaryUnit()): @@ -912,6 +906,16 @@ class RisingFactorial(DefinedFunction): else: return 1/reduce(lambda r, i: r*(x-i), xrange(1, abs(int(k))+1), 1) +class ApplyRisingFactorial(Apply): + + def tostr(self, level=0): + r = 'rf(%s)' % ', '.join([a.tostr() for a in self.args]) + + if self.precedence <= level: + return '(%s)' % (r) + else: + return r + class FallingFactorial(DefinedFunction): """Falling factorial (related to rising factorial) is a double valued function arising in concrete mathematics, hypergeometric functions @@ -977,6 +981,16 @@ class FallingFactorial(DefinedFunction): else: return 1/reduce(lambda r, i: r*(x+i), xrange(1, abs(int(k))+1), 1) +class ApplyFallingFactorial(Apply): + + def tostr(self, level=0): + r = 'ff(%s)' % ', '.join([a.tostr() for a in self.args]) + + if self.precedence <= level: + return '(%s)' % (r) + else: + return r + Basic.singleton['pochhammer'] = RisingFactorial Basic.singleton['rf'] = RisingFactorial @@ -986,6 +1000,8 @@ Basic.singleton['ff'] = FallingFactorial ############## REAL and IMAGINARY PARTS, CONJUNGATION, ARGUMENT ############### ############################################################################### +# TBD : re, im, arg + ############################################################################### ########################## TRIGONOMETRIC FUNCTIONS ############################ ############################################################################### diff --git a/sympy/core/integer_sequences.py b/sympy/core/integer_sequences.py dissimilarity index 71% index fbc1be5..fff0ec6 100644 --- a/sympy/core/integer_sequences.py +++ b/sympy/core/integer_sequences.py @@ -1,46 +1,130 @@ - -from basic import Basic, Atom -from function import DefinedFunction, Lambda, Function, Apply - -class IntegerSequence(DefinedFunction): - """ - http://en.wikipedia.org/wiki/Category:Integer_sequences - """ - - pass - -#class Factorial(IntegerSequence): - -# nofargs = 1 - -# def _eval_apply(self, arg): -# if isinstance(arg, Basic.Zero): return Basic.One() -# if isinstance(arg, Basic.Integer) and arg.is_positive: -# r = arg.p -# m = 1 -# while r: -# m *= r -# r -= 1 -# return Basic.Integer(m) - # return Gamma(arg) - -class Binomial(IntegerSequence): - """ - http://en.wikipedia.org/wiki/Binomial_coefficient - """ - - nofargs = 2 - - def _eval_apply(self, z, k): - from sympy.modules.specfun.factorials import Factorial - if isinstance(k, Basic.Integer): - assert not k.is_negative,`k` - l = [] - zk = z - k - for n in xrange(1, k.p+1): - l.append(zk+n) - return (Basic.Mul(*l)/Factorial()(k)).expand() - -#Basic.singleton['factorial'] = Factorial -Basic.singleton['binomial'] = Binomial - + +from basic import Basic, S, Memoizer +from numbers import Integer, Rational +from function import DefinedFunction, Lambda, Function, Apply + +class IntegerSequence(DefinedFunction): + """ + http://en.wikipedia.org/wiki/Category:Integer_sequences + """ + + pass + +############################################################################### +########################### FACTORIALS and BINOMIAL ########################### +############################################################################### + +class Factorial(IntegerSequence): + + nofargs = 1 + + def _eval_apply(self, arg): + arg = Basic.sympify(arg) + + if isinstance(arg, Basic.Number): + if isinstance(arg, Basic.Zero): + return S.One + elif isinstance(arg, Basic.Integer): + if arg.is_negative: + return S.Zero + else: + k, n = arg.p, 1 + + while k: + n *= k + k -= 1 + + return Integer(n) + +class Binomial(IntegerSequence): + + nofargs = 2 + + def _eval_apply(self, r, k): + r, k = map(Basic.sympify, (r, k)) + + if isinstance(k, Basic.Integer): + if k.is_negative: + return S.Zero + else: + rk, result = r - k, [] + + for i in xrange(1, k.p+1): + result.append(rk+i) + + numer = Basic.Mul(*result) + denom = Factorial()(k) + + return (numer/denom).expand() + +Basic.singleton['factorial'] = Factorial +Basic.singleton['binomial'] = Binomial + +############################################################################### +###################### BERNOULLI NUMBERS and POLYNOMIALS ###################### +############################################################################### + +def _bernoulli_sum(m, bound): + r, a = S.Zero, int(Binomial()(m+3, m-6)) + + for j in xrange(1, bound+1): + r += a * Bernoulli()(m-6*j) + + a *= reduce(lambda r, k: r*k, range(m-6*j-5, m-6*j+1), 1) + a //= reduce(lambda r, k: r*k, range(6*j+4, 6*j+10), 1) + + return r + +@Memoizer((int, long)) +def _b0mod6(m): + return Rational(m+3, 3) - _bernoulli_sum(m, m//6) + +@Memoizer((int, long)) +def _b2mod6(m): + return Rational(m+3, 3) - _bernoulli_sum(m, (m-2)//6) + +@Memoizer((int, long)) +def _b4mod6(m): + return -Rational(m+3, 6) - _bernoulli_sum(m, (m-4)//6) + +class Bernoulli(DefinedFunction): + + def _eval_apply(self, arg, sym=None): + arg = Basic.sympify(arg) + + if isinstance(arg, Basic.Number): + if isinstance(arg, Basic.Zero): + return S.One + elif isinstance(arg, Basic.One): + if sym is None: + return -S.Half + else: + return sym - S.Half + elif isinstance(arg, Basic.Integer): + if arg.is_nonnegative: + if sym is None: + m = int(arg) + + val_table = { + 0 : _b0mod6, + 2 : _b2mod6, + 4 : _b4mod6 + } + + try: + value = val_table[m % 6](m) + return value / Binomial()(m+3, m) + except KeyError: + return S.Zero + else: + n, result = int(arg), [] + + for k in xrange(n + 1): + result.append(Binomial()(n, k)*self(k)*sym**(n-k)) + + return Basic.Add(*result) + else: + raise ValueError("Bernoulli numbers are defined only" + " for nonnegative integer indices.") + +Basic.singleton['bernoulli'] = Bernoulli diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py index c730fb4..a1e6a47 100644 --- a/sympy/core/numbers.py +++ b/sympy/core/numbers.py @@ -60,7 +60,7 @@ class Number(Atom, RelMeths, ArithMeths): Floating point numbers are represented by the Real class. Integer numbers (of any size), together with rational numbers (again, there - is no limit on their size) are represented by the Rational class. + is no limit on their size) are represented by the Rational class. If you want to represent for example 1+sqrt(2), then you need to do: @@ -673,7 +673,7 @@ class Integer(Rational): def __rfloordiv__(self, other): return Integer(Integer(other).p // self.p) - + class Zero(Singleton, Integer): p = 0 @@ -757,7 +757,7 @@ class Infinity(Singleton, Rational): is_bounded = False is_finite = None is_odd = None - + def tostr(self, level=0): return 'oo' @@ -793,7 +793,7 @@ class NegativeInfinity(Singleton, Rational): is_positive = False is_bounded = False is_finite = False - + precedence = 40 # same as Add def tostr(self, level=0): @@ -924,20 +924,77 @@ class Pi(NumberSymbol): is_real = True is_positive = True - is_negative = False # XXX Forces is_negative/is_nonnegative + is_negative = False is_irrational = True + def _eval_evalf(self): + return Real(decimal_math.pi()) + def approximation_interval(self, number_cls): - if issubclass(number_cls,Integer): - return (Integer(3),Integer(4)) - elif issubclass(number_cls,Rational): - pass + if issubclass(number_cls, Integer): + return (Integer(3), Integer(4)) + elif issubclass(number_cls, Rational): + return (Rational(223,71), Rational(22,7)) def tostr(self, level=0): return 'Pi' +class GoldenRatio(NumberSymbol): + + is_real = True + is_positive = True + is_negative = False + is_irrational = True + def _eval_evalf(self): - return Real(decimal_math.pi()) + return Real(decimal_math.golden_ratio()) + + def approximation_interval(self, number_cls): + if issubclass(number_cls, Integer): + return (S.One, Rational(2)) + elif issubclass(number_cls, Rational): + pass + + def tostr(self, level=0): + return 'GoldenRatio' + +class EulerGamma(NumberSymbol): + + is_real = True + is_positive = True + is_negative = False + is_irrational = None + + def _eval_evalf(self): + return + + def approximation_interval(self, number_cls): + if issubclass(number_cls, Integer): + return (S.Zero, S.One) + elif issubclass(number_cls, Rational): + return (S.Half, Rational(3, 5)) + + def tostr(self, level=0): + return 'EulerGamma' + +class Catalan(NumberSymbol): + + is_real = True + is_positive = True + is_negative = False + is_irrational = None + + def _eval_evalf(self): + return + + def approximation_interval(self, number_cls): + if issubclass(number_cls, Integer): + return (S.Zero, S.One) + elif issubclass(number_cls, Rational): + return (Rational(9, 10), S.One) + + def tostr(self, level=0): + return 'Catalan' class ImaginaryUnit(Singleton, Atom, RelMeths, ArithMeths): @@ -990,3 +1047,7 @@ Basic.singleton['Pi'] = Pi Basic.singleton['I'] = ImaginaryUnit Basic.singleton['oo'] = Infinity Basic.singleton['nan'] = NaN + +Basic.singleton['GoldenRatio'] = GoldenRatio +Basic.singleton['EulerGamma'] = EulerGamma +Basic.singleton['Catalan'] = Catalan \ No newline at end of file diff --git a/sympy/modules/__init__.py b/sympy/modules/__init__.py index 12485e3..87694d7 100644 --- a/sympy/modules/__init__.py +++ b/sympy/modules/__init__.py @@ -10,5 +10,5 @@ from matrices import * from geometry import * from polynomials import * from utilities import * -from specfun import * +#from specfun import * from integrals import * diff --git a/sympy/modules/specfun/tests/test_orthogonal_polynomials.py b/sympy/modules/specfun/tests/test_orthogonal_polynomials.py index 06c715c..863b21e 100644 --- a/sympy/modules/specfun/tests/test_orthogonal_polynomials.py +++ b/sympy/modules/specfun/tests/test_orthogonal_polynomials.py @@ -1,6 +1,6 @@ from sympy import * from sympy.core.orthogonal_polynomials import * -from sympy.modules.specfun.orthogonal_polynomials import Chebyshev3 +from sympy.modules.specfun.orthogonal_polynomials import Chebyshev3, legendre, legendre_zero, chebyshev_zero x = Symbol('x') -- 2.11.4.GIT