Symbolic Chebyshev polynomials

[sagetrac-maldun]

Apply:

• attachment: trac_9706_chebyshev_collection.patch
• attachment: 9706_review.patch
• attachment: trac_9706_new_base_classes.patch
• attachment: trac_9706-review-ts.patch

Depends on #864
Depends on #9640
Depends on #10018
Depends on #11868
Depends on #15422

CC: @fredrik-johansson @sagetrac-fstan @kcrisman

Component: symbolics

Keywords: orthogonal polynomials, symbolics

Author: Stefan Reiterer

Branch/Commit: u/vbraun/chebyshev @ bb227b6

Reviewer: Burcin Erocal, Travis Scrimshaw, Stefan Reiterer, Jeroen Demeyer

Issue created by migration from https://trac.sagemath.org/ticket/9706

[sagetrac-maldun]

Attachment: orthogonal_polys.py.gz

A new version of the orthogonal_polys.py file.

[sagetrac-maldun]

Newer version, with legendre_P, and faster evaluation of symbolic expressions

[sagetrac-maldun]

Attachment: orthogonal_polys.2.py.gz

Attachment: orthogonal_polys.3.py.gz

Version from 10. August 2010

[sagetrac-maldun]

Latest version. It holds classes of all polys (but not all completed yet)

[sagetrac-maldun]

comment:1

Attachment: orthogonal_polys.4.py.gz

All Polys now have their own class.
Much faster evaluation is added.
Numerical evaluation is provided.
Except for legendre_Q, gen_legendre_P, and gen_legendre_Q these aren't ready yet

[sagetrac-maldun]

comment:2

Replying to @sagetrac-maldun:

All Polys now have their own class.
Much faster evaluation is added.
Numerical evaluation is provided.
Except for legendre_Q, gen_legendre_P, and gen_legendre_Q these aren't ready yet

orthogonal_polys4.py hold all changes but is not a patch yet, because it holds old code fragments,
which I have to clean up...

[sagetrac-maldun]

comment:4

I added in the latest patch (and orthogonal_polys.4.py contains these changes also) a new symbolic evaluation method for the orthogonal polynomials: Instead of call Maxima or use of the recursion, the polynomial is evaluated just using explicit formulas from Abramowitz and Stegun. This is an O(n) algorithm of course.

a little comparison on my machine:
old version:

sage: time chebyshev_T(10,x);
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.04 s
sage: time chebyshev_T(100,x);
CPU times: user 0.13 s, sys: 0.01 s, total: 0.14 s
Wall time: 0.23 s
sage: time chebyshev_T(1000,x);
CPU times: user 5.01 s, sys: 0.01 s, total: 5.02 s
Wall time: 6.98 s
sage time chebyshev_T(5000,x);
??? (I got no output her after 2min)

sage: time gegenbauer(10,5,x);
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.05 s
sage: time gegenbauer(100,5,x);
CPU times: user 0.19 s, sys: 0.00 s, total: 0.19 s
Wall time: 0.29 s
sage: time gegenbauer(1000,5,x);
CPU times: user 5.46 s, sys: 0.02 s, total: 5.48 s
Wall time: 7.79 s

New Version
sage: time chebyshev_T(10,x);
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01 s
sage: time chebyshev_T(100,x);
CPU times: user 0.06 s, sys: 0.00 s, total: 0.06 s
Wall time: 0.08 s
sage: time chebyshev_T(1000,x);
CPU times: user 1.22 s, sys: 0.00 s, total: 1.22 s
Wall time: 1.22 s
sage: time chebyshev_T(5000,x);
CPU times: user 27.17 s, sys: 0.15 s, total: 27.32 s
Wall time: 27.46 s

sage: time gegenbauer(10,5,x);
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01 s
sage: time gegenbauer(100,5,x);
CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
Wall time: 0.04 s
sage: time gegenbauer(1000,5,x);
CPU times: user 1.08 s, sys: 0.01 s, total: 1.09 s
Wall time: 1.11 s

A little bit faster :) I also don't need to spawn an instance of maxima which makes the initialisation faster.

And now also wider symbolic evaluation is possible:

old version:
sage: var('a')
a
sage: gegenbauer(3,a,x)
...
NameError: name 'a' is not defined

new version:
sage: var('a')
a
sage: gegenbauer(3,a,x)
4/3x^3gamma(a + 3) - 2xgamma(a + 2)

The code needs now some cleanup, especially the documentations.
The complete versions for legendre_Q, gen_legendre_P, and gen_legendre_Q will not be finished
soon since the mpmath functions, don't seem to work correctly...
I only provide a call function for maxima for them now.

[fredrik-johansson]

comment:5

The complete versions for legendre_Q, gen_legendre_P, and gen_legendre_Q will not be finished soon since the mpmath functions, don't seem to work correctly...

Care to elaborate?

[sagetrac-maldun]

Attachment: orthogonal_polys.5.py.gz

Latest version from 12. August 2010 (with bugfix in legendre_P)

[sagetrac-maldun]

comment:6

Killed bug in legendre_P

[sagetrac-maldun]

comment:7

Replying to @fredrik-johansson:

The complete versions for legendre_Q, gen_legendre_P, and gen_legendre_Q will not be finished soon since the mpmath functions, don't seem to work correctly...

Care to elaborate?

Sorry for the late answer, I was on holidays.

In mpmath I have probs with the legenp and legenq functions. For some inputs I get this error:

sage: mpmath.call(mpmath.legenp,5,1,2)
---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)

/home/maldun/prog/sage/ortho/<ipython console> in <module>()

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/sage/libs/mpmath/utils.so in sage.libs.mpmath.utils.call (sage/libs/mpmath/utils.c:5021)()

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/functions/hypergeometric.pyc in legenp(ctx, n, m, z, type, **kwargs)
   1481             T = [1+z, 1-z], [g, -g], [], [1-m], [-n, n+1], [1-m], 0.5*(1-z)
   1482             return (T,)
-> 1483         return ctx.hypercomb(h, [n,m], **kwargs)
   1484     if type == 3:
   1485         def h(n,m):

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/functions/hypergeometric.pyc in hypercomb(ctx, function, params, discard_known_zeros, **kwargs)
    125                     [ctx.gamma(a) for a in alpha_s] + \
    126                     [ctx.rgamma(b) for b in beta_s] + \
--> 127                     [ctx.power(w,c) for (w,c) in zip(w_s,c_s)])
    128                 if verbose:
    129                     print "    Value:", v

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/ctx_base.pyc in power(ctx, x, y)
    417             3.16470269330255923143453723949e+12978188
    418         """
--> 419         return ctx.convert(x) ** ctx.convert(y)
    420 
    421     def _zeta_int(ctx, n):

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/sage/libs/mpmath/ext_main.so in sage.libs.mpmath.ext_main.mpnumber.__pow__ (sage/libs/mpmath/ext_main.c:13946)()

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/sage/libs/mpmath/ext_main.so in sage.libs.mpmath.ext_main.binop (sage/libs/mpmath/ext_main.c:4588)()

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in mpf_pow(s, t, prec, rnd)
    340     # General formula: s**t = exp(t*log(s))
    341     # TODO: handle rnd direction of the logarithm carefully
--> 342     c = mpf_log(s, prec+10, rnd)
    343     return mpf_exp(mpf_mul(t, c), prec, rnd)
    344 

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in mpf_log(x, prec, rnd)
    725     # optimal between 1000 and 100,000 digits.
    726     if wp <= LOG_TAYLOR_PREC:
--> 727         m = log_taylor_cached(lshift(man, wp-bc), wp)
    728         if mag:
    729             m += mag*ln2_fixed(wp)

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in log_taylor_cached(x, prec)
    643     else:
    644         a = n << (cached_prec - LOG_TAYLOR_SHIFT)
--> 645         log_a = log_taylor(a, cached_prec, 8)
    646         log_taylor_cache[n, cached_prec] = (a, log_a)
    647     a >>= dprec

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in log_taylor(x, prec, r)
    607     """
    608     for i in xrange(r):
--> 609         x = isqrt_fast(x<<prec)
    610     one = MPZ_ONE << prec
    611     v = ((x-one)<<prec)//(x+one)

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libintmath.pyc in isqrt_fast_python(x)
    240                     y = (y + x//y) >> 1
    241         return y
--> 242     bc = bitcount(x)
    243     guard_bits = 10
    244     x <<= 2*guard_bits

/home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libintmath.pyc in python_bitcount(n)
     78     if bc != 300:
     79         return bc
---> 80     bc = int(math.log(n, 2)) - 4
     81     return bc + bctable[n>>bc]
     82 

OverflowError: cannot convert float infinity to integer

[fredrik-johansson]

comment:8

That looks strange. I get:

sage: import sage.libs.mpmath.all as mpmath
sage: mpmath.call(mpmath.legenp, 5,1,2)
-2.96434298694874e-22 - 912.574269237852*I
sage: mpmath.call(mpmath.legenp, 5,1,2, prec=100)
-2.1062923756778274648015607872e-36 - 912.57426923785222402727329118*I

[sagetrac-maldun]

comment:9

Replying to @fredrik-johansson:

That looks strange. I get:

sage: import sage.libs.mpmath.all as mpmath
sage: mpmath.call(mpmath.legenp, 5,1,2)
-2.96434298694874e-22 - 912.574269237852*I
sage: mpmath.call(mpmath.legenp, 5,1,2, prec=100)
-2.1062923756778274648015607872e-36 - 912.57426923785222402727329118*I

Hm strange. Today I install the new Sage version, perhaps it will then work again

[sagetrac-maldun]

comment:10

Replying to @fredrik-johansson:

That looks strange. I get:

sage: import sage.libs.mpmath.all as mpmath
sage: mpmath.call(mpmath.legenp, 5,1,2)
-2.96434298694874e-22 - 912.574269237852*I
sage: mpmath.call(mpmath.legenp, 5,1,2, prec=100)
-2.1062923756778274648015607872e-36 - 912.57426923785222402727329118*I

It was the old version!a Thanx for pointing that out, I will continue soon =)

[sagetrac-maldun]

Version from 19. August 2010

[sagetrac-maldun]

comment:11

Attachment: orthogonal_polys.6.py.gz

So now a "beta" is ready with full support of all classes.

Only the Legendre functions are still using Maxima.

some advances for the future:

-Zernike polys (this should be done in the next time, since explicit formulas are available)
-support for numpy_eval. (But this will be done, when the scipy package is updated to 0.8, else it has no sense, because the current version of scipy does not support ortho polys well, but the newer can handle them)

Now I need some people for testing this out =)

[sagetrac-maldun]

comment:12

And there was an interisting bug:

the import of mpmath at the beginning of the file caused the whole trouble I had with the numeric evaluation of the legendre functions....

I think I should report this..

[sagetrac-maldun]

Added numpy support, eliminated some bugs (19.08.2010)

[sagetrac-maldun]

comment:14

Attachment: orthogonal_polys.7.py.gz

-support for numpy_eval. (But this will be done, when the scipy package is updated to 0.8, else it has no sense, because the current version of scipy does not support ortho polys well, but the newer can handle them)

I decided to give at least some numpy support for compability reasons.
But this is a bad hack...when scipy 0.8 comes I use scipy itself, I change this to a better version :)

[sagetrac-maldun]

comment:17

Some of the old doctests fail.
But it is not my fault, it seem's that it is a bug in the SymbolicFunction class.

see: #9769

[jdemeyer]

comment:117

Something like that looks right indeed. maldun: what do you think?

Perhaps the only code so far that could be truly generic is the __call__ method.

[sagetrac-maldun]

comment:118

after some thinking, I guess you are right. A General OrthogonalPolynomial is sophisticated, but it needs too much tweaking, and too much exception cathing, which makes the code unsafe.

@ __call__ : I'm not even sure about that, since we check for negative integers, but some ortho polys get only an analytical expression with non negative integers, and no algebraic meaning.

I propose the following

• OrthogonalPolynomial: Naming Conventions and __init__ method
• ChebyshevPolynomial: base Class for Cheby_t and Cheby_u (Current Orthogonal Polynomials)
• LegendrePolynomials
  .... etc.

[sagetrac-maldun]

Attachment: trac_9706_new_base_classes.patch.gz

proposed new structure

[sagetrac-maldun]

comment:119

I added an intermediate class between OrthogonalPolynomials and the Chebyshev polynomals namely ChebyshevPolynomial.

Any comments on that?

[tscrim]

comment:120

Sorry for lagging behind on reviewing this, but I think the class structure is good. I'll write a small review patch tomorrow afternoon (in the US Pacific timezone) on a few tweaks (and pet peeves of mine).

[jdemeyer]

comment:121

I am currently working on finishing the Sage 5.13 release, so I don't feel like reviewing this now.

[tscrim]

comment:122

Okay, I've made the tweaks I wanted to. So it's a positive review from me, and someone will just need to review my changes.

[sagetrac-maldun]

comment:123

Okay thanks. I set it on positive review then.
And good idea to replace the old % format string operator, since it will be deprecated in Python 3 (although I will miss it ...)

[nbruin]

comment:124

Replying to @sagetrac-maldun:

And good idea to replace the old % format string operator, since it will be deprecated in Python 3

This has been a rumour that didn't pan out. I don't think a time line has been set towards actual deprecation of %. The possibility of deprecation is mentioned in
http://www.python.org/dev/peps/pep-3101/
but it's careful to argue that deprecation is not required at all. At some point it was planned to deprecate % in 3.0 and remove it in 3.1 but that hasn't happened. I suspect it's so ingrained by now that deprecation and removal is impractical.

[jdemeyer]

comment:125

Travis: are you sure you applied the dependency #15422? I am suspicious because you added

-            (2^7 + O(2^8))*t^5 + (O(2^8))*t^4 + (2^6 + O(2^8))*t^3 + (O(2^8))*t^2 + (1 + 2^6 + O(2^8))*t + (O(2^8))
+            (2^7 + O(2^8))*t^5 + (2^6 + O(2^8))*t^3 + (1 + 2^6 + O(2^8))*t

which seems to be a doctest failure.

[jdemeyer]

comment:126

I also don't agree with all changes of

if condition:
    return A
else:
    return B

to

if condition:
    return A
return B

but I guess that's one of your pet peeves.

I personally prefer

if condition:
    return A
else:
    return B

if there are clearly two cases and the code might as well have been written as

if not condition:
    return B
else:
    return A

So I personally would keep the if/else structure for the if n % 2 == 0 test. And for the n >= 0 test, I would say that

if n < 0:
    return B
return A

would be a lot better that what you did:

if n >= 0:
    return A
return B

The first feels much better to me because you put the normal case outside the if block and the special cases if n == 0 and if n < 0 would be inside ifs.

(Of course these are all personal preferences, I'm not asking you to change this, maybe just think about it.)

[tscrim]

Attachment: trac_9706-review-ts.patch.gz

[tscrim]

comment:128

I had missed the #15422 dependency. Fixed.

As for the if/else, I've seen (sometimes long) else blocks that constitute the normal cases and I try to be uniform about it. Plus with the extra indents, I sometimes have difficulties knowing what the correct indent level should be (although not in this case). I do agree that the cases should be swapped -- I had only removed the else:. I added back in the else blocks for the if n % 2 == 0.

[sagetrac-maldun]

comment:130

First stage done. Next step: Legendre Polys.

I will keep up the new design + incorporate more doc tests next time. I guess it is important to look at the algebraic aspects too concerning doc tests.

Thank all reviewers for their hard work and the good input!

[vbraun]

New commits:

bb227b6	#9706: review patch.
51fcd6b	trac 9706: Propose new class structure
4677a15	Symbolic Chebyshev polynomials: reviewer patch
75d1e43	trac 9706: Collective patch. Bugfixes, extensions, optimizations, documentation, doctests for chebyshev_T, chebyshev_U and base class for ortho polys

[vbraun]

Commit: bb227b6

[vbraun]

Branch: u/vbraun/chebyshev