Comment #11 on issue 3226 by smi...@gmail.com: high-order derivatives should be cse-simplified
http://code.google.com/p/sympy/issues/detail?id=3226

cse(sqrt((a+c*t)**2+(x+t)*(a+c*t)**2+a*(x+t))
... )
([(x0, t + x), (x1, (a + c*t)**2)], [sqrt(a*x0 + x0*x1 + x1)])


Without preprocessing:

_t=time();ok=diff(sqrt((a+c*t)**2+(x+t)*(a+c*t)**2+a*(x+t)),t,10);time()-_t
6.296999931335449
ok.count_ops()
32211
_t=time();ok=factor_terms(signsimp(diff(sqrt((a+c*t)**2+(x+t)
... *(a+c*t)**2+a*(x+t)),t,10)));time()-_t;ok.count_ops()
15.574000120162964
530
print filldedent(ok)

14175*(-40960*c**6*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
c*t)**2) - 61440*c**6*(2*a + 2*c*t + c*(t + x) + c)**2 +
215040*c**5*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
c*t)**2)**2*(2*a + 2*c*t + c*(t + x) + c)/(a*(t + x) + (a + c*t)**2*(t
+ x) + (a + c*t)**2) + 143360*c**5*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a
+ c*t) + (a + c*t)**2)*(2*a + 2*c*t + c*(t + x) + c)**3/(a*(t + x) +
(a + c*t)**2*(t + x) + (a + c*t)**2) + 7168*c**5*(2*a + 2*c*t + c*(t +
x) + c)**5/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2) -
80640*c**4*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
c*t)**2)**4/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**2 -
322560*c**4*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
c*t)**2)**3*(2*a + 2*c*t + c*(t + x) + c)**2/(a*(t + x) + (a +
c*t)**2*(t + x) + (a + c*t)**2)**2 - 80640*c**4*(a + 2*c*(a + c*t)*(t
+ x) + 2*c*(a + c*t) + (a + c*t)**2)**2*(2*a + 2*c*t + c*(t + x) +
c)**4/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**2 +
177408*c**3*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
c*t)**2)**5*(2*a + 2*c*t + c*(t + x) + c)/(a*(t + x) + (a + c*t)**2*(t
+ x) + (a + c*t)**2)**3 + 147840*c**3*(a + 2*c*(a + c*t)*(t + x) +
2*c*(a + c*t) + (a + c*t)**2)**4*(2*a + 2*c*t + c*(t + x) +
c)**3/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**3 -
27456*c**2*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
c*t)**2)**7/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**4 -
96096*c**2*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
c*t)**2)**6*(2*a + 2*c*t + c*(t + x) + c)**2/(a*(t + x) + (a +
c*t)**2*(t + x) + (a + c*t)**2)**4 + 25740*c*(a + 2*c*(a + c*t)*(t +
x) + 2*c*(a + c*t) + (a + c*t)**2)**8*(2*a + 2*c*t + c*(t + x) +
c)/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**5 - 2431*(a +
2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a + c*t)**2)**10/(a*(t + x) +
(a + c*t)**2*(t + x) + (a + c*t)**2)**6)/(1024*(a*(t + x) + (a +
c*t)**2*(t + x) + (a + c*t)**2)**(7/2))


Start new session with

def di(eq, *syms):
    # find common subexpressions
    r, e = cse(eq)
    e = e[0]
    # last must be first
    r = list(reversed(r))
    # for Xi, f(xi), replace x0 -> x0(xi) and diff and simplify
    s2f = dict([(ri[0], ri[0](*ri[1].free_symbols)) for ri in r])
    rv = e.xreplace(s2f).diff(*syms)
    # replace generic derivatives w/ true cse expression derivatives and
    # replace generic functions w/ true cse functions
    f = [b for _, b in s2f.iteritems()]
    F = [B for _, B in r]
    f2F = dict(list(zip(f, F)))
reps = df2dF = dict([(f.diff(*syms), F.diff(*syms)) for f, F in f2F.iteritems()])
    reps.update(f2F)
    rv = rv.xreplace(reps).doit()
    return factor_terms(signsimp(rv))


_t=time();ok=di(sqrt((a+c*t)**2+(x+t)*(a+c*t)**2+a*(x+t)),t,10).doit();time(
)-_t;ok.count_ops()
38.84599995613098
552
print filldedent(ok)

14175*(-61440*c**6*(a*c + 2*a + 2*c*t + c*(t + x))**2 -
40960*c**6*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) + (a + c*t)**2 +
1) + 7168*c**5*(a*c + 2*a + 2*c*t + c*(t + x))**5/(a*(a + c*t)**2 + t
+ x + (a + c*t)**2*(t + x)) + 143360*c**5*(a*c + 2*a + 2*c*t + c*(t +
x))**3*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) + (a + c*t)**2 +
1)/(a*(a + c*t)**2 + t + x + (a + c*t)**2*(t + x)) + 215040*c**5*(a*c
+ 2*a + 2*c*t + c*(t + x))*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) +
(a + c*t)**2 + 1)**2/(a*(a + c*t)**2 + t + x + (a + c*t)**2*(t + x)) -
80640*c**4*(a*c + 2*a + 2*c*t + c*(t + x))**4*(2*a*c*(a + c*t) +
2*c*(a + c*t)*(t + x) + (a + c*t)**2 + 1)**2/(a*(a + c*t)**2 + t + x +
(a + c*t)**2*(t + x))**2 - 322560*c**4*(a*c + 2*a + 2*c*t + c*(t +
x))**2*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) + (a + c*t)**2 +
1)**3/(a*(a + c*t)**2 + t + x + (a + c*t)**2*(t + x))**2 -
80640*c**4*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) + (a + c*t)**2 +
1)**4/(a*(a + c*t)**2 + t + x + (a + c*t)**2*(t + x))**2 +
147840*c**3*(a*c + 2*a + 2*c*t + c*(t + x))**3*(2*a*c*(a + c*t) +
2*c*(a + c*t)*(t + x) + (a + c*t)**2 + 1)**4/(a*(a + c*t)**2 + t + x +
(a + c*t)**2*(t + x))**3 + 177408*c**3*(a*c + 2*a + 2*c*t + c*(t +
x))*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) + (a + c*t)**2 +
1)**5/(a*(a + c*t)**2 + t + x + (a + c*t)**2*(t + x))**3 -
96096*c**2*(a*c + 2*a + 2*c*t + c*(t + x))**2*(2*a*c*(a + c*t) +
2*c*(a + c*t)*(t + x) + (a + c*t)**2 + 1)**6/(a*(a + c*t)**2 + t + x +
(a + c*t)**2*(t + x))**4 - 27456*c**2*(2*a*c*(a + c*t) + 2*c*(a +
c*t)*(t + x) + (a + c*t)**2 + 1)**7/(a*(a + c*t)**2 + t + x + (a +
c*t)**2*(t + x))**4 + 25740*c*(a*c + 2*a + 2*c*t + c*(t +
x))*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) + (a + c*t)**2 +
1)**8/(a*(a + c*t)**2 + t + x + (a + c*t)**2*(t + x))**5 -
2431*(2*a*c*(a + c*t) + 2*c*(a + c*t)*(t + x) + (a + c*t)**2 +
1)**10/(a*(a + c*t)**2 + t + x + (a + c*t)**2*(t + x))**6)/(1024*(a*(a
+ c*t)**2 + t + x + (a + c*t)**2*(t + x))**(7/2))

Compare with non-preprocessed result
test_numerically(ok,S('''14175*(-40960*c**6*(a + 2*c*(a + c*t)*(t + x) + 2*c
*(a + c*t) + (a +
... c*t)**2) - 61440*c**6*(2*a + 2*c*t + c*(t + x) + c)**2 +
... 215040*c**5*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
... c*t)**2)**2*(2*a + 2*c*t + c*(t + x) + c)/(a*(t + x) + (a + c*t)**2*(t
... + x) + (a + c*t)**2) + 143360*c**5*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a
... + c*t) + (a + c*t)**2)*(2*a + 2*c*t + c*(t + x) + c)**3/(a*(t + x) +
... (a + c*t)**2*(t + x) + (a + c*t)**2) + 7168*c**5*(2*a + 2*c*t + c*(t +
... x) + c)**5/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2) -
... 80640*c**4*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
... c*t)**2)**4/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**2 -
... 322560*c**4*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
... c*t)**2)**3*(2*a + 2*c*t + c*(t + x) + c)**2/(a*(t + x) + (a +
... c*t)**2*(t + x) + (a + c*t)**2)**2 - 80640*c**4*(a + 2*c*(a + c*t)*(t
... + x) + 2*c*(a + c*t) + (a + c*t)**2)**2*(2*a + 2*c*t + c*(t + x) +
... c)**4/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**2 +
... 177408*c**3*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
... c*t)**2)**5*(2*a + 2*c*t + c*(t + x) + c)/(a*(t + x) + (a + c*t)**2*(t
... + x) + (a + c*t)**2)**3 + 147840*c**3*(a + 2*c*(a + c*t)*(t + x) +
... 2*c*(a + c*t) + (a + c*t)**2)**4*(2*a + 2*c*t + c*(t + x) +
... c)**3/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**3 -
... 27456*c**2*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
... c*t)**2)**7/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**4 -
... 96096*c**2*(a + 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a +
... c*t)**2)**6*(2*a + 2*c*t + c*(t + x) + c)**2/(a*(t + x) + (a +
... c*t)**2*(t + x) + (a + c*t)**2)**4 + 25740*c*(a + 2*c*(a + c*t)*(t +
... x) + 2*c*(a + c*t) + (a + c*t)**2)**8*(2*a + 2*c*t + c*(t + x) +
... c)/(a*(t + x) + (a + c*t)**2*(t + x) + (a + c*t)**2)**5 - 2431*(a +
... 2*c*(a + c*t)*(t + x) + 2*c*(a + c*t) + (a + c*t)**2)**10/(a*(t + x) +
... (a + c*t)**2*(t + x) + (a + c*t)**2)**6)/(1024*(a*(t + x) + (a +
... c*t)**2*(t + x) + (a + c*t)**2)**(7/2))
... '''))
False

So I must be making a mistake someplace. And the preprocessing took longer. Does anyone see my mistake?

--
You received this message because this project is configured to send all issue notifications to this address.
You may adjust your notification preferences at:
https://code.google.com/hosting/settings

--
You received this message because you are subscribed to the Google Groups 
"sympy-issues" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sympy-issues+unsubscr...@googlegroups.com.
To post to this group, send email to sympy-issues@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy-issues?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.


Reply via email to