Sage: Ticket #9706: New Version of orthogonal Polynomials (Part I: Base Class and Chebyshev Polynomials)
https://trac.sagemath.org/ticket/9706
<p>
The current implementation of orthogonal polynomials is just a wrapper around maxima. (see <a class="ext-link" href="http://wiki.sagemath.org/symbolics/"><span class="icon"></span>http://wiki.sagemath.org/symbolics/</a>)
This update holds the following changes:
</p>
<ul><li>using of the pynac class for symbolic functions.
</li><li>faster evaluation in general
</li><li>evaluation of special values
</li><li>mpmath for numeric evaluation
</li></ul><p>
<strong>Remarks:</strong>
</p>
<ul><li>The current patch needs scipy-0.8. One has to install it before testing (see <a class="closed ticket" href="https://trac.sagemath.org/ticket/9808" title="#9808: task: Upgrade numpy to 1.5.0 and scipy to 0.8 (closed: fixed)">#9808</a> for spkg's and installation instructions)
</li><li>Some of the old doctests in the old file don't work any more, due to coercion problems with pynac (see <a class="closed ticket" href="https://trac.sagemath.org/ticket/9769" title="#9769: defect: symbolic function do not work with numpy.int64 arguments (closed: duplicate)">#9769</a>)
</li><li>Some doctests in Sage change, due to the fact that new <a class="missing wiki">BuiltIn?</a> functions are added. symbolic.random_test.py had output changes since the random expression creation changed of course. The tests in pynac.pyx also changed, but this has a strange behavior (see below).
</li></ul><p>
Apply:
</p>
<ul><li><a class="attachment" href="https://trac.sagemath.org/attachment/ticket/9706/trac_9706_chebyshev.patch" title="Attachment 'trac_9706_chebyshev.patch' in Ticket #9706">trac_9706_chebyshev.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/9706/trac_9706_chebyshev.patch" title="Download"></a>
</li><li><a class="attachment" href="https://trac.sagemath.org/attachment/ticket/9706/trac_9706-review-ts.patch" title="Attachment 'trac_9706-review-ts.patch' in Ticket #9706">trac_9706-review-ts.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/9706/trac_9706-review-ts.patch" title="Download"></a>
</li><li><a class="attachment" href="https://trac.sagemath.org/attachment/ticket/9706/trac_9706_new_clenshaw.patch" title="Attachment 'trac_9706_new_clenshaw.patch' in Ticket #9706">trac_9706_new_clenshaw.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/9706/trac_9706_new_clenshaw.patch" title="Download"></a>
</li><li><a class="attachment" href="https://trac.sagemath.org/attachment/ticket/9706/trac_9706_chebyshev_patch_new.patch" title="Attachment 'trac_9706_chebyshev_patch_new.patch' in Ticket #9706">trac_9706_chebyshev_patch_new.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/9706/trac_9706_chebyshev_patch_new.patch" title="Download"></a>
</li></ul>en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/9706
Trac 1.2Stefan ReitererMon, 09 Aug 2010 00:13:02 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.py</em>
</li>
</ul>
<p>
A new version of the orthogonal_polys.py file.
</p>
TicketStefan ReitererMon, 09 Aug 2010 00:13:36 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.2.py</em>
</li>
</ul>
<p>
Newer version, with legendre_P, and faster evaluation of symbolic expressions
</p>
TicketStefan ReitererTue, 10 Aug 2010 16:57:34 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.3.py</em>
</li>
</ul>
<p>
Version from 10. August 2010
</p>
TicketStefan ReitererWed, 11 Aug 2010 00:07:02 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.4.py</em>
</li>
</ul>
<p>
Latest version. It holds classes of all polys (but not all completed yet)
</p>
TicketStefan ReitererWed, 11 Aug 2010 00:09:36 GMT
https://trac.sagemath.org/ticket/9706#comment:1
https://trac.sagemath.org/ticket/9706#comment:1
<p>
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
</p>
TicketStefan ReitererWed, 11 Aug 2010 00:11:08 GMT
https://trac.sagemath.org/ticket/9706#comment:2
https://trac.sagemath.org/ticket/9706#comment:2
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:1" title="Comment 1">maldun</a>:
</p>
<blockquote class="citation">
<p>
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
</p>
</blockquote>
<p>
orthogonal_polys4.py hold all changes but is not a patch yet, because it holds old code fragments,
which I have to clean up...
</p>
TicketFredrik JohanssonWed, 11 Aug 2010 13:02:19 GMTcc set
https://trac.sagemath.org/ticket/9706#comment:3
https://trac.sagemath.org/ticket/9706#comment:3
<ul>
<li><strong>cc</strong>
<em>Fredrik Johansson</em> added
</li>
</ul>
TicketStefan ReitererThu, 12 Aug 2010 10:25:41 GMT
https://trac.sagemath.org/ticket/9706#comment:4
https://trac.sagemath.org/ticket/9706#comment:4
<p>
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.
</p>
<p>
a little comparison on my machine:
old version:
</p>
<p>
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)
</p>
<p>
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
</p>
<p>
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
</p>
<p>
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
</p>
<p>
A little bit faster :) I also don't need to spawn an instance of maxima which makes the initialisation faster.
</p>
<p>
And now also wider symbolic evaluation is possible:
</p>
<p>
old version:
sage: var('a')
a
sage: gegenbauer(3,a,x)
...
<a class="missing wiki">NameError?</a>: name 'a' is not defined
</p>
<p>
new version:
sage: var('a')
a
sage: gegenbauer(3,a,x)
4/3*x<sup>3*gamma(a + 3) - 2*x*gamma(a + 2)
</sup></p>
<p>
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.
</p>
TicketFredrik JohanssonThu, 12 Aug 2010 11:12:10 GMT
https://trac.sagemath.org/ticket/9706#comment:5
https://trac.sagemath.org/ticket/9706#comment:5
<blockquote class="citation">
<p>
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...
</p>
</blockquote>
<p>
Care to elaborate?
</p>
TicketStefan ReitererThu, 12 Aug 2010 11:16:04 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.5.py</em>
</li>
</ul>
<p>
Latest version from 12. August 2010 (with bugfix in legendre_P)
</p>
TicketStefan ReitererThu, 12 Aug 2010 11:16:40 GMT
https://trac.sagemath.org/ticket/9706#comment:6
https://trac.sagemath.org/ticket/9706#comment:6
<p>
Killed bug in legendre_P
</p>
TicketStefan ReitererMon, 16 Aug 2010 11:51:22 GMT
https://trac.sagemath.org/ticket/9706#comment:7
https://trac.sagemath.org/ticket/9706#comment:7
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:5" title="Comment 5">fredrik.johansson</a>:
</p>
<blockquote class="citation">
<blockquote class="citation">
<p>
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...
</p>
</blockquote>
<p>
Care to elaborate?
</p>
</blockquote>
<p>
Sorry for the late answer, I was on holidays.
</p>
<p>
In mpmath I have probs with the legenp and legenq functions. For some inputs I get this error:
</p>
<pre class="wiki">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
</pre>
TicketFredrik JohanssonMon, 16 Aug 2010 15:55:46 GMT
https://trac.sagemath.org/ticket/9706#comment:8
https://trac.sagemath.org/ticket/9706#comment:8
<p>
That looks strange. I get:
</p>
<pre class="wiki">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
</pre>
TicketStefan ReitererTue, 17 Aug 2010 08:54:50 GMT
https://trac.sagemath.org/ticket/9706#comment:9
https://trac.sagemath.org/ticket/9706#comment:9
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:8" title="Comment 8">fredrik.johansson</a>:
</p>
<blockquote class="citation">
<p>
That looks strange. I get:
</p>
<pre class="wiki">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
</pre></blockquote>
<p>
Hm strange. Today I install the new Sage version, perhaps it will then work again
</p>
TicketStefan ReitererWed, 18 Aug 2010 16:31:58 GMT
https://trac.sagemath.org/ticket/9706#comment:10
https://trac.sagemath.org/ticket/9706#comment:10
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:8" title="Comment 8">fredrik.johansson</a>:
</p>
<blockquote class="citation">
<p>
That looks strange. I get:
</p>
<pre class="wiki">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
</pre></blockquote>
<p>
It was the old version!a Thanx for pointing that out, I will continue soon =)
</p>
TicketStefan ReitererThu, 19 Aug 2010 16:24:16 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.6.py</em>
</li>
</ul>
<p>
Version from 19. August 2010
</p>
TicketStefan ReitererThu, 19 Aug 2010 16:30:43 GMT
https://trac.sagemath.org/ticket/9706#comment:11
https://trac.sagemath.org/ticket/9706#comment:11
<p>
So now a "beta" is ready with full support of all classes.
</p>
<p>
Only the Legendre functions are still using Maxima.
</p>
<p>
some advances for the future:
</p>
<p>
-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)
</p>
<p>
Now I need some people for testing this out =)
</p>
TicketStefan ReitererThu, 19 Aug 2010 16:34:50 GMT
https://trac.sagemath.org/ticket/9706#comment:12
https://trac.sagemath.org/ticket/9706#comment:12
<p>
And there was an interisting bug:
</p>
<p>
the import of mpmath at the beginning of the file caused the whole trouble I had with the numeric evaluation of the legendre functions....
</p>
<p>
I think I should report this..
</p>
TicketStefan ReitererThu, 19 Aug 2010 16:35:08 GMTtype changed
https://trac.sagemath.org/ticket/9706#comment:13
https://trac.sagemath.org/ticket/9706#comment:13
<ul>
<li><strong>type</strong>
changed from <em>defect</em> to <em>enhancement</em>
</li>
</ul>
TicketStefan ReitererThu, 19 Aug 2010 19:54:21 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.7.py</em>
</li>
</ul>
<p>
Added numpy support, eliminated some bugs (19.08.2010)
</p>
TicketStefan ReitererThu, 19 Aug 2010 20:09:59 GMT
https://trac.sagemath.org/ticket/9706#comment:14
https://trac.sagemath.org/ticket/9706#comment:14
<blockquote class="citation">
<p>
-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)
</p>
</blockquote>
<p>
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 :)
</p>
TicketStefan ReitererThu, 19 Aug 2010 23:56:01 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:15
https://trac.sagemath.org/ticket/9706#comment:15
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
</ul>
TicketStefan ReitererThu, 19 Aug 2010 23:57:41 GMTmilestone set
https://trac.sagemath.org/ticket/9706#comment:16
https://trac.sagemath.org/ticket/9706#comment:16
<ul>
<li><strong>milestone</strong>
set to <em>sage-5.0</em>
</li>
</ul>
TicketStefan ReitererFri, 20 Aug 2010 00:08:44 GMT
https://trac.sagemath.org/ticket/9706#comment:17
https://trac.sagemath.org/ticket/9706#comment:17
<p>
Some of the old doctests fail.
But it is not my fault, it seem's that it is a bug in the <a class="missing wiki">SymbolicFunction?</a> class.
</p>
<p>
see: <a class="ext-link" href="http://trac.sagemath.org/sage_trac/ticket/9769"><span class="icon"></span>http://trac.sagemath.org/sage_trac/ticket/9769</a>
</p>
TicketStefan ReitererFri, 20 Aug 2010 00:17:36 GMTmilestone changed
https://trac.sagemath.org/ticket/9706#comment:18
https://trac.sagemath.org/ticket/9706#comment:18
<ul>
<li><strong>milestone</strong>
changed from <em>sage-5.0</em> to <em>sage-4.5.3</em>
</li>
</ul>
TicketStefan ReitererTue, 24 Aug 2010 21:31:01 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.8.py</em>
</li>
</ul>
<p>
Latest version with some code cleanup (no program changes)
</p>
TicketStefan ReitererThu, 26 Aug 2010 14:50:26 GMTowner changed
https://trac.sagemath.org/ticket/9706#comment:19
https://trac.sagemath.org/ticket/9706#comment:19
<ul>
<li><strong>owner</strong>
changed from <em>Burcin Erocal</em> to <em>burcin, maldun</em>
</li>
</ul>
TicketBurcin ErocalSat, 28 Aug 2010 12:03:30 GMT
https://trac.sagemath.org/ticket/9706#comment:20
https://trac.sagemath.org/ticket/9706#comment:20
<p>
Hi Stefan,
</p>
<p>
can you post a patch corresponding to <a class="attachment" href="https://trac.sagemath.org/attachment/ticket/9706/orthogonal_polys.8.py" title="Attachment 'orthogonal_polys.8.py' in Ticket #9706">attachment:orthogonal_polys.8.py</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/9706/orthogonal_polys.8.py" title="Download"></a> for review?
</p>
<p>
Thanks,<br />
Burcin
</p>
TicketStefan ReitererSat, 28 Aug 2010 16:49:03 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_ortho_polys.patch</em>
</li>
</ul>
<p>
Patch for latest version with some code cleanup (no program changes)
</p>
TicketStefan ReitererSat, 28 Aug 2010 16:50:33 GMT
https://trac.sagemath.org/ticket/9706#comment:21
https://trac.sagemath.org/ticket/9706#comment:21
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:20" title="Comment 20">burcin</a>:
</p>
<blockquote class="citation">
<p>
Hi Stefan,
</p>
<p>
can you post a patch corresponding to <a class="attachment" href="https://trac.sagemath.org/attachment/ticket/9706/orthogonal_polys.8.py" title="Attachment 'orthogonal_polys.8.py' in Ticket #9706">attachment:orthogonal_polys.8.py</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/9706/orthogonal_polys.8.py" title="Download"></a> for review?
</p>
<p>
Thanks,<br />
Burcin
</p>
</blockquote>
<p>
Done!
</p>
TicketFredrik JohanssonFri, 03 Sep 2010 13:07:05 GMT
https://trac.sagemath.org/ticket/9706#comment:22
https://trac.sagemath.org/ticket/9706#comment:22
<p>
Why is mpmath's precision used by default? Shouldn't the default be RR / CC precision? Actually, does _evalf_ ever get called without this information?
</p>
<p>
Some complex tests would be nice.
</p>
TicketStefan ReitererFri, 03 Sep 2010 19:59:28 GMT
https://trac.sagemath.org/ticket/9706#comment:23
https://trac.sagemath.org/ticket/9706#comment:23
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:22" title="Comment 22">fredrik.johansson</a>:
</p>
<blockquote class="citation">
<p>
Why is mpmath's precision used by default? Shouldn't the default be RR / CC precision? Actually, does _evalf_ ever get called without this information?
</p>
<p>
Some complex tests would be nice.
</p>
</blockquote>
<p>
This is a good point, and it shouldn't be a problem to change that.
But I don't think it's a big deal, because the function takes the "parents" precision, which means, if my input is RR it evals it with RR's precision.
</p>
<p>
Of course can you call _evalf_ just with (), and then the default value is used.
</p>
<p>
I just sticked to the old's version tests, and expanded it. Of course it's possible to expand the tests. I hope I will find some time for it soon, since I have some other more urgent things todo also.
</p>
TicketStefan ReitererFri, 03 Sep 2010 20:05:52 GMT
https://trac.sagemath.org/ticket/9706#comment:24
https://trac.sagemath.org/ticket/9706#comment:24
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:23" title="Comment 23">maldun</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:22" title="Comment 22">fredrik.johansson</a>:
</p>
<blockquote class="citation">
<p>
Why is mpmath's precision used by default? Shouldn't the default be RR / CC precision? Actually, does _evalf_ ever get called without this information?
</p>
<p>
Some complex tests would be nice.
</p>
</blockquote>
<p>
This is a good point, and it shouldn't be a problem to change that.
But I don't think it's a big deal, because the function takes the "parents" precision, which means, if my input is RR it evals it with RR's precision.
</p>
<p>
Of course can you call _evalf_ just with (), and then the default value is used.
</p>
</blockquote>
<p>
Ok sorry, wrong explination: when your input are exact data types like ZZ ore QQ then the parent has no precision, then you need a default value
</p>
TicketStefan ReitererSun, 05 Sep 2010 10:27:48 GMT
https://trac.sagemath.org/ticket/9706#comment:25
https://trac.sagemath.org/ticket/9706#comment:25
<p>
Since it seems that numpy-1.4.1, and scipy 0.8 should work now (see <a class="closed ticket" href="https://trac.sagemath.org/ticket/9808" title="#9808: task: Upgrade numpy to 1.5.0 and scipy to 0.8 (closed: fixed)">#9808</a>) I programmed a version which uses scipy itself to evaluate the orthogonal polys for numpy arrays.
When the new versions of numpy/scipy become merged into sage I will provide a patch for these.
</p>
<p>
Another thing I have to mention are these 2 failde doctests:
</p>
<ul><li>sage -t -long "devel/sage/sage/symbolic/random_tests.py"
</li><li>sage -t -long "devel/sage/sage/symbolic/pynac.pyx"
</li></ul><pre class="wiki">sage -t -long "devel/sage/sage/symbolic/random_tests.py"
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/random_tests.py", line 17:
sage: [f for (one,f,arity) in _mk_full_functions()]
Expected:
[Ei, abs, arccos, arccosh, arccot, arccoth, arccsc, arccsch,
arcsec, arcsech, arcsin, arcsinh, arctan, arctan2, arctanh,
binomial, ceil, conjugate, cos, cosh, cot, coth, csc, csch,
dickman_rho, dilog, dirac_delta, elliptic_e, elliptic_ec,
elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi, erf, exp,
factorial, floor, heaviside, imag_part, integrate,
kronecker_delta, log, polylog, real_part, sec, sech, sgn, sin,
sinh, tan, tanh, unit_step, zeta, zetaderiv]
Got:
[Ei, abs, arccos, arccosh, arccot, arccoth, arccsc, arccsch, arcsec, arcsech, arcsin, arcsinh, arctan, arctan2, arctanh, binomial, ceil, chebyshev_T, chebyshev_U, conjugate, cos, cosh, cot, coth, csc, csch, dickman_rho, dilog, dirac_delta, elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi, erf, exp, factorial, floor, gegenbauer, gen_laguerre, gen_legendre_P, gen_legendre_Q, heaviside, hermite, imag_part, integrate, jacobi_P, kronecker_delta, laguerre, legendre_P, legendre_Q, log, polylog, real_part, sec, sech, sgn, sin, sinh, tan, tanh, unit_step, zeta, zetaderiv]
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/random_tests.py", line 237:
sage: random_expr(50, nvars=3, coeff_generator=CDF.random_element)
Expected:
(euler_gamma - v3^(-e) + (v2 - factorial(-e/v2))^(((2.85879036573 - 1.18163393202*I)*v2 + (2.85879036573 - 1.18163393202*I)*v3)*pi - 0.247786879678 + 0.931826724898*I)*arccsc((0.891138386848 - 0.0936820840629*I)/v1) + (-0.553423153995 + 0.5481180572*I)*v3 + 0.149683576515 - 0.155746451854*I)*v1 + arccsch(pi + e)*elliptic_f(khinchin*v2, 1.4656989704 + 0.863754357069*I)
Got:
-v1*e^((0.0666829501658 + 0.206976992303*I)/(v3 + e))/v3 + hermite(-(v3^(-0.48519994364 - 0.485764091302*I) - log((1.21734510331 - 1.22580558833*I)*pi*v1 + zeta((0.781366128261 + 0.957400336147*I)*v1*e + (-1.8919687109 + 0.753422167447*I)*elliptic_f(v1, v1))*arccsch(v3)))*v1, (-0.647983235144 + 1.20665952957*I)*v1 + (0.0909404921682 + 0.281538203756*I)/v3)
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/random_tests.py", line 239:
sage: random_expr(5, verbose=True)
Exception raised:
Traceback (most recent call last):
File "/home/maldun/sage/sage-4.5.2/local/bin/ncadoctest.py", line 1231, in run_one_test
self.run_one_example(test, example, filename, compileflags)
File "/home/maldun/sage/sage-4.5.2/local/bin/sagedoctest.py", line 38, in run_one_example
OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
File "/home/maldun/sage/sage-4.5.2/local/bin/ncadoctest.py", line 1172, in run_one_example
compileflags, 1) in test.globs
File "<doctest __main__.example_5[5]>", line 1, in <module>
random_expr(Integer(5), verbose=True)###line 239:
sage: random_expr(5, verbose=True)
File "/home/maldun/sage/sage-4.5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 254, in random_expr
return random_expr_helper(size, internal, leaves, verbose)
File "/home/maldun/sage/sage-4.5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 210, in random_expr_helper
return r[1](*children)
File "element.pyx", line 1529, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11992)
File "coerce.pyx", line 713, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:6126)
File "element.pyx", line 1527, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11973)
File "expression.pyx", line 2269, in sage.symbolic.expression.Expression._div_ (sage/symbolic/expression.cpp:11444)
ZeroDivisionError: Symbolic division by zero
**********************************************************************
2 items had failures:
1 of 4 in __main__.example_0
2 of 6 in __main__.example_5
***Test Failed*** 3 failures.
For whitespace errors, see the file /home/maldun/.sage//tmp/.doctest_random_tests.py
[7.7 s]
----------------------------------------------------------------------
The following tests failed:
sage -t -long "devel/sage/sage/symbolic/random_tests.py"
Total time for all tests: 7.8 seconds
</pre><p>
I quite understand these, because we have introduced new functions, but I don't understand the exception in the last one
</p>
<pre class="wiki">sage -t -long "devel/sage/sage/symbolic/pynac.pyx"
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 386:
sage: get_sfunction_from_serial(i) == foo
Expected:
True
Got:
False
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 388:
sage: py_latex_function_pystring(i, (x,y^z))
Expected:
'my args are: x, y^z'
Got:
'\\mathrm{bar}\\left(x, y^{z}\\right)'
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 478:
sage: get_sfunction_from_serial(i) == foo
Expected:
True
Got:
False
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 480:
sage: py_print_fderivative(i, (0, 1, 0, 1), (x, y^z))
Expected:
D[0, 1, 0, 1]func_with_args(x, y^z)
Got:
D[0, 1, 0, 1](foo)(x, y^z)
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 540:
sage: get_sfunction_from_serial(i) == foo
Expected:
True
Got:
False
**********************************************************************
File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 542:
sage: py_latex_fderivative(i, (0, 1, 0, 1), (x, y^z))
Expected:
D[0, 1, 0, 1]func_with_args(x, y^z)
Got:
D[0, 1, 0, 1]\left(\mathrm{bar}\right)\left(x, y^{z}\right)
**********************************************************************
3 items had failures:
2 of 19 in __main__.example_14
2 of 14 in __main__.example_16
2 of 18 in __main__.example_18
***Test Failed*** 6 failures.
For whitespace errors, see the file /home/maldun/.sage//tmp/.doctest_pynac.py
[7.3 s]
----------------------------------------------------------------------
The following tests failed:
sage -t -long "devel/sage/sage/symbolic/pynac.pyx"
Total time for all tests: 7.3 seconds
</pre><p>
And these are really strange, because when I type then into sage by hand everything works. wtf??
Can anyone have a look at these?
</p>
TicketStefan ReitererSun, 05 Sep 2010 10:30:07 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>orthogonal_polys.9.py</em>
</li>
</ul>
<p>
ortho polys with scipy support
</p>
TicketStefan ReitererSun, 03 Oct 2010 00:49:32 GMTstatus, milestone changed
https://trac.sagemath.org/ticket/9706#comment:26
https://trac.sagemath.org/ticket/9706#comment:26
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
<li><strong>milestone</strong>
changed from <em>sage-4.6</em> to <em>sage-5.0</em>
</li>
</ul>
TicketKarl-Dieter CrismanFri, 15 Oct 2010 18:42:08 GMT
https://trac.sagemath.org/ticket/9706#comment:27
https://trac.sagemath.org/ticket/9706#comment:27
<p>
Just cc:ing myself by commenting.
</p>
<p>
Also, there seems to be a lot of stuff in the latest Python file that is the same as the original one (in terms of explanation, not code). Maybe posting an updated patch (once the numpy/scipy-fest is over, which is hopefully the case) would help some of us figure this out. Thanks for working on this - there is still a lot of overhauling that symbolics could use, but this is a great step.
</p>
TicketStefan ReitererFri, 15 Oct 2010 23:31:17 GMTstatus, description changed
https://trac.sagemath.org/ticket/9706#comment:28
https://trac.sagemath.org/ticket/9706#comment:28
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=28">diff</a>)
</li>
</ul>
<p>
@kcrisman thanks for paying attention. I added now an updated patch and extended instructions.
</p>
<p>
the doctest changes in <strong>symbolic.random_tests.py</strong> are easy to explain: new functions are involved -> new random expressions. But I had to change
<code>random_expr(50, nvars=3, coeff_generator=CDF.random_element) </code> to <code>random_expr(60, nvars=3, coeff_generator=CDF.random_element) </code> or else one gets an expression generated where a division through zero occours.
</p>
<p>
As mentioned on sage-devel I repaired the doctests in <strong>symbolic.pynac.pyx</strong>, the trick is to enlarge the range of the for loop:
<code>for i in range(get_ginac_serial(), get_ginac_serial()+50):</code>
changed to
<code>for i in range(get_ginac_serial(), get_ginac_serial()+100):</code>
now it works. My explaination: since we have new functions we have longer to search, and then we reach our goal. What I can not explain is, that it works, when I type it in by hand.
</p>
<p>
All doctests pass now, so I think a review would be nice.
</p>
<p>
-maldun
</p>
TicketStefan ReitererFri, 15 Oct 2010 23:55:00 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:29
https://trac.sagemath.org/ticket/9706#comment:29
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=29">diff</a>)
</li>
</ul>
<p>
Cleaned up discription of the ticket and some comments in the ortho polys file.
</p>
TicketKarl-Dieter CrismanSat, 16 Oct 2010 03:27:13 GMT
https://trac.sagemath.org/ticket/9706#comment:30
https://trac.sagemath.org/ticket/9706#comment:30
<p>
I don't have time to review this for a while, but did take a quick look - thanks for polishing that patch! I don't think we are allowed to import numpy or scipy like that anymore, but rather have to do it in an individual function (lest startup times get huge). I don't quite understand exactly how that works, but anyway such a blanket import statement probably isn't appropriate, the way I understand what others have said.
</p>
TicketStefan ReitererSat, 16 Oct 2010 12:53:19 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_orthogonal_polys.patch</em>
</li>
</ul>
<p>
Latest version of orthogonal polys with scipy support, and changed doctests. Tested in sage-4.6.alpha3
</p>
TicketStefan ReitererSat, 16 Oct 2010 13:04:50 GMT
https://trac.sagemath.org/ticket/9706#comment:31
https://trac.sagemath.org/ticket/9706#comment:31
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:30" title="Comment 30">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
I don't have time to review this for a while, but did take a quick look - thanks for polishing that patch! I don't think we are allowed to import numpy or scipy like that anymore, but rather have to do it in an individual function (lest startup times get huge). I don't quite understand exactly how that works, but anyway such a blanket import statement probably isn't appropriate, the way I understand what others have said.
</p>
</blockquote>
<p>
But thanks for giving feedback! I know that this patch isn't easy for review because the code grew from 650 to about 2300 lines of code. But I'm happy to get at least some info.
</p>
<p>
You are right the imports didn't change since I started this ticket and importing the whole numpy and scipy packages is to much. This isn't a very good Idea if one thinks about performance either. I changed that now so that only functions that are really needed are importet. I did this also for mpmath but the problem with the global import remains.
(see above). Also changed some errors in the discription I missed and repaired a wrong doctest.
</p>
<p>
PS: If diffs or more changelogs are needed let me know. I'm keeping track with git on my machine of the changes.
</p>
TicketBurcin ErocalSat, 16 Oct 2010 15:38:16 GMTstatus, cc changed; reviewer set
https://trac.sagemath.org/ticket/9706#comment:32
https://trac.sagemath.org/ticket/9706#comment:32
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
<li><strong>cc</strong>
<em>Flavia Stan</em> added
</li>
<li><strong>reviewer</strong>
set to <em>Burcin Erocal</em>
</li>
</ul>
<p>
Great work Stefan. Your patch looks good overall, but it needs a lot of polish. Thank you very much for this.
</p>
<p>
Here are some quick comments after reading <a class="attachment" href="https://trac.sagemath.org/attachment/ticket/9706/trac_9706_orthogonal_polys.patch" title="Attachment 'trac_9706_orthogonal_polys.patch' in Ticket #9706">attachment:trac_9706_orthogonal_polys.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/9706/trac_9706_orthogonal_polys.patch" title="Download"></a>. I didn't try to apply and run the code yet. It would be better if other people try this as well since I am really short on time these days.
</p>
<ul><li>I suggest you use your real name in the HG headers. This information is used for copyright/license issues as well. In the future it might cause a lot of trouble if people have to chase down <code>maldun</code> for copyright questions.
</li><li>You shouldn't import any part of <code>numpy</code> at the module level. This slows down startup too much. See <a class="closed ticket" href="https://trac.sagemath.org/ticket/3561" title="#3561: defect: [with patch; positive review] make it so sage does *not* import numpy ... (closed: fixed)">#3561</a> for example. I'd say the same holds for <code>mpmath</code> and <code>scipy</code>.
</li><li>line 385-386 has this:
<pre class="wiki">Then after using one of these functions, it changes:: (The value is now
False for chebyshev_T because chebyshev_T uses clenshaw method instead...)
</pre>I don't think this is valid Sphinx.
</li><li>delete line 412
<pre class="wiki">#load /home/maldun/sage/sage-4.5.2/devel/sage-ortho/sage/functions/orthogonal_polys.py
</pre></li><li>line 419: he -> the
</li><li>There are no doctests for the <code>OrthogonalPolynomial</code> class, make sure your file passes <code>sage -coverage</code>
</li><li>The commented timings in the docstring of <code>OrthogonalPolynomial._clenshaw_method_()</code> are confusing. It would be better if you provide a function in the same file that does these timings automatically and prints out the results. You should at least delete this from the documentation though.
</li><li>In the docstring of <code>OrthogonalPolynomial._eval_()</code>
<ul><li>remove the empty first line (line 494) of
</li><li>remove the commented out timings as well
</li><li>you need an empty line after <code>EXAMPLES::</code>
</li><li>the empty last line should be removed
</li></ul></li><li>add some comments to the <code>OrthogonalPolynomial._eval_()</code> method to indicate what you're trying to do with these tests.
<ul><li>lines 583-593 have a confusing comment and a bug
<pre class="wiki">try:
#s = maxima(self._maxima_init_evaled_(*args))
#This above is very inefficient! The older
#methods were much faster...
return self._maxima_init_evaled_(*args)
except TypeError:
return None
if self._maxima_name in repr(s):
return None
else:
return s.sage()
</pre></li></ul></li><li>You don't need to state "Class for" on line 598, "The Chebyshev ..." is enough.
</li><li>Why do you delete the <code>chebyshev_T(2,x)</code> test on line 371? You can just add the new ones after that.
</li><li>line 626, <code>EXAMPLES:</code> -> <code>EXAMPLES::</code>
</li><li>Don't use <code>*args</code> or <code>**kwds</code> when you don't need them. Name the arguments and be explicit. Remember the "Zen of Python", "Explicit is better than implicit."
</li><li>OK, generally, fix the docstrings to conform to Sphinx standards. This should be documented somewhere in the developers guide.
</li><li>line 673, <code>_maxima_init_evaled_()</code> doesn't have doctests.
</li><li>line 678 - , <code>_clenshaw_method_()</code>
<ul><li>docstring is not indented properly.
</li><li>It would be better to put the recursion formula in the docstring.
</li></ul></li><li>line 790 <code>_clenshaw_method_()</code> doesn't have doctests.
</li><li>There is something wrong with the <code>_maxima_init_evaled_()</code> on line 821. Are you sure this function shouldn't just return a string to be run in maxima? How do we know that doctest actually calls this function? In any case, the right way to convert a maxima object to sage is to run <code>.sage()</code> on it. Never use <code>sage_eval()</code> on a string in the Sage library.
</li><li>Calls to mpmath should be able to use the precision directly from the type of the argument now. Are you sure all this is necessary:
<pre class="wiki">try:
step_parent = kwds['parent']
except KeyError:
step_parent = parent(args[-1])
try:
precision = step_parent.prec()
except AttributeError:
precision = RR.prec()
</pre>See <a class="closed ticket" href="https://trac.sagemath.org/ticket/9566" title="#9566: enhancement: Allow sage.libs.mpmath.call(..., parent=something) (closed: fixed)">#9566</a>.
</li><li>line 924, change the error message to something more professional. "Derivative w.r.t. to the index is not supported, yet, and perhaps never will be..." is not acceptable. "Derivatives with respect to the index is not supported." would be enough.
</li><li>Document the derivative formula in the docstring, using proper math notation
</li><li>What needs to be discussed from the comments on line 968-974?
</li><li>Same for lines 1058-1060?
</li><li>no doctests for <code>_clenshaw_method_()</code> on line 1156.
</li><li>no doctests for <code>_maxima_init_evaled_()</code> on line 1189.
</li></ul><p>
I give up at this point. It seems that there are similar issues in the rest of the file as well.
</p>
<p>
After you clean up the code according to the comments above, perhaps a native English speaker like Karl-Dieter or Minh can help with the documentation.
</p>
<p>
Thanks again for all your work.
</p>
TicketKarl-Dieter CrismanThu, 10 Nov 2011 14:11:28 GMTcc changed
https://trac.sagemath.org/ticket/9706#comment:33
https://trac.sagemath.org/ticket/9706#comment:33
<ul>
<li><strong>cc</strong>
<em>Karl-Dieter Crisman</em> added
</li>
</ul>
TicketJeroen DemeyerTue, 13 Aug 2013 15:35:53 GMTmilestone changed
https://trac.sagemath.org/ticket/9706#comment:34
https://trac.sagemath.org/ticket/9706#comment:34
<ul>
<li><strong>milestone</strong>
changed from <em>sage-5.11</em> to <em>sage-5.12</em>
</li>
</ul>
TicketStefan ReitererSat, 09 Nov 2013 22:29:50 GMT
https://trac.sagemath.org/ticket/9706#comment:35
https://trac.sagemath.org/ticket/9706#comment:35
<p>
Hi!
</p>
<p>
I will now retry to build the new orthogonal polynomials.
The last time I ran out of time due to my phd studies/theses
this time i will split the changes up into several patches. So it will be easier to
apply the changes step by step, and the review process gets simpler.
</p>
<p>
Hope this time everything will work out!
</p>
TicketStefan ReitererSat, 09 Nov 2013 22:30:36 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:36
https://trac.sagemath.org/ticket/9706#comment:36
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
</ul>
TicketTravis ScrimshawSat, 23 Nov 2013 07:58:00 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:37
https://trac.sagemath.org/ticket/9706#comment:37
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=37">diff</a>)
</li>
</ul>
<p>
Here is a review patch which does a bunch of documentation formatting tweaks. There are probably one or two other things that will need to be addressed, but I'd like to get the ball rolling on this again (and I need some sleep right now).
</p>
<p>
Best,<br />
Travis
</p>
<p>
For patchbot:
</p>
<p>
Apply: trac_9706_chebyshev.patch, trac_9706-review-ts.patch
</p>
TicketStefan ReitererFri, 29 Nov 2013 16:44:09 GMTsummary changed
https://trac.sagemath.org/ticket/9706#comment:38
https://trac.sagemath.org/ticket/9706#comment:38
<ul>
<li><strong>summary</strong>
changed from <em>New Version of orthogonal Polynomials</em> to <em>New Version of orthogonal Polynomials (Part I: Base Class and Chebyshev Polynomials)</em>
</li>
</ul>
<p>
Thanks for reviewing. It would be great if the new Chebyshev Polynomials could be added. If this ticket is done I will open the next issue and start working on the Legendre Polynomials
</p>
TicketTravis ScrimshawMon, 02 Dec 2013 23:45:29 GMT
https://trac.sagemath.org/ticket/9706#comment:39
https://trac.sagemath.org/ticket/9706#comment:39
<p>
Okay, I've done a few other tweaks and I'd okay with it. If you're happy with my changes, then go ahead and set this to positive review.
</p>
<p>
For patchbot:
</p>
<p>
Apply: trac_9706_chebyshev.patch, trac_9706-review-ts.patch
</p>
TicketStefan ReitererTue, 03 Dec 2013 12:15:58 GMTstatus, reviewer changed
https://trac.sagemath.org/ticket/9706#comment:40
https://trac.sagemath.org/ticket/9706#comment:40
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
<li><strong>reviewer</strong>
changed from <em>Burcin Erocal</em> to <em>Burcin Erocal, Travis Scrimshaw</em>
</li>
</ul>
<p>
Thanks for your hard work in correcting all those small mistakes!
</p>
<p>
I'm very happy, that finally the new ortho polys are going into sage!
</p>
TicketKarl-Dieter CrismanTue, 03 Dec 2013 16:07:49 GMTreviewer, author changed
https://trac.sagemath.org/ticket/9706#comment:41
https://trac.sagemath.org/ticket/9706#comment:41
<ul>
<li><strong>reviewer</strong>
changed from <em>Burcin Erocal, Travis Scrimshaw</em> to <em>Burcin Erocal, Travis Scrimshaw, Stefan Reiterer</em>
</li>
<li><strong>author</strong>
changed from <em>Stefan Reiterer</em> to <em>Stefan Reiterer, Travis Scrimshaw</em>
</li>
</ul>
<p>
I'm elevating Travis to an author because these are quite substantial review changes - thanks for being so meticulous on the formatting etc!
</p>
<p>
I'd love one final check from either of you. There are a lot of imports added; I think most should be okay but the Maxima-related ones scare me, so if you can just check that startup time hasn't budged more than a couple milliseconds, that would be helpful. I don't think this should import numpy, at least!
</p>
TicketJeroen DemeyerTue, 03 Dec 2013 17:12:07 GMT
https://trac.sagemath.org/ticket/9706#comment:42
https://trac.sagemath.org/ticket/9706#comment:42
<p>
Sorry to spoil the party, but this is a regression:
</p>
<pre class="wiki">sage: K.<a> = NumberField(x^3-x-1)
sage: chebyshev_T(10^3,a)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-17-aa97c56dd147> in <module>()
----> 1 chebyshev_T(Integer(10)**Integer(3),a)
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5279)()
TypeError: cannot coerce arguments: no canonical coercion from Number Field in a with defining polynomial x^3 - x - 1 to Symbolic Ring
</pre><p>
(and yes: Chebyshev polynomials are important in number theory)
</p>
TicketJeroen DemeyerTue, 03 Dec 2013 17:12:15 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:43
https://trac.sagemath.org/ticket/9706#comment:43
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>needs_work</em>
</li>
</ul>
TicketJeroen DemeyerTue, 03 Dec 2013 17:18:56 GMT
https://trac.sagemath.org/ticket/9706#comment:44
https://trac.sagemath.org/ticket/9706#comment:44
<p>
I would also like to point out that PARI is faster at evaluating Chebyshev polynomials:
</p>
<pre class="wiki">sage: timeit('''chebyshev_T(10^5,2)''')
5 loops, best of 3: 270 ms per loop
sage: timeit('''pari('polchebyshev (10^5,1,2)')''')
625 loops, best of 3: 447 µs per loop
</pre><pre class="wiki">sage: timeit('''chebyshev_T(10^5, Mod(2,1009))''')
5 loops, best of 3: 208 ms per loop
sage: timeit('''pari('polchebyshev(10^5, 1, Mod(2,1009))')''')
625 loops, best of 3: 11.5 µs per loop
</pre><p>
We should definately use PARI to evaluate Chebyshev polynomials for certain types of input.
</p>
TicketJeroen DemeyerTue, 03 Dec 2013 17:19:11 GMT
https://trac.sagemath.org/ticket/9706#comment:45
https://trac.sagemath.org/ticket/9706#comment:45
<p>
Another regression:
</p>
<pre class="wiki">sage: parent(chebyshev_T(10^2, RIF(2)))
Real Field with 53 bits of precision
</pre>
TicketJeroen DemeyerTue, 03 Dec 2013 17:29:00 GMT
https://trac.sagemath.org/ticket/9706#comment:46
https://trac.sagemath.org/ticket/9706#comment:46
<p>
The "Clenshaw method" uses a very naive method of evaluating the recursion which needs <code>O(n)</code> steps, while there is a much faster method (which compute <code>T_2n</code> and <code>U_2n</code> in function of <code>T_n</code> and <code>U_n</code>) which only needs <code>O(log(n))</code> steps.
</p>
<p>
Even this is totally feasible:
</p>
<pre class="wiki">sage: timeit('''pari('polchebyshev(10^10, 1, Mod(2,1009))')''')
625 loops, best of 3: 16.3 µs per loop
</pre>
TicketJeroen DemeyerTue, 03 Dec 2013 17:32:45 GMT
https://trac.sagemath.org/ticket/9706#comment:47
https://trac.sagemath.org/ticket/9706#comment:47
<p>
This is also bad:
</p>
<pre class="wiki">sage: R.<x> = QQ[]
sage: parent(chebyshev_T(5, x))
Symbolic Ring
</pre>
TicketJeroen DemeyerTue, 03 Dec 2013 17:35:24 GMT
https://trac.sagemath.org/ticket/9706#comment:48
https://trac.sagemath.org/ticket/9706#comment:48
<p>
Another suggestion: why not make the base class more general and call it <code>SymbolicPolynomial</code>? I think that's a more natural class of functions, which could (in the future) also include cyclotomic polynomials for example.
</p>
TicketJeroen DemeyerTue, 03 Dec 2013 17:42:30 GMT
https://trac.sagemath.org/ticket/9706#comment:49
https://trac.sagemath.org/ticket/9706#comment:49
<p>
I propose the logic for evaluating <code>chebyshev_T(n, x)</code> should be:
</p>
<ol><li>if <code>x</code> is symbolic, then use the method of the current patch.
</li><li>if <code>x</code> is not symbolic, try evaluation using PARI.
</li><li>if conversion to PARI fails (for example for <code>RIF</code>), use an efficient <code>O(log(n))</code> recursion algorithm.
</li></ol>
TicketFredrik JohanssonTue, 03 Dec 2013 18:08:57 GMT
https://trac.sagemath.org/ticket/9706#comment:50
https://trac.sagemath.org/ticket/9706#comment:50
<p>
Jeroen: do you know a reference for the recursion pari uses?
</p>
TicketJeroen DemeyerTue, 03 Dec 2013 20:30:50 GMT
https://trac.sagemath.org/ticket/9706#comment:51
https://trac.sagemath.org/ticket/9706#comment:51
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:50" title="Comment 50">fredrik.johansson</a>:
</p>
<blockquote class="citation">
<p>
Jeroen: do you know a reference for the recursion pari uses?
</p>
</blockquote>
<p>
No, but it's pretty straight-forward (think doubling formulas for <code>cos</code> and <code>sin</code>).
</p>
TicketStefan ReitererWed, 04 Dec 2013 07:47:43 GMT
https://trac.sagemath.org/ticket/9706#comment:52
https://trac.sagemath.org/ticket/9706#comment:52
<p>
Thank you all for the input! I think it is still a good idea that I only implement chebychev polys for now, since there are a lot of improvements out there.
</p>
<p>
@jdemeyer
</p>
<p>
@Bugs I will look into this. And yes I'm aware that Chebyshev Polynomials are important to number theory since there are quite interesting generalizations on general fields.
</p>
<p>
@Clenshaw: PARI is a good hint, I will look into this. And I already think I know how to generalize the clenshaw method, to get O(log N).
The reason for this quite naive choice, was that the method can be applied to all ortho polys. Nevertheless I will adapt it on Chebychev Polys since there are more possibilities since we can use trigonometric formulas. I think the benefit of implementing it directly in sage is that there is less trouble if one wants to use more general data types, since there are no type casts. I will try to find an optimal way for this. (Maybe an additional switch)
</p>
<p>
@<a class="missing wiki">SymbolicPolynomial?</a>: I don't think this is a good idea because ortho polys are quite special even among the polynomials. But If you really would like to have a <a class="missing wiki">SymbolicPolynomial?</a> class I would propose to introduce the <a class="missing wiki">SymbolicPolynomial?</a> class, and derive the <a class="missing wiki">OrthogonalPolynomials?</a> from that class.
I make the following suggestion: I will finish the <a class="missing wiki">OrthogonalPolynomials?</a> with the current design. And then open an new ticket where we discuss the design of a general polynomial parent class. Fortunately, such design changes are very easy to implement in Python, and I don't see any big problem in introducing an intermediate class.
But if you want to introduce such a class there are some major decisions to make:
-Where do we put this class? (such a general class should not belong to orthogonoal_polys.py )
-What should all <a class="missing wiki">SymbolicPolynomials?</a> have in common?
-What should they have concerning other general polynomials?
</p>
<p>
But it would be really good to add the ortho polys, and I really want to finish this task.
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 08:38:33 GMT
https://trac.sagemath.org/ticket/9706#comment:53
https://trac.sagemath.org/ticket/9706#comment:53
<blockquote class="citation">
<p>
<code>SymbolicPolynomial</code>: I don't think this is a good idea because ortho polys are quite special even among the polynomials.
</p>
</blockquote>
<p>
What's special about orthogonal polynomials from a <strong>computer algebra</strong> point of view? I can tell you that "symbolic polynomials" are special because you generally want to be able to evaluate them for any ring element (as opposed to other symbolic functions, which often only make sense in real or complex fields).
</p>
<blockquote class="citation">
<p>
I make the following suggestion: I will finish the <a class="missing wiki">OrthogonalPolynomials?</a> with the current design.
</p>
</blockquote>
<p>
Sure...
</p>
TicketStefan ReitererWed, 04 Dec 2013 09:08:49 GMT
https://trac.sagemath.org/ticket/9706#comment:54
https://trac.sagemath.org/ticket/9706#comment:54
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:53" title="Comment 53">jdemeyer</a>:
</p>
<blockquote class="citation">
<blockquote class="citation">
<p>
<code>SymbolicPolynomial</code>: I don't think this is a good idea because ortho polys are quite special even among the polynomials.
</p>
</blockquote>
<p>
What's special about orthogonal polynomials from a <strong>computer algebra</strong> point of view? I can tell you that "symbolic polynomials" are special because you generally want to be able to evaluate them for any ring element (as opposed to other symbolic functions, which often only make sense in real or complex fields).
</p>
</blockquote>
<p>
It makes a difference in an <strong>OO-Design</strong> point of view, because you can apply a whole bunch of evaluation techniques and tricks due to the three term recursion (e.g. clenshaw method or the eval_numpy method are not needed for symbolic polynomials, since there are no methods for that). Thats the reason why I suggested to derive them from a base class on top to avoid redundant methods. That would be a clean solution. And I'm already thinking on some features, which I want to implement in future version, which are very specific to ortho polys. And if you want to introduce such a general class, you should not put it into orthogonal_polys.py, but somewhere else, were it fits better into the class hierarchy.
</p>
<p>
Please don't get me wrong, I'm not stating, that I think it is a bad idea to introduce such an abstract base class, but If you want a clean OO Design, it is not done by a simple renaming.
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 09:14:39 GMT
https://trac.sagemath.org/ticket/9706#comment:55
https://trac.sagemath.org/ticket/9706#comment:55
<p>
Evaluating <code>chebyshev_T(n,x)</code> can be done as
</p>
<pre class="wiki">(Matrix(2,2,[x,x^2-1,1,x])^n)[0,0]
</pre><p>
While <code>chebyshev_U(n-1,x)</code> equals
</p>
<pre class="wiki">(Matrix(2,2,[x,x^2-1,1,x])^n)[1,0]
</pre><p>
These can be evaluated with <code>O(log(n))</code> operations.
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 09:16:03 GMT
https://trac.sagemath.org/ticket/9706#comment:56
https://trac.sagemath.org/ticket/9706#comment:56
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:54" title="Comment 54">maldun</a>:
</p>
<blockquote class="citation">
<p>
And I'm already thinking on some features, which I want to implement in future version, which are very specific to ortho polys.
</p>
</blockquote>
<p>
That alone is a very good to keep the <code>OrthogonalPolynomial</code> class.
</p>
TicketStefan ReitererWed, 04 Dec 2013 09:24:30 GMT
https://trac.sagemath.org/ticket/9706#comment:57
https://trac.sagemath.org/ticket/9706#comment:57
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:55" title="Comment 55">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Evaluating <code>chebyshev_T(n,x)</code> can be done as
</p>
<pre class="wiki">(Matrix(2,2,[x,x^2-1,1,x])^n)[0,0]
</pre><p>
While <code>chebyshev_U(n-1,x)</code> equals
</p>
<pre class="wiki">(Matrix(2,2,[x,x^2-1,1,x])^n)[1,0]
</pre><p>
These can be evaluated with <code>O(log(n))</code> operations.
</p>
</blockquote>
<p>
Thanks for the hint. I have also an own idea to implement this. My implementation should be optimal with respect to the flop count, but yours could be faster since matrix multiplication and powers are well optimized. I will compare both methods and use the faster one.
</p>
<p>
For reference: The method which I mean is based on the generalized recursion formula (originating from the cosine addition theorem):
T_{n+m} + T{n-m} = 2 T_n T_m
</p>
<p>
For T_N one can now use the binary representation of N to recursively build T_N in O(log N) time
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 11:59:36 GMT
https://trac.sagemath.org/ticket/9706#comment:58
https://trac.sagemath.org/ticket/9706#comment:58
<p>
This should be a <code>ValueError</code>:
</p>
<pre class="wiki">sage: chebyshev_T(1/2,0)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-7-830f13ad2f0d> in <module>()
----> 1 chebyshev_T(Integer(1)/Integer(2),Integer(0))
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5531)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _eval_(self, *args)
489 if not is_Expression(args[-1]):
490 try:
--> 491 return self._evalf_(*args)
492 except AttributeError:
493 pass
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _evalf_(self, *args, **kwds)
606 precision = step_parent.prec()
607 except AttributeError:
--> 608 precision = RR.prec()
609
610 from sage.libs.mpmath.all import call as mpcall
NameError: global name 'RR' is not defined
</pre>
TicketJeroen DemeyerWed, 04 Dec 2013 12:01:50 GMT
https://trac.sagemath.org/ticket/9706#comment:59
https://trac.sagemath.org/ticket/9706#comment:59
<p>
This should be <code>ArithmeticError</code> (I guess), since deriving w.r.t. the index simply isn't defined:
</p>
<pre class="wiki">sage: var('n,x')
(n, x)
sage: chebyshev_T(n,x).diff(n)
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
<ipython-input-14-a23f5209eb49> in <module>()
----> 1 chebyshev_T(n,x).diff(n)
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.derivative (sage/symbolic/expression.cpp:16561)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/misc/derivative.so in sage.misc.derivative.multi_derivative (sage/misc/derivative.c:2715)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._derivative (sage/symbolic/expression.cpp:16951)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _derivative_(self, *args, **kwds)
714 diff_param = kwds['diff_param']
715 if diff_param == 0:
--> 716 raise NotImplementedError("derivative w.r.t. to the index is not supported yet")
717
718 return args[0]*chebyshev_U(args[0]-1,args[1])
NotImplementedError: derivative w.r.t. to the index is not supported yet
</pre><p>
Also: doctest your exceptions.
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 12:07:52 GMT
https://trac.sagemath.org/ticket/9706#comment:60
https://trac.sagemath.org/ticket/9706#comment:60
<p>
I personally like to define <code>chebyshev_T(-n,x) = chebyshev_T(n,x)</code> and <code>chebyshev_U(-n,x) = -chebyshev_U(n-2,x)</code>. This makes sense from a number-theoretic point of view and is also consistent with <code>chebyshev_T(n,cos(x)) = cos(n*x)</code> and <code>chebyshev_U(n-1,cos(x)) = sin(n*x)/sin(x)</code>.
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 14:42:23 GMT
https://trac.sagemath.org/ticket/9706#comment:61
https://trac.sagemath.org/ticket/9706#comment:61
<p>
For real/complex fields, you should use
</p>
<pre class="wiki">chebyshev_T(n,x) = ((x+sqrt(x^2-1))^n + (x-sqrt(x^2-1))^n)/2
chebyshev_U(n-1,x) = ((x+sqrt(x^2-1))^n - (x-sqrt(x^2-1))^n)/(2*sqrt(x^2-1))
</pre>
TicketFredrik JohanssonWed, 04 Dec 2013 15:06:49 GMT
https://trac.sagemath.org/ticket/9706#comment:62
https://trac.sagemath.org/ticket/9706#comment:62
<p>
The doubling recursion formulas should be better also for real/complex fields, I think: computing nth powers is basically the same amount of work, and it's best to avoid the square roots (especially for real |x| < 1).
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 15:16:43 GMT
https://trac.sagemath.org/ticket/9706#comment:63
https://trac.sagemath.org/ticket/9706#comment:63
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:62" title="Comment 62">fredrik.johansson</a>:
</p>
<blockquote class="citation">
<p>
The doubling recursion formulas should be better also for real/complex fields, I think: computing nth powers is basically the same amount of work, and it's best to avoid the square roots (especially for real |x| < 1).
</p>
</blockquote>
<p>
The matrix algorithm does seem more numerically stable (checked by using both algorithms inside <code>RIF</code>). So it's easy then, if there is one algorithm which is obviously the best.
</p>
TicketStefan ReitererWed, 04 Dec 2013 15:33:48 GMT
https://trac.sagemath.org/ticket/9706#comment:64
https://trac.sagemath.org/ticket/9706#comment:64
<p>
I found an interesting paper on this topic: <a class="ext-link" href="http://www.mathematik.uni-kassel.de/~koepf/cheby.pdf"><span class="icon"></span>http://www.mathematik.uni-kassel.de/~koepf/cheby.pdf</a>
</p>
<p>
Maybe this could give some input to the discussion.
</p>
<p>
It states that for an expanded representation the approach I have
initally chosen (series expansion) should be the best for the symbolic
evaluation, since this is somehow that what the user expects from ohter CAS.
However, the proposed recursive/symbolic method would be interesting too.
(p16f), since it gives a more compact form for large n. Maybe a switch for
n >= 100?
</p>
<p>
Considering rational numbers as input, the recursion formula works best, due to this paper. This should be considered too.
</p>
<p>
@Matrix Multiplication: It's also my favorite, but I will compare first.
</p>
TicketFredrik JohanssonWed, 04 Dec 2013 16:45:02 GMT
https://trac.sagemath.org/ticket/9706#comment:65
https://trac.sagemath.org/ticket/9706#comment:65
<p>
Here is an implementation of the divide-and-conquer algorithm that doesn't require caching (it only makes one recursive call). It might look even nicer if one rewrites it in iterative form. I think it's also equivalent to the code in pari. It should be faster than matrix powering by a constant factor, just like the analogous Fibonacci number algorithms.
</p>
<pre class="wiki">def chebyshev_t(n, x):
# returns (T(n,x), T(n-1,x)), or (T(n,x), _) if both=False
def recur(n, x, both=False):
if n == 0:
return 1, x
if n == 1:
return x, 1
a, b = recur((n+1)//2, x, both or n % 2)
if n % 2 == 0:
return 2*a^2 - 1, both and 2*a*b - x
else:
return 2*a*b - x, both and 2*b^2 - 1
return recur(n, x, False)[0]
</pre><p>
Come to think of it, it might even be useful to publicly expose a method that returns both T(n,x) and T(n-1,x) simultaneously.
</p>
<p>
Similar code for U (using same algorithm as Pari):
</p>
<pre class="wiki">def chebyshev_u(n, x):
def recur(n, x, both=False):
if n == 0:
return 1, both and 2*x
if n == 1:
return 2*x, both and 4*x^2-1
a, b = recur((n-1)//2, x, True)
if n % 2 == 0:
return (b+a)*(b-a), both and 2*b*(x*b-a)
else:
return 2*a*(b-x*a), both and (b+a)*(b-a)
return recur(n, x, False)[0]
</pre><p>
Edit: streamlined the code.
</p>
TicketStefan ReitererWed, 04 Dec 2013 19:46:44 GMT
https://trac.sagemath.org/ticket/9706#comment:66
https://trac.sagemath.org/ticket/9706#comment:66
<p>
I think your recursive implementation is very good. If you try to implement it iteratively, you have to consider some cases (current in row even/odd and next in row even/odd), and the code gets quite ugly in my opinion. And I think Knuth is right. One should prefer readable code over over optimzed "faster" code...
</p>
TicketStefan ReitererWed, 04 Dec 2013 20:39:45 GMT
https://trac.sagemath.org/ticket/9706#comment:67
https://trac.sagemath.org/ticket/9706#comment:67
<p>
this would be a functioning iterative algorithm:
</p>
<pre class="wiki">def chebyshev_t(n,x):
if n == 0:
return 1
elif n == 1:
return x
elif n == 2:
return 2*x**2-1
else:
T_c = x
T_p = 2*x**2 -1
for k in range(floor(log(n,2)),0,-1):
T_p_old = T_p
T_c_old = T_c
if (n//2**(k-1)) % 2 == 0: # next is even
T_p = 2*T_p_old*T_c_old - x
T_c = 2*T_c_old**2 - 1
elif (n//2**(k-1)) % 2 == 1: # next is odd
T_p = 2*T_p_old**2 - 1
T_c = 2*T_p_old*T_c_old - x
# Cases for output
if log(n - 1,2) in ZZ or log(n-2,2) in ZZ:
return T_c
elif n % 2 == 0:
if n//2 % 2 == 0:
return T_c
else:
return T_p
elif n % 2 == 1:
if n//2 % 2 == 0:
return T_p
else:
return T_c
</pre><p>
I made it shorter. What is preferable? recursive or iterative? Normaly iterative, but in this case it is longer due to the indices battles ...
</p>
<p>
Edit: I also measured time: recursive is faster, even if I do some optimization (e.g. xrange)
</p>
TicketStefan ReitererWed, 04 Dec 2013 22:01:05 GMT
https://trac.sagemath.org/ticket/9706#comment:68
https://trac.sagemath.org/ticket/9706#comment:68
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:42" title="Comment 42">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Sorry to spoil the party, but this is a regression:
</p>
<pre class="wiki">sage: K.<a> = NumberField(x^3-x-1)
sage: chebyshev_T(10^3,a)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-17-aa97c56dd147> in <module>()
----> 1 chebyshev_T(Integer(10)**Integer(3),a)
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5279)()
TypeError: cannot coerce arguments: no canonical coercion from Number Field in a with defining polynomial x^3 - x - 1 to Symbolic Ring
</pre><p>
(and yes: Chebyshev polynomials are important in number theory)
</p>
</blockquote>
<p>
I had now a look on these regressions: The problem originates from the fact, tat the ortho polys are now symbolic functions instead of simple function call in maxima.
the symbolic functions have a lot more mechanisms concerning type coercions.
</p>
<p>
The only way around this would be a hack to chatch this coercion troubles, e.g.
</p>
<pre class="wiki">def __call__(self,n,x):
try:
super(Func_chebyshev_T,self).__call__(n,x)
except TypeError:
return self._clenshaw_method_(n,x)
</pre><p>
would work for that problem and some others. But how to catch such things correctly?
Use a try catch, or make several is_instance checks? This could get a quite long list on coercions...
</p>
TicketStefan ReitererWed, 04 Dec 2013 23:30:50 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_new_clenshaw.patch</em>
</li>
</ul>
<p>
new fixes and patches
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 23:51:34 GMT
https://trac.sagemath.org/ticket/9706#comment:69
https://trac.sagemath.org/ticket/9706#comment:69
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:68" title="Comment 68">maldun</a>:
</p>
<blockquote class="citation">
<p>
would work for that problem and some others. But how to catch such things correctly?
</p>
</blockquote>
<p>
I think there are only two cases: the "symbolic" case and the "algebraic" case. The latter means that we really consider the polynomial as a polynomial, not a symbolic function. In <code>chebyshev(n,x)</code>, if either <code>n</code> or <code>x</code> is symbolic, we are in the symbolic case, otherwise we're in the algebraic case. In the algebraic case, the index <code>n</code> must be a concrete integer and we use the iterative algorithm.
</p>
TicketJeroen DemeyerWed, 04 Dec 2013 23:53:48 GMT
https://trac.sagemath.org/ticket/9706#comment:70
https://trac.sagemath.org/ticket/9706#comment:70
<p>
maldun: I don't like all the log's in your approach (I don't think you need them), but otherwise I'm happy with either the recursive or iterative approach.
</p>
TicketStefan ReitererThu, 05 Dec 2013 08:17:51 GMT
https://trac.sagemath.org/ticket/9706#comment:71
https://trac.sagemath.org/ticket/9706#comment:71
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:69" title="Comment 69">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:68" title="Comment 68">maldun</a>:
</p>
<blockquote class="citation">
<p>
would work for that problem and some others. But how to catch such things correctly?
</p>
</blockquote>
<p>
I think there are only two cases: the "symbolic" case and the "algebraic" case. The latter means that we really consider the polynomial as a polynomial, not a symbolic function. In <code>chebyshev(n,x)</code>, if either <code>n</code> or <code>x</code> is symbolic, we are in the symbolic case, otherwise we're in the algebraic case. In the algebraic case, the index <code>n</code> must be a concrete integer and we use the iterative algorithm.
</p>
</blockquote>
<p>
Thank you for the input, that's a great idea!
</p>
<p>
Based on that I propose the following switch:
</p>
<pre class="wiki">def __call__(self,n,x):
if n in ZZ: # check if n is integer -> consider polynomial as algebraic structure
self._eval_(n,x) # Let eval methode decide which is best
else:
super(OrthogonalPolynomial,self).__call__(n,x)
</pre>
TicketJeroen DemeyerThu, 05 Dec 2013 08:43:18 GMT
https://trac.sagemath.org/ticket/9706#comment:72
https://trac.sagemath.org/ticket/9706#comment:72
<p>
I think that <code>chebyshev_T(1/2, 2)</code> should raise a <code>ValueError</code> (or can we make sense of this?). So, in your code there should really be 3 cases: integer, symbolic and "something else" which is always an error.
</p>
<p>
So, I would do something like
</p>
<pre class="wiki">def __call__(self,n,x):
if is_Expression(n):
return super(OrthogonalPolynomial, self).__call__(n,x)
# We consider the polynomial really as a polynomial,
# not a symbolic expression.
try:
n = ZZ(n)
except StandardError:
raise ValueError("Index for symbolic polynomials must be an integer")
return self._eval_polynomial(n, x)
</pre>
TicketJeroen DemeyerThu, 05 Dec 2013 08:48:18 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:73
https://trac.sagemath.org/ticket/9706#comment:73
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=73">diff</a>)
</li>
</ul>
TicketStefan ReitererThu, 05 Dec 2013 08:54:31 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:74
https://trac.sagemath.org/ticket/9706#comment:74
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=74">diff</a>)
</li>
</ul>
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:70" title="Comment 70">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
maldun: I don't like all the log's in your approach (I don't think you need them), but otherwise I'm happy with either the recursive or iterative approach.
</p>
</blockquote>
<p>
I removed the logs, it is much faster now, but still slower than the recursive implementation about a constant factor of 2 (which is not much considering that we are talking about µs, but still)
</p>
TicketStefan ReitererThu, 05 Dec 2013 09:06:11 GMT
https://trac.sagemath.org/ticket/9706#comment:75
https://trac.sagemath.org/ticket/9706#comment:75
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:72" title="Comment 72">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
I think that <code>chebyshev_T(1/2, 2)</code> should raise a <code>ValueError</code> (or can we make sense of this?). So, in your code there should really be 3 cases: integer, symbolic and "something else" which is always an error.
</p>
<p>
So, I would do something like
</p>
<pre class="wiki">def __call__(self,n,x):
if is_Expression(n):
return super(OrthogonalPolynomial, self).__call__(n,x)
# We consider the polynomial really as a polynomial,
# not a symbolic expression.
try:
n = ZZ(n)
except StandardError:
raise ValueError("Index for symbolic polynomials must be an integer")
return self._eval_polynomial(n, x)
</pre></blockquote>
<p>
No there are several other cases to consider: You can also have complex and real input values if you consider a chebyshev polynomial as extension of the Hypergeometric function 1F2 like mpmath does: <a class="ext-link" href="http://mpmath.googlecode.com/svn/trunk/doc/build/functions/orthogonal.html#chebyt"><span class="icon"></span>http://mpmath.googlecode.com/svn/trunk/doc/build/functions/orthogonal.html#chebyt</a>
That's the reason for the message in the _derive_ method, since it would be theoretically possible to differentiate a chebyshev polynomial with respect to the index.
</p>
<p>
So the way I proposed makes sense, since this could be important for symbolic computation purposes where hypergeometric functions play an important role (eg. Zeilberger algorithm), and of course for analytical considerations.
</p>
TicketJeroen DemeyerThu, 05 Dec 2013 09:27:15 GMT
https://trac.sagemath.org/ticket/9706#comment:76
https://trac.sagemath.org/ticket/9706#comment:76
<p>
OK, I agree.
</p>
<p>
What remains to do:
</p>
<ul><li>support <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:60" title="Comment 60">negative indices</a>
</li><li>add many doctests (essentially, all the examples I mentioned on this ticket should become doctests)
</li><li>for ease of reviewing, fold all the patches into one patch.
</li></ul>
TicketStefan ReitererThu, 05 Dec 2013 19:12:45 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_chebyshev_patch_new.patch</em>
</li>
</ul>
<p>
Patch in a whole with corrections etc.
</p>
TicketStefan ReitererThu, 05 Dec 2013 19:14:23 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:77
https://trac.sagemath.org/ticket/9706#comment:77
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
</ul>
<p>
New patch attached. I incorporated all changes discussed.
</p>
<p>
@Pari I tried the evaluation with pari, but with Fredericks recursion, there is no gain in speed, due to the type checks beforehand. And the recursion in sage avoid type conversions.
</p>
TicketJeroen DemeyerThu, 05 Dec 2013 19:40:55 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:78
https://trac.sagemath.org/ticket/9706#comment:78
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=78">diff</a>)
</li>
</ul>
<p>
meldun: please adjust the "apply" section in the ticket description so it's clear which patch(es) should be applied.
</p>
TicketJeroen DemeyerThu, 05 Dec 2013 19:44:49 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:79
https://trac.sagemath.org/ticket/9706#comment:79
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=79">diff</a>)
</li>
</ul>
TicketJeroen DemeyerThu, 05 Dec 2013 19:46:19 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_chebyshev.patch</em>
</li>
</ul>
<p>
cheby changes
</p>
TicketJeroen DemeyerThu, 05 Dec 2013 19:46:35 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:80
https://trac.sagemath.org/ticket/9706#comment:80
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=80">diff</a>)
</li>
</ul>
<p>
Folded the 4 patches.
</p>
TicketJeroen DemeyerThu, 05 Dec 2013 19:47:04 GMTsummary changed
https://trac.sagemath.org/ticket/9706#comment:81
https://trac.sagemath.org/ticket/9706#comment:81
<ul>
<li><strong>summary</strong>
changed from <em>New Version of orthogonal Polynomials (Part I: Base Class and Chebyshev Polynomials)</em> to <em>Symbolic Chebyshev polynomials</em>
</li>
</ul>
TicketJeroen DemeyerThu, 05 Dec 2013 19:50:06 GMT
https://trac.sagemath.org/ticket/9706#comment:82
https://trac.sagemath.org/ticket/9706#comment:82
<p>
What's the advantage of the <code>_clenshaw_method_</code> over the recursive method? I see no need for the two different implementations and suggest to remove <code>_clenshaw_method_</code>.
</p>
TicketJeroen DemeyerThu, 05 Dec 2013 19:52:17 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:83
https://trac.sagemath.org/ticket/9706#comment:83
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
<p>
More doctests are needed, this still doesn't work:
</p>
<pre class="wiki">sage: chebyshev_T(1/2, 0)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-2-c142cc68c50b> in <module>()
----> 1 chebyshev_T(Integer(1)/Integer(2), Integer(0))
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in __call__(self, n, x)
554 return self._eval_(n,x) # Let eval methode decide which is best
555 else: # Consider OrthogonalPolynomial as symbol
--> 556 return super(OrthogonalPolynomial,self).__call__(n,x)
557
558
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (sage/symbolic/function.cpp:8126)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:5531)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _eval_(self, *args)
493 if not is_Expression(args[-1]):
494 try:
--> 495 return self._evalf_(*args)
496 except AttributeError:
497 pass
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _evalf_(self, *args, **kwds)
641 precision = step_parent.prec()
642 except AttributeError:
--> 643 precision = RR.prec()
644
645 from sage.libs.mpmath.all import call as mpcall
NameError: global name 'RR' is not defined
</pre>
TicketJeroen DemeyerThu, 05 Dec 2013 20:00:04 GMT
https://trac.sagemath.org/ticket/9706#comment:84
https://trac.sagemath.org/ticket/9706#comment:84
<p>
Why do we have <code>__call__</code> and <code>_eval_</code> given they do essentially the same thing? Please keep in mind <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:69" title="Comment 69">69</a> also in the <code>_eval_()</code> function (I really don't understand why it needs to be so complicated).
</p>
<p>
This is also still broken:
</p>
<pre class="wiki">sage: parent(chebyshev_T(4, RIF(5)))
Real Field with 53 bits of precision
</pre><p>
Please fix all issues that I mentioned <strong>and turn them into doctests</strong>.
</p>
TicketJeroen DemeyerThu, 05 Dec 2013 22:35:59 GMT
https://trac.sagemath.org/ticket/9706#comment:85
https://trac.sagemath.org/ticket/9706#comment:85
<p>
For the <code>REFERENCES</code>, see <a class="ext-link" href="http://sagemath.org/doc/developer/conventions.html#docstring-markup-with-rest-and-sphinx"><span class="icon"></span>http://sagemath.org/doc/developer/conventions.html#docstring-markup-with-rest-and-sphinx</a> for the right syntax.
</p>
TicketStefan ReitererThu, 05 Dec 2013 23:14:40 GMT
https://trac.sagemath.org/ticket/9706#comment:86
https://trac.sagemath.org/ticket/9706#comment:86
<p>
@clenshaw_method: there is a difference. clenshaw method also applies a direct formula for small n and calls the recursive method else. The difference is that the recursive evaluation does not give an expanded representation of the polynomial, which is wanted for small n, because that was the standard till now and people expect this, especially if you are used to mathematica or maxima. Expanding huge expressions costs a lot of time, and this approach is much faster in that situation. Of course it is a matter of naming. But the reason why I have 2 methods, is to avoid too long code segments, and splitting them apart is better for readability. It also is important concerning other orthogonal polynomials.
</p>
<p>
@<span class="underline">call</span> & _eval_ : This convention is part of the <a class="missing wiki">BuiltinFunction?</a> structure.
<span class="underline">call</span> does all the stuff like coercions, transforming into a symbolic expression (e.g. if n is a symbolic value don't return a polynomial but hold the closed form.)
_eval cares about evaluating the polynomial (e.g return a number if x is a number etc.)
Look into the symbolic.function module for more details
And eval is so complicated because there are several cases to consider: correct evaluation of symbolic expressions, numerical expressions and numpy arrays etc.
This is also part of the <a class="missing wiki">BuiltinFunction?</a> structure. And you also have to keep in mind that this method should work for all ortho polys.
</p>
<p>
@bugs Sorry, during the patch merging process I had forgotten to apply a patch, which I'm now missing, since I work on different machines ... I will correct this tomorrow.
It's annoying since I already had fixed it ...
</p>
TicketStefan ReitererThu, 05 Dec 2013 23:55:09 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_bugfixes_and_doctests.patch</em>
</li>
</ul>
<p>
incorporated things that were already been done ...
</p>
TicketStefan ReitererFri, 06 Dec 2013 00:00:16 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:87
https://trac.sagemath.org/ticket/9706#comment:87
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=87">diff</a>)
</li>
</ul>
<p>
Ok I incorporated the bugfixes and doctests <strong>again</strong> ...
</p>
<p>
there are still some minor changes (formatting and new doctests) todo. Please let me know if you find more bugs.
</p>
TicketJeroen DemeyerFri, 06 Dec 2013 10:11:44 GMT
https://trac.sagemath.org/ticket/9706#comment:88
https://trac.sagemath.org/ticket/9706#comment:88
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:86" title="Comment 86">maldun</a>:
</p>
<blockquote class="citation">
<p>
@clenshaw_method: there is a difference. clenshaw method also applies a direct formula for small n and calls the recursive method else. The difference is that the recursive evaluation does not give an expanded representation of the polynomial, which is wanted for small n
</p>
</blockquote>
<p>
OK, fine. But for simplicity, you could simply call <code>_cheb_recur_(...).expand()</code> instead which would achieve the same thing without an additional method.
</p>
<blockquote class="citation">
<p>
_eval cares about evaluating the polynomial (e.g return a number if x is a number etc.)
Look into the symbolic.function module for more details
And eval is so complicated because there are several cases to consider: correct evaluation of symbolic expressions, numerical expressions and numpy arrays etc.
</p>
</blockquote>
<p>
If you really think the complexity is justified (I have a hard time believing that), you should add comments in the code to describe the various cases, because I'm having a hard time understanding the logic. A comment like <code># A faster check would be nice...</code> doesn't mean much to me because I don't understand what you're trying to do.
</p>
TicketJeroen DemeyerFri, 06 Dec 2013 10:14:52 GMT
https://trac.sagemath.org/ticket/9706#comment:89
https://trac.sagemath.org/ticket/9706#comment:89
<p>
I think <code>mpmath</code> should only be used for the "pure" <code>RealField()</code> and <code>ComplexField()</code> and <code>RDF</code> and <code>CDF</code>, nothing else.
</p>
<p>
This is bad:
</p>
<pre class="wiki">sage: chebyshev_T(5,Qp(3)(2))
...
TypeError: unable to coerce to a ComplexNumber: <type 'sage.rings.padics.padic_capped_relative_element.pAdicCappedRelativeElement'>
</pre><p>
and the way you use <code>RIF</code> is also kind of stupid since the computation should really be done using the recursive formula.
</p>
TicketJeroen DemeyerFri, 06 Dec 2013 10:17:52 GMT
https://trac.sagemath.org/ticket/9706#comment:90
https://trac.sagemath.org/ticket/9706#comment:90
<p>
Also, you should never simply <code>print</code> stuff. Either delete those prints or use a <a class="ext-link" href="http://docs.python.org/2/library/warnings.html#available-functions"><span class="icon"></span>Python warning</a>.
</p>
<pre class="wiki">sage: chebyshev_T(10^6,RealField(256)(2))
Warning: mpmath returns NoConvergence!
Switching to clenshaw_method, but it may not be stable!
1.764019427245793968639371137094247875315949668035854027376792158135922267593e571947
</pre><p>
The message is also misleading, since for some inputs it retries <code>mpmath</code> anyway:
</p>
<pre class="wiki">sage: chebyshev_T(100001/2, 2)
---------------------------------------------------------------------------
NoConvergence Traceback (most recent call last)
<ipython-input-34-9c95a5ff4276> in <module>()
----> 1 chebyshev_T(Integer(100001)/Integer(2), Integer(2))
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/ort
557 return self._eval_(n,x) # Let eval methode decide which is best
558 else: # Consider OrthogonalPolynomial as symbol
--> 559 return super(OrthogonalPolynomial,self).__call__(n,x)
560
561
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFp:8126)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orth
496 if not is_Expression(args[-1]):
497 try:
--> 498 return self._evalf_(*args)
499 except AttributeError:
500 pass
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.pyc in _evalf_(self, *args,
651 from sage.libs.mpmath.all import chebyt as mpchebyt
652
--> 653 return mpcall(mpchebyt,args[0],args[-1],prec = precision)
654
655 def _maxima_init_evaled_(self, *args):
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/libs/mpmath/utils.so in sage.libs.mpmath.utils.call (sage/libs/
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/libs/mpmath/ext_main.so in sage.libs.mpmath.ext_main.wrapped_specfun.__call__ (sage/lit_main.c:17447)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/mpmath/functions/orthogonal.pyc in chebyt(ctx, n, x, **kwarg
444 if (not x) and ctx.isint(n) and int(ctx._re(n)) % 2 == 1:
445 return x * 0
--> 446 return ctx.hyp2f1(-n,n,(1,2),(1-x)/2, **kwargs)
447
448 @defun_wrapped
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/mpmath/functions/hypergeometric.pyc in hyp2f1(ctx, a, b, c, z, **kwargs)
248 @defun
249 def hyp2f1(ctx,a,b,c,z,**kwargs):
--> 250 return ctx.hyper([a,b],[c],z,**kwargs)
251
252 @defun
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/mpmath/functions/hypergeometric.pyc in hyper(ctx, a_s, b_s, z, **kwargs)
223 elif q == 0: return ctx._hyp1f0(a_s[0][0], z)
224 elif p == 2:
--> 225 if q == 1: return ctx._hyp2f1(a_s, b_s, z, **kwargs)
226 elif q == 2: return ctx._hyp2f2(a_s, b_s, z, **kwargs)
227 elif q == 3: return ctx._hyp2f3(a_s, b_s, z, **kwargs)
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/mpmath/functions/hypergeometric.pyc in _hyp2f1(ctx, a_s, b_s, z, **kwargs)
442 if absz <= 0.8 or (ctx.isint(a) and a <= 0 and a >= -1000) or \
443 (ctx.isint(b) and b <= 0 and b >= -1000):
--> 444 return ctx.hypsum(2, 1, (atype, btype, ctype), [a, b, c], z, **kwargs)
445
446 orig = ctx.prec
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/mpmath/ctx_mp.pyc in hypsum(ctx, p, q, flags, coeffs, z, accurate_small, **kwargs)
640 mag_dict = {}
641 zv, have_complex, magnitude = summator(coeffs, v, prec, wp, \
--> 642 epsshift, mag_dict, **kwargs)
643 cancel = -magnitude
644 jumps_resolved = True
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/mpmath/libmp/libhyper.pyc in _hypsum(coeffs, z, prec, wp, epsshift, magnitude_check, **kwargs)
319 def _hypsum(coeffs, z, prec, wp, epsshift, magnitude_check, **kwargs):
320 return hypsum_internal(p, q, param_types, ztype, coeffs, z,
--> 321 prec, wp, epsshift, magnitude_check, kwargs)
322
323 return "(none)", _hypsum
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/libs/mpmath/ext_main.so in sage.libs.mpmath.ext_main.hypsum_internal (sage/libs/mpmath/ext_main.c:24850)()
/usr/local/src/sage-5.13.beta1/local/lib/python2.7/site-packages/sage/libs/mpmath/ext_impl.so in sage.libs.mpmath.ext_impl.MPF_hypsum (sage/libs/mpmath/ext_impl.c:20614)()
NoConvergence: Hypergeometric series converges too slowly. Try increasing maxterms.
</pre>
TicketJeroen DemeyerFri, 06 Dec 2013 10:24:09 GMT
https://trac.sagemath.org/ticket/9706#comment:91
https://trac.sagemath.org/ticket/9706#comment:91
<p>
Don't forget to fix the <code>REFERENCES</code>.
</p>
TicketJeroen DemeyerFri, 06 Dec 2013 13:27:04 GMT
https://trac.sagemath.org/ticket/9706#comment:92
https://trac.sagemath.org/ticket/9706#comment:92
<pre class="wiki">sage: chebyshev_U(-1, Mod(2,5))
...
RuntimeError: maximum recursion depth exceeded
</pre>
TicketJeroen DemeyerFri, 06 Dec 2013 13:28:40 GMT
https://trac.sagemath.org/ticket/9706#comment:93
https://trac.sagemath.org/ticket/9706#comment:93
<p>
This is totally meaningless (this should be a <code>ValueError</code>):
</p>
<pre class="wiki">sage: chebyshev_U(1.5, Mod(8,9))
63.8753125671848
</pre>
TicketJeroen DemeyerFri, 06 Dec 2013 13:42:35 GMT
https://trac.sagemath.org/ticket/9706#comment:94
https://trac.sagemath.org/ticket/9706#comment:94
<p>
<code>mpmath</code> is slower than <code>_cheb_recur_()</code>, so I think <code>mpmath</code> should be only used in cases where the index is not an integer.
</p>
<p>
Perhaps the <code>_eval_</code> logic should be
</p>
<ol><li><code>n</code> in <code>ZZ</code> => use recursion. If <code>x</code> is symbolic and <code>abs(n) <= 50</code>, <code>expand()</code> the result
</li><li><code>n</code> or <code>x</code> symbolic => use symbolic evaluation
</li><li><code>n</code> and <code>x</code> real/complex => use <code>mpmath</code>
</li><li>anything else => <code>raise ValueError</code>
</li></ol>
TicketJeroen DemeyerFri, 06 Dec 2013 13:46:24 GMT
https://trac.sagemath.org/ticket/9706#comment:95
https://trac.sagemath.org/ticket/9706#comment:95
<p>
One should be able to convert to maxima (otherwise, what's the point of introducing symbolic Chebyshev polynomials):
</p>
<pre class="wiki">sage: var('n,x')
(n, x)
sage: maxima( chebyshev_T(n, cos(x)) )
...
TypeError: unable to convert x (=n) to an integer
</pre>
TicketStefan ReitererSun, 08 Dec 2013 15:13:42 GMT
https://trac.sagemath.org/ticket/9706#comment:96
https://trac.sagemath.org/ticket/9706#comment:96
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:93" title="Comment 93">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
This is totally meaningless (this should be a <code>ValueError</code>):
</p>
<pre class="wiki">sage: chebyshev_U(1.5, Mod(8,9))
63.8753125671848
</pre></blockquote>
<p>
Indeed, but it's not my fault. It appears, that the <a class="missing wiki">BuiltinFunction?</a> calls the _eval_numpy_ method. E.g.
</p>
<pre class="wiki">sage: csc(Mod(8,9))
1.01075621840010
</pre><p>
makes no sense either but works. Maybe we should open a ticket on this?
</p>
<p>
I have worked out now a new patch, where _eval_ returns None, like expected, but still returns this numerical value. But it should be patched in the Symbolic function classes and not here.
</p>
TicketStefan ReitererSun, 08 Dec 2013 15:20:11 GMT
https://trac.sagemath.org/ticket/9706#comment:97
https://trac.sagemath.org/ticket/9706#comment:97
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:95" title="Comment 95">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
One should be able to convert to maxima (otherwise, what's the point of introducing symbolic Chebyshev polynomials):
</p>
<pre class="wiki">sage: var('n,x')
(n, x)
sage: maxima( chebyshev_T(n, cos(x)) )
...
TypeError: unable to convert x (=n) to an integer
</pre></blockquote>
<p>
One of the reasons I started this rewriting buisness, was the fact, that Maxima is very limited in respect to Orthogonal polynomials. And if we consider advanced use for symbolic methods like 'Creative Telescoping' it makes perfect sense to allow a non
Maxima conform representation, since Sage symbolic capabilities don't relie on Maxima allone.
</p>
TicketJeroen DemeyerSun, 08 Dec 2013 17:17:05 GMT
https://trac.sagemath.org/ticket/9706#comment:98
https://trac.sagemath.org/ticket/9706#comment:98
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:97" title="Comment 97">maldun</a>:
</p>
<blockquote class="citation">
<p>
One of the reasons I started this rewriting buisness, was the fact, that Maxima is very limited in respect to Orthogonal polynomials. And if we consider advanced use for symbolic methods like 'Creative Telescoping' it makes perfect sense to allow a non
Maxima conform representation, since Sage symbolic capabilities don't relie on Maxima allone.
</p>
</blockquote>
<p>
This looks like a poor excuse to me, since Maxima can deal with Chebyshev polynomials just fine. The following works:
</p>
<pre class="wiki">sage: maxima('chebysheb_t(n,x)')
</pre><p>
The fact that the conversion to Maxima doesn't work is a fault of your patch, don't blame Maxima for that.
</p>
TicketJeroen DemeyerSun, 08 Dec 2013 17:17:45 GMT
https://trac.sagemath.org/ticket/9706#comment:99
https://trac.sagemath.org/ticket/9706#comment:99
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:96" title="Comment 96">maldun</a>:
</p>
<blockquote class="citation">
<pre class="wiki">sage: csc(Mod(8,9))
1.01075621840010
</pre><p>
makes no sense either but works. Maybe we should open a ticket on this?
</p>
</blockquote>
<p>
I agree, this is a problem.
</p>
TicketStefan ReitererSun, 08 Dec 2013 18:14:30 GMT
https://trac.sagemath.org/ticket/9706#comment:100
https://trac.sagemath.org/ticket/9706#comment:100
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:98" title="Comment 98">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:97" title="Comment 97">maldun</a>:
</p>
<blockquote class="citation">
<p>
One of the reasons I started this rewriting buisness, was the fact, that Maxima is very limited in respect to Orthogonal polynomials. And if we consider advanced use for symbolic methods like 'Creative Telescoping' it makes perfect sense to allow a non
Maxima conform representation, since Sage symbolic capabilities don't relie on Maxima allone.
</p>
</blockquote>
<p>
This looks like a poor excuse to me, since Maxima can deal with Chebyshev polynomials just fine. The following works:
</p>
<pre class="wiki">sage: maxima('chebysheb_t(n,x)')
</pre><p>
The fact that the conversion to Maxima doesn't work is a fault of your patch, don't blame Maxima for that.
</p>
</blockquote>
<p>
I'm not blaming maxima, but this never worked on the main branch in the first place:
</p>
<pre class="wiki">maxima( chebyshev_T(n, cos(x)) )
...
TypeError: unable to convert x (=n) to an integer
</pre><p>
so I did not break anything ...
</p>
<p>
You can argue that this is an open improvement, but it's definitely no regression.
Nevertheless it will be corrected in the next patch.
</p>
TicketStefan ReitererSun, 08 Dec 2013 21:04:46 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9707_chebyshev_corrections.patch</em>
</li>
</ul>
<p>
Bunch of corrections and new doctests
</p>
TicketStefan ReitererSun, 08 Dec 2013 21:06:41 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:101
https://trac.sagemath.org/ticket/9706#comment:101
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=101">diff</a>)
</li>
</ul>
<p>
I incorporated now the things you wished for, except things for ducumentation.
I also added many new doctests.
</p>
<p>
I also cleaned up the evaluation methods.
</p>
TicketStefan ReitererSun, 08 Dec 2013 22:22:18 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_chebyshev_corrections.patch</em>
</li>
</ul>
<p>
corrected name of file
</p>
TicketStefan ReitererSun, 08 Dec 2013 22:22:39 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_chebyshev_docu.patch</em>
</li>
</ul>
<p>
docu changes
</p>
TicketStefan ReitererSun, 08 Dec 2013 22:22:43 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_chebyshev_docu.2.patch</em>
</li>
</ul>
<p>
docu changes
</p>
TicketStefan ReitererSun, 08 Dec 2013 22:24:27 GMTstatus, description changed
https://trac.sagemath.org/ticket/9706#comment:102
https://trac.sagemath.org/ticket/9706#comment:102
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=102">diff</a>)
</li>
</ul>
<p>
added changes in docu.
</p>
<p>
Still one bug to fix
</p>
TicketStefan ReitererSun, 08 Dec 2013 22:27:03 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:103
https://trac.sagemath.org/ticket/9706#comment:103
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
TicketStefan ReitererSun, 08 Dec 2013 22:41:00 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_chebyshev_U-1.patch</em>
</li>
</ul>
<p>
fixed bugs with chebyshev_U(-1,...)
</p>
TicketStefan ReitererSun, 08 Dec 2013 22:41:58 GMTstatus, description changed
https://trac.sagemath.org/ticket/9706#comment:104
https://trac.sagemath.org/ticket/9706#comment:104
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=104">diff</a>)
</li>
</ul>
<p>
Finally fixed the chebyshev_U(-1,...) issue
</p>
TicketStefan ReitererSun, 08 Dec 2013 23:06:11 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_chebyshev_collection.patch</em>
</li>
</ul>
<p>
collection of all patches
</p>
TicketStefan ReitererSun, 08 Dec 2013 23:06:46 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:105
https://trac.sagemath.org/ticket/9706#comment:105
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=105">diff</a>)
</li>
</ul>
TicketJeroen DemeyerMon, 09 Dec 2013 09:55:44 GMTreviewer, description changed
https://trac.sagemath.org/ticket/9706#comment:106
https://trac.sagemath.org/ticket/9706#comment:106
<ul>
<li><strong>reviewer</strong>
changed from <em>Burcin Erocal, Travis Scrimshaw, Stefan Reiterer</em> to <em>Burcin Erocal, Travis Scrimshaw, Stefan Reiterer, Jeroen Demeyer</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=106">diff</a>)
</li>
</ul>
TicketJeroen DemeyerMon, 09 Dec 2013 10:59:42 GMT
https://trac.sagemath.org/ticket/9706#comment:107
https://trac.sagemath.org/ticket/9706#comment:107
<p>
Working on reviewer patch...
</p>
TicketJeroen DemeyerMon, 09 Dec 2013 14:31:06 GMTpriority, description changed; dependencies set
https://trac.sagemath.org/ticket/9706#comment:108
https://trac.sagemath.org/ticket/9706#comment:108
<ul>
<li><strong>priority</strong>
changed from <em>minor</em> to <em>major</em>
</li>
<li><strong>dependencies</strong>
set to <em>#864, #9640, #10018, #11868, #15422</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=108">diff</a>)
</li>
</ul>
TicketJeroen DemeyerMon, 09 Dec 2013 14:34:31 GMT
https://trac.sagemath.org/ticket/9706#comment:109
https://trac.sagemath.org/ticket/9706#comment:109
<p>
Review patch, mainly simplifies the code: less special cases and use actual arguments instead of <code>args</code> and <code>kwds</code>. Also add some more doctests.
</p>
TicketStefan ReitererMon, 09 Dec 2013 17:09:37 GMT
https://trac.sagemath.org/ticket/9706#comment:110
https://trac.sagemath.org/ticket/9706#comment:110
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:109" title="Comment 109">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Review patch, mainly simplifies the code: less special cases and use actual arguments instead of <code>args</code> and <code>kwds</code>. Also add some more doctests.
</p>
</blockquote>
<p>
Is it really a good idea to replace *args and <strong>kwds in the <a class="missing wiki">OrthogonalPolynomial?</a> class?
</strong></p>
<p>
Since not all ortho polys only have 2 arguments, e.g. Gegenbauer polynomials <a class="ext-link" href="http://en.wikipedia.org/wiki/Gegenbauer_polynomials"><span class="icon"></span>http://en.wikipedia.org/wiki/Gegenbauer_polynomials</a> which have an additional parameter alpha
It makes perfect sense for the chebyshev polys though
</p>
<p>
I know that the base class looks strange from the point of the Chebyshev polys,
but there is some reasoning behind this =)
</p>
TicketJeroen DemeyerMon, 09 Dec 2013 18:41:31 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>9706_review.patch</em>
</li>
</ul>
TicketJeroen DemeyerMon, 09 Dec 2013 18:43:51 GMT
https://trac.sagemath.org/ticket/9706#comment:111
https://trac.sagemath.org/ticket/9706#comment:111
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:110" title="Comment 110">maldun</a>:
</p>
<blockquote class="citation">
<p>
Is it really a good idea to replace *args and <strong>kwds in the <a class="missing wiki">OrthogonalPolynomial?</a> class?
</strong></p>
<p>
Since not all ortho polys only have 2 arguments, e.g. Gegenbauer polynomials <a class="ext-link" href="http://en.wikipedia.org/wiki/Gegenbauer_polynomials"><span class="icon"></span>http://en.wikipedia.org/wiki/Gegenbauer_polynomials</a> which have an additional parameter alpha
It makes perfect sense for the chebyshev polys though
</p>
</blockquote>
<p>
Ok, it's always difficult to design an interface for something you don't have an implementation for, but I made the change such that the Gegenbauer polynomials should (in theory) work.
</p>
TicketStefan ReitererMon, 09 Dec 2013 20:10:55 GMT
https://trac.sagemath.org/ticket/9706#comment:112
https://trac.sagemath.org/ticket/9706#comment:112
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:111" title="Comment 111">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:110" title="Comment 110">maldun</a>:
</p>
<blockquote class="citation">
<p>
Is it really a good idea to replace *args and <strong>kwds in the <a class="missing wiki">OrthogonalPolynomial?</a> class?
</strong></p>
<p>
Since not all ortho polys only have 2 arguments, e.g. Gegenbauer polynomials <a class="ext-link" href="http://en.wikipedia.org/wiki/Gegenbauer_polynomials"><span class="icon"></span>http://en.wikipedia.org/wiki/Gegenbauer_polynomials</a> which have an additional parameter alpha
It makes perfect sense for the chebyshev polys though
</p>
</blockquote>
<p>
Ok, it's always difficult to design an interface for something you don't have an implementation for, but I made the change such that the Gegenbauer polynomials should (in theory) work.
</p>
</blockquote>
<p>
Maybe you should have a short look ino this patch: <a class="ext-link" href="http://trac.sagemath.org/attachment/ticket/9706/trac_9706_orthogonal_polys.patch"><span class="icon"></span>http://trac.sagemath.org/attachment/ticket/9706/trac_9706_orthogonal_polys.patch</a>
it contains already prototypes of the most ortho polys.
</p>
<p>
I have some remarks concerning your patch:
</p>
<pre class="wiki">sage: chebyshev_T(n,Mod(0,8))
1/2*(-1)^(1/2*n)*((-1)^n + 1)
</pre><p>
but this makes no sense since 1/2 is not defined in Mod(8). This was the reason for the 0 in CC check at this point.
</p>
<p>
You evaluate numerical expressions for n in NN with recursion. this is favorable for chebyshev polynomials, but not for all ortho polys you can evaluate numeric values in O(log n). You have already problems with the legendre polynomials, since the coefficients depend on n, and the recursion is not stable. Thus other evaluation methods should be used. Thats the reason why _evalf_ with mpmath should come first.
In case of chebyshev I catched this with an explicit call of the recursion in _evalf_.
</p>
<p>
The _old_maxima_ method is used for some oddballs, where the only useful implementation is in maxima, and for some special cases. So removing is probably not a good idea.
</p>
<p>
The reason why negative values are checked in special values, is that in later versions, or for other polys non integral negative values can be treated analogously, as in the case of their algebraic counterpart. E.g negative legendre polynomials would return associated legendre functions, but have no algebraic sense.
</p>
<p>
Another reason to allow the special values check for non symbolic input, is that e.g. (-1)<sup>100000000 is evaluated faster, than applying the recursion, or other special points.
I use this feature sometimes to evaluate certain polynomials at the end points of an intervall e.g. [-1,1]
</sup></p>
<p>
<strong>Question:</strong> Since the Legendre Polynomials already conflicts with the design in the ortho poly class, should we add legendre_P too, to get a better overview, what actually should be kept in the orthogonal poly class? then we could also add more reasonable doctests, to understand the evaluation methods better.
</p>
<p>
I have to admit, that about 3 years have passed since I wrote the initial version, and I'm not remembering why I took some design desicisions back then, but slowly my memories are coming back.
</p>
TicketJeroen DemeyerMon, 09 Dec 2013 20:26:11 GMT
https://trac.sagemath.org/ticket/9706#comment:113
https://trac.sagemath.org/ticket/9706#comment:113
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:112" title="Comment 112">maldun</a>:
</p>
<blockquote class="citation">
<p>
<strong>Question:</strong> Since the Legendre Polynomials already conflicts with the design in the ortho poly class, should we add legendre_P too, to get a better overview, what actually should be kept in the orthogonal poly class?
</p>
</blockquote>
<p>
Absolutely not. The patch is already big enough now.
</p>
TicketJeroen DemeyerMon, 09 Dec 2013 20:30:04 GMT
https://trac.sagemath.org/ticket/9706#comment:114
https://trac.sagemath.org/ticket/9706#comment:114
<p>
If the various orthogonal polynomials are so different, then perhaps the simple answer is that we shouldn't try to force a generic <code>_eval_</code> which will work for all orthogonal polynomials.
</p>
<p>
We could have a common superclass for both kinds of Chebyshev polynomials and implement <code>_eval_()</code> there. For Legendre polynomials, we could implement a different <code>_eval_()</code>. That would be a much simpler solution than making an overly complicated generic <code>_eval_()</code>.
</p>
TicketTravis ScrimshawMon, 09 Dec 2013 20:47:52 GMTauthor changed
https://trac.sagemath.org/ticket/9706#comment:115
https://trac.sagemath.org/ticket/9706#comment:115
<ul>
<li><strong>author</strong>
changed from <em>Stefan Reiterer, Travis Scrimshaw</em> to <em>Stefan Reiterer</em>
</li>
</ul>
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:114" title="Comment 114">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
If the various orthogonal polynomials are so different, then perhaps the simple answer is that we shouldn't try to force a generic <code>_eval_</code> which will work for all orthogonal polynomials.
</p>
<p>
We could have a common superclass for both kinds of Chebyshev polynomials and implement <code>_eval_()</code> there. For Legendre polynomials, we could implement a different <code>_eval_()</code>. That would be a much simpler solution than making an overly complicated generic <code>_eval_()</code>.
</p>
</blockquote>
<p>
I think this would be a good course of action, and that we should put other orthogonal polynomials in other tickets. However, I think it might be worthwhile to at least diagram/pseudocode the overall class structure we want at the end of the day. Currently I believe the proposal is something like:
</p>
<pre class="wiki">* Orthogonal polynomials
* Chebyshev polynomials
- general _evel_(x, n) method
* Chebyshev T
- specific code (ex. _evalf_() method), recursions, ...
* Chebyshev U
- specific code (ex. _evalf_() method), recursions, ...
* Legendre polynomials
- general _evel_(x, n) method
* Legendre P
* Legendre Q
* Gegenbauer polynomials
- an _evel_(x, n, alpha) method
</pre><p>
Hopefully my notation is clear
</p>
TicketJeroen DemeyerMon, 09 Dec 2013 21:04:25 GMT
https://trac.sagemath.org/ticket/9706#comment:116
https://trac.sagemath.org/ticket/9706#comment:116
<p>
Something like that looks right indeed. maldun: what do you think?
</p>
<p>
Perhaps the only code so far that could be truly generic is the <code>__call__</code> method.
</p>
TicketJeroen DemeyerMon, 09 Dec 2013 21:04:34 GMT
https://trac.sagemath.org/ticket/9706#comment:117
https://trac.sagemath.org/ticket/9706#comment:117
<p>
Something like that looks right indeed. maldun: what do you think?
</p>
<p>
Perhaps the only code so far that could be truly generic is the <code>__call__</code> method.
</p>
TicketStefan ReitererMon, 09 Dec 2013 21:21:56 GMT
https://trac.sagemath.org/ticket/9706#comment:118
https://trac.sagemath.org/ticket/9706#comment:118
<p>
after some thinking, I guess you are right. A General <a class="missing wiki">OrthogonalPolynomial?</a> is sophisticated, but it needs too much tweaking, and too much exception cathing, which makes the code unsafe.
</p>
<p>
@ <span class="underline">call</span> : 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.
</p>
<p>
I propose the following
</p>
<ul><li><a class="missing wiki">OrthogonalPolynomial?</a>: Naming Conventions and <span class="underline">init</span> method
</li><li><a class="missing wiki">ChebyshevPolynomial?</a>: base Class for Cheby_t and Cheby_u (Current Orthogonal Polynomials)
</li><li><a class="missing wiki">LegendrePolynomials?</a>
</li></ul><p>
.... etc.
</p>
TicketStefan ReitererMon, 09 Dec 2013 22:02:25 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706_new_base_classes.patch</em>
</li>
</ul>
<p>
proposed new structure
</p>
TicketStefan ReitererWed, 11 Dec 2013 06:01:54 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:119
https://trac.sagemath.org/ticket/9706#comment:119
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=119">diff</a>)
</li>
</ul>
<p>
I added an intermediate class between <a class="missing wiki">OrthogonalPolynomials?</a> and the Chebyshev polynomals namely <a class="missing wiki">ChebyshevPolynomial?</a>.
</p>
<p>
Any comments on that?
</p>
TicketTravis ScrimshawWed, 11 Dec 2013 06:38:40 GMT
https://trac.sagemath.org/ticket/9706#comment:120
https://trac.sagemath.org/ticket/9706#comment:120
<p>
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).
</p>
TicketJeroen DemeyerWed, 11 Dec 2013 07:22:50 GMT
https://trac.sagemath.org/ticket/9706#comment:121
https://trac.sagemath.org/ticket/9706#comment:121
<p>
I am currently working on finishing the Sage 5.13 release, so I don't feel like reviewing this now.
</p>
TicketTravis ScrimshawWed, 11 Dec 2013 22:57:33 GMTdescription changed
https://trac.sagemath.org/ticket/9706#comment:122
https://trac.sagemath.org/ticket/9706#comment:122
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/9706?action=diff&version=122">diff</a>)
</li>
</ul>
<p>
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.
</p>
TicketStefan ReitererThu, 12 Dec 2013 19:28:30 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:123
https://trac.sagemath.org/ticket/9706#comment:123
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
</ul>
<p>
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 ...)
</p>
TicketNils BruinThu, 12 Dec 2013 21:04:05 GMT
https://trac.sagemath.org/ticket/9706#comment:124
https://trac.sagemath.org/ticket/9706#comment:124
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/9706#comment:123" title="Comment 123">maldun</a>:
</p>
<blockquote class="citation">
<p>
And good idea to replace the old % format string operator, since it will be deprecated in Python 3
</p>
</blockquote>
<p>
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
<a class="ext-link" href="http://www.python.org/dev/peps/pep-3101/"><span class="icon"></span>http://www.python.org/dev/peps/pep-3101/</a>
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.
</p>
TicketJeroen DemeyerMon, 16 Dec 2013 09:30:48 GMT
https://trac.sagemath.org/ticket/9706#comment:125
https://trac.sagemath.org/ticket/9706#comment:125
<p>
Travis: are you sure you applied the dependency <a class="closed ticket" href="https://trac.sagemath.org/ticket/15422" title="#15422: defect: factorization of non-squarefree polynomials over the p-adics (closed: fixed)">#15422</a>? I am suspicious because you added
</p>
<div class="wiki-code"><div class="code"><pre><span class="gd">- (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))
</span><span class="gi">+ (2^7 + O(2^8))*t^5 + (2^6 + O(2^8))*t^3 + (1 + 2^6 + O(2^8))*t
</span></pre></div></div><p>
which seems to be a doctest failure.
</p>
TicketJeroen DemeyerMon, 16 Dec 2013 09:44:08 GMT
https://trac.sagemath.org/ticket/9706#comment:126
https://trac.sagemath.org/ticket/9706#comment:126
<p>
I also don't agree with all changes of
</p>
<pre class="wiki">if condition:
return A
else:
return B
</pre><p>
to
</p>
<pre class="wiki">if condition:
return A
return B
</pre><p>
but I guess that's one of your pet peeves.
</p>
<p>
I personally prefer
</p>
<pre class="wiki">if condition:
return A
else:
return B
</pre><p>
if there are clearly two cases and the code might as well have been written as
</p>
<pre class="wiki">if not condition:
return B
else:
return A
</pre><p>
So I personally would keep the <code>if/else</code> structure for the <code>if n % 2 == 0</code> test. And for the <code>n >= 0</code> test, I would say that
</p>
<pre class="wiki">if n < 0:
return B
return A
</pre><p>
would be a lot better that what you did:
</p>
<pre class="wiki">if n >= 0:
return A
return B
</pre><p>
The first feels much better to me because you put the normal case outside the <code>if</code> block and the special cases <code>if n == 0</code> and <code>if n < 0</code> would be inside <code>if</code>s.
</p>
<p>
(Of course these are all personal preferences, I'm not asking you to change this, maybe just think about it.)
</p>
TicketJeroen DemeyerMon, 16 Dec 2013 11:01:06 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:127
https://trac.sagemath.org/ticket/9706#comment:127
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>needs_info</em>
</li>
</ul>
TicketTravis ScrimshawMon, 16 Dec 2013 15:07:19 GMTattachment set
https://trac.sagemath.org/ticket/9706
https://trac.sagemath.org/ticket/9706
<ul>
<li><strong>attachment</strong>
set to <em>trac_9706-review-ts.patch</em>
</li>
</ul>
TicketTravis ScrimshawMon, 16 Dec 2013 15:15:04 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:128
https://trac.sagemath.org/ticket/9706#comment:128
<ul>
<li><strong>status</strong>
changed from <em>needs_info</em> to <em>needs_review</em>
</li>
</ul>
<p>
I had missed the <a class="closed ticket" href="https://trac.sagemath.org/ticket/15422" title="#15422: defect: factorization of non-squarefree polynomials over the p-adics (closed: fixed)">#15422</a> dependency. Fixed.
</p>
<p>
As for the <code>if/else</code>, I've seen (sometimes long) <code>else</code> 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 <code>else:</code>. I added back in the <code>else</code> blocks for the <code>if n % 2 == 0</code>.
</p>
TicketJeroen DemeyerMon, 16 Dec 2013 15:46:38 GMTstatus changed
https://trac.sagemath.org/ticket/9706#comment:129
https://trac.sagemath.org/ticket/9706#comment:129
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
</ul>
TicketStefan ReitererTue, 17 Dec 2013 12:01:13 GMT
https://trac.sagemath.org/ticket/9706#comment:130
https://trac.sagemath.org/ticket/9706#comment:130
<p>
First stage done. Next step: Legendre Polys.
</p>
<p>
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.
</p>
<p>
Thank all reviewers for their hard work and the good input!
</p>
TicketVolker BraunSat, 21 Dec 2013 14:30:42 GMTcommit, branch set
https://trac.sagemath.org/ticket/9706#comment:131
https://trac.sagemath.org/ticket/9706#comment:131
<ul>
<li><strong>commit</strong>
set to <em>bb227b60e93b7d8806b32b00874351b9201123cd</em>
</li>
<li><strong>branch</strong>
set to <em>u/vbraun/chebyshev</em>
</li>
</ul>
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=bb227b6"><span class="icon"></span>bb227b6</a></td><td><code>#9706: review patch.</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=51fcd6b"><span class="icon"></span>51fcd6b</a></td><td><code>trac 9706: Propose new class structure</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=4677a15"><span class="icon"></span>4677a15</a></td><td><code>Symbolic Chebyshev polynomials: reviewer patch</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=75d1e43"><span class="icon"></span>75d1e43</a></td><td><code>trac 9706: Collective patch. Bugfixes, extensions, optimizations, documentation, doctests for chebyshev_T, chebyshev_U and base class for ortho polys</code>
</td></tr></table>
TicketVolker BraunSat, 21 Dec 2013 18:43:54 GMTstatus changed; resolution set
https://trac.sagemath.org/ticket/9706#comment:132
https://trac.sagemath.org/ticket/9706#comment:132
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>closed</em>
</li>
<li><strong>resolution</strong>
set to <em>fixed</em>
</li>
</ul>
Ticket