# Why can Sage not compute this efficiently?

Consider this function:

 def A261042_list(len):
f = lambda k: x*2*(k+1)*(k+2)
g = [0]*len
g[len-1] = 1
for k in range(len-2,-1,-1):
g[k] = 1 - f(k)/(f(k)-1/g[k+1])
return taylor(g[0], x, 0, len-1).list()

A261042_list(6)


But taking larger arguments, for example A261042_list(16), it takes longer than I have patience to wait. In contrast Maple and Mathematica return immediately with small arguments like this.

Provided I have not made a mistake I do not ask for a workaround, I ask if this behavior is acknowledged as a bug.

Edit: By the comments of vdelecroix we can write

def A261042_list(len):
f = lambda k: x*2*(k+1)*(k+2)
g = 1
for k in range(len-2,-1,-1):
g = (1-f(k)/(f(k)-1/g)).simplify_rational()
return taylor(g, x, 0, len-1).list()


On the other hand I am still not satisfied that we are forced to use 'simplify_rational' here.

edit retag close merge delete

.

sage: oeis(261042)
A261042: allocated for Peter Luschny


:P

( 2015-08-08 15:29:26 -0600 )edit

Clear the browser cache and try again :)

( 2015-08-08 16:14:22 -0600 )edit

Sort by » oldest newest most voted

One reason is because there is no simplification performed if you do not explicitely ask for them. If for example you print g[0] when len=5 you will see

-4*x/(4*x + 1/(12*x/(12*x + 1/(24*x/(24*x + 1/(40*x/(40*x - 1) - 1)) - 1)) - 1)) + 1


In order to avoid these nested fractions that are complicated to compute with, you might just add the line

g[k] = g[k].simplify_rational()


Another way is to not use the symbolic ring and use rational fractions defined over QQ

def A261042_list_bis(l):
x = polygen(QQ)
f = lambda k: x*2*(k+1)*(k+2)
g = [0]*l
g[l-1] = 1
for k in range(l-2,-1,-1):
g[k] = 1 - f(k)/(f(k)-1/g[k+1])
return PowerSeriesRing(QQ,'x',default_prec=l)(g[0]).list()


The last line is not very nice but I did not find out something better to compute the Taylor expansion.

more

You: "..there is no simplification performed if you do not explicitly ask for them." This is the point! I did neither ask Maple nor Mathematica explicitly to perform any simplifications and they compute fast. That's what I'm used to and I am sure many other user would appreciate to get it this way with Sage also. (In fact I know respected mathematicians which precisely for this reason prefer Maple or Mathematica over Sage.)

( 2015-08-08 16:05:31 -0600 )edit

Since I used simplify_rational() in the OEIS version I choose your answer.

( 2015-08-08 16:27:34 -0600 )edit

Actually, let me also say that you can simplify this function. There is no need to create the whole array. You can use only one function g initialized with 0 and update it with g = 1 - f(k)/(f(k)-1/g).

( 2015-08-09 04:31:02 -0600 )edit

In this form it would lead immediately to a division by zero. However I appreciate your suggestion (works with initial g=1).

( 2015-08-09 06:42:21 -0600 )edit

The reason is the following: add a print g[k] after the definition of g[k], you will see the following expressions (for len=10):

-180*x/(180*x - 1) + 1
-144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) + 1
-112*x/(112*x + 1/(144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) - 1)) + 1
-84*x/(84*x + 1/(112*x/(112*x + 1/(144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) - 1)) - 1)) + 1
-60*x/(60*x + 1/(84*x/(84*x + 1/(112*x/(112*x + 1/(144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) - 1)) - 1)) - 1)) + 1
-40*x/(40*x + 1/(60*x/(60*x + 1/(84*x/(84*x + 1/(112*x/(112*x + 1/(144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) - 1)) - 1)) - 1)) - 1)) + 1
-24*x/(24*x + 1/(40*x/(40*x + 1/(60*x/(60*x + 1/(84*x/(84*x + 1/(112*x/(112*x + 1/(144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) - 1)) - 1)) - 1)) - 1)) - 1)) + 1
-12*x/(12*x + 1/(24*x/(24*x + 1/(40*x/(40*x + 1/(60*x/(60*x + 1/(84*x/(84*x + 1/(112*x/(112*x + 1/(144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) - 1)) - 1)) - 1)) - 1)) - 1)) - 1)) + 1
-4*x/(4*x + 1/(12*x/(12*x + 1/(24*x/(24*x + 1/(40*x/(40*x + 1/(60*x/(60*x + 1/(84*x/(84*x + 1/(112*x/(112*x + 1/(144*x/(144*x + 1/(180*x/(180*x - 1) - 1)) - 1)) - 1)) - 1)) - 1)) - 1)) - 1)) - 1)) + 1


So, you can see that they become quite large, so that the computation become very long, which explains why your computation is slow.

So, to check our hypothesis (and solve the problem), what you have to do is to simplify the expressions to keep something small along the computation. You can work on fraction field (very fast), or, if you want to stay in the symbolic ring, then you should replace:

g[k] = 1 - f(k)/(f(k)-1/g[k+1])


with

g[k] = (1 - f(k)/(f(k)-1/g[k+1])).full_simplify()


Now, with len=20, you get without waiting:

[1,
4,
64,
2176,
126976,
11321344,
1431568384,
243680935936,
53725527801856,
14893509177769984,
5070334006399074304,
2079588119566033616896,
1011390382859091900891136,
575501120339508919401447424,
378784713733072451034702413824,
285539131625477547496925147693056,
244410463271028625971446390840098816,
235752874763994894181564461756568305664,
254537317950438535955585357041409084882944,
305766504814246648318888059404547150542012416]

more