Why can Sage not compute this efficiently?

asked 2015-08-08 22:19:01 +0100

Peter Luschny
548 ●14 ●37 ●58

updated 2015-08-09 14:00:46 +0100

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 flag offensive close merge delete

Comments

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

tmonteil ( 2015-08-08 22:29:26 +0100 )edit

Clear the browser cache and try again :)

Peter Luschny ( 2015-08-08 23:14:22 +0100 )edit

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2015-08-08 22:41:31 +0100

vdelecroix
7645 ●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.

edit flag offensive delete 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 ( 2015-08-08 23:05:31 +0100 )edit

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

Peter Luschny ( 2015-08-08 23:27:34 +0100 )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).

vdelecroix ( 2015-08-09 11:31:02 +0100 )edit

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

Peter Luschny ( 2015-08-09 13:42:21 +0100 )edit

add a comment

answered 2015-08-08 22:40:03 +0100

tmonteil
27423 ●31 ●204 ●517 http://wiki.sagemath.o...

updated 2015-08-08 22:43:55 +0100

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]

edit flag offensive delete 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? edit

Comments

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Why can Sage not compute this efficiently?