Why can Sage not compute this efficiently?

Peter Luschny
548 ●14 ●37 ●58

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.

Comments

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

tmonteil ( 9 years ago )

Clear the browser cache and try again :)

Peter Luschny ( 9 years ago )

add a comment

2 Answers

Sort by » oldest newest most voted

answered 9 years ago

vdelecroix
7640 ●20 ●85 ●168 http://www.labri.fr/pe...

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

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

link

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

Peter Luschny ( 9 years ago )

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

Peter Luschny ( 9 years ago )

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

vdelecroix ( 9 years ago )

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

Peter Luschny ( 9 years ago )

add a comment

answered 9 years ago

tmonteil
27483 ●32 ●205 ●519 http://wiki.sagemath.o...

updated 9 years ago

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]

link

add a comment

Your Answer

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

Add Answer

Why can Sage not compute this efficiently?

Comments

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Why can Sage not compute this efficiently? savecancel

Comments

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Why can Sage not compute this efficiently?