Ask Your Question

# Why can Sage not compute this efficiently?

Consider this function:

 def A261042_list(len):
f = lambda k: x*2*(k+1)*(k+2)
g = *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, 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

## Comments

.

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


:P

## 2 Answers

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


to your for loop.

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 = *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).list()


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

more

## Comments

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

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

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

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

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

## Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

## Stats

Asked: 2015-08-08 22:19:01 +0200

Seen: 164 times

Last updated: Aug 09 '15