Opened 11 years ago

# New Version of orthogonal Polynomials (Part I: Base Class and Chebyshev Polynomials) — at Version 74

Reported by: Owned by: maldun burcin, maldun major sage-6.1 symbolics orthogonal polynomials, symbolics fredrik.johansson, fstan, kcrisman Stefan Reiterer, Travis Scrimshaw Burcin Erocal, Travis Scrimshaw, Stefan Reiterer N/A

The current implementation of orthogonal polynomials is just a wrapper around maxima. (see http://wiki.sagemath.org/symbolics/) This update holds the following changes:

• using of the pynac class for symbolic functions.
• faster evaluation in general
• evaluation of special values
• mpmath for numeric evaluation

Remarks:

• The current patch needs scipy-0.8. One has to install it before testing (see #9808 for spkg's and installation instructions)
• Some of the old doctests in the old file don't work any more, due to coercion problems with pynac (see #9769)
• Some doctests in Sage change, due to the fact that new BuiltIn? 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).

Apply:

### Changed 11 years ago by maldun

A new version of the orthogonal_polys.py file.

### Changed 11 years ago by maldun

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

### Changed 11 years ago by maldun

Version from 10. August 2010

### Changed 11 years ago by maldun

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

### comment:1 follow-up: ↓ 2 Changed 11 years ago by maldun

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

### comment:2 in reply to: ↑ 1 Changed 11 years ago by maldun

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

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

### comment:4 Changed 11 years ago by maldun

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

a little comparison on my machine: old version:

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

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

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

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

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

And now also wider symbolic evaluation is possible:

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

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

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

### comment:5 follow-up: ↓ 7 Changed 11 years ago by fredrik.johansson

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

Care to elaborate?

### comment:6 Changed 11 years ago by maldun

Killed bug in legendre_P

### comment:7 in reply to: ↑ 5 Changed 11 years ago by maldun

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

Care to elaborate?

Sorry for the late answer, I was on holidays.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

OverflowError: cannot convert float infinity to integer

Last edited 8 years ago by jdemeyer (previous) (diff)

### comment:8 follow-ups: ↓ 9 ↓ 10 Changed 11 years ago by fredrik.johansson

That looks strange. I get:

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


### comment:9 in reply to: ↑ 8 Changed 11 years ago by maldun

That looks strange. I get:

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


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

### comment:10 in reply to: ↑ 8 Changed 11 years ago by maldun

That looks strange. I get:

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


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

### Changed 11 years ago by maldun

Version from 19. August 2010

### comment:11 follow-ups: ↓ 12 ↓ 14 Changed 11 years ago by maldun

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

Only the Legendre functions are still using Maxima.

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

Now I need some people for testing this out =)

### comment:12 in reply to: ↑ 11 Changed 11 years ago by maldun

And there was an interisting bug:

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

I think I should report this..

### comment:13 Changed 11 years ago by maldun

• Type changed from defect to enhancement

### Changed 11 years ago by maldun

Added numpy support, eliminated some bugs (19.08.2010)

### comment:14 in reply to: ↑ 11 Changed 11 years ago by maldun

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

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

### comment:15 Changed 11 years ago by maldun

• Status changed from new to needs_review

### comment:16 Changed 11 years ago by maldun

• Milestone set to sage-5.0

### comment:17 Changed 11 years ago by maldun

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

### comment:18 Changed 11 years ago by maldun

• Milestone changed from sage-5.0 to sage-4.5.3

### comment:19 Changed 11 years ago by maldun

• Owner changed from burcin to burcin, maldun

### comment:20 follow-up: ↓ 21 Changed 11 years ago by burcin

Hi Stefan,

can you post a patch corresponding to attachment:orthogonal_polys.8.py for review?

Thanks,
Burcin

### comment:21 in reply to: ↑ 20 Changed 11 years ago by maldun

Hi Stefan,

can you post a patch corresponding to attachment:orthogonal_polys.8.py for review?

Thanks,
Burcin

Done!

### comment:22 follow-up: ↓ 23 Changed 11 years ago by fredrik.johansson

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?

Some complex tests would be nice.

### comment:23 in reply to: ↑ 22 ; follow-up: ↓ 24 Changed 11 years ago by maldun

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?

Some complex tests would be nice.

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.

Of course can you call _evalf_ just with (), and then the default value is used.

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.

### comment:24 in reply to: ↑ 23 Changed 11 years ago by maldun

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?

Some complex tests would be nice.

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.

Of course can you call _evalf_ just with (), and then the default value is used.

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

### comment:25 Changed 11 years ago by maldun

Since it seems that numpy-1.4.1, and scipy 0.8 should work now (see #9808) 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.

Another thing I have to mention are these 2 failde doctests:

• sage -t -long "devel/sage/sage/symbolic/random_tests.py"
• sage -t -long "devel/sage/sage/symbolic/pynac.pyx"
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
**********************************************************************
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


I quite understand these, because we have introduced new functions, but I don't understand the exception in the last one

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)
**********************************************************************
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


And these are really strange, because when I type then into sage by hand everything works. wtf?? Can anyone have a look at these?

### Changed 11 years ago by maldun

ortho polys with scipy support

### comment:26 Changed 11 years ago by maldun

• Milestone changed from sage-4.6 to sage-5.0
• Status changed from needs_review to needs_work

### comment:27 Changed 11 years ago by kcrisman

Just cc:ing myself by commenting.

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.

### comment:28 Changed 11 years ago by maldun

• Description modified (diff)
• Status changed from needs_work to needs_review

@kcrisman thanks for paying attention. I added now an updated patch and extended instructions.

the doctest changes in symbolic.random_tests.py are easy to explain: new functions are involved -> new random expressions. But I had to change random_expr(50, nvars=3, coeff_generator=CDF.random_element)  to random_expr(60, nvars=3, coeff_generator=CDF.random_element)  or else one gets an expression generated where a division through zero occours.

As mentioned on sage-devel I repaired the doctests in symbolic.pynac.pyx, the trick is to enlarge the range of the for loop: for i in range(get_ginac_serial(), get_ginac_serial()+50): changed to for i in range(get_ginac_serial(), get_ginac_serial()+100): 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.

All doctests pass now, so I think a review would be nice.

-maldun

### comment:29 Changed 11 years ago by maldun

• Description modified (diff)

Cleaned up discription of the ticket and some comments in the ortho polys file.

### comment:30 follow-up: ↓ 31 Changed 11 years ago by kcrisman

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.

### Changed 11 years ago by maldun

Latest version of orthogonal polys with scipy support, and changed doctests. Tested in sage-4.6.alpha3

### comment:31 in reply to: ↑ 30 Changed 11 years ago by maldun

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.

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.

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.

PS: If diffs or more changelogs are needed let me know. I'm keeping track with git on my machine of the changes.

### comment:32 Changed 11 years ago by burcin

• Reviewers set to Burcin Erocal
• Status changed from needs_review to needs_work

Great work Stefan. Your patch looks good overall, but it needs a lot of polish. Thank you very much for this.

Here are some quick comments after reading attachment:trac_9706_orthogonal_polys.patch. 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.

• 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 maldun for copyright questions.
• You shouldn't import any part of numpy at the module level. This slows down startup too much. See #3561 for example. I'd say the same holds for mpmath and scipy.
• line 385-386 has this:
Then after using one of these functions, it changes:: (The value is now
False for chebyshev_T because chebyshev_T uses clenshaw method instead...)

I don't think this is valid Sphinx.
• delete line 412
#load /home/maldun/sage/sage-4.5.2/devel/sage-ortho/sage/functions/orthogonal_polys.py

• line 419: he -> the
• There are no doctests for the OrthogonalPolynomial class, make sure your file passes sage -coverage
• The commented timings in the docstring of OrthogonalPolynomial._clenshaw_method_() 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.
• In the docstring of OrthogonalPolynomial._eval_()
• remove the empty first line (line 494) of
• remove the commented out timings as well
• you need an empty line after EXAMPLES::
• the empty last line should be removed
• add some comments to the OrthogonalPolynomial._eval_() method to indicate what you're trying to do with these tests.
• lines 583-593 have a confusing comment and a bug
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()

• You don't need to state "Class for" on line 598, "The Chebyshev ..." is enough.
• Why do you delete the chebyshev_T(2,x) test on line 371? You can just add the new ones after that.
• line 626, EXAMPLES: -> EXAMPLES::
• Don't use *args or **kwds when you don't need them. Name the arguments and be explicit. Remember the "Zen of Python", "Explicit is better than implicit."
• OK, generally, fix the docstrings to conform to Sphinx standards. This should be documented somewhere in the developers guide.
• line 673, _maxima_init_evaled_() doesn't have doctests.
• line 678 - , _clenshaw_method_()
• docstring is not indented properly.
• It would be better to put the recursion formula in the docstring.
• line 790 _clenshaw_method_() doesn't have doctests.
• There is something wrong with the _maxima_init_evaled_() 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 .sage() on it. Never use sage_eval() on a string in the Sage library.
• 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:
try:
step_parent = kwds['parent']
except KeyError:
step_parent = parent(args[-1])

try:
precision = step_parent.prec()
except AttributeError:
precision = RR.prec()

See #9566.
• 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.
• Document the derivative formula in the docstring, using proper math notation
• What needs to be discussed from the comments on line 968-974?
• Same for lines 1058-1060?
• no doctests for _clenshaw_method_() on line 1156.
• no doctests for _maxima_init_evaled_() on line 1189.

I give up at this point. It seems that there are similar issues in the rest of the file as well.

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.

Thanks again for all your work.

### comment:34 Changed 8 years ago by jdemeyer

• Milestone changed from sage-5.11 to sage-5.12

### comment:35 Changed 8 years ago by maldun

Hi!

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.

Hope this time everything will work out!

### comment:36 Changed 8 years ago by maldun

• Status changed from needs_work to needs_review

### comment:37 Changed 8 years ago by tscrim

• Description modified (diff)

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).

Best,
Travis

For patchbot:

Apply: trac_9706_chebyshev.patch​, trac_9706-review-ts.patch

Last edited 8 years ago by tscrim (previous) (diff)

### comment:38 Changed 8 years ago by maldun

• Summary changed from New Version of orthogonal Polynomials to New Version of orthogonal Polynomials (Part I: Base Class and Chebyshev Polynomials)

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

### comment:39 Changed 8 years ago by tscrim

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.

For patchbot:

Apply: trac_9706_chebyshev.patch​, trac_9706-review-ts.patch

### comment:40 Changed 8 years ago by maldun

• Reviewers changed from Burcin Erocal to Burcin Erocal, Travis Scrimshaw
• Status changed from needs_review to positive_review

Thanks for your hard work in correcting all those small mistakes!

I'm very happy, that finally the new ortho polys are going into sage!

### comment:41 Changed 8 years ago by kcrisman

• Authors changed from Stefan Reiterer to Stefan Reiterer, Travis Scrimshaw
• Reviewers changed from Burcin Erocal, Travis Scrimshaw to Burcin Erocal, Travis Scrimshaw, Stefan Reiterer

I'm elevating Travis to an author because these are quite substantial review changes - thanks for being so meticulous on the formatting etc!

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!

### comment:42 follow-up: ↓ 68 Changed 8 years ago by jdemeyer

Sorry to spoil the party, but this is a regression:

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


(and yes: Chebyshev polynomials are important in number theory)

### comment:43 Changed 8 years ago by jdemeyer

• Status changed from positive_review to needs_work

### comment:44 Changed 8 years ago by jdemeyer

I would also like to point out that PARI is faster at evaluating Chebyshev polynomials:

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

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


We should definately use PARI to evaluate Chebyshev polynomials for certain types of input.

### comment:45 Changed 8 years ago by jdemeyer

Another regression:

sage: parent(chebyshev_T(10^2, RIF(2)))
Real Field with 53 bits of precision


### comment:46 Changed 8 years ago by jdemeyer

The "Clenshaw method" uses a very naive method of evaluating the recursion which needs O(n) steps, while there is a much faster method (which compute T_2n and U_2n in function of T_n and U_n) which only needs O(log(n)) steps.

Even this is totally feasible:

sage: timeit('''pari('polchebyshev(10^10, 1, Mod(2,1009))')''')
625 loops, best of 3: 16.3 µs per loop


### comment:47 Changed 8 years ago by jdemeyer

sage: R.<x> = QQ[]
sage: parent(chebyshev_T(5, x))
Symbolic Ring

Last edited 8 years ago by jdemeyer (previous) (diff)

### comment:48 Changed 8 years ago by jdemeyer

Another suggestion: why not make the base class more general and call it SymbolicPolynomial? I think that's a more natural class of functions, which could (in the future) also include cyclotomic polynomials for example.

### comment:49 Changed 8 years ago by jdemeyer

I propose the logic for evaluating chebyshev_T(n, x) should be:

1. if x is symbolic, then use the method of the current patch.
2. if x is not symbolic, try evaluation using PARI.
3. if conversion to PARI fails (for example for RIF), use an efficient O(log(n)) recursion algorithm.

### comment:50 follow-up: ↓ 51 Changed 8 years ago by fredrik.johansson

Jeroen: do you know a reference for the recursion pari uses?

### comment:51 in reply to: ↑ 50 Changed 8 years ago by jdemeyer

Jeroen: do you know a reference for the recursion pari uses?

No, but it's pretty straight-forward (think doubling formulas for cos and sin).

### comment:52 follow-up: ↓ 53 Changed 8 years ago by maldun

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.

@jdemeyer

@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.

@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)

@SymbolicPolynomial?: 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 SymbolicPolynomial? class I would propose to introduce the SymbolicPolynomial? class, and derive the OrthogonalPolynomials? from that class. I make the following suggestion: I will finish the OrthogonalPolynomials? 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 SymbolicPolynomials? have in common? -What should they have concerning other general polynomials?

But it would be really good to add the ortho polys, and I really want to finish this task.

Last edited 8 years ago by maldun (previous) (diff)

### comment:53 in reply to: ↑ 52 ; follow-up: ↓ 54 Changed 8 years ago by jdemeyer

SymbolicPolynomial: I don't think this is a good idea because ortho polys are quite special even among the polynomials.

What's special about orthogonal polynomials from a computer algebra 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).

I make the following suggestion: I will finish the OrthogonalPolynomials? with the current design.

Sure...

### comment:54 in reply to: ↑ 53 ; follow-up: ↓ 56 Changed 8 years ago by maldun

SymbolicPolynomial: I don't think this is a good idea because ortho polys are quite special even among the polynomials.

What's special about orthogonal polynomials from a computer algebra 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).

It makes a difference in an OO-Design 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.

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.

Last edited 8 years ago by maldun (previous) (diff)

### comment:55 follow-up: ↓ 57 Changed 8 years ago by jdemeyer

Evaluating chebyshev_T(n,x) can be done as

(Matrix(2,2,[x,x^2-1,1,x])^n)[0,0]


While chebyshev_U(n-1,x) equals

(Matrix(2,2,[x,x^2-1,1,x])^n)[1,0]


These can be evaluated with O(log(n)) operations.

### comment:56 in reply to: ↑ 54 Changed 8 years ago by jdemeyer

And I'm already thinking on some features, which I want to implement in future version, which are very specific to ortho polys.

That alone is a very good to keep the OrthogonalPolynomial class.

### comment:57 in reply to: ↑ 55 Changed 8 years ago by maldun

Evaluating chebyshev_T(n,x) can be done as

(Matrix(2,2,[x,x^2-1,1,x])^n)[0,0]


While chebyshev_U(n-1,x) equals

(Matrix(2,2,[x,x^2-1,1,x])^n)[1,0]


These can be evaluated with O(log(n)) operations.

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.

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

For T_N one can now use the binary representation of N to recursively build T_N in O(log N) time

Last edited 8 years ago by maldun (previous) (diff)

### comment:58 Changed 8 years ago by jdemeyer

This should be a ValueError:

sage: chebyshev_T(1/2,0)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
----> 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


### comment:59 Changed 8 years ago by jdemeyer

This should be ArithmeticError (I guess), since deriving w.r.t. the index simply isn't defined:

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


### comment:60 Changed 8 years ago by jdemeyer

I personally like to define chebyshev_T(-n,x) = chebyshev_T(n,x) and chebyshev_U(-n,x) = -chebyshev_U(n-2,x). This makes sense from a number-theoretic point of view and is also consistent with chebyshev_T(n,cos(x)) = cos(n*x) and chebyshev_U(n-1,cos(x)) = sin(n*x)/sin(x).

### comment:61 Changed 8 years ago by jdemeyer

For real/complex fields, you should use

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))


### comment:62 follow-up: ↓ 63 Changed 8 years ago by fredrik.johansson

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).

### comment:63 in reply to: ↑ 62 Changed 8 years ago by jdemeyer

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).

The matrix algorithm does seem more numerically stable (checked by using both algorithms inside RIF). So it's easy then, if there is one algorithm which is obviously the best.

### comment:64 Changed 8 years ago by maldun

I found an interesting paper on this topic: http://www.mathematik.uni-kassel.de/~koepf/cheby.pdf

Maybe this could give some input to the discussion.

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?

Considering rational numbers as input, the recursion formula works best, due to this paper. This should be considered too.

@Matrix Multiplication: It's also my favorite, but I will compare first.

Last edited 8 years ago by maldun (previous) (diff)

### comment:65 Changed 8 years ago by fredrik.johansson

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.

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]


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.

Similar code for U (using same algorithm as Pari):

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]


Edit: streamlined the code.

Last edited 8 years ago by fredrik.johansson (previous) (diff)

### comment:66 Changed 8 years ago by maldun

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...

### comment:67 Changed 8 years ago by maldun

this would be a functioning iterative algorithm:

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


I made it shorter. What is preferable? recursive or iterative? Normaly iterative, but in this case it is longer due to the indices battles ...

Edit: I also measured time: recursive is faster, even if I do some optimization (e.g. xrange)

Last edited 8 years ago by maldun (previous) (diff)

### comment:68 in reply to: ↑ 42 ; follow-up: ↓ 69 Changed 8 years ago by maldun

Sorry to spoil the party, but this is a regression:

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


(and yes: Chebyshev polynomials are important in number theory)

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.

The only way around this would be a hack to chatch this coercion troubles, e.g.

def __call__(self,n,x):

try:
super(Func_chebyshev_T,self).__call__(n,x)
except TypeError:
return self._clenshaw_method_(n,x)


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...

### Changed 8 years ago by maldun

new fixes and patches

### comment:69 in reply to: ↑ 68 ; follow-up: ↓ 71 Changed 8 years ago by jdemeyer

would work for that problem and some others. But how to catch such things correctly?

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 chebyshev(n,x), if either n or x is symbolic, we are in the symbolic case, otherwise we're in the algebraic case. In the algebraic case, the index n must be a concrete integer and we use the iterative algorithm.

### comment:70 follow-up: ↓ 74 Changed 8 years ago by jdemeyer

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.

### comment:71 in reply to: ↑ 69 ; follow-up: ↓ 72 Changed 8 years ago by maldun

would work for that problem and some others. But how to catch such things correctly?

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 chebyshev(n,x), if either n or x is symbolic, we are in the symbolic case, otherwise we're in the algebraic case. In the algebraic case, the index n must be a concrete integer and we use the iterative algorithm.

Thank you for the input, that's a great idea!

Based on that I propose the following switch:

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)


### comment:72 in reply to: ↑ 71 Changed 8 years ago by jdemeyer

I think that chebyshev_T(1/2, 2) should raise a ValueError (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.

So, I would do something like

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)

Last edited 8 years ago by jdemeyer (previous) (diff)

### comment:73 Changed 8 years ago by jdemeyer

• Description modified (diff)

### comment:74 in reply to: ↑ 70 Changed 8 years ago by maldun

• Description modified (diff)