Ask Your Question
0

Why can Sage not compute this efficiently?

asked 2015-08-08 15:19:01 -0500

Peter Luschny gravatar image

updated 2015-08-09 07:00:46 -0500

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

:P

tmonteil gravatar imagetmonteil ( 2015-08-08 15:29:26 -0500 )edit

Clear the browser cache and try again :)

Peter Luschny gravatar imagePeter Luschny ( 2015-08-08 16:14:22 -0500 )edit

2 answers

Sort by ยป oldest newest most voted
2

answered 2015-08-08 15:41:31 -0500

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

Peter Luschny gravatar imagePeter Luschny ( 2015-08-08 16:05:31 -0500 )edit

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

Peter Luschny gravatar imagePeter Luschny ( 2015-08-08 16:27:34 -0500 )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 gravatar imagevdelecroix ( 2015-08-09 04:31:02 -0500 )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 gravatar imagePeter Luschny ( 2015-08-09 06:42:21 -0500 )edit
1

answered 2015-08-08 15:40:03 -0500

tmonteil gravatar image

updated 2015-08-08 15:43:55 -0500

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 more

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

Asked: 2015-08-08 15:19:01 -0500

Seen: 94 times

Last updated: Aug 09 '15