Emmanuel Charpentier's profile - activity

I'm not well versed in this kind of problems, but I thing that chapter 17 of this book should give you a head start.

Well...

sage: r.library('fortunes')
sage: r.fortune("'TFM'")

This is all documented in TFM. Those who WTFM don't want to have to WTFM again
on the mailing list. RTFM.
   -- Barry Rowlingson
      R-help (October 2003)

NOTE : This has replaced, (not complemented) my original answer concluding that this was a (Python) bug. Apologies for the confusion...

As pointed out by Serge Lelièvre, this seems a limitation of Python 2 floats, which does not exist in Python 3, where one can do :

sage: def F(u): return arg((u-3)^(1/3))
sage: F(1)
arg((-2)^(1/3))
sage: F(1).n()
1.04719755119660
sage: F(2).n()
1.04719755119660
sage: plot(F,(1,2))
Launched png viewer for Graphics object consisting of 1 graphics primitive
sage: integrate(F(x),x,1,2)
1/3*pi

See this sage-devel thread for further illumination...

In the original example:

sage: cauchy(z)=y.substitute(solve(z*y^3 +y^2 - 2*z*y+2, y)[0])
sage: def G(u): return arg(cauchy(u))
sage: numerical_integral(lambda u:G(u), 1, 2)
(0.6116545934235456, 6.790730127365591e-15)

This question deserves a new tag, such as Nine-billion-names-of-god, but I don't know how to retag a question...

I do not understand what you're trying to accomplish (computing the nine billions names of God ?).

I note that you are doing nothing with the result of do_calculations : it's neither used (to what ends ?), printed (where ?) nor stored (where ?).

For argument's sake, let's have a look at the size of problem. your elements are all the nine-long lists whose elements belong to [0, 1, 2, 3]. There are 4^9=262144 of them. Your (potential) table of results has therefore (4^9)^2=68719476736 elements (about seventy billions).

How much time is necessary to compute one element ? Let's see : pick two elements at random and operate on them. Using your original definitions in a Python3-based Sage :

sage: A=tuple(r.sample([0, 1, 2, 3], 9, replace=True).sage())
sage: B=tuple(r.sample([0, 1, 2, 3], 9, replace=True).sage())
sage: time(do_calculations(A,B))
CPU times: user 1.21 ms, sys: 0 ns, total: 1.21 ms
Wall time: 1.22 ms
[0, 0, 0, 0, 0, 1, 1, 1, 0]

You would need about (((4^9)^2)*1.22e-3).n(digits=3)=8.38e7 seconds (i. about (((4^9)^2)*1.22e-3/(24*60*60*365.25)).n(digits=3)=2.66 years) to perform the computations.

Storing the result (using an efficient encoding, i. e. two bits per element) but no previsible gain from string compressin algorithm (since your computation doesn't seem to have any special properties), you would need about 17.5 GB of storage.

Therefore, absent any "interesting" properties of your addition and multiply functions, inducing "interesting" symmetries (i. e. a structure) on the set of your nine-letter words, the "brute-force" approach is doomed.

May I suggest that using (a part of) 2.66 years to study the possible algebraic properties of your elements and operations might be more fruitful ?

sage: X = [1, 5, 7, 8, 8.7, 10, 13, 15] sage: Y = [2, 8, 13, 10, 9, 6.3, 2, -1] sage: points(zip(X,Y)) Launched png viewer for Graphics object consisting of 1 graphics primitive

sage: line(zip(X,Y))

?timeit would enlighten you...

Homework ?

Homework ?

The Web server software seems to have munched my answer irretrievably. Sorry...

That was the point of the old Sage notebook, now considered as unmaintainable.

Attempts are done to interface Sagemath and Jupyterlab, which would allow such multiple interactions. Search for that in the sage-devel archives. IIRC, this will be much easier with a Python3-based Sagemath.

You can also develop a dedicated Web server using Sagemath as a backend. The real difficulty is to choose between a "single-shot" service (i. e. any Sagemath use is self-contained) and a session-based service (i.e. queries at time t can use definitions and statements made at t'<t). in="" the="" first="" case,="" you="" use="" sagemath="" almost="" as="" a="" library,="" in="" the="" second="" case,="" your="" application="" must="" maintain="" its="" sessions="" (or,="" equivalently,="" launch="" one="" sagemath="" session="" per="" user.<="" p="">

HTH,

A bug, IMHO...

Homework ?

Okay. Time to think a bit.

Let's try to solve this in SR, the symbolic ring. Let's get rid of the goddamn numeric constants :

sage: var("k_1","k_2")
(k_1, k_2)
sage: E1="-J_d^2*k_yy*k_zz*omega^4 + 0.382*J_d^2*k_yy*omega^6 + 0.382*J_d^2*k_zz
....: *omega^6 - 0.145924*J_d^2*omega^8 + J_d*k_phiphi*k_yy*k_zz*omega^2 - 0.382
....: *J_d*k_phiphi*k_yy*omega^4 - 0.382*J_d*k_phiphi*k_zz*omega^4 + 0.145924*J_
....: d*k_phiphi*omega^6 + J_d*k_thetatheta*k_yy*k_zz*omega^2 - 0.382*J_d*k_thet
....: atheta*k_yy*omega^4 - 0.382*J_d*k_thetatheta*k_zz*omega^4 + 0.145924*J_d*k
....: _thetatheta*omega^6 - J_d*k_yphi^2*k_zz*omega^2 + 0.382*J_d*k_yphi^2*omega
....: ^4 - J_d*k_yy*k_ztheta^2*omega^2 + 0.382*J_d*k_ztheta^2*omega^4 + J_p^2*Om
....: ega^2*k_yy*k_zz*omega^2 - 0.382*J_p^2*Omega^2*k_yy*omega^4 - 0.382*J_p^2*O
....: mega^2*k_zz*omega^4 + 0.145924*J_p^2*Omega^2*omega^6 - k_phiphi*k_thetathe
....: ta*k_yy*k_zz + 0.382*k_phiphi*k_thetatheta*k_yy*omega^2 + 0.382*k_phiphi*k
....: _thetatheta*k_zz*omega^2 - 0.145924*k_phiphi*k_thetatheta*omega^4 + k_phip
....: hi*k_yy*k_ztheta^2 - 0.382*k_phiphi*k_ztheta^2*omega^2 + k_thetatheta*k_yp
....: hi^2*k_zz - 0.382*k_thetatheta*k_yphi^2*omega^2 - k_yphi^2*k_ztheta^2"
sage: E2=SR(E1.replace("0.382","k_1").replace("0.145924","k_2"))

What's this ?

sage: E2.expand()
J_p^2*Omega^2*k_2*omega^6 - J_d^2*k_2*omega^8 - J_p^2*Omega^2*k_1*k_yy*omega^4 - J_p^2*Omega^2*k_1*k_zz*omega^4 + J_d^2*k_1*k_yy*omega^6 + J_d^2*k_1*k_zz*omega^6 + J_d*k_2*k_phiphi*omega^6 + J_d*k_2*k_thetatheta*omega^6 + J_p^2*Omega^2*k_yy*k_zz*omega^2 + J_d*k_1*k_yphi^2*omega^4 - J_d*k_1*k_phiphi*k_yy*omega^4 - J_d*k_1*k_thetatheta*k_yy*omega^4 + J_d*k_1*k_ztheta^2*omega^4 - J_d*k_1*k_phiphi*k_zz*omega^4 - J_d*k_1*k_thetatheta*k_zz*omega^4 - J_d^2*k_yy*k_zz*omega^4 - k_2*k_phiphi*k_thetatheta*omega^4 - k_1*k_thetatheta*k_yphi^2*omega^2 + k_1*k_phiphi*k_thetatheta*k_yy*omega^2 - k_1*k_phiphi*k_ztheta^2*omega^2 - J_d*k_yy*k_ztheta^2*omega^2 + k_1*k_phiphi*k_thetatheta*k_zz*omega^2 - J_d*k_yphi^2*k_zz*omega^2 + J_d*k_phiphi*k_yy*k_zz*omega^2 + J_d*k_thetatheta*k_yy*k_zz*omega^2 - k_yphi^2*k_ztheta^2 + k_phiphi*k_yy*k_ztheta^2 + k_thetatheta*k_yphi^2*k_zz - k_phiphi*k_thetatheta*k_yy*k_zz

Looks like a polynomial in omega. Let's check this : sage: LK=E2.expand().coefficients(omega); LK

[[-k_yphi^2*k_ztheta^2 + k_phiphi*k_yy*k_ztheta^2 + k_thetatheta*k_yphi^2*k_zz - k_phiphi*k_thetatheta*k_yy*k_zz,
  0],
 [J_p^2*Omega^2*k_yy*k_zz - k_1*k_thetatheta*k_yphi^2 + k_1*k_phiphi*k_thetatheta*k_yy - k_1*k_phiphi*k_ztheta^2 - J_d*k_yy*k_ztheta^2 + k_1*k_phiphi*k_thetatheta*k_zz - J_d*k_yphi^2*k_zz + J_d*k_phiphi*k_yy*k_zz + J_d*k_thetatheta*k_yy*k_zz,
  2],
 [-J_p^2*Omega^2*k_1*k_yy - J_p^2*Omega^2*k_1*k_zz + J_d*k_1*k_yphi^2 - J_d*k_1*k_phiphi*k_yy - J_d*k_1*k_thetatheta*k_yy + J_d*k_1*k_ztheta^2 - J_d*k_1*k_phiphi*k_zz - J_d*k_1*k_thetatheta*k_zz - J_d^2*k_yy*k_zz - k_2*k_phiphi*k_thetatheta,
  4],
 [J_p^2*Omega^2*k_2 + J_d^2*k_1*k_yy + J_d^2*k_1*k_zz + J_d*k_2*k_phiphi + J_d*k_2*k_thetatheta,
  6],
 [-J_d^2*k_2, 8]]

Hey ! it's indeed a polynomial of degree 4 in omega^2... Let's give its coefficients a more sy(n|mpa)thetic name (i. e. create a simpler-to-write polynomial having the same value given the right substitution) and check it :

sage: SK=[var("K_{}".format(u[1]))==u[0] for u in LK]
sage: E3=sum([SK[u].lhs()*omega^(2*u) for u in range(len(SK))])
sage: bool(E3.subs(SK)==E2)
True

Transform this in a polynomial in omega^2 :

sage: E4=E3.subs([omega^(2*u)==OmSq^u for u in range(len(SK))]);E4
K_8*OmSq^4 + K_6*OmSq^3 + K_4*OmSq^2 + K_2*OmSq + K_0

Now, this fourth-degree polynomial has four (possibly equal, possibly complex) roots:

sage: Sol4=[u.subs(SK) for u in solve(E4, OmSq)]
sage: len(Sol4)
4

Each of these solutions s for OmSq gives two possible solutions for omega (namely $-\sqrt{s}$ and $\sqrt{s}$) :

sage: Sol3=flatten([[omega==-sqrt(s.rhs()),omega==sqrt(s.rhs())] for s in Sol4])
sage: len(Sol3)
8

If one insists to reincorporate the goddamn numeric constants in the final solution, this is easy :

sage: Sol=[u.subs([k_1==0.382, k_2==0.145924]) for u in Sol3]

Printing the solutions is left as an exercise for the (masochistic) reader (fitted with an A0 printer and lots of time). Similarly, checking them (by substitution in E2) is quite slow...

An alternative (=more general (=better)) solution is to create a multivariate polynomial ring with the right unknowns (please choose easier names !) and compute the solutions by searching the "right" ideal. This would also allow you to create the set of roots without having to compute an explicit solution for them (which do not, in general, exists for degrees >=5).

HTH,

Indeed :

sage: var("y,z")
(y, z)
sage: cauchy(z)=solve(z*y^3 +y^2 - 2*z*y+2 ,y)[0].rhs()
sage: def foo(u):return(arg(cauchy(u)).n())
sage: [foo(t) for t in (1,1.1..2)]
[0.828222717321238,
## Snip...
 0.377665318352514]

Therefore, this function can be evaluated. But :

sage: plot(foo,(1,2))
verbose 0 (3635: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 200 points.
verbose 0 (3635: plot.py, generate_plot_points) Last error message: 'negative number cannot be raised to a fractional power'
Launched png viewer for Graphics object consisting of 0 graphics primitives

... not in plot. Trying to plot symbolically yelds :

TypeError: Cannot evaluate symbolic expression to a numeric value.

A bug, IMHO

POST-SCRIPTUM : It turns out that vectors have a list` method, that you can use to lighten the notations :

sage: solve((A*X-Y).list(),X.list())
[[x_1 == -(a_2_2*y_1 - a_1_2*y_2)/(a_1_2*a_2_1 - a_1_1*a_2_2), x_2 == (a_2_1*y_1 - a_1_1*y_2)/(a_1_2*a_2_1 - a_1_1*a_2_2)]]

Lighter and cleaner, indeed...

A*x returns a list

Nope :

sage: (A*X).parent()
Vector space of dimension 2 over Symbolic Ring

But A*x is an iterable, as well as Y, and map iterates through those, producing at each step one equation, and returns the list of those equations. Similarly, I do not use X as the second argument of solve, but the list of its components.

?map should enlighten you :-).

Try this :

sage: A=matrix([[var("a_{}_{}".format(u,v), latex_name="a_{{{},{}}}".format(u,v)) for v in (1..2)] for u in (1..2)]);A
[a_1_1 a_1_2]
[a_2_1 a_2_2]
sage: X=vector([var("x_{}".format(u)) for u in (1..2)]);X
(x_1, x_2)
sage: Y=vector([var("y_{}".format(u)) for u in (1..2)]);Y
(y_1, y_2)
sage: Sol=solve(map(lambda u,v:u==v, A*X, Y),[u for u in X]);Sol
[[x_1 == -(a_2_2*y_1 - a_1_2*y_2)/(a_1_2*a_2_1 - a_1_1*a_2_2), x_2 == (a_2_1*y_1 - a_1_1*y_2)/(a_1_2*a_2_1 - a_1_1*a_2_2)]]

i. e. $$\left[\left[x_{1} = -\frac{{a_{2,2}} y_{1} - {a_{1,2}} y_{2}}{{a_{1,2}} {a_{2,1}} - {a_{1,1}} {a_{2,2}}}, x_{2} = \frac{{a_{2,1}} y_{1} - {a_{1,1}} y_{2}}{{a_{1,2}} {a_{2,1}} - {a_{1,1}} {a_{2,2}}}\right]\right]$$.

which looks entirely reasonable...

HTH,

DanieleC : I agree. I misread/misunderstood the whishes of the anonymous OP. I'm not sure it makes sense...

You might be interested by this tutorial : the quaternions are a good example of a non-commutative ring. Matrices of such elements would have the properties you seek.

You may find inspiration in this free textbook and its Sagemath supplement...

Also, ISTR that Maxima allows you to define non-commutating sets of variables ; however, I have ni way to dive in its doc now, and I can't remember if it allows to define matrices of such elements...

Same difference. RTFM again...

No time to re-run your (heavy) program. But a cursory look makes me think that the problem might be with

show(Unum.substitute((varsol[0])=0,(varsol[1])=1)

That I'd write

show(Unum.substitute([varsol[0])==0, varsol[1]==1])

RTFM...

Okay. After 12 days, your (hypothetic) homework is past due. Let's look at a couple of solutions.

"when you're lost, draw what you know." Undergrads' wisdom. Let's see:

Without loss of generality, set up a framework placing the fence center in $O=(0,0)$ and the tether anchor in $A=(r_0,0)$. In this framework a point $(x,y)$ is edible by the goat if it is :

• inside the fence, implying $x^2+y^2<=r_0^2$, and
• closer to the tether's anchor than $r_1$, implying $(r-x)^2+y^2<=r1^2$.

In other words, the edible region is the intersection of two disks $D_0$, centered in $O$ with radius $r_0$ and $D_1$, centered in $A$ with radius $r_1$. This region can be exactly partitioned by the straight line segment joining the intersections of the two circles boundaries of these disks:

Each of these regions is the intersection of a disk of radius $r$ and a half-plane cutting the disk at a distance $X<=r$ of its center.

Therefore, a first, (almost) purely geometrical, solution to your problem is to derive the area of such a surface, and apply it to both the "yellow" and the "red" surfaces. Let's consider the "yellow" one, delimited bu th arc $B'AB$ and the straight segment $BB'$. This surface is symmetric about $OA$ ; therefore, its area is twice of the surface bounded by the arc $AB$ and the straight segments $BX$ and $XA$.

This latter surface can be in turn considered as the difference between

• the circular sector delimited by$OA$, $OB$, and
• the arc $AB$ and the triangle $BOX$.

Let's call $\theta$ the angle $\widehat{AOB}$. It can be determined from the position of $B$. We have of course $\cos\theta=\frac{OX}{r}=\frac{B_x}{r}$ and $\tan\theta={B_y}{B_x}$. At the high school level, this "determination" is a bit fuzzy, but will be clarified in first year of college...

The circular sector has area $\frac{r^2\theta}{2}$, the triangle has area $\frac{OX*XB}{2}=\frac{B_x B_y}{2}$. The yellow surface has area twice their difference, therefore $r_0^2 \theta - B_x\cdot B_y$.

For similar reasons (left to the reader as an exercise ; beware the signs !), with $\cos\theta_1=\frac{B_y}{r_1}$, the "red" sector has area $r_1^2\theta_1-(r_0-B_x) B_y$.

All that is lacking are the coordinates of $B$ which are determined by the fact that it belongs to both the boundary circles (the fence and the tether). Since there are two solution points, we'll select the one with positive ordinate :

x,y,r_0,r_1=var("x,y,r_0,r_1")
assume(x,"real",y,"real",r_0>0,r_1>0,r_1<2*r_0)
O=vector([0,0])
A=vector([r_0,0])
b=vector([x,y])
E0=(b-O).norm()^2==r_0^2
E1=(b-A).norm()^2==r_1^2
Sol=solve([E0, E1],[x,y], solution_dict=True)
if bool(Sol[0].get(y)>0):
    B=vector([Sol[0].get(x),Sol[0].get(y)])
else:
    B=vector([Sol[1].get(x),Sol[1].get(y)])
X=B[0]
Y=B[1]

That gives us $B=\left(\frac{2 \, r_{0}^{2} - r_{1}^{2}}{2 \, r_{0}},\ \frac{\sqrt{4 \, r_{0}^{2} - r_{1}^{2}} r_{1}}{2 \, r_{0}}\right)$.

This also determines the angle $\theta=\widehat{AOB}$ by either theta=arccos(X/r_0) or theta=arctan2(Y,X) (Note : the third "obvious" choice, theta=arcsin(Y./r_0 is unsuitable here, given the choice of branches of Sage's implementation of arcsin). We determine $\theta_1$ in a similar way (again, mind the signs...).

Summing the partial areas and finding the root (value of $r_1$ leaving an edible area of 10 when $r_0=10$) is left as an exercise to the reader...

This first answer will be edited by comparing it to an analytical solution using calculus (much easier, but with a surprise Sage bug...).

EDIT : Now, an analytical solution using all-singin', all-dancin' calculus.

Let's consider an element of the "yellow" surface : an (infinitely) small vertical strip delimited bu the lines of equation $x=t$ and $x=t+dt$. The surface of this strip will be $(y_u(t)-y_l(t))dt$, where $(t, y_l(t)$ and $(t, y_u(t)$ satisfy $E0$. The total "yellow" area will be the summation of such elemental areas from $t=X$ to $t=r_0$, i.e. $\displaystyle{\int_X^{r_0} (y_u(t)-y_l(t)) dt}$.

Our previous work already gave us $X$. We need to compute explicitly the bounds $y_l$ and $y_u$, and the integral:

Sol0=solve(E0,y)
Y0l(x)=Sol0[0].rhs()
Y0u(x)=Sol0[1].rhs()
g0(x)=(Y0u-Y0l)(x).integrate(x)
Yellow=(g0(r_0)-g0(X)).expand().factor()

A couple notes:

• We should be able to obtain directly the definite integral integrate((Y0u-Y0l)(t), t, X, r_0) ; but we aren't : Trac#27816 will crash Sage...

• The integral, $\frac{2 \, \pi r_{0}^{3} + 4 \, r_{0}^{3} \arcsin\left(-\frac{2 \, r_{0}^{2} - r_{1}^{2}}{2 \, r_{0}^{2}}\right) - 2 \, \sqrt{4 \, r_{0}^{2} - r_{1}^{2}} r_{0} r_{1} + \frac{\sqrt{4 \, r_{0}^{2} - r_{1}^{2}} r_{1}^{3}}{r_{0}}}{4 \, r_{0}}$, bears some similarity with the geometrical solution seen above. This not random...

For the same reasons, we can get the area of the "red" surface :

Sol1=solve(E1,y)
Y1l(x)=Sol1[0].rhs()
Y1u(x)=Sol1[1].rhs()
g1(x)=(Y1u-Y1l)(x).integrate(x)
Red=(g1(X)-g1(r_0-r_1)).expand().factor()

which allows us to compute the expression of the edible fraction as a finction of the fence radius and the tether length :

EdibleFraction(r_0,r_1)=((Yellow+Red)/(pi*r_0^2)).expand().factor()

That is

$$ \frac{\frac{{\left(2 \, \pi r_{0} - 4 \, r_{0} arcsin\left(\frac{r_{1}}{2 \, r_{0}}\right) - \frac{\sqrt{4 \, r_{0}^{2} - r_{1}^{2}} r_{1}}{r_{0}}\right)} r_{1}^{2}}{r_{0}} + \frac{2 \, \pi r_{0}^{3} + 4 \, r_{0}^{3} arcsin\left(-\frac{2 \, r_{0}^{2} - r_{1}^{2}}{2 \, r_{0}^{2}}\right) - 2 \, \sqrt{4 \, r_{0}^{2} - r_{1}^{2}} r_{0} r_{1} + \frac{\sqrt{4 \, r_{0}^{2} - r_{1}^{2}} r_{1}^{3}}{r_{0}}}{r_{0}}}{400 \, \pi} $$

which looks reasonable:

The numerical answer is immediate :

sage: find_root(lambda u:EdibleFraction(10,u)-1/2, 0, 20)
11.587284730181219

One can check that trying to get the definite integral directly gives curious results :

• sympy can't give us an usable answer ;

• fricas gives a different answer, using arctan instead of arcsin (we saw that this was another possibility) ;

• giac can do the integration, but this result can't be used numerically ; this is discussed onsage-devel, but not yet filed as a bug.

• mathematica needs a bit of rewriting (it doesn't tolerate underscore in variable names...), and gives a solution close to the one given by fricas.

Next installment : numerical integration :

• Why you shouldn't.

• Why you sometimes have to.

• How not to do it.

• How to speed it

  • Use of Cython.

  • Trading speed for precision with stochastic integration.

Stay tuned...

Those variables are generated by Maxima, which has non problem generating new symbols for (e. g.) integration constants or integers denoting the multiplicity of solutions. Those variables denoting arbitrary quantities, there is no reason for them to share names...

One might note that other solvers (e. g. Sympy) may be less wasteful.

A workaround is to collect the names of the variables appearing in the solution(s) and declare them as vars if they have to be used outside the call context. Or substitute them with more significant names.

No answer to my previous question ? My hunch might have been right... Nevertheless, I'll throw a couple bones at you :

• your problem has an analytic/geometric solution (high school level).
• You can also write a numerical integral solution much more simply. (think double integration)...

As you write them, your $p_k$ are first-degrees polynomials of $n$ variables at most. Is that really what you want ?