CoCalc -- test.ipynb

⁴⁹⁵ views
License: GPL3
ubuntu2204

Kernel: SageMath 10.6

In [1]:

Out[1]:

(3.14159292035398, 3.14159265358979)

In [0]:

In [0]:

load("mc_interact.sage")

Source of the problem below

#STARTTEXT

Source of the problem below

#ENDTEXT

#STARTTEXT

This problem suggests lots of other problems. Here's one: Given two intersecting circles find the set of points $\cal{P}$ equidistant from each circle.

Coordinatize the circles so that the larger one (or either one if they are the same size) is the unit circle and the other one has center $(c,0)$ with $0 \lt c$ and radius $r \le 1$ . Since the circles intersect, we have $-1 \lt c-r \lt 1$ and $c+r \gt 1$ . From rough sketches, we can see the set $\cal{P}$ is going to be an oval shape passing through the two points of intersection of the circles. We use Sagemath's symbolic calculator to work out the equation for the set $\cal{P}$ .

#ENDTEXT

In [66]:


#STARTPROB

var('t,s,r,c')
A=vector([0,0])
B=lambda c:vector([c,0])
P=lambda s,t:(1-s)*vector([cos(t*pi/180),sin(t*pi/180)])
#The vector equation
#
Veqn2=(1-s)^2-2*P(s,t)*B(c)+c^2-(r+s)^2
show(Veqn2)
soln=solve(Veqn2,s)
show(soln)
print(soln)
sval=lambda t,c,r:-1/2*(c^2 - r^2 - 2*c*cos(1/180*pi*t) + 1)/(c*cos(1/180*pi*t) - r - 1)
def drawit(c,r):
    pic=circle((0,0),1,linestyle='dotted')
    pic+=circle((c,0),r)+point([c,0])+point([P(sval(i,c,r),i) for i in range(1,360)],color='red')
    return pic
drawit(.7,1)

#ENDPROB

Out[66]:

$\displaystyle 2 \, c {\left(s - 1\right)} \cos\left(\frac{1}{180} \, \pi t\right) + c^{2} - {\left(r + s\right)}^{2} + {\left(s - 1\right)}^{2}$

$\displaystyle \left[s = -\frac{c^{2} - r^{2} - 2 \, c \cos\left(\frac{1}{180} \, \pi t\right) + 1}{2 \, {\left(c \cos\left(\frac{1}{180} \, \pi t\right) - r - 1\right)}}\right]$

[
s == -1/2*(c^2 - r^2 - 2*c*cos(1/180*pi*t) + 1)/(c*cos(1/180*pi*t) - r - 1)
]

In [5]:

int(3./2)

Out[5]:

1

In [0]:

In [0]:

In [0]:

In [6]:

#First, we draw the picture 
G=Graphics()
sqr=[(0,0),(1,0),(1,1),(0,1)]
sqr=[vector(s) for s in sqr]
sqrpic=polygon(sqr,fill=false,color='black')
arc1=arc((.5,0),.5,sector=(0,pi.n()),color='black')
arc2=arc((1,.5),.5,sector=(pi/2,3*pi/2),color='grey')
arc3=arc((0,0),1,sector=(0,pi),color='black')
arc4=arc((.25,0),.75,sector=(0,pi))
G+=sqrpic+arc3
show(G)

Out[6]:

In [20]:

var('t,s,r')
G=circle((0,0),1);show(G)
A=vector([0,0])
B=lambda r:vector([r,0])
P=lambda s,t:(1-s)*vector([cos(t*pi/180),sin(t*pi/180)])
Veqn2=(expand((P(s,t)-B(r))*(P(s,t)-B(r))-(r+s)^2)).simplify_full()
show(Veqn2)
soln=solve(Veqn2,s)
print(soln)
sol=lambda t,r:1/2*(2*r*cos(1/180*pi*t) - 1)/(r*cos(1/180*pi*t) - r - 1)
P1=lambda t,r:P(sol(t,r),t)
#sol =lambda t,r: -1/2*r + 1/2*sqrt(r^2 - 4*r*cos(1/180*pi*t) + 2)

show(sol(60,.5))
show(sol(t,r))

Out[20]:

$\displaystyle -2 \, {\left(r + 1\right)} s + 2 \, {\left(r s - r\right)} \cos\left(\frac{1}{180} \, \pi t\right) + 1$

[
s == 1/2*(2*r*cos(1/180*pi*t) - 1)/(r*cos(1/180*pi*t) - r - 1)
]

$\displaystyle 0.200000000000000$

$\displaystyle \frac{2 \, r \cos\left(\frac{1}{180} \, \pi t\right) - 1}{2 \, {\left(r \cos\left(\frac{1}{180} \, \pi t\right) - r - 1\right)}}$

In [0]:

In [40]:

var('t,s,r')
A=vector([0,0])
B=lambda r:vector([1-r,0])
P=lambda s,t:(1-s)*vector([cos(t*pi/180),sin(t*pi/180)])
Veqn2=(expand((P(s,t)-B(r))*(P(s,t)-B(r))-(r+s)^2)).simplify_full()
show(Veqn2)
soln=solve(Veqn2,s)
print(soln)
sol=lambda t,r: ((r - 1)*cos(1/180*pi*t) - r + 1)/((r - 1)*cos(1/180*pi*t) + r + 1)
#sol=lambda t,r:1/2*(2*r*cos(1/180*pi*t) - 1)/(r*cos(1/180*pi*t) - r - 1)
#sol =lambda t,r: -1/2*r + 1/2*sqrt(r^2 - 4*r*cos(1/180*pi*t) + 2)
arcr=lambda c,r:arc(B(r),r,sector=(0,pi.n()),color='black')
show(sol(t,r))

pic=circle((0,0),1)
circs=lambda t,r:circle((1-r,0),.6)+point([1-r,0])+line([P(sol(i,r),i) for i in range(1,360)],color='red')+circle((.2,0),.8)
lns=lambda t,r:line([B(r),P(sol(t,r),t)])+line([(0,0),P(sol(t,r),t)/norm(P(sol(t,r),t))])
pic=lambda t,r:circle((0,0),1)+circs(t,r)+lns(t,r)

Out[40]:

$\displaystyle -2 \, {\left(r + 1\right)} s - 2 \, {\left({\left(r - 1\right)} s - r + 1\right)} \cos\left(\frac{1}{180} \, \pi t\right) - 2 \, r + 2$

[
s == ((r - 1)*cos(1/180*pi*t) - r + 1)/((r - 1)*cos(1/180*pi*t) + r + 1)
]

$\displaystyle \frac{{\left(r - 1\right)} \cos\left(\frac{1}{180} \, \pi t\right) - r + 1}{{\left(r - 1\right)} \cos\left(\frac{1}{180} \, \pi t\right) + r + 1}$

In [43]:

pic(70,.3)

Out[43]:

In [0]:

In [16]:

var('t,s,r,c')
A=vector([0,0])
B=lambda c,r:vector([c-r,0])
P=lambda s,t:(1-s)*vector([cos(t*pi/180),sin(t*pi/180)])
Veqn2=(expand((P(s,t)-B(c,r))*(P(s,t)-B(c,r))-(r+s)^2)).simplify_full()
show(Veqn2)
soln=solve(Veqn2,s)
show(soln)
print(soln)
sol=lambda t,c,r:-1/2*(c^2 - 2*c*r - 2*(c - r)*cos(1/180*pi*t) + 1)/((c - r)*cos(1/180*pi*t) - r - 1)
#sol=lambda t,r: ((r - 1)*cos(1/180*pi*t) - r + 1)/((r - 1)*cos(1/180*pi*t) + r + 1)
#sol=lambda t,r:1/2*(2*r*cos(1/180*pi*t) - 1)/(r*cos(1/180*pi*t) - r - 1)
#sol =lambda t,r: -1/2*r + 1/2*sqrt(r^2 - 4*r*cos(1/180*pi*t) + 2)
#arcr=lambda r:arc(B(r),r,sector=(0,pi.n()),color='black')
#show(sol(60,.5))
#show(sol(t,r))

G+arcr(.6)+point([.4,0])+line([P(sol(i,.6),i) for i in range(1,1790)],color='red')+circle((.2,0),.8)

Out[16]:

$\displaystyle c^{2} - 2 \, c r - 2 \, {\left(r + 1\right)} s + 2 \, {\left({\left(c - r\right)} s - c + r\right)} \cos\left(\frac{1}{180} \, \pi t\right) + 1$

$\displaystyle \left[s = -\frac{c^{2} - 2 \, c r - 2 \, {\left(c - r\right)} \cos\left(\frac{1}{180} \, \pi t\right) + 1}{2 \, {\left({\left(c - r\right)} \cos\left(\frac{1}{180} \, \pi t\right) - r - 1\right)}}\right]$

[
s == -1/2*(c^2 - 2*c*r - 2*(c - r)*cos(1/180*pi*t) + 1)/((c - r)*cos(1/180*pi*t) - r - 1)
]

In [0]:

In [72]:

#Checking that the radius of the small circle is 4/33
(sol(43.6028189727036,.5)-4/33.).n()

Out[72]:

-1.38777878078145e-16

In [76]:

(41/59).n(digits=100)

Out[76]:

0.6949152542372881355932203389830508474576271186440677966101694915254237288135593220338983050847457627

In [0]:

In [28]:

sol=(2 + sin(1/2*pi*t)^2  - 2*cos(1/2*pi*t))/(8+sin(1/2*pi*t)^2)
show(sol)
sl=lambda t:(2 + sin(pi*t/180)^2  - 2*cos(pi*t/180))/(8+sin(pi*t/180)^2)
show(sl(t))
sl1=(cos(pi/180*t)^2+2*cos(pi/180*t)-3)/(cos(pi/180*t)^2-9)
show(sl1)
sl1=lambda t:(cos(pi/180*t)-1)/(cos(pi/180*t)-3)
show(sl1(t))
P1=lambda t:P(sl1(t),t)
# Checking that sl==sl1 (it is)
print([sl(i).n()-sl1(i).n() for i in range(1,17)])

Out[28]:

$\displaystyle \frac{\sin\left(\frac{1}{2} \, \pi t\right)^{2} - 2 \, \cos\left(\frac{1}{2} \, \pi t\right) + 2}{\sin\left(\frac{1}{2} \, \pi t\right)^{2} + 8}$

$\displaystyle \frac{\sin\left(\frac{1}{180} \, \pi t\right)^{2} - 2 \, \cos\left(\frac{1}{180} \, \pi t\right) + 2}{\sin\left(\frac{1}{180} \, \pi t\right)^{2} + 8}$

$\displaystyle \frac{\cos\left(\frac{1}{180} \, \pi t\right)^{2} + 2 \, \cos\left(\frac{1}{180} \, \pi t\right) - 3}{\cos\left(\frac{1}{180} \, \pi t\right)^{2} - 9}$

$\displaystyle \frac{\cos\left(\frac{1}{180} \, \pi t\right) - 1}{\cos\left(\frac{1}{180} \, \pi t\right) - 3}$

[-1.32814766129474e-18, -2.71050543121376e-18, -2.10335221462188e-17, -1.51788304147971e-18, -8.67361737988404e-19, -1.30104260698261e-17, -3.42607886505419e-17, 9.54097911787244e-18, 6.93889390390723e-18, 1.56125112837913e-17, 5.20417042793042e-18, -6.93889390390723e-18, 1.38777878078145e-17, -1.73472347597681e-18, -6.93889390390723e-18, -1.38777878078145e-17]

In [29]:

sl1(43.6028189727036).n()-4/33.

Out[29]:

-1.38777878078145e-16

In [1]:

Out[1]:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 P
NameError: name 'P' is not defined

In [145]:

trac=line([(P(sl.subs(t=3*i),3*i)) for i in range(1,61)],color='red')
G+trac+arc4+point(P1(43.6028189727036),size=20)+circle(P1(43.6028189727036),sl1(43.6028189727036))

Out[145]:

In [4]:

var('h')
H=lambda h:vector([1,h])
EQ=(P(sl1(t),t)-H(h))*(P(sl1(t),t)-H(h))-(h-sl1(t))^2
sol2=solve(EQ,h)
show(EQ)
show(sol2)
EQ1=h-sol2[0].rhs()
EQ1

Out[4]:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 3
      1 var('h')
      2 H=lambda h:vector([Integer(1),h])
----> 3 EQ=(P(sl1(t),t)-H(h))*(P(sl1(t),t)-H(h))-(h-sl1(t))**Integer(2)
      4 sol2=solve(EQ,h)
      5 show(EQ)
NameError: name 'P' is not defined

In [52]:

c=lambda h:((h+2)/(4*h^2))^2
x=lambda h:(c(h)-1)/(c(h)+1) 
tsol=lambda h:acos(x(h))*180/pi.n()
tsol(.5)

Out[52]:

43.6028189727036

In [3]:

var('x,y,h,rd')
C=vector([x,y])
C1=vector([0,0])
C2=vector([1/2,0])
C3=vector([1,h])
eq1=(C-C1)*(C-C1)-(1-rd)^2
eq2=(C-C2)*(C-C2)-(1/2+rd)^2
eq3=(C-C3)*(C-C3)-(1/2-rd)^2
print(eq1)
print(eq2)
print(eq3)

Out[3]:

-(rd - 1)^2 + x^2 + y^2
-1/4*(2*rd + 1)^2 + 1/4*(2*x - 1)^2 + y^2
(h - y)^2 - 1/4*(2*rd - 1)^2 + (x - 1)^2

In [222]:

sol1=solve(eq1-eq2,rd)
show(sol1)
eq4=eq1.subs(rd=sol1[0].rhs());show(eq4)

solxy=solve(eq4,y);print(solxy)
eq5=eq3.subs(rd=sol1[0].rhs(),y=solxy[1].rhs());show((eq5.expand()));print(eq5.expand())
eq6=(h^2 - 5/3*x + 17/12)^2-( 4/3*sqrt(-2*x^2 + x + 1)*h)^2
solve(eq6,x)

xsol=lambda h: 1/4*(92*h^2 - 12*sqrt(-32*h^4 + 18)*h + 85)/(32*h^2 + 25)
ysol=lambda h:2/3*sqrt(-2*xsol(h)^2 + xsol(h) + 1)
show(xsol(1/2))
show(ysol(1/2))

Out[222]:

$\displaystyle \left[\mathit{rd} = -\frac{1}{3} \, x + \frac{1}{3}\right]$

$\displaystyle -\frac{1}{9} \, {\left(x + 2\right)}^{2} + x^{2} + y^{2}$

[
y == -2/3*sqrt(-2*x^2 + x + 1),
y == 2/3*sqrt(-2*x^2 + x + 1)
]

$\displaystyle h^{2} - \frac{4}{3} \, \sqrt{-2 \, x^{2} + x + 1} h - \frac{5}{3} \, x + \frac{17}{12}$

h^2 - 4/3*sqrt(-2*x^2 + x + 1)*h - 5/3*x + 17/12

$\displaystyle \frac{7}{11}$

$\displaystyle \frac{20}{33}$

In [220]:

G+point([(xsol(.5),ysol(.5)),(xsol(.8),ysol(.8)),(xsol(.85),ysol(.85)),(xsol(.4),ysol(.4)),(xsol(.3),ysol(.3))])

Out[220]:

In [0]:

In [127]:

sol2=h + 8*(1 - cos(1/180*pi*t) )/(cos(1/180*pi*t)^2 + 2*(cos(1/180*pi*t) - 3)*sqrt(1-cos(1/180*pi*t)^2) - 4*cos(1/180*pi*t) + 3)
show(sol2)
eq1=h + 8*(1 - x )/(x^2 + 2*(x - 3)*sqrt(1-x^2) - 4*x + 3)
eq2=h*(x^2  - 4*x + 3)+ 8*(1 - x )+ 2*h*(x - 3)*sqrt(1-x^2) 
eq3=(h*(x^2  - 4*x + 3)+ 8*(1 - x ))^2/( 2*h*(x - 3))^2-(1-x^2)
show(eq3)
solve(eq3,x)
SOLX=[x == 1,
 x == -1/15*(1/2)^(2/3)*(-I*sqrt(3) + 1)*((27*h + 16)^2/h^2 - 15*(27*h^2 + 64*h + 64)/h^2)/(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) - 1/30*(1/2)^(1/3)*(I*sqrt(3) + 1)*(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) + 1/15*(27*h + 16)/h,
 x == -1/15*(1/2)^(2/3)*(I*sqrt(3) + 1)*((27*h + 16)^2/h^2 - 15*(27*h^2 + 64*h + 64)/h^2)/(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) - 1/30*(1/2)^(1/3)*(-I*sqrt(3) + 1)*(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) + 1/15*(27*h + 16)/h,
 x == 2/15*(1/2)^(2/3)*((27*h + 16)^2/h^2 - 15*(27*h^2 + 64*h + 64)/h^2)/(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) + 1/15*(1/2)^(1/3)*(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) + 1/15*(27*h + 16)/h]

Out[127]:

$\displaystyle h - \frac{8 \, {\left(\cos\left(\frac{1}{180} \, \pi t\right) - 1\right)}}{\cos\left(\frac{1}{180} \, \pi t\right)^{2} + 2 \, \sqrt{-\cos\left(\frac{1}{180} \, \pi t\right)^{2} + 1} {\left(\cos\left(\frac{1}{180} \, \pi t\right) - 3\right)} - 4 \, \cos\left(\frac{1}{180} \, \pi t\right) + 3}$

$\displaystyle x^{2} + \frac{{\left({\left(x^{2} - 4 \, x + 3\right)} h - 8 \, x + 8\right)}^{2}}{4 \, h^{2} {\left(x - 3\right)}^{2}} - 1$

In [129]:

SOLX[3].rhs()
solx=2/15*(1/2)^(2/3)*((27*h + 16)^2/h^2 - 15*(27*h^2 + 64*h + 64)/h^2)/(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) + 1/15*(1/2)^(1/3)*(2*(27*h + 16)^3/h^3 - 45*(27*h^2 + 64*h + 64)*(27*h + 16)/h^3 - 675*(27*h^2 - 48*h - 64)/h^2 + 5760*sqrt(54*h^4 - 4/3*h^3 - 64/3*h^2 + 320/3*h + 256/3)/h^3)^(1/3) + 1/15*(27*h + 16)/h

In [132]:

x1=solx.subs(h=1/2).n();y1=sqrt(1-x1^2)
G+point([x1,y1])

Out[132]:

In [119]:

sol2

Out[119]:

[h == -2*(2*cos(1/180*pi*t)^2 + sin(1/180*pi*t)^2 - 4*cos(1/180*pi*t) + 2)/(cos(1/180*pi*t)^2 + 2*(cos(1/180*pi*t) - 3)*sin(1/180*pi*t) - 4*cos(1/180*pi*t) + 3)]

In [117]:

[((num-nu).subs(t=10*i)).n() for i in range(1,17)]

Out[117]:

[-0.122461193376488,
 -0.497006953663421,
 -1.14359353944898,
 -2.09058526543021,
 -3.36810188968151,
 -5.00000000000000,
 -6.99558882055134,
 -9.34224391575730,
 -12.0000000000000,
 -14.8989856010991,
 -17.9402334069727,
 -21.0000000000000,
 -23.9373053996508,
 -26.6040074452375,
 -28.8564064605510,
 -30.5671708188125]

In [22]:

tsol=lambda t:(1-cos(pi/180*t))/3
Q=lambda t: (1-tsol(t))*vector([cos(pi/180*t),sin(pi/180*t)])
G+line([P1(i) for i in range(1,179)],color='red')+point(P1(60))

Out[22]:

In [36]:

t=90.;((1-norm(P1(t)) )-(norm(P1(t)-vector([1/2,0])) -1/2)).n()

Out[36]:

-1.11022302462516e-16

In [0]:

In [94]:

var('x')
EQ3=(2*(x-1)*(x-3))+h*((x^2+2*x-6)*sqrt(1-x^2)-4*x+3)
EQ4=(h*(4*x-3-2)*(x-1)*(x-3))/(h*(x^2+2*x-6))^2-(1-x^2)
sl4=solve(EQ4.subs(h=1/2),x)
show(sl4)
EQ5=sl4[1].rhs()
show(EQ5)

Out[94]:

$\displaystyle \left[x = 1, 0 = x^{5} + 5 \, x^{4} - 4 \, x^{3} - 24 \, x^{2} - 22 \, x + 66\right]$

$\displaystyle x^{5} + 5 \, x^{4} - 4 \, x^{3} - 24 \, x^{2} - 22 \, x + 66$

In [0]:

In [87]:

Out[87]:

[x == -(sqrt(-x^2 + 1)*h - 2*h + sqrt(-7*h^2*x^2 + 11*h^2 - (7*h^2 + 2*h)*sqrt(-x^2 + 1) + 10*h + 4) - 4)/(sqrt(-x^2 + 1)*h + 2),
 x == -(sqrt(-x^2 + 1)*h - 2*h - sqrt(-7*h^2*x^2 + 11*h^2 - (7*h^2 + 2*h)*sqrt(-x^2 + 1) + 10*h + 4) - 4)/(sqrt(-x^2 + 1)*h + 2)]

In [0]:

In [0]:

1/2*(2*cos(1/2*pi*t)^2 + 2*sin(1/2*pi*t)^2 - sqrt(4*cos(1/2*pi*t)^2 + 3*sin(1/2*pi*t)^2 - 4*cos(1/2*pi*t) + 1) - cos(1/2*pi*t))/(cos(1/2*pi*t)^2 + sin(1/2*pi*t)^2 - 1)

In [34]:

(sqrt( cos(1/2*pi*t)^2/4+  - cos(1/2*pi*t) + 1) )

Out[34]:

3*cos(1/2*pi*t)^2 + 3*sin(1/2*pi*t)^2 - 3

In [0]:

In [0]:

In [7]:

def invert(P, C, r):
    """
    Inverts point P about a circle with center C and radius r.
    P, C: vectors of length 2 (or tuples)
    r: radius (float or symbolic)
    Returns: vector
    """
    P = vector(P)
    C = vector(C)
    v = P - C
    d2 = v.norm()^2
    if d2 == 0:
        raise ValueError("Inversion is undefined for P = C (the center of the circle)")
    return C + (r^2/d2) * v

def Invert(P,C,r):
    return [invert(p,C,r) for p in P]

In [4]:

load('quaddefs.sage')

In [7]:

var('x y a b')
A=vector((a,sqrt(1-a^2)))
B=vector((b,sqrt(1-b^2)))
eq1=x*A-y*B
eq2=((1-x)^2*(A*A))-((y-1/2)^2*(B*B))
show(eq1)
show(eq2)
eq3=eq1.subs(b=a*x/y)
show(eq3)

Out[7]:

$\displaystyle \left(a x - b y,\,\sqrt{-a^{2} + 1} x - \sqrt{-b^{2} + 1} y\right)$

$\displaystyle {\left(x - 1\right)}^{2} - \frac{1}{4} \, {\left(2 \, y - 1\right)}^{2}$

$\displaystyle \left(0,\,\sqrt{-a^{2} + 1} x - \sqrt{-\frac{a^{2} x^{2}}{y^{2}} + 1} y\right)$

In [32]:

Out[32]:

3.53846153846154

In [3]:

load('orthogtraps19.sage')

In [20]:

mod(12,5)

Out[20]:

2

In [36]:

pol=[(0,0),(1,0),(1,1)]
P=dpoly(pol)
P.perimeter

Out[36]:

3.41421356237309

In [32]:

P.perimeter

Out[32]:

<bound method dpoly.perimeter of >

In [22]:

P.perimeter=sum([norm(P.verts[i]-P.verts[mod(i+1,len(P.verts)) ]) for i in range(len(P.verts))])
P.perimeter

Out[22]:

4.00000000000000

In [0]:

In [82]:

var('x y s')
eqs=[x^2+y^2-(30*s)^2,x^2+(y-10)^2-(40*s-10)^2,(x-10)^2+y^2-(10-35*s)^2]
soln=solve(eqs,s)
eq1=eqs[0].expand()
eq2=eqs[1].expand()-eq1
eq3=eqs[2].expand()-eq1
show(eq1,eq2,eq3)

Out[82]:

$\displaystyle -900 \, s^{2} + x^{2} + y^{2} -700 \, s^{2} + 800 \, s - 20 \, y -325 \, s^{2} + 700 \, s - 20 \, x$

In [8]:

var('x')
expr=sqrt(2.)*x+sin(x/5.+1/7.)
show(expr)

Out[8]:

$\displaystyle 1.41421356237310 \, x + \sin\left(0.200000000000000 \, x + 0.142857142857143\right)$

In [0]:

In [55]:

x1=solve(eq3,x)[0].rhs()
y1=solve(eq2,y)[0].rhs()

In [85]:

eq4=(eq1.subs(x=x1,y=y1)).expand().simplify()
show(eq4)
eq4=(eq4/s^2).expand()
show(eq4)
soln=solve(eq4,s)
sol1=soln[0].rhs()
sol2=soln[1].rhs()
show([sol1,sol2])
show([sol1.n(),sol2.n()])

Out[85]:

$\displaystyle \frac{23825}{16} \, s^{4} - \frac{7875}{2} \, s^{3} + 1925 \, s^{2}$

$\displaystyle \frac{23825}{16} \, s^{2} - \frac{7875}{2} \, s + 1925$

$\displaystyle \left[-\frac{8}{953} \, \sqrt{6461} + \frac{1260}{953}, \frac{8}{953} \, \sqrt{6461} + \frac{1260}{953}\right]$

$\displaystyle \left[0.647384295014192, 1.99689692219462\right]$

In [9]:

def bash_help():
    lst=['Various useful bash commands in a jupiter worksheet:',
         "!pwd  to see the current directory."
         "!ls to list the files in the current directory.",
         "!grep 'class' 'file'   to list the classes defined in file.",
         "!grep 'class \\|def' file  to list the classes defined in file.",
         "!ls --help  to get a list of the options for ls.  Works for grep too."]
    for l in lst:
        print(l)
bash_help()

Out[9]:

Various useful bash commands in a jupiter worksheet:
!pwd  to see the current directory.!ls to list the files in the current directory.
!grep 'class' 'file'   to list the classes defined in file.
!grep 'class \|def' file  to list the classes defined in file.
!ls --help  to get a list of the options for ls.  Works for grep too.

In [3]:

!ls

Out[3]:

 2018-06-20-181627.zip
 2018-11-03-175954.x11
 2019-01-13-120733.sage
 2020-03-07-191811.term
 2020-08-13-163609.ipynb
 2020-08-13-163609.sagews
 2023-01-04-105231.ipynb
 2023-02-07-132638.term
 2023-03-30-232356.sagews
 2024-04-17-190210.ipynb
 2024-05-25-170940.ipynb
 2024-06-02-122435.term
 An_interesting_sagelet.sagews
 Bicycles.html
 Circumscribing_Ellipses.html
 CnMd.ipynb
 CnMd.sagews
 MRCA_Howard2.html
 Notebooks
 Orthogonal_Trapezoids-1.html
 Orthogonal_Trapezoids.html
 Pappus.png
'Problem 4.sagews'
 Problem_4.html
 Projective_Geometry_with_Sagemath.html
 Projective_Geometry_with_Sagemath.ipynb
 QuadSpace-backup.html
 QuadSpace.html
 QuadSpace.pdf
 QuadSpace.sagews
 QuadSpace.tex
 QuadSpace_graphs-backup.html
 QuadSpace_graphs.html
 QuadSpace_graphs.ipynb
 QuadSpace_graphs.sage
 QuadSpace_graphs2.ipynb
 QuadSpacedefs.sage
 Quadrisection.html
 Question.html
'Showing discontinuities in 2d functions.sagews'
 Symmetric_orthogonal_traps2.html
 Tangential_quads.html
 Trials
 Untitled-checkpoint.ipynb
'algebraic trig functions and their derivatives.sagews'
 april_try.sagews
 bill.png
 bn2.png
 fibo.ipynb
 file1.sage
 file2.sage
 graphs.sagews
 heat_index.html
 heat_index.ipynb
 images
 index.html
 interpolation.sagews
 levelngrad-bu.html
 levelngrad2.html
 matrix_factors.html
 may_try.sagews
 misc
 newquad.html
 newquad.ipynb
 nonconvexsolution1.png
 nonconvexsolution2.png
 o3.png
 o5.png
 o7.png
 o9.png
 openapi_test.ipynb
 orthog_experiments.sagews
 orthog_experiments2.sagews
 orthog_traps.html
 orthogonal_traps.html
 orthogtraps19.sage
 pappus.ipynb
 pappus1.html
 parallelograms.html
 pde_plots.html
 plot_tensor.ipynb
 projective_geom.ipynb
 projective_geom.pdf
 projective_geom.sagews
 projective_geom_27_0.png
 projective_geom_4_0.png
 projective_geom_7_1.png
 projective_geom_8_1.png
 quad_index.html
 quadapps12_23.sage
 quaddefs.sage
 quadrilateral_prob.html
 quadrisection20.sage
 quadtri-1.html
 rational_trig
 sagecellblank.html
 sagelets.sage
 space_of_quads.pdf
 stuff
 sudoku.ipynb
 sudoku.sagews
 symmetric_abc.html
 symotrap.sage
 symotrap2.sage
 template.html
 test-no-output.ipynb
 test.ipynb
 testing.html
 testquad.html
 testquad.ipynb
 testquadapps12_23.ipynb
 tmp1.sage
 tmp2.sage
 trapezoid_experiments.html
 twitter_probs
 uksagelets
 unihyper.ipynb
 unit_circle.ipynb

In [7]:

!grep 'class \| def' 'CnMd.sagews'

Out[7]:

class elip(object):
    def __init__(self,h=0,k=0,a=2,b=1):
class elip2(object):
    def __init__(self,h=0,k=0,a=2,b=1):
︡e93af470-3a3c-42a1-b361-53bf55b441ff︡{"stderr":"Error in lines 1-1\nTraceback (most recent call last):\n  File \"/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py\", line 1191, in execute\n    flags=compile_flags), namespace, locals)\n  File \"\", line 1, in <module>\nNameError: name 'self' is not defined\n"}︡{"done":true}
︡9b984208-1792-470c-b134-3ff214ebb796︡{"stderr":"Error in lines 1-1\nTraceback (most recent call last):\n  File \"/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py\", line 1191, in execute\n    flags=compile_flags), namespace, locals)\n  File \"\", line 1, in <module>\nNameError: name 'self' is not defined\n"}︡{"done":true}
def doit(v=input_grid(1, 4, default=[0,0,2,1], label='$[h_1,k_1,r_1,b_1]$'),u=input_grid(1, 3,default=[[2,3,1]], label='$[h_2,k_2,r_2]$',to_value=lambda x:flatten(x)),th=input_box(45,width=6,label='$\\alpha$'), th1=input_box(180,width=6,label='$\\theta$'),movie=checkbox(false),auto_update=true):
    def f(X=[1,5,1]):
    def g(X=[1,5,1]):
    def f1(X=[1,5]):
    def g1(X=[1,5]):
class hyper_traps(object):
    def __init__(self,H1=[(0,2),(0,0),60],H2=[(1,1),(-.5,1),80]):
class hyperold(object):
    def __init__(self,Inpts=[[(-1,0),(1,0)],60,70]):
class hyper(object):
    def __init__(self,Inpts=[[(-1,0),(1,0)],60,70]):
class hypgorth_trap(object):
    def __init__(self,AC1=[(1,0),(-1,0)],AC2=[(0,.8),(0,-1)],alp1=30,bet1=90,thet1=5,alp2=50, bet2=70, bp=[0,0]):
class dpoly(object):
    def __init__(self,size,decimal=true):
    def __repr__(self):
    def fit(self,P,Q):
    def opvert(self,k,Q):
    def getvops(self,Q):
    def perpvert(self,k,Q):
    def perpvert(self,k,Q):
    def get_perp(self,P,Q,T):
    def perpvert(self,k,Q):
    def getperps(self,Q):
    def get_all(self,verts):
    def getpol(self,X,A,B):
class cpoly8(object):
    def __init__(self,h,numtraps=0,frs=12,L=[0,0],Rt=([0,0],0)):
        def AR(t):
        def upcls(eps=.5,cen=P2):
class angtrap8(object):
    def __init__(self,A1=[2,0],C1=[-1,0],alp=80,bet=50,thet=20,trns=[0,0],mode=0):
        def s(t):
        def ds(t):
        def Alp(t):
class angorth8_trap(object):
    def __init__(self,A1=[1,0],C1=[-.9,0],A2=[0,.8],C2=[0,-.85],alp1=30,bet1=90,thet1=30,alp2=120,bet2=60,trs=[0,0]):
class quadra(object):
    def __init__(self,pol=[(0,0),(8,0),(12,3),(8,6)]):
        def edgy(n=0,colr='gray'):
        def arcy(n=0,colr='gray'):
class hyper2(object):
    def __init__(self,Inpts=[[(-1,0),(1,0)],60,70]):
class hyp2gorth_trap(object):
    def __init__(self,AC1=[(1,0),(-1,0)],AC2=[(0,.8),(0,-1)],alp1=30,bet1=90,thet1=5,alp2=50, bet2=70, bp=[0,0]):
class hyptrap(object):
    def __init__(self,A1=[2,0],C1=[-1,0],alp=80,bet=100,thet=20,trns=[0,0],mode=0):
        def s(t):
        def ds(t):
        def Alp(t):

In [0]:

showhelp()

In [68]:

!grep 'class\def.' .sage

In [0]:

In [60]:

!grep 'class\|def\|".' symotrap.sage

Out[60]:

def meps(n=10^2):
    """meps(n=10^2) returns n*(machine epsilon)"""
def betw(A,B,X,Y):
def chkcon(pol):
    """returns true if pol is a convex polygon oriented counterclockwise."""
def Refll(Pol=[[0,0],[1,0],[1,3]],ln=[[0,0],[1,0]]):
    """Refll(Pol=[[0,0],[1,0],[1,3]],ln=[[0,0],[1,0]]) returns the reflection of Pol about the line ln"""
def Rot(Pol=[[0,0],[1,0],[1,3]],cen=[0,0],th=30,deg=true,dec=true):
    """Rot(Pol=[[0,0],[1,0],[1,3]],cen=[0,0],th=30,deg=true) returns the rotation of Pol by th degrees about cen."""
def Trans(Pol=[[0,0],[1,0],[1,3]],v=[1,2]):
    """Trans(Pol=[[0,0],[1,0],[1,3]],v=[1,2]) returns the list [p+v for p in Pol]. That is, it translates the polygon by v."""
def Roll(Pol=[[0,0],[1,0],[1,3]],num=1,left=true):
    """Roll(Pol=[[0,0],[1,0],[1,3]],num=1,left=true) return the list rolled to the left num times"""
def area(A,B,C):
    """area(A,B,C) returns the area of the triangle with vertices A, B, and C."""
def polarea(P,warn=false):
    """polarea(P) returns (a, ars) where a is the area of the polygon with vertices P (a list of ordered pairs),
    That is, triangles P[0] P[i] P[i+1] and P[0] P[i+1] P[i+2] have disjoint interiors."""
def Refll(Pol=[[0,0],[1,0],[1,3]],ln=[[0,0],[1,0]]):
    """Refll(Pol=[[0,0],[1,0],[1,3]],ln=[[0,0],[1,0]]) returns the reflection of Pol about the line ln"""
def intsct(*args):
    """intsct(A,An,B,Bn) returns the point of intersection of the segments AAn and BBn"""
def proj(u,v):
def offset(eps,dl):
def ang(u,v,deg=true,dec=true):
    """return the angle between the vectors u and v, measured counterclockwise in degrees."""
def dir(th,degs=true):
    """returns the direction [cos(th),sin(th)]"""
def getang(v=[-2,.10],degs=true,dec=true,clkws=false):
    """if v is a vector, returns the angle the vector v makes with the positive x-axis in decimal degrees measured counterclockwise.
    if v is a list [A,B,C] of points, returns the angle ABC in degrees measured counterclockwise."""
def getdir(A=[1,0],B=[2,-1],th=-20,degs=true,dec=true):
    """Returns the direction from A towards B rotated th degrees counterclockwise if th > 0, rotated clockwise otherwise."""
class cpoly5(object):
    """c=cpoly5(h,frs=6) takes h, a dpoly object and creates a cpoly5 object with lots of attributes, the most important one of which is
    c.trap a list of the angorth2_trap objects constructed from h."""
    def __init__(self,h,frs=12,L=[0,0],Rt=([0,0],0)):
        def upcls(eps=.5,cen=P2):
def toreal(x):
def concyclic(m=3,lst=[]):
    """concyclic(n) returns the the vertices P of the regular n-gon centered at (0,0) listed counterclockwise starting at P[0]=(1,0)"""
class Xm(object):
    def __init__(self,m=1,n=1,mode=0):
        def info():
class XM(object):
    def __init__(self,m=1,n=1,mode=0):
        def info():
def upcls(pic,cen=[0,0],ep=1,axes=false,bord=false):
def view(pic,ep=.1,cen=[0,0],axes=true,win=false,bord=false,figsize=''):
class angtrap4(object):
    def __init__(self,A1=[2,0],C1=[-1,0],alp=80,bet=50,thet=20,trns=[0,0],mode=0):
        def s(t):
        def ds(t):
        def Alp(t):
class angorth5_trap(object):
    def __init__(self,A1=[1,0],C1=[-.9,0],A2=[0,.8],C2=[0,-.85],alp1=30,bet1=90,thet1=30,alp2=120,bet2=60,trs=[0,0]):
        def chk():
            def chk_view(eps=.5,t0=.9,dt=.05,sz=15,fz=10):
            def info():
def Unit(ang=45,deg=true):
def getsym2_trap(A1=[1,0],C1=[-1,0],alp1=100,bet1=60,thet1=29,rp=1/2,up=true):

In [0]:

In [15]:

!ls

Out[15]:

 2018-06-20-181627.zip
 2018-11-03-175954.x11
 2019-01-13-120733.sage
 2020-03-07-191811.term
 2020-08-13-163609.ipynb
 2020-08-13-163609.sagews
 2023-01-04-105231.ipynb
 2023-02-07-132638.term
 2023-03-30-232356.sagews
 2024-04-17-190210.ipynb
 2024-05-25-170940.ipynb
 2024-06-02-122435.term
 An_interesting_sagelet.sagews
 Bicycles.html
 Circumscribing_Ellipses.html
 CnMd.ipynb
 CnMd.sagews
 MRCA_Howard2.html
 Notebooks
 Orthogonal_Trapezoids-1.html
 Orthogonal_Trapezoids.html
 Pappus.png
'Problem 4.sagews'
 Problem_4.html
 Projective_Geometry_with_Sagemath.html
 Projective_Geometry_with_Sagemath.ipynb
 QuadSpace-backup.html
 QuadSpace.html
 QuadSpace.pdf
 QuadSpace.sagews
 QuadSpace.tex
 QuadSpace_graphs-backup.html
 QuadSpace_graphs.html
 QuadSpace_graphs.ipynb
 QuadSpace_graphs.sage
 QuadSpace_graphs2.ipynb
 QuadSpacedefs.sage
 Quadrisection.html
 Question.html
'Showing discontinuities in 2d functions.sagews'
 Symmetric_orthogonal_traps2.html
 Tangential_quads.html
 Trials
 Untitled-checkpoint.ipynb
'algebraic trig functions and their derivatives.sagews'
 april_try.sagews
 bill.png
 bn2.png
 fibo.ipynb
 file1.sage
 file2.sage
 graphs.sagews
 heat_index.html
 heat_index.ipynb
 images
 index.html
 interpolation.sagews
 levelngrad-bu.html
 levelngrad2.html
 matrix_factors.html
 may_try.sagews
 misc
 newquad.html
 newquad.ipynb
 nonconvexsolution1.png
 nonconvexsolution2.png
 o3.png
 o5.png
 o7.png
 o9.png
 openapi_test.ipynb
 orthog_experiments.sagews
 orthog_experiments2.sagews
 orthog_traps.html
 orthogonal_traps.html
 orthogtraps19.sage
 pappus.ipynb
 pappus1.html
 parallelograms.html
 pde_plots.html
 plot_tensor.ipynb
 projective_geom.ipynb
 projective_geom.pdf
 projective_geom.sagews
 projective_geom_27_0.png
 projective_geom_4_0.png
 projective_geom_7_1.png
 projective_geom_8_1.png
 quad_index.html
 quadapps12_23.sage
 quaddefs.sage
 quadrilateral_prob.html
 quadrisection20.sage
 quadtri-1.html
 rational_trig
 sagecellblank.html
 sagelets.sage
 space_of_quads.pdf
 stuff
 sudoku.ipynb
 sudoku.sagews
 symmetric_abc.html
 symotrap.sage
 symotrap2.sage
 template.html
 test-no-output.ipynb
 test.ipynb
 testing.html
 testquad.html
 testquad.ipynb
 testquadapps12_23.ipynb
 tmp1.sage
 tmp2.sage
 trapezoid_experiments.html
 twitter_probs
 uksagelets
 unihyper.ipynb
 unit_circle.ipynb

In [8]:

!grep 'class\|def\|"\|self.' QuadSpace_graphs.sage >stuff

In [6]:

!ls

Out[6]:

 2018-06-20-181627.zip
 2018-11-03-175954.x11
 2019-01-13-120733.sage
 2020-03-07-191811.term
 2020-08-13-163609.ipynb
 2020-08-13-163609.sagews
 2023-01-04-105231.ipynb
 2023-02-07-132638.term
 2023-03-30-232356.sagews
 2024-04-17-190210.ipynb
 2024-05-25-170940.ipynb
 2024-06-02-122435.term
 An_interesting_sagelet.sagews
 Bicycles.html
 Circumscribing_Ellipses.html
 CnMd.ipynb
 CnMd.sagews
 MRCA_Howard2.html
 Notebooks
 Orthogonal_Trapezoids-1.html
 Orthogonal_Trapezoids.html
 Pappus.png
'Problem 4.sagews'
 Problem_4.html
 Projective_Geometry_with_Sagemath.html
 Projective_Geometry_with_Sagemath.ipynb
 QuadSpace-backup.html
 QuadSpace.html
 QuadSpace.pdf
 QuadSpace.sagews
 QuadSpace.tex
 QuadSpace_graphs-backup.html
 QuadSpace_graphs.html
 QuadSpace_graphs.ipynb
 QuadSpace_graphs.sage
 QuadSpace_graphs2.ipynb
 QuadSpacedefs.sage
 Quadrisection.html
 Question.html
'Showing discontinuities in 2d functions.sagews'
 Symmetric_orthogonal_traps2.html
 Tangential_quads.html
 Trials
 Untitled-checkpoint.ipynb
'algebraic trig functions and their derivatives.sagews'
 april_try.sagews
 bill.png
 bn2.png
 fibo.ipynb
 file1.sage
 file2.sage
 graphs.sagews
 heat_index.html
 heat_index.ipynb
 images
 index.html
 interpolation.sagews
 levelngrad-bu.html
 levelngrad2.html
 matrix_factors.html
 may_try.sagews
 misc
 newquad.html
 newquad.ipynb
 nonconvexsolution1.png
 nonconvexsolution2.png
 o3.png
 o5.png
 o7.png
 o9.png
 openapi_test.ipynb
 orthog_experiments.sagews
 orthog_experiments2.sagews
 orthog_traps.html
 orthogonal_traps.html
 orthogtraps19.sage
 pappus.ipynb
 pappus1.html
 parallelograms.html
 pde_plots.html
 plot_tensor.ipynb
 projective_geom.ipynb
 projective_geom.pdf
 projective_geom.sagews
 projective_geom_27_0.png
 projective_geom_4_0.png
 projective_geom_7_1.png
 projective_geom_8_1.png
 quad_index.html
 quadapps12_23.sage
 quaddefs.sage
 quadrilateral_prob.html
 quadrisection20.sage
 quadtri-1.html
 rational_trig
 sagecellblank.html
 sagelets.sage
 space_of_quads.pdf
 sudoku.ipynb
 sudoku.sagews
 symmetric_abc.html
 symotrap.sage
 symotrap2.sage
 template.html
 test-no-output.ipynb
 test.ipynb
 testing.html
 testquad.html
 testquad.ipynb
 testquadapps12_23.ipynb
 tmp1.sage
 tmp2.sage
 trapezoid_experiments.html
 twitter_probs
 uksagelets
 unihyper.ipynb
 unit_circle.ipynb

In [0]:

In [19]:

!grep -rni 'def ' ~/sagelets/Trials/*.sagews

Out[19]:

/home/user/sagelets/Trials/Brilliant_November_Problem.sagews:50:def prob(n=30,ns=1,s=.1):
/home/user/sagelets/Trials/Brilliant_November_Problem.sagews:73:def prob(n=30,ns=1,s=.1):
/home/user/sagelets/Trials/Brilliant_November_Problem.sagews:77:def maxprob(n=30,rng=[1,3],s=.1):
/home/user/sagelets/Trials/Brilliant_November_Problem.sagews:90:def _(n=input_box(default=30,width=10),c=input_box(default=3,width=10),p=input_box(default=1/10,width=10),auto_update=false):
/home/user/sagelets/Trials/Brilliant_November_Problem.sagews:102:︡07c16dcb-c52e-4adb-921a-b0458308821a︡{"done":true,"md":"# STARTPROB\n\ndef prob(n=30,ns=1,s=.1):\n    \"\"\"prob(n,ns,s) returns the probability of exactly ns successes in n trials if the probability of success is s.\"\"\"\n    return binomial(n,ns)*s^ns*(1-s)^(n-ns)\n\ndef maxprob(n=30,rng=[1,3],s=.1):\n    sp=[0,0]\n    rng=range(rng[0],rng[1]+1)\n    if n<rng[-1]:\n        return 1\n    for k in range(rng[0],n+1):\n        spn=sum([prob(k,i,s) for i in rng])\n        if spn<sp[1]:\n            return sp\n        else:\n            sp=(k,spn)\n    return sp\n@interact(layout=dict(top=[['n','c','p']]))\ndef _(n=input_box(default=30,width=10),c=input_box(default=3,width=10),p=input_box(default=1/10,width=10),auto_update=false):\n    if c==0:\n        show('c needs to be between 1 and n.')\n        return\n    ans=maxprob(n,[1,c],p)\n    stg='If '+str(n)+' people are invited to a party and the probability of a person bringing a can of cranberries is '+str(p.n(digits=5))+' if asked,\\n'\n    stg+='then how many people should be asked to maximize the probability that at least 1 and no more than '+str(c)+' cans will be brought?\\n'\n    answer='Answer: Ask '+str(ans[0])+' people for a maximum probability of '+str(ans[1].n(digits=5))+'.'\n    print(stg)\n    show(answer)\n\n#ENDPROB"}
/home/user/sagelets/Trials/Brilliant_November_Problem.sagews:104:def prob(n=30,ns=1,s=.1):
/home/user/sagelets/Trials/Brilliant_November_Problem.sagews:108:def maxprob(n=30,rng=[1,3],s=.1):
/home/user/sagelets/Trials/PopularMechanicsProblem.sagews:65:def pol(u=-8,v=13,clr='black',symbol=false):
/home/user/sagelets/Trials/PopularMechanicsProblem.sagews:130:def pol(u=-8,v=13,clr='black',symbol=false):
/home/user/sagelets/Trials/PopularMechanicsProblem.sagews:213:def pol(u=-8,v=13,clr='black',symbol=false):
/home/user/sagelets/Trials/PopularMechanicsProblem.sagews:224:def _(t1=input_box(30,width=10),auto_update=false):
/home/user/sagelets/Trials/PopularMechanicsProblem.sagews:279:def pol(u=-8,v=13,t1=30,b=10,clr='black',symbol=false):
/home/user/sagelets/Trials/PopularMechanicsProblem.sagews:348:def pol(u=-8,v=13,b=10,clr='black',symbol=false):
/home/user/sagelets/Trials/PopularMechanicsProblem.sagews:360:def _(b=input_box(10,width=10),t1=input_box(30,width=10),auto_update=false):
/home/user/sagelets/Trials/Regular.sagews:95:def reglr(n=5,rot=0,clr='black'):
/home/user/sagelets/Trials/Regular.sagews:101:def comm(crd=[6,9],rot=0,clr=['blue','red','green','black','turquoise','magenta'] ):
/home/user/sagelets/Trials/Regular.sagews:129:def comm(crd=[6,9],rot=[0,0],clr=['blue','red','green','black','turquoise','magenta'] ):
/home/user/sagelets/Trials/Regular.sagews:170:def cs
/home/user/sagelets/Trials/Regular.sagews:224:def reglr(n=5,rot=0,clr='black'):
/home/user/sagelets/Trials/Regular.sagews:230:def comm(crd=[6,9],rot=0,clr=['blue','red','green','black','turquoise','magenta'] ):
/home/user/sagelets/Trials/Regular.sagews:258:def cscale(ps=[3,5],rot=2/3,clr=['blue','red','green','black','turquoise','magenta']):
/home/user/sagelets/Trials/Regular.sagews:277:def _(crd=input_box(default=[6,9,12,10],width=30),rot=input_box(default=2/3,width=20),auto_update=false):
/home/user/sagelets/Trials/Regular.sagews:288:def reglr(n=5,clr='black',ang=0,fill=False):
/home/user/sagelets/Trials/Regular.sagews:292:def comm(crd=[6,9],clr=['blue','red','green','black','turquoise'],angs=[0,0,0],fill=False):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:53:def stateprob(a=10,b=7):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:56:def answer(a=10,b=7):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:61:def sol():
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:65:def __(a=input_box(default=10,width=10),b=input_box(default=7,width=10),answer=input_box(width=10),solution=false,auto_update=false):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:100:def stateprob(b=56,p=10,d=6,c=5):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:103:def answer(b=56,p=10,d=6,c=5):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:106:def sol():
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:111:def  __(b=input_box(default=56,width=10),p=input_box(default=10,width=10), d=input_box(default=6, width=10),c=input_box(default=5,width=10), dogs=input_box(width=10), cats=input_box(width=10), solution=false,auto_update=false):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:164:def lst2num(input=[2,3,4]):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:170:def strip3s(lst=[3,3,3,3,2]):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:181:def curious(input=405):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:211:def _(input=input_box(default=326785,width=40)):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:247:def lst2num(input=[2,3,4]):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:253:def strip3s(lst=[3,3,3,3,2]):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:264:def curious(input=405):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:294:def _(input=input_box(default=326785,width=40)):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:313:def lst2num(input=[2,3,4]):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:319:def strip3s(lst=[3,3,3,3,2]):
/home/user/sagelets/Trials/Smullyan_puzzles.sagews:330:def curious(input=405):
/home/user/sagelets/Trials/activate-1.sagews:151:def set_axes_labels(graph, xlabel, ylabel, zlabel,minpt=true, **kwds):
/home/user/sagelets/Trials/activate-1.sagews:197:def _(t3=input_box(default=6/10,width=10),b1=input_box(default=3/10,width=10),c1=input_box(default=9/10,width=10),d1=input_box(default=1/4,width=10),e1=input_box(default=1/2,width=10)):
/home/user/sagelets/Trials/testing-1.sagews:113:def game(p=[.8,.15,.05],pv=[1,3,10],ns=[8,6]):
/home/user/sagelets/Trials/testing-1.sagews:129:def trials(p=[.8,.15,.05],pv=[1,3,10],ns=[8,3],nt=100):
/home/user/sagelets/Trials/testing-1.sagews:141:def expval(p=[.8,.15,.05],v=[1,3,10],ns=[4,4]):
/home/user/sagelets/Trials/testing-1.sagews:147:def tabval(p=[.7,.25,.05],v=[1,3,10],n=8):
/home/user/sagelets/Trials/testing-1.sagews:163:def _(p=input_box(default=[8/10,3/20,1/20],width=30),pv=input_box(default=[1,3,10],width=20),ns=input_box(default=8,width=8), ntrl=input_box(default=1000,width=10)):
/home/user/sagelets/Trials/testing-1.sagews:197:def game(p=[.8,.15,.05],pv=[1,3,10],ns=[8,6]):
/home/user/sagelets/Trials/testing-1.sagews:214:def trials(p=[.8,.15,.05],pv=[1,3,10],ns=[8,3],nt=100):
/home/user/sagelets/Trials/testing-1.sagews:228:def expval(p=[.8,.15,.05],v=[1,3,10],ns=[4,4]):
/home/user/sagelets/Trials/testing-1.sagews:234:def tabval(p=[.7,.25,.05],v=[1,3,10],n=8):
/home/user/sagelets/Trials/testing-1.sagews:262:def game(p=[.8,.15,.05],pv=[1,3,10],ns=[8,6]):
/home/user/sagelets/Trials/testing-1.sagews:278:def trials(p=[.8,.15,.05],pv=[1,3,10],ns=[8,3],nt=100):
/home/user/sagelets/Trials/testing-1.sagews:291:def expval(p=[.8,.15,.05],v=[1,3,10],ns=[4,4]):
/home/user/sagelets/Trials/testing-1.sagews:297:def tabval(p=[.7,.25,.05],v=[1,3,10],n=8):
/home/user/sagelets/Trials/testing-1.sagews:361:def maxv(p=[.8,.15,.05],v=[1,3,10],n=10):
/home/user/sagelets/Trials/testing-1.sagews:393:def f():
/home/user/sagelets/Trials/testing-1.sagews:419:def game2(p=[.8,.15,.05],pv=[1,3,10],ns=[10,6]):
/home/user/sagelets/Trials/testing-1.sagews:456:def trials2(p=[.05,.15,.8],pv=[10,3,1],ns=[8,4],nt=100):
/home/user/sagelets/Trials/testing.sagews:113:def game(p=[.8,.15,.05],pv=[1,3,10],ns=[8,6]):
/home/user/sagelets/Trials/testing.sagews:129:def trials(p=[.8,.15,.05],pv=[1,3,10],ns=[8,3],nt=100):
/home/user/sagelets/Trials/testing.sagews:141:def expval(p=[.8,.15,.05],v=[1,3,10],ns=[4,4]):
/home/user/sagelets/Trials/testing.sagews:147:def tabval(p=[.7,.25,.05],v=[1,3,10],n=8):
/home/user/sagelets/Trials/testing.sagews:163:def _(p=input_box(default=[8/10,3/20,1/20],width=30),pv=input_box(default=[1,3,10],width=20),ns=input_box(default=8,width=8), ntrl=input_box(default=1000,width=10)):
/home/user/sagelets/Trials/testing.sagews:197:def game(p=[.8,.15,.05],pv=[1,3,10],ns=[8,6]):
/home/user/sagelets/Trials/testing.sagews:214:def trials(p=[.8,.15,.05],pv=[1,3,10],ns=[8,3],nt=100):
/home/user/sagelets/Trials/testing.sagews:228:def expval(p=[.8,.15,.05],v=[1,3,10],ns=[4,4]):
/home/user/sagelets/Trials/testing.sagews:234:def tabval(p=[.7,.25,.05],v=[1,3,10],n=8):
/home/user/sagelets/Trials/testing.sagews:262:def game(p=[.8,.15,.05],pv=[1,3,10],ns=[8,6]):
/home/user/sagelets/Trials/testing.sagews:278:def trials(p=[.8,.15,.05],pv=[1,3,10],ns=[8,3],nt=100):
/home/user/sagelets/Trials/testing.sagews:291:def expval(p=[.8,.15,.05],v=[1,3,10],ns=[4,4]):
/home/user/sagelets/Trials/testing.sagews:297:def tabval(p=[.7,.25,.05],v=[1,3,10],n=8):
/home/user/sagelets/Trials/testing.sagews:361:def maxv(p=[.8,.15,.05],v=[1,3,10],n=10):
/home/user/sagelets/Trials/testing.sagews:393:def f():
/home/user/sagelets/Trials/testing.sagews:419:def game2(p=[.8,.15,.05],pv=[1,3,10],ns=[10,6]):
/home/user/sagelets/Trials/testing.sagews:456:def trials2(p=[.05,.15,.8],pv=[10,3,1],ns=[8,4],nt=100):

In [16]:

!grep --help

Out[16]:

Usage: grep [OPTION]... PATTERNS [FILE]...
Search for PATTERNS in each FILE.
Example: grep -i 'hello world' menu.h main.c
PATTERNS can contain multiple patterns separated by newlines.

Pattern selection and interpretation:
  -E, --extended-regexp     PATTERNS are extended regular expressions
  -F, --fixed-strings       PATTERNS are strings
  -G, --basic-regexp        PATTERNS are basic regular expressions
  -P, --perl-regexp         PATTERNS are Perl regular expressions
  -e, --regexp=PATTERNS     use PATTERNS for matching
  -f, --file=FILE           take PATTERNS from FILE
  -i, --ignore-case         ignore case distinctions in patterns and data
      --no-ignore-case      do not ignore case distinctions (default)
  -w, --word-regexp         match only whole words
  -x, --line-regexp         match only whole lines
  -z, --null-data           a data line ends in 0 byte, not newline

Miscellaneous:
  -s, --no-messages         suppress error messages
  -v, --invert-match        select non-matching lines
  -V, --version             display version information and exit
      --help                display this help text and exit

Output control:
  -m, --max-count=NUM       stop after NUM selected lines
  -b, --byte-offset         print the byte offset with output lines
  -n, --line-number         print line number with output lines
      --line-buffered       flush output on every line
  -H, --with-filename       print file name with output lines
  -h, --no-filename         suppress the file name prefix on output
      --label=LABEL         use LABEL as the standard input file name prefix
  -o, --only-matching       show only nonempty parts of lines that match
  -q, --quiet, --silent     suppress all normal output
      --binary-files=TYPE   assume that binary files are TYPE;
                            TYPE is 'binary', 'text', or 'without-match'
  -a, --text                equivalent to --binary-files=text
  -I                        equivalent to --binary-files=without-match
  -d, --directories=ACTION  how to handle directories;
                            ACTION is 'read', 'recurse', or 'skip'
  -D, --devices=ACTION      how to handle devices, FIFOs and sockets;
                            ACTION is 'read' or 'skip'
  -r, --recursive           like --directories=recurse
  -R, --dereference-recursive  likewise, but follow all symlinks
      --include=GLOB        search only files that match GLOB (a file pattern)
      --exclude=GLOB        skip files that match GLOB
      --exclude-from=FILE   skip files that match any file pattern from FILE
      --exclude-dir=GLOB    skip directories that match GLOB
  -L, --files-without-match  print only names of FILEs with no selected lines
  -l, --files-with-matches  print only names of FILEs with selected lines
  -c, --count               print only a count of selected lines per FILE
  -T, --initial-tab         make tabs line up (if needed)
  -Z, --null                print 0 byte after FILE name

Context control:
  -B, --before-context=NUM  print NUM lines of leading context
  -A, --after-context=NUM   print NUM lines of trailing context
  -C, --context=NUM         print NUM lines of output context
  -NUM                      same as --context=NUM
      --group-separator=SEP  print SEP on line between matches with context
      --no-group-separator  do not print separator for matches with context
      --color[=WHEN],
      --colour[=WHEN]       use markers to highlight the matching strings;
                            WHEN is 'always', 'never', or 'auto'
  -U, --binary              do not strip CR characters at EOL (MSDOS/Windows)

When FILE is '-', read standard input.  With no FILE, read '.' if
recursive, '-' otherwise.  With fewer than two FILEs, assume -h.
Exit status is 0 if any line is selected, 1 otherwise;
if any error occurs and -q is not given, the exit status is 2.

Report bugs to: [email protected]
GNU grep home page: <https://www.gnu.org/software/grep/>
General help using GNU software: <https://www.gnu.org/gethelp/>

In [0]:

In [87]:

A1=vector([x1.subs(s=sol1),y1.subs(s=sol1)])

A2=vector([x1.subs(s=sol2),y1.subs(s=sol2)])
B=vector([10,0]);C=vector([0,0]);D=vector([0,10])
pol1=polygon([A1,B,C,D],fill=false)+text('('+str(A1[0].n(digits=4))+','+str(A1[1].n(digits=4))+')',(A1.n()+vector([.5,.5]).n(digits=4)))
pol2=polygon([A2,B,C,D],fill=false,color='red')
show(pol1)

Out[87]:

In [12]:

def round_constants(expr, digits=4):
    """Rounds constants in a symbolic expression to a specified number of decimal digits.

    Args:
        expr: The symbolic expression.
        digits: The number of decimal digits to round to.

    Returns:
        The expression with constants rounded.
    """
    return expr.n(digits=digits)

# Example usage (add this to demonstrate the function)
var('x, y')
expr = sqrt(2)*x + pi*y + e
rounded_expr = round_constants(expr)
print(expr)
print(rounded_expr)

Out[12]:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[12], line 16
     14 var('x, y')
     15 expr = sqrt(Integer(2))*x + pi*y + e
---> 16 rounded_expr = round_constants(expr)
     17 print(expr)
     18 print(rounded_expr)
Cell In[12], line 11, in round_constants(expr, digits)
      1 def round_constants(expr, digits=Integer(4)):
      2     """Rounds constants in a symbolic expression to a specified number of decimal digits.
      3 
      4     Args:
   (...)
      9         The expression with constants rounded.
     10     """
---> 11     return expr.n(digits=digits)
File /ext/sage/10.6/src/sage/structure/element.pyx:903, in sage.structure.element.Element.n()
    901         0.666666666666667
    902     """
--> 903     return self.numerical_approx(prec, digits, algorithm)
    904 
    905 def _mpmath_(self, prec=53, rounding=None):
File /ext/sage/10.6/src/sage/symbolic/expression.pyx:6674, in sage.symbolic.expression.Expression.numerical_approx()
   6672     res = x.pyobject()
   6673 else:
-> 6674     raise TypeError("cannot evaluate symbolic expression numerically")
   6675 
   6676 # Important -- the  we get might not be a valid output for numerical_approx in
TypeError: cannot evaluate symbolic expression numerically

In [0]:

In [11]:

var('x y')
expr=1/3.*x+sin(y/7.)-y/10.^(-17)
show(expr.n(digits=4))

Out[11]:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[11], line 3
      1 var('x y')
      2 expr=Integer(1)/RealNumber('3.')*x+sin(y/RealNumber('7.'))-y/RealNumber('10.')**(-Integer(17))
----> 3 show(expr.n(digits=Integer(4)))
File /ext/sage/10.6/src/sage/structure/element.pyx:903, in sage.structure.element.Element.n()
    901         0.666666666666667
    902     """
--> 903     return self.numerical_approx(prec, digits, algorithm)
    904 
    905 def _mpmath_(self, prec=53, rounding=None):
File /ext/sage/10.6/src/sage/symbolic/expression.pyx:6674, in sage.symbolic.expression.Expression.numerical_approx()
   6672     res = x.pyobject()
   6673 else:
-> 6674     raise TypeError("cannot evaluate symbolic expression numerically")
   6675 
   6676 # Important -- the  we get might not be a valid output for numerical_approx in
TypeError: cannot evaluate symbolic expression numerically

In [0]:

In [7]:

load('Trials/mc_interact.sage')

In [9]:

g = sphere()
ga = set_axes_labels(g, 'X', 'Y', 'Z',color='black',fontsize=20,axes=true)
show(ga)

Out[9]:

-1.0 -1.0 -1.0

Graphics3d Object

In [4]:


pol=mapv([(29/90, 61/40),  (-29/90, 61/40),(-1, 0),(1,0)])

In [1]:

print('hi')

Out[1]:

hi

In [1]:

var('x t3 b c d e')
S=2*(1-t3)*(b*x+c*(1-x))+3*t3*(d*(1-x)+e*x)
expand(S)

Out[1]:

-2*b*t3*x + 2*c*t3*x - 3*d*t3*x + 3*e*t3*x - 2*c*t3 + 3*d*t3 + 2*b*x - 2*c*x + 2*c

In [0]:

In [24]:

A = np.array([
  [[1,2,3],    [4,5,6],    [7,8,9]],
  [[11,12,13], [14,15,16], [17,18,19]],
  [[21,22,23], [24,25,26], [27,28,29]],
  ])
B = np.array([
  [[1,2,3],    [4,5,6],    [7,8,9]],
  [[11,12,13], [14,15,16], [17,18,19]],
  [[21,22,23], [24,25,26], [27,28,29]],
  ])
C = A + B
print(C)

Out[24]:

[[[ 2  4  6]
  [ 8 10 12]
  [14 16 18]]

 [[22 24 26]
  [28 30 32]
  [34 36 38]]

 [[42 44 46]
  [48 50 52]
  [54 56 58]]]

In [2]:

from numpy import array
from numpy import tensordot
A = array([1,2])
B = array([3,4])
C = tensordot(A, B, axes=0)
print(C)

Out[2]:

[[3 4]
 [6 8]]

In [4]:

import tensorly as tly
import numpy as np

Out[4]:

/ext/sage/10.3/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/scikits/__init__.py:1: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__("pkg_resources").declare_namespace(__name__)
/ext/sage/10.3/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/scikits/__init__.py:1: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('scikits')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  __import__("pkg_resources").declare_namespace(__name__)

In [35]:

dea=tly.decomposition.CP(a,1)

In [36]:

show(a,dea)

Out[36]:

$\displaystyle \begin{array}{l} \verb|[[[0|\verb| |\verb|1]|\\ \verb| |\verb|[2|\verb| |\verb|3]]|\\ \\ \verb| |\verb|[[4|\verb| |\verb|5]|\\ \verb| |\verb|[6|\verb| |\verb|7]]]| \end{array} \begin{array}{l} \verb|Rank-[[[0|\verb| |\verb|1]|\\ \verb| |\verb|[2|\verb| |\verb|3]]|\\ \\ \verb| |\verb|[[4|\verb| |\verb|5]|\\ \verb| |\verb|[6|\verb| |\verb|7]]]|\verb| |\verb|CP|\verb| |\verb|decomposition.| \end{array}$

In [24]:

apf=tly.decomposition.CP.fit(a)
show(apf[0],apf[1])

Out[24]:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[24], line 1
----> 1 apf=tly.decomposition.CP.fit(a)
      2 show(apf[Integer(0)],apf[Integer(1)])
TypeError: DecompositionMixin.fit() missing 1 required positional argument: 'tensor'

In [5]:

a = np.array([range(8)])

a.shape = (2, 2, 2)

A = np.array(('a', 'b', 'c', 'd'))

A.shape = (2, 2)

show(a, A)

Out[5]:

$\displaystyle \begin{array}{l} \verb|[[[0|\verb| |\verb|1]|\\ \verb| |\verb|[2|\verb| |\verb|3]]|\\ \\ \verb| |\verb|[[4|\verb| |\verb|5]|\\ \verb| |\verb|[6|\verb| |\verb|7]]]| \end{array} \begin{array}{l} \verb|[['a'|\verb| |\verb|'b']|\\ \verb| |\verb|['c'|\verb| |\verb|'d']]| \end{array}$

In [16]:

a[0,1,1]

Out[16]:

3

In [49]:

pol=[(-1, 0),(1, 0),(1,2),(-1,2)]

In [50]:

P=getpoly(pol)

Out[50]:

trap  0
rt0 =  0
rt1 =  0.3333333332593333

In [9]:

#save(P.quad_pic,"fivequads.png")

In [60]:

plot(P.C.AR,(0,pi.n()/2),aspect_ratio=.020)

Out[60]:

In [18]:

P.C.QAr+plot(0,(0,pi.n()),title="% deviation from 1/4 area")

Out[18]:

In [3]:

pol,pic=poly4(0.48739496, 0.9091, 1, 0.24369748)
pol=mapv([(29/90, 61/40),  (-29/90, 61/40),(-1, 0),(1,0)])
A=pol[0];B=pol[1];C=pol[2];D=pol[3]
X=intsct([A,C],[B,D]);X
AX=norm(A-X)
BX=norm(B-X)
DX=norm(D-X)
angBXA=getang([A-C,(0,0),B-D])
angBXA,AX,BX,DX

Out[3]:

(81.8526016111049, 0.491876458806038, 0.491876458806038, 1.52651314801874)

In [1]:

load('~/Forgons/quadapps.sage')

In [20]:

T=angorth5_trap(A1=[1,0],C1=[-.9,0],A2=[0,.8],C2=[0,-.85],alp1=30,bet1=90,thet1=30,alp2=120,bet2=60)

In [30]:

pol=T.T2.T1+T.T1.T1

Out[30]:

[(0.000000000000000, 0.800000000000000),
 (-0.591887121551798, 1.14172618895780),
 (0.000000000000000, -0.850000000000000),
 (0.418527397346263, -0.608363094478902)]

In [31]:

T.T1.T1

Out[31]:

[(1, 0),
 (0.443502884254440, 0.321293759578949),
 (-0.900000000000000, 0.000000000000000),
 (-0.900000000000000, -0.454377992302390)]

In [101]:

def getparams(pol,checkall=false):
    """returns the parameters for the member of M space which is similar to pol"""
    def getp(pol):
        pol=mapv(pol)
        A=pol[0];B=pol[1];C=pol[2];D=pol[3]
        X=intsct([A,C],[B,D])
        ax=norm(X-A);cx=norm(X-C);bx=norm(X-B);dx=norm(X-D)
        th=getang([A-C,(0,0),B-D])
        t1=2*ax/(ax+cx);t2=th/90;t3=(bx+dx)/(ax+cx);t4=bx/(bx+dx)
        if  t1<=1+meps(10^5) and t2<= 1+meps(10^5)  and t3>0 and t3<=1+meps(10^5)  and t4<1:
            if ((rclose(t1,1) or rclose(t2,1)) and t4 > .5+meps(10^5)) or (rclose(t3,1) and ax > dx):
                return 'nope'
            return [t1,t2,t3,t4]
        return 'nope'
    allst=[]
    for i in range(4):
        x=getp(Roll(pol,i))
        if type(x)==list and not checkall:
            return x
        else:
            allst+=[x];
    #pol.reverse()
    pol=Refll(pol,[pol[0],pol[2]],sense=true)
    for i in range(4):
        x=getp(Roll(pol,i))
        x=getp(Roll(pol,i))
        if type(x)==list and not checkall:
            return x
        else:
            allst+=[x];
    return allst

def classify_quad(t1,t2,t3,t4,full=5):
    '''takes a set of parameters or vertices for a quadrilateral and classifies
    it.'''
    eps=meps(10^5)
    cls=[]
    if max([abs(t1-1),abs(t2-1),abs(t3-1),abs(t4-1/2)])<eps:
        cls+=['the square, dimension 0']
        if len(cls)==full:
            return ": ".join(cls)
    if max([abs(t3-1),abs(t1-1),abs(t4-1/2)])<eps :
        cls+=['rectangle, dimension 1']
        if len(cls)==full:
            return ": ".join(cls)
    if max([abs(t1-1),abs(t2-1),abs(t4-1/2)])<eps:
        cls+=['rhombus, dimension 1']
        if len(cls)==full:
            return ": ".join(cls)
    if max([abs(t4-1/2),abs(t2-1),abs(t3-sqrt(t1*(2-t1)))])<eps:
        cls+=['right kite, dimension 1']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t3-1)<eps and abs(t1-1)<eps:
        cls+=['equidiagonal bimajor quadrilateral, dimension 2']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t1-1)<eps and abs(t4-1/2)<eps:
        cls+=['parallelogram, dimension 2']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t4-1+t1/2)<eps and abs(t3-1)<eps:
        cls+=['isosceles short trapezoid, dimension 2']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t4-t1/2)<eps and abs(t3-1)<eps:
        cls+=['isosceles tall trapezoid, dimension 2']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t4-1+t1/2)<eps:
        cls+=['short trapezoid, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t4-t1/2)<eps:
        cls+=['tall trapezoid, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t1-1+sqrt(1-4*t3^2*t4*(1-t4)))<eps:
        cls+=['cyclic quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t2-1)<eps:
        cls+=['orthodiagonal quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t3-1)<eps:
        cls+=['equidiagonal quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t1-1)<eps:
        cls+=['bimajor quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(t4-1/2)<eps:
        cls+=['biminor quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
# Now check tangential or extangential
    pol=poly4(t1,t2,t3,t4)[0]
    ab=norm(pol[1]-pol[0]);bc=norm(pol[1]-pol[2]);cd=norm(pol[2]-pol[3]);da=norm(pol[3]-pol[0])
    if abs(ab+cd-(bc+da))<eps:
        cls+=['tangential quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(ab+bc-(cd+da))<eps:
        cls+=['major extangential quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if abs(da+ab-(cd+bc))<eps:
        cls+=['minor extangential quadrilateral, dimension 3']
        if len(cls)==full:
            return ": ".join(cls)
    if t4<1/2*t1:
        cls+=['tall quadrilateral, dimension 4']
        return ": ".join(cls)
    if t4> 1-1/2*t1:
        cls+=['short quadrilateral, dimension 4']
        return ": ".join(cls)
    if t4> 1/2*t1 and t4 < 1-1/2*t1:
        cls+=['middle quadrilateral, dimension 4']
        return ": ".join(cls)
    return ": ".join(cls)

In [22]:

help(poly4)

Out[22]:

Help on function poly4 in module __main__:

poly4(t1, t2, t3, t4, A=(1, 0), C=(-1, 0), labels=False, model=True)
    returns the quadrilateral ABCD defined by the parameters t1 t2 t3 t4.

In [5]:

pol,pic=poly4(.5,.7,.9,.3,labels=true)
show(pic,figsize=4)

Out[5]:

In [43]:

T.picT.T1.D1

Out[43]:

(-0.900000000000000, -0.454377992302390)

In [36]:

help(window)

Out[36]:

Help on function window in module __main__:

window(L=0, U=0, ep=0.100000000000000, sz=0, cen=(0, 0), figsize=5, color='blue', alpha=0, axes=False, aspect_ratio=1)

In [32]:

T=angorth5_trap(A1=[1,0],C1=[-.9,0],A2=[0,.8],C2=[0,-.85],alp1=30,bet1=90,thet1=30,alp2=120,bet2=60,trs=[0,0])

In [48]:

(T.pic1(0)+point([T.T1.D1])+window(ep=.1))

Out[48]:

In [0]:

In [13]:

param=mapn(getparams(pol),digs=10);param

Out[13]:

[0.4873949580, 0.9094733512, 1.000000000, 0.2436974790]

In [6]:

@interact(layout=dict(top=[['parms',  'classify','quads','position']]))
def _(parms=input_box(default=[1,1,1,1/2],width=20), quads=checkbox(default=false,label='Quadrisection(s)'), classify=checkbox(default=false,label='Classsify'),auto_update=false):
    global pol,D
    t1=parms[0];t2=parms[1];t3=parms[2];t4=parms[3]
    pol,pic=poly4(t1,t2,t3,t4)

    Pic=pic
    if quads  and t2>10^(-12):
        D=dpoly(pol)
        Pic+=getquads2(D)
    if classify:
        tis=classify_quad(t1,t2,t3,t4)
        Pic+=text(tis,(0,1.5))
    show(Pic,figsize=8,axes=false,aspect_ratio=1)

Out[6]:

Manual interactive function <function _ at 0x7f73b8e17740> with 3 widgets
  parms: EvalText(value='[1, 1, 1, 1…

In [2]:

pic.show(title='bill was here')

Out[2]:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 1
----> 1 pic.show(title='bill was here')
NameError: name 'pic' is not defined

In [18]:

def classify_quad(t1,t2,t3,t4,full=1):
    '''takes a set of parameters or vertices for a quadrilateral and classifies
    it.'''
    eps=meps(10^4)
    cls=[]
    if max([abs(t1-1),abs(t2-1),abs(t3-1),abs(t4-1/2)])<eps:
        cls+=['the square, dimension 0']
        if len(cls)==full:
            return cls
    if max([abs(t3-1),abs(t1-1),abs(t4-1/2)])<eps :
        cls+=['rectangle, dimension 1']
        if len(cls)==full:
            return cls        
    if max([abs(t1-1),abs(t2-1),abs(t4-1/2)])<eps:
        cls+=['rhombus, dimension 1']
        if len(cls)==full:
            return cls        
    if max([abs(t4-1/2),abs(t2-1),abs(t3-sqrt(t1*(2-t1)))])<eps:
        cls+=['right kite, dimension 1']
        if len(cls)==full:
            return cls
    if abs(t3-1)<eps and abs(t1-1)<eps:
        cls+=['equidiagonal bimajor quadrilateral, dimension 2']
        if len(cls)==full:
            return cls
    if abs(t1-1)<eps and abs(t4-1/2)<eps:
        cls+=['parallelogram, dimension 2']
        if len(cls)==full:
            return cls
    if abs(t4-1+t1/2)<eps and abs(t3-1)<eps:
        cls+=['isosceles short trapezoid, dimension 2']
        if len(cls)==full:
            return cls
    if abs(t4-t1/2)<eps and abs(t3-1)<eps:
        cls+=['isosceles tall trapezoid, dimension 2']
        if len(cls)==full:
            return cls
    if abs(t4-1+t1/2)<eps:
        cls+=['short trapezoid, dimension 3']
        if len(cls)==full:
            return cls
    if abs(t4-t1/2)<eps:
        cls+=['tall trapezoid, dimension 3']
        if len(cls)==full:
            return cls
    if abs(t1-1+sqrt(1-4*t3^2*t4*(1-t4)))<eps:
        cls+=['cyclic quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
    if abs(t2-1)<eps:
        cls+=['orthodiagonal quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
    if abs(t3-1)<eps:
        cls+=['equidiagonal quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
    if abs(t1-1)<eps:
        cls+=['bimajor quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
    if abs(t4-1/2)<eps:
        cls+=['biminor quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
# Now check tangential or extangential
    pol=poly4(t1,t2,t3,t4)[0]
    ab=norm(pol[1]-pol[0]);bc=norm(pol[1]-pol[2]);cd=norm(pol[2]-pol[3]);da=norm(pol[3]-pol[0])
    if abs(ab+cd-(bc+da))<eps:
        cls+=['tangential quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
    if abs(ab+bc-(cd+da))<eps:
        cls+=['major extangential quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
    if abs(da+ab-(cd+bc))<eps:
        cls+=['minor extangential quadrilateral, dimension 3']
        if len(cls)==full:
            return cls
    if t4<1/2*t1:
        cls+=['tall quadrilateral, dimension 4']
        return ": ".join(cls)
    if t4> 1-1/2*t1:
        cls+=['short quadrilateral, dimension 4']
        return cls
    if t4> 1/2*t1 and t4 < 1-1/2*t1:
        cls+=['middle quadrilateral, dimension 4']
        return cls
    return cls

In [20]:

cls=classify_quad(.487400000000000, 0.910000000000000, 1, 0.2437,full=10)
for c in cls:
    show(c)

Out[20]:

In [0]:

In [227]:

Eqn1=eqn.subs(k=sol[1][0].rhs(),a=sol[1][1].rhs(),b=sol[1][2].rhs())
Eqn2=eqn.subs(k=sol[2][0].rhs(),a=sol[2][1].rhs(),b=sol[2][2].rhs())
Eqn3=eqn.subs(k=sol[3][0].rhs(),a=sol[3][1].rhs(),b=sol[3][2].rhs())
Eqn4=eqn.subs(k=sol[4][0].rhs(),a=sol[4][1].rhs(),b=sol[4][2].rhs())

In [244]:

var('x1,y1,z1,w1')
eqn=(k - y)^2/b^2 + x^2/a^2 - 1

In [208]:

Eqn1=eqn.subs( k= (-440751549542831198201161/5533220865273519395330948), a = 1/92686284036895753342881653781412080*sqrt(22419555699070124546480457610617294529)*sqrt(386192850153732305595340224089217), b = 1054057/5533220865273519395330948*sqrt(22419555699070124546480457610617294529))

In [213]:

Eqn1=22244708168854980802291596907538899200/22419555699070124546480457610617294529.*x^2 + 1/24908937055963900416841725604336190477515580448721.*(5533220865273519395330948.*y + 440751549542831198201161.)^2 - 1.00000000000000 
Eqn1

Out[213]:

0.9922011152869369424814976667934448505*x^2 + (4.014623336809837343425993029288696660772513783868e-50)*(5.53322086527351939533095e24*y + 4.4075154954283119820116e23)^2 - 1.00000000000000

In [231]:

ell=implicit_plot(Eqn1,(-1.1,1.1),(-1,1))

In [0]:

In [0]:

In [246]:

t1,t2,t3,t4=.5,.5,.5,.45
pol,pic=poly4(t1,t2,t3,t4)

quads=true;classify=true
if classify:
    cls=classify_quad(t1,t2,t3,t4,full=5)
    print(cls)
clrs=['black','magenta','blue','green','red']
D=dpoly(pol)
Pic=D.pic
C=cpoly5(D)
for i in range(len(C.quads)):
    q=C.quads[i]
    Pic+=line(q[0],color=clrs[i])+line(q[1],color=clrs[i])
Pic+=Pic+line(C.oval,color='red')
for i in range(len(C.trap)):
    tr=C.trap[i]
    Pic+=line([tr.T1.X(0),tr.T1.Y(0)],linestyle='dotted',color=clrs[i])+line([tr.T2.X(0),tr.T2.Y(0)],linestyle='dotted',color=clrs[i])

Gth=Graphics()
Ov=Graphics()
Ar=C.AR;
for i in range(len(C.trap)):
    tr=C.trap[i];TH=tr.THET;Th=tr.Thet
    f=lambda t:Ar((t-TH)/Th+i)
    Gth+=plot(f,(TH,TH+Th),color=clrs[i])+point([(TH,f(TH))],color=clrs[i])
    Ov+=line([tr.Cn(j/10) for j in range(11)],color=clrs[i])
show(Pic+Ov,figsize=5,axes=false,aspect_ratio=1)
show(Gth,figsize=3)
show(Ov,figsize=3)

Out[246]:

middle quadrilateral, dimension 4
trap  0
rt1 =  0.17968012009352616
trap  1
done
trap  2
done
trap  3
done

In [0]:

In [158]:

var('x,y,z,w,t1,t2,t3,t4')
A=vector([1,0]);B=vector([x,y])
C=vector([-1,0]);D=vector([z,w])
P=intsct([A,B],[C,D],'nodec')
Q=intsct([A,D],[C,B],'nodec')
show(P);show(Q)

Out[158]:

$\displaystyle \left(-\frac{{\left(w {\left(z + 1\right)} - w {\left(z - 1\right)}\right)} {\left(x - 1\right)}}{w {\left(x - 1\right)} - y {\left(z + 1\right)}} + 1,\,-\frac{{\left(w {\left(z + 1\right)} - w {\left(z - 1\right)}\right)} y}{w {\left(x - 1\right)} - y {\left(z + 1\right)}}\right)$

$\displaystyle \left(\frac{{\left({\left(x + 1\right)} y - {\left(x - 1\right)} y\right)} {\left(z - 1\right)}}{w {\left(x + 1\right)} - y {\left(z - 1\right)}} + 1,\,\frac{{\left({\left(x + 1\right)} y - {\left(x - 1\right)} y\right)} w}{w {\left(x + 1\right)} - y {\left(z - 1\right)}}\right)$

In [166]:

var('t1 t2 t3 t4')
B=vector([t1,0])+2*t3*t4*vector([t2,sqrt(1-t2^2)])
D=vector([t1,0])-2*t3*(1-t4)*vector([t2,sqrt(1-t2^2)])
P1=intsct([A,B],[C,D],'nodec')
p=P1[0];q=P1[1]
P2=vector([p.numerator()/p.denominator(),q.numerator().factor()/q.denominator()])
Q1=intsct([A,D],[C,B],'nodec')
r=Q1[0];s=Q1[1]
Q2=vector([r.numerator()/r.denominator(),s.numerator().factor()/s.denominator()])
show(P2)
show(Q2)

Out[166]:

$\displaystyle \left(\frac{4 \, t_{2} t_{3} t_{4}^{2} - 4 \, t_{2} t_{3} t_{4} + 2 \, t_{1} t_{4} - t_{1} + 1}{t_{1} + 2 \, t_{4} - 1},\,\frac{4 \, \sqrt{-t_{2}^{2} + 1} t_{3} {\left(t_{4} - 1\right)} t_{4}}{t_{1} + 2 \, t_{4} - 1}\right)$

$\displaystyle \left(-\frac{4 \, t_{2} t_{3} t_{4}^{2} - 4 \, t_{2} t_{3} t_{4} + 2 \, t_{1} t_{4} - t_{1} - 1}{t_{1} - 2 \, t_{4} + 1},\,-\frac{4 \, \sqrt{-t_{2}^{2} + 1} t_{3} {\left(t_{4} - 1\right)} t_{4}}{t_{1} - 2 \, t_{4} + 1}\right)$

In [163]:

P2=vector([p.numerator()/p.denominator(),q.numerator().factor()/q.denominator()]);show(P2)

Out[163]:

In [44]:

tr=C.trap[0];x1=.15;x2=.75
plot(tr.Ar1,(x1,x2),color='green')

Out[44]:

In [0]:

In [0]:

In [93]:

def backf(f,a,b):
    c=b-a
    def bf(t):
        n=floor(t/c)
        r=(t-n*c)
        if mod(n,2)==0:
            return f(a+r)
        else:
            return f(b-r)
    return bf

In [141]:

f=lambda x:x*(1-x)*(x-2)
a=0;b=2
g=backf(f,a,b)
plot(g,(0,6),figsize=4)

Out[141]:

In [62]:

In [53]:

Out[53]:

1

In [0]:

In [0]:

In [0]:

In [0]:

#TceGraph=Graphics()
#for i in range(6):
#    TceGraph+=implicit_plot3d(tce.subs(t2=i/5.),(t4,.01,1),(t3,.001,1),(t1,.001,1),color='magenta',frame=false,alpha=.3,mesh=true)
#TceGraph

In [0]:

#CycGraph=Graphics()
#for i in range(6):
#    CycGraph+=plot3d(cyc+i/5.,(t4,.0,1),(t3,0,1),color='turquoise',frame=false,alpha=1,mesh=true)+plot3d(cyc+1,(t4,.0,1),(t3,0,1),color='turquoise',frame=false,alpha=1,mesh=true)
#CycGraph

In [0]:

ln(220.),ln(406.),ln(3660.),ln(3685.)

In [0]:

def gparms(pol):
    X=intsct([pol[0],pol[2]],[pol[1],pol[3]]);scl=2/norm(pol[0]-pol[2])
    return scl*(norm(pol[0]-X)),getang([pol[1],X,pol[0]])/90,scl/2*(norm(pol[1]-pol[3])),norm(pol[1]-X)/norm(pol[1]-pol[3])
    
def tancirc(alp=90,bet=30,gam=40,r=1):
    cnvt=pi.n()/180
    A=vector([1,0]);C=vector([-1,0])
    if alp+bet+gam>180:
        tq=true
    else:
        tq=false
    if tq:
        br=-alp/2*cnvt;tr=alp/2*cnvt;tl=(alp/2+bet)*cnvt;bl=(-alp/2-gam)*cnvt
    else:
        br=(180+alp/2+gam)*cnvt;tr=(180-bet-alp/2)*cnvt;tl=(180-alp/2)*cnvt;bl=(180+alp/2)*cnvt
    degs=[br,tr,tl,bl]
    tpts=[sqrt(2.)/2*vector([cos(d),sin(d)]) for d in degs];
    tlns=[[vector([t[1],-t[0]])+t,t] for t in tpts];
    if tq:
        pol=[intsct(tlns[0],tlns[1]),intsct(tlns[1],tlns[2]),intsct(tlns[2],tlns[3]),intsct(tlns[3],tlns[0])]
    else:
        pol=[intsct(tlns[2],tlns[3]),intsct(tlns[3],tlns[1]),intsct(tlns[1],tlns[0]),intsct(tlns[0],tlns[2])]
    scl=2/norm(pol[0]-pol[2])
    rad=scl*sqrt(2.)/2
    pol1=[scl*(p-pol[0])+A for p in pol]
    tpts1=[scl*(t-pol[0])+A for t in tpts]
    cen1=-scl*pol[0]+A
    ang=180-getang(pol1[2]-pol1[0]);ang=abs(ang)
    pol2=Rot(pol1,pol1[0],sign(pol1[2][1])*ang)
    tpts2=Rot(tpts1,pol1[0],sign(pol1[2][1])*ang);
    cen2=Rot([cen1],pol1[0],sign(pol1[2][1])*ang)[0]
    t1,t2,t3,t4=gparms(pol2)
    #parms=[t1,t2,t3,t4];print(parms)
    if t3>1:
        print('t 3> 1')
        scl=1/t3;rad=scl*rad
        parms=[t1,t2,t3,t4]
        Vct=vector([cos(t2/2*pi.n()/2),sin(t2/2*pi.n()/2)])
        X=intsct([pol[0],pol[2]],[pol[1],pol[3]])
        pol2=Scale(Roll(Reverse(Refll(pol2,[X,X+Vct])),-2),scl)
        trns=A-pol2[0]
        pol2=Trans(pol2,trns)
        tpts2=Trans(Scale(Roll(Reverse(Refll(tpts2,[X,X+Vct])),-1),scl),trns)
        tpts2[1],tpts2[2]=tpts2[2],tpts2[1]
        cen2=Trans(Scale(Refll([cen2],[X,X+Vct]),scl),trns)[0]
        t1,t2,t3,t4=gparms(pol2)
    if t2>1:
        print('t2 > 1')
        pol2=Refll(pol2,[(0,0),(1,0)])
          
        return pol2,poly4(t1,t2,t3,t4)[1],[t1,t2,t3,t4,cen2,rad,tpts2]
    B=pol2[1];D=pol2[3];X=intsct([A,C],[B,D])
    D=A+r*(D-A)
    pol2[3]=D
    a=norm(A-B);b=norm(B-C);c=norm(C-D);d=norm(D-A)
    t1=norm(A-X);t2=getang([A-C,(0,0),B-D])/90;t3=norm(B-D)/2;t4=norm(B-X)/norm(B-D)
    clr='black'
    parms=[t1,t2,t3,t4]
    if not test_parms(t1,t2,t3,t4):
        clr='red'
    pic=circle(cen2,rad,color=clr,linestyle='dotted')+point(tpts2+[cen2],color=clr)+line([tpts2[0],cen2],color=clr,linestyle='dotted')
    pic+=line([tpts2[1],cen2],color=clr,linestyle='dotted')
    pic+=line([A,C],color=clr,linestyle='dotted')+line([pol2[1],pol2[3]],color=clr,linestyle='dotted')+line([tpts2[2],cen2],color=clr,linestyle='dotted')
    pic+=line([tpts2[3],cen2],color=clr,linestyle='dotted')+point(pol2+[intsct([pol2[0],pol2[2]],[pol2[1],pol2[3]])],color=clr)
    if not tq:
        var('k x1 x2')
        pic+=line([pol2[0],tpts2[2]],color=clr,linestyle='dotted')+line([pol2[1],tpts2[1]],color=clr,linestyle='dotted')
        pic+=line([pol2[2],tpts2[0]],color=clr,linestyle='dotted')+line([pol2[0],tpts2[3]],color=clr,linestyle='dotted')
    pic+=polygon(pol2,fill=false,color=clr)
    return pol2,pic,parms

In [0]:

pol2,pic,parms=tancirc(150,30,70,1)
show(pic)
t1,t2,t3,t4=gparms(pol2)
pol3,pic1=poly4(2-t1,t2,t3,1-t4)
show(pic1)
print(parms);print(2-t1,t2,t3,1-t4)
show(classify_quad(2-t1,t2,t3,1-t4))

In [0]:

pol,pic=poly4(1/3,1/2,1,1-1/6)
pol1,pic1=poly4(1/3,1/2,1,1/6)
show(gtsidesangs(pol))
show(gtsidesangs(pol1))
show(pic);show(pic1)

In [0]:

def good_parms(t1,t2,t3,t4):
    pol,pic=poly4(t1,t2,t3,t4)

    if t3>1: #reflect about bisector of angle BXA
        print('t3>1')
        t1,t2,t3,t4=2*t4,t2,1/t3,t1/2
        pol1 ,pic1 =poly4(t1,t2,t3,t4)
        pic+=pic1
        show([t1,t2,t3,t4])
        show(gtsidesangs(pol1))
        pol=pol1
    if t2>1 and t1<=1: #reflect about x-axis
        print('t2>1 and t1<=1')
        t1,t2,t3,t4=t1,2-t2,t3,1-t4#;show(poly4(t1,t2,t3,t4)[1]);return
        pol2 ,pic2 =poly4(t1,t2,t3,t4) 
        pic+=pic2
        show([t1,t2,t3,t4])
        show(gtsidesangs(pol2))
        pol=pol2
    if t1>1+meps(10^5): #reflect about y-axis
        print('t1>1')
        t1,t2,t3,t4=2-t1,2-t2,t3,t4
        pol3,pic3=poly4(t1,t2,t3,t4)
        pic+=pic3
        show([t1,t2,t3,t4])
        show(gtsidesangs(pol3))
        pol=pol3

    show(pic,axes=false,aspect_ratio=1)
    return pol,[t1,t2,t3,t4]

def transform(t1,t2,t3,t4,trans=''):
    pol,pic=poly4(t1,t2,t3,t4)
    #show(pic) 
    show([t1,t2,t3,t4])
    show(gtsidesangs(pol))
    if trans=='RBXA': #reflect about bisector of angle BXA
        t1,t2,t3,t4=2*t4,t2,1/t3,t1/2
        show('Reflect about bisector of angle BXA: t1,t4=2*t4,t1/2')
        pol1 ,pic1 =poly4(t1,t2,t3,t4)
        #pic+=(polygon(pol1,fill=false,color='purple')+line([pol1[0],pol1[2]])+line([pol1[1],pol1[3]]))
        pic+=pic1
        #show(pic)
        show([t1,t2,t3,t4])
        show(gtsidesangs(pol1))
        return t1,t2,t3,t4
    if trans=='Ry': #reflect about y-axis
        t1,t2,t3,t4=2-t1,2-t2,t3,t4
        show('Reflect about y-axes: t1,t2=2-t1,2-t2')
        pol1  ,pic1  =poly4(t1,t2,t3,t4)
        #pic+=(polygon(pol1,fill=false,color='brown')+line([pol1[0],pol1[2]])+line([pol1[1],pol1[3]]))
        pic+=pic1
        #show(pic)
        show([t1,t2,t3,t4])
        show(gtsidesangs(pol1))
        return t1,t2,t3,t4

    if trans=='Rx': #reflect about x-axis
        t1,t2,t3,t4=t1,2-t2,t3,1-t4 
        show('Reflect about x-axes: t2,t4=2-t2,1-t4')
        pol2 ,pic2 =poly4(t1,t2,t3,t4)
        #pic+=(polygon(pol2,fill=false,color='green')+line([pol2[0],pol2[2]])+line([pol2[1],pol2[3]]))
        pic+=pic2
        #show(pic)
        show([t1,t2,t3,t4])
        show(gtsidesangs(pol2))
        return t1,t2,t3,t4
    if trans=='RtO': #rotate 180 about origin
        print('Rotate 180 about origin: t1=2-t1')
        t1,t2,t3,t4=2-t1,t2,t3,1-t4
        pol1,pic1=poly4(t1,t2,t3,t4)
        #pic+=(polygon(pol1,fill=false,color='turquoise')+line([pol1[0],pol1[2]])+line([pol1[1],pol1[3]]))
        pic+=pic1
        show(pic)
        show([t1,t2,t3,t4])
        show(gtsidesangs(pol1))
        pol=pol3
        return t1,t2,t3,t4



def showpols(pol,digs=3):
    if type(pol)==list:
        pol=gparms(pol)
    t1=pol[0];t2=pol[1];t3=pol[2];t4=pol[3]
    G=Graphics();parms=[];sidangs=[]
    pic1=poly4(t1,t2,t3,t4)[1]+text(mapnd([t1,t2,t3,t4],digs=digs),(0,1)) #original
    pic2=poly4(2-t1,2-t2,t3,t4)[1]+text(mapnd([2-t1,2-t2,t3,t4],digs=digs),(0,1)) #reflect about y-axis
    pic3=poly4(t1,2-t2,t3,1-t4 )[1]+text(mapnd([t1,2-t2,t3,1-t4],digs=digs),(0,1)) #reflect about x-axis
    pic4=poly4(2-t1,t2,t3,1-t4 )[1]+text(mapnd([2-t1,t2,t3,1-t4],digs=digs),(0,1)) #rotate 180 about (0,0)
    disp=graphics_array([[pic1,pic2],[pic3,pic4]])
    show(disp,figsize=6)
    if t3>1-meps(10^5):
        t1,t2,t3,t4=2*t4,t2,1/t3,t1/2
        pic1=poly4(t1,t2,t3,t4)[1]+text(mapnd([t1,t2,t3,t4],digs=digs),(0,1)) #original
        pic2=poly4(2-t1,2-t2,t3,t4)[1]+text(mapnd([2-t1,2-t2,t3,t4],digs=digs),(0,1)) #reflect about y-axis
        pic3=poly4(t1,2-t2,t3,1-t4 )[1]+text(mapnd([t1,2-t2,t3,1-t4 ],digs=digs),(0,1)) #reflect about x-axis
        pic4=poly4(2-t1,t2,t3,1-t4 )[1]+text(mapnd([2-t1,t2,t3,1-t4 ],digs=digs),(0,1)) #rotate 180 about (0,0)
        disp=graphics_array([[pic1,pic2],[pic3,pic4]])
        show(disp,figsize=6)

In [0]:

def mapnd(parms,digs=3):
    return [Rational(p).n(digits=digs) for p in parms]
mapnd([sqrt(2.)])

In [0]:

pol,pic=poly4(2/3,2/3,1,5/6)
pol1,pic1=poly4(1/3,2/3,1,2/3)
show(graphics_array([pic,pic1]),aspect_ratio=1,figsize=7)

In [0]:

def mapnd(parms,digs=3):
    if digs==0:
        return parms
    par=[]
    for p in parms:
        if type(p)==sage.rings.real_mpfr.RealNumber:
            par+=[Rational(p).n(digits=digs)]
        else:
            par+=[p.n(digits=digs)]
    return par

In [0]:

#showpols(mapv([(0,1),(1,2),(2,4),(-3,3)]))
showpols((1/4,1/2,1,1/8))

In [0]:

def goodparms(t1,t2,t3,t4):
    if test_parms(t1,t2,t3,t4):
        return t1,t2,t3,t4
    if test_parms(2-t1,2-t2,t3,t4):
        return 2-t1,2-t2,t3,t4
    if test_parms(t1,2-t2,t3,1-t4):
        return t1,2-t2,t3,1-t4
    if test_parms(2-t1,t2,t3,1-t4):
        return 2-t1,t2,t3,1-t4
    if t3>1-meps(10^5):
        t1,t2,t3,t4=2*t4,t2,1/t3,t1/2   
    if test_parms(t1,t2,t3,t4):
        return t1,t2,t3,t4
    if test_parms(2-t1,2-t2,t3,t4):
        return 2-t1,2-t2,t3,t4
    if test_parms(t1,2-t2,t3,1-t4):
        return t1,2-t2,t3,1-t4
    if test_parms(2-t1,t2,t3,1-t4):
        return 2-t1,t2,t3,1-t4

In [0]:

def test_parms(t1=.8,t2=1,t3=.91,t4=.5):
    if t4<=meps(10^5) or t4>=1-meps(10^5):
        return false
    if min(t1,t2,t3)<meps(10^5) or max(t1,t2,t3)>1+meps(10^5):
        return false
    if min(abs(t1-1),abs(t2-1))<meps(10^5) and t4>1/2-meps(10^5):
        return false
    if abs(t3-1)<meps(10^5) and (t4 > t1/2+meps(10^5) and t4<1-t1/2 ):
        return false
    return true
test_parms()

In [0]:

load('~/Forgons/')

In [0]:

@interact(layout=dict(top=[['t1','t2','t3','t4'],['classify','quads']]))
def _(t1=slider(default=1,vmin=0,vmax=1,step_size=1/100),t2=slider(default=1,vmin=0,vmax=1,step_size=1/100), t3=slider(default=1,vmin=0,vmax=1,step_size=1/100),t4=slider(default=1/2,vmin=0,vmax=1,step_size=1/100), quads=checkbox(default=false,label='Quadrisection(s)'),classify=checkbox(default=false,label='Classsify'),auto_update=false):
    global pol,D
    pol,pic=poly4(t1,t2,t3,t4)
    T1,T2,T3,T4=goodparms(t1,t2,t3,t4)
    pol1,pic1=poly4(T1,T2,T3,T4)
    Pic=Base+pic+pic1
    sidsangs(pol1)
    show('parms: '[T1,T2,T3,T4])
    show('sid')
    if quads  and t2>10^(-12):
        D=dpoly(pol)
        Pic+=getquads2(D)
    if classify:
        tis=classify_quad([t1,t2,t3,t4])
        print(tis)
    Pic.show(figsize=6,axes=false,aspect_ratio=1)

In [0]:

In [8]:

[(cos(i*2*pi.n()/5),sin(i*2*pi.n()/5)) for i in range(5)]

Out[8]:

[(1.00000000000000, 0.000000000000000),
 (0.309016994374947, 0.951056516295154),
 (-0.809016994374947, 0.587785252292473),
 (-0.809016994374948, -0.587785252292473),
 (0.309016994374947, -0.951056516295154)]

In [0]:

poly4(1,1/2,1/2,1/2)[1]

In [0]:

show(poly4(3/4,1,1 ,1/4)[1],figsize=10,aspect_ratio=1)

In [0]:

load('~/Forgons/quadapps.sage')

In [35]:

A=matrix([(1439384,12200),(0,.00001)])
Ai=A.inverse().n(digits=6)
b=vector([5,.02])
c=Ai*b
A*c

Out[35]:

(8.00000, 0.0200000)

In [0]:

In [0]:

In [0]:

In [0]:

In [32]:

In [37]:

np.linalg.cond(A)

Out[37]:

143948740534.56198

In [0]:

In [0]:

def mypoly(pol,clr='grey',lbls=false):
    pol1=pol
    pic=polygon(pol1,color=clr,alpha=.1)+line(pol1+[pol1[0]],color='black')+line([pol1[0],pol1[2]],color='black',linestyle='dotted',axes=false)
    pic+=line([pol1[1],pol1[3]],color='black',linestyle='dotted')+point(pol1,color='black')
    X=intsct([pol1[1],pol1[3]],[pol1[0],pol1[2]])
    M1=.5*(pol1[0]+pol1[2]);M2=.5*(pol1[1]+pol1[3])
    pic+=point([X,M1,M2],color='black')
    if lbls:
        eps=.06;dl=.06
        pic+=text("A",pol1[0]+offset(eps,0),color='black')
        pic+=text("C",pol1[2]+offset(-eps,0),color='black')
        pic+=text("B",pol1[1]+offset(0,dl),color='black')
        pic+=text("D",pol1[3]+offset(0,-dl),color='black')
        pic+=text("X",X+offset(-eps,dl),color='black')
    return pic

In [0]:

def gtsidesangs(pol):
    sids=[norm(pol[0]-pol[1]),norm(pol[2]-pol[1]),norm(pol[2]-pol[3]),norm(pol[0]-pol[3])]
    angs=[getang([pol[3],pol[0],pol[1]]),getang([ pol[0],pol[1],pol[2]]),getang([pol[1],pol[2],pol[3]]),getang([pol[2],pol[3],pol[0]])]
    return sids,angs

In [0]:

var('x1 y1 x2 y2 t s)
A=vector([1,0])
B=vector([x1,y1])
D=vecttor([x2,y2])
P1=A+t*(A-B);P2=A+s*(A-D)

In [0]:

load('~/Forgons/quadapps.sage')

In [0]:

In [0]:

integrate(cos(ln(x))/x,x,1,2)

In [0]:

integrate(cos(x),x,0,log(2))

In [7]:

27*75

Out[7]:

2025

In [0]:

def showpols3(parms,digs=3):
    t1=parms[0];t2=parms[1];t3=parms[2];t4=parms[3]
    G=Graphics();parms=[];sidangs=[]
    pic1=poly4(t1,t2,t3,t4)[1]+text([t1,t2,t3,t4],(0,1)) #original
    pic2=poly4(2-t1,2-t2,t3,t4)[1]+text([2-t1,2-t2,t3,t4],(0,1)) #reflect about y-axis
    pic3=poly4(t1,2-t2,t3,1-t4 )[1]+text([t1,2-t2,t3,1-t4],(0,1)) #reflect about x-axis
    pic4=poly4(2-t1,t2,t3,1-t4 )[1]+text([2-t1,t2,t3,1-t4],(0,1)) #rotate 180 about (0,0)
    disp=graphics_array([[pic1,pic2],[pic3,pic4]])
    show(disp,figsize=6)
    if t3>1-meps(10^5):
        t1,t2,t3,t4=2*t4,t2,1/t3,t1/2
        pic1=poly4(t1,t2,t3,t4)[1]+text([t1,t2,t3,t4],(0,1)) #original
        pic2=poly4(2-t1,2-t2,t3,t4)[1]+text([2-t1,2-t2,t3,t4],(0,1)) #reflect about y-axis
        pic3=poly4(t1,2-t2,t3,1-t4 )[1]+text([t1,2-t2,t3,1-t4],(0,1)) #reflect about x-axis
        pic4=poly4(2-t1,t2,t3,1-t4 )[1]+text([2-t1,t2,t3,1-t4],(0,1)) #rotate 180 about (0,0)
        disp=graphics_array([[pic1,pic2],[pic3,pic4]]) 
        show(disp,figsize=6)
        
def showpols2(parms):
    t1=parms[0];t2=parms[1];t3=parms[2];t4=parms[3]
    G=Graphics();parms=[];sidangs=[];pics=[]
    eps=.1
    t1t4sqr=line([(0,0),(1,0),(1,1),(0,1),(1,1/2),(0,0)],color='black')+line([(0,0),(0,1),(0,1/2),(1,1/2)],color='grey')
    t1t4sqr+=text('t1',(1/2,-eps),color='black')+text('t4',(1+eps,1/23),color='black')
    pics+=[poly4(t1,t2,t3,t4)[1]+text([t1,t2,t3,t4],(0,1))] #original
    t1t4sqr+=point((t1,t4))
    if 2-t1<=1 and 1-t4>0:
        pics+=[poly4(2-t1,t2,t3,1-t4 )[1]+text([2-t1,t2,t3,1-t4],(0,1))] #rotate 180 about (0,0)
        t1t4sqr+=point([(2-t1,1-t4)],color='red')
    if t3>1-meps(10^5):
        t1,t2,t3,t4=2*t4,t2,1/t3,t1/2
        if t1<=1:
            pics+=[poly4(t1,t2,t3,t4)[1]+text([t1,t2,t3,t4],(0,1))] #original
            t1t4sqr+=point([( t1, t4)])
        if 2-t1<=1 and 1-t4>0:
            pics+=[poly4(2-t1,t2,t3,1-t4 )[1]+text([2-t1,t2,t3,1-t4],(0,1))] #rotate 180 about (0,0)       
            t1t4sqr+=point([(2-t1,1-t4)])
    show(graphics_array(pics))
    show(t1t4sqr,aspect_ratio=1)
@interact(layout=dict(top=[['parameters']]))
def _(parms=input_box((1/2,1,1,1/2),width=20),auto_update=false):
    showpols2(parms)

In [0]:

showpols3(( 6/28, 9/10, 1, 25/28))

In [0]:

gtsidesangs(poly4( 3/14,9/10,1,3/28)[0] )

In [0]:

gtsidesangs(poly4( 3/14,9/10,1,25/28)[0] )

In [0]:

var('x y z')
eq = x^2 + y^2 + z^2 - 1
implicit_plot3d(eq, (x,-1,1), (y,-1,1), (z,-1,1), plot_points=10,mesh=true, boundary_offset=0.0,opacity=.4,color='yellow')

In [0]:

!hellno.sh ~/Forgons/stuff.ipynb
!hellno.sh ~/sagelets/stuff.ipynb

In [0]:

show(pic,ymax=2.2,ymin=-2.5,xmin=-2,xmax=4)

In [0]:

def threed_setup():
    var('t1 t2 t3 t4')
    C=Cube(names=['t4','t3','t1'],clrs=['grey','grey','grey'])
    pic=C.pic+point((1,1,1),viewpoint=[[-0.8,-0.3,-0.5],200])
    g1=polygon3d([(0,0,0),(1,0,1/2),(1,1,1/2),(0,1,0)],color='yellow',transparency=.1)
    p=pi.n()
    cyc=1-sqrt(1-4*t3^2*t4*(1-t4))
    #cyc=t3-1/2*sqrt((t1*(2-t1))/(1-t4))
    a=sqrt(4*t3^2*t4^2*sin(1.57079632679490*t2)^2 + (2*t3*t4*cos(1.57079632679490*t2) - t1)^2)
    b=sqrt(4*t3^2*t4^2*sin(1.57079632679490*t2)^2 + (2*t3*t4*cos(1.57079632679490*t2) - t1 + 2)^2)
    c=sqrt(4*t3^2*(t4 - 1)^2*sin(1.57079632679490*t2)^2 + (2*t3*(t4 - 1)*cos(1.57079632679490*t2) - t1 + 2)^2)
    d=sqrt(4*t3^2*(t4 - 1)^2*sin(1.57079632679490*t2)^2 + (2*t3*(t4 - 1)*cos(1.57079632679490*t2) - t1)^2)
    tce=a+c-b-d
    extceB=a+d-c-b
    extceA=a+b-c-d
    pic3=pic+text3d('(1,1,1)',(1,1,1.04),color='black')+text3d('(0,0,0)',(0,0,-.04),color='black')
    pic3+=text3d('(1,0,0)',(1.04, .04,-.04),color='black')+text3d('(0,0,1)',(0,0,1.04),color='black')+text3d('(0,1,0)',(0, 1.04,-.04),color='black')
    return pic3,cyc,tce,extceA,extceB

In [0]:

pic3,cyc,tce,extceA,extceB=threed_setup()
pic3,cyc,tce,extceA,extceB

In [0]:

prompt="""What do you look like?"""

In [0]:

print(prompt)

In [0]:


def levcurv(fun=x^2+y^2,x0=-1.,y0=0.,dt=.01, dist=.5, clr='black',grad=false,grid=100,t3d=false):
    """levcurv is used by patch, which used by full_patch:
    levcurv(fun=x^2+y^2,x0=-1,y0=0,dt=.01, dist=.5, clr='black',grad=false,grid=100) returns an approximate piece of the level curve of fun starting at (x0,y0) and moving in the direction of the gradient (rotated 90 degs counterclockwise) for an approximate distance dist by steps approximately dt units long. If dt is negative, it moves the opposite direction. If grad=true, it moves in the direction of the gradient (the direction of greatest ascent), and if also dt < 0, it moves in the direction of greatest descent.  levcurv returns two lists of points: the first is   the points along the curve; the second is a list of points on the curve approximately grid*dt units apart on the curve."""
    steps=floor(dist/abs(dt))
    if not type(fun)==list:
        f(x,y)=fun
        gradf=vector([derivative(f,x),derivative(f,y)])
    else:
        f1(x,y)=fun[0];f2(x,y)=fun[1]
        gradf=vector([f1,f2])
    if not grad:
        tanf=vector([-gradf[1],gradf[0]])
    else:
        tanf=gradf
    p=[vector([x0,y0])]
    gpts=[p[0]]

    for i in range(steps):
        dst=0
        while dst<abs(dt):
            lev=tanf(p[-1][0],p[-1][1]).n()
            d=min(1,1/norm(lev))
            p+=[p[-1]+d*lev/norm(lev)*dt]
            dst+=abs(d*dt)
        if mod(i,grid)==0:
            gpts+=[p[-1]]
    if t3d:
        p=[[q[0],q[1],f(q[0],q[1])] for q in p]
        gpts=[[g[0],g[1],f(g[0],g[1])] for g in gpts]
    return p,gpts

In [0]:

var('x y')
t1a=.5;t2a=.5;
f=(extceA-cyc+t1a).subs(t1=t1a,t2=t2a)
p,lev=levcurv(f.subs(t4=x,t3=y),0,0)
#pic3+line([1],color='black')
lev

In [0]:

var('x tth');sin(pi/2-tth).simplify_full()

In [0]:

t10=.4
strt=vector([0,0])
lev=[strt];msh=.01

f.subs(t4=lev[0][0],t3=lev[0][1])

In [0]:

r=.5
pol=mapv([(1,0),(r^2,2*r),(-r^2,2*r),(-1,0)])
cen,rad=circumcirc([pol[0],pol[1],pol[2]],pic=false);print(cen,rad)
pic=polygon(pol,fill=false)
parms=gparms(pol);print(parms)
pol1,pic1=poly4(parms[0],parms[1],parms[2],parms[3])
pic+=circle((0,r),r)+circle(cen,rad)
D=dpoly(pol)
C=cpoly5(D)
show(pic)
show(pic1)
showpols(pol)

In [0]:


pol1,pic1=poly4(1,.7,.5,1/2)

D=getpoly(pol1)
D.C.QAr

In [0]:

D.quad_pic+pic1
D.C.QAr

In [0]:

var('A B C s t x y')
A=vector([1,0]);B=vector([x,y]);C=vector([-1,0]);
s=(t*B*(A-B)-B*A)/((C-B)*(A+t*(B-A)))
X=A+t*(t*(B-A));Y=B+s*(C-B)
eqn=det(matrix([B,Y-X]))-1/2*det(matrix([B-A,C-A]))
eqn
solt=solve(eqn,t)
solt1=solt[0].rhs()
solt2=solt[1].rhs()
solt1

In [0]:

help(poly4)

In [0]:

t3=.6;
h0=t3;x0=.6;y0=h0*(1-x0);s0=1/(y0+x0/h0)
x1=-s0*y0;y1=s0*x0
pol,pic=poly4(1,1,t3,.5)
ep=.1
pic1=pic+text('(1,0)',(1+ep,0),color='black')+text('(0,h)',(0,h0+ep/2),color='black')+text('(-1,0)',(-1-ep,0),color='black',aspect_ratio=1)
pic2=pic1+point((x0,y0),color='black')+text('(x,y)',(x0+ep,y0+ep/2),color='black')+line([(0,0),(x0,y0)],color='black',linestyle='dashed')
pic3=pic2+point((x1,y1),color='black')+text('(-sy,sx)',(x1-ep,y1+ep/2),color='black')+line([(0,0),(x1,y1)],color='black',linestyle='dashed')
pic4=pic3+point((-x0,-y0),color='black')+text('(-x,-y)',(-x0-ep,-y0-ep/2),color='black')+line([(0,0),(-x0,-y0)],color='black',linestyle='dashed')
pic5=pic4+point((-x1,-y1),color='black')+text('( sy,-sx)',(-x1+ep,-y1-ep/2),color='black')+line([(0,0),(-x1,-y1)],color='black',linestyle='dashed')
pic5+polygon([(0,0),(x0,y0),(0,h0),(x1,y1)],color='lightgrey')

In [0]:

t1=solt1.subs(x=pol[1][0],y=pol[1][1]).real().n()
X1=X.subs(t=t1,x=pol[1][0],y=pol[1][1])
Y1=Y.subs(s=)

In [0]:

pic

In [0]:

pol=Roll(mapv([(-1, 0), (1, 0), (1/3, 61/40), (-1/3, 61/40)]),2)
H=getpoly(pol)

In [0]:

parms=gparms(mapv(pol))
parms=[p.n() for p in parms]

parms

In [0]:

H.quad_pic

In [0]:

#pol1,pic1=poly4(0.487394957983193, 0.909473351234499, 1.00000000000000, 0.243697478991597);
pol1,pic1=poly4(0.48739, 0.90947, 1.00000000000000, 0.24369);
G=getpoly(pol1)
G.quad_pic.show()

In [0]:

G.quad_pic+point(G.D.alverts)+line([pol1[0]])

In [0]:

#kite
var('x h')
y=(1-x)*h;s=1/(x/h+y)
A=(x+s*y)*h/2
pretty_print((h/2-A).numerator().collect(x))
pretty_print((h/2-A).denominator().collect(x))
#dA=derivative(A,x)
#pretty_print(dA.numerator().collect(x))
#pretty_print(dA.denominator().collect(x))

In [0]:

In [0]:

#parallelogram
var('x h a')
y=(x-1)/(a-1)*h
s=h/(h*y-x*(a+1))
A=det(matrix([[a,h],[-s*y-x,s*x-y]]))/2
pretty_print((A/h-1).numerator().collect(x))

solnx=solve(A-h/2,x) 
sol0=solnx[0].rhs();sol1=solnx[1].rhs()
def f(a,h):
    if a^2+h^2<1:
        return rtsave.subs(a=a,h=h)
    return 0
rtsave=((sol0+sol1)/2).simplify_full()
pretty_print(rtsave)
dA=derivative(A,x)/h
pretty_print(dA.numerator().collect(x))

In [0]:

plot3d(f,(a,0,1),(h,0,1))

In [0]:

pretty_print((solx[0]-a).numerator().collect(a))
pretty_print((solx[0]-a).denominator().collect(a))

In [0]:

eq=()

In [0]:

a0=.6;h0=.5
plot(A.subs(a==a0,h=h0)-h0,(x,.1,.9))

In [0]:

solnx[0].rhs().subs(a=.6,h=.5)

In [0]:

load('~/Forgons/quadapps.sage')

In [0]:

#pol1,pic1=poly4(0.48739, 0.90947, 1.00000000000000, 0.24369); #5 quadrisections
#pol1,pic1=poly4(0.48739, 0.90947, 1.00000000000000, 0.24369); #3 quadrisections
#pol,pic=poly4(.49,.89,1,.21) #3 quadrisections
#pol,pic=poly4(.47,.89,1,.21) #3 quads 1 in each of 3 orthogonal traps
#pol,pic=poly4(.434,.89,1,.217) #5 quadrisections
#pol,pic=poly4(.43,.89,1,.215) #5 quadrisections
#pol,pic=poly4(.42,.89,1,.21) #5 quadrisections
#pol,pic=poly4(.4,.88,1,.2) #5 quadrisections (2/5,22/25,1,1/5)
pol,pic=poly4(.4,.88,1,.2)
T=getpoly(pol)
T.C.QAr.show(figsize=4),(T.quad_pic+line([pol[1],pol[3]],linestyle='dotted',color='lightgrey')+T.C.trap[0].pic(0)).show(axes=true)

In [0]:

Tr=T.C.trap[0]
(Tr.pic(.5)+T.D.pic).show(axes=false)

plot(Tr.Ar1(t)+Tr.Ar2(t),(t,0,1),color='red')

In [0]:

da1=norm(Tr.T1.X(t)-Tr.Cn(t))*norm(Tr.T1.Y(t)-Tr.Cn(t));
da2=norm(Tr.T2.X(t)-Tr.Cn(t))*norm(Tr.T2.Y(t)-Tr.Cn(t));
plot(-da1(t=t)+da2(t=t),(0,1))

In [0]:

plot(Tr.DAr,(0,1))

In [0]:

((pic)+window(ep=.3,cen=Tr.Cn(.5))).show(axes=false)

In [0]:

T.C.QAr

In [0]:

var('t1 t2 t3 t4 T1 T4')
A=vector([1,0]);C=vector([-1,0]);X=vector([T1,0]);V=vector([t2,sqrt(1-t2^2)]);B=X+2*t3*t4*V;D=X-2*t3*T4*V
a =(B-A)*(B-A);a.simplify_full()
a=4*(T1 - 1)*t2*t3*t4 + 4*t3^2*t4^2 + t1^2

In [0]:

var('x y r s');l=-1;u=1;ar=.6;r0=-1/2;s0=-1/2;ll=-.1;uu=.1
pic1=plot3d(x^3+y^3 +3*x*y, (x,l,u),(y,l,u),aspect_ratio=[1,1,ar],mesh=true)+polygon3d([(-2,-2,0),(2,-2,0),(2,2,0),(-2,2,0)],alpha=.2)
tpl=r^3+s^3+3*r*s+3*(r^2+s)*(x-r)+3*(s^2+r)*(y-s)
tq=tpl+3*(x-r)^2+(y-s)^2+(x-r)*(y-s)
tpic=plot3d(tpl.subs(r=r0,s=s0),(x,r0+ll,r0+uu),(y,s0+ll,s0+uu),aspect_ratio=[1,1,ar],color='green')
tqpic=plot3d(tq.subs(r=r0,s=s0),(x,r0+ll,r0+uu),(y,s0+ll,s0+uu),aspect_ratio=[1,1,ar],color='red')
pic1+tpic+tqpic+point3d([(r0,s0,tpl.subs(x==0,y==0,r==r0,s==s0))],color='black',pointsize=100)

In [5]:

2^12

Out[5]:

4096

In [4]:

x=((1+sqrt(5))/2)^12;x.expand()

Out[4]:

72*sqrt(5) + 161

In [0]:

#ExtceAGraph=Graphics()
#for i in range(6):
#    ExtceAGraph+=implicit_plot3d(extceA.subs(t2=i/5.+.002),(t4,.01,1),(t3,.001,1),(t1,.001,1),color='green',frame=false,alpha=2,mesh=true)
#ExtceAGraph

In [17]:

ln(300.)/ln(3.)

Out[17]:

5.19180654857877

In [21]:

1.+3.*log(5.)/log(3.)
1.+3*log(5.,3.)

Out[21]:

5.39492056215378

In [2]:

var('y,h,k,a,b,x,x1,y1,z1,w1')
eqn=(x-h)^2/a^2+(y-k)^2/b^2-1
eqn

Out[2]:

(h - x)^2/a^2 + (k - y)^2/b^2 - 1

In [548]:

var('h,k,a2,b2,x,x1,y1,z1,w1')
#h=0
#assume(y1>0,w1>0,a>0,b>=a)
#assume(a2>0,b2>0)
eqn=x^2/a2+(y-k)^2/b2-1

eq1=eqn.subs(x=1,y=0)
eq2=eqn.subs(x=-1,y=0)
eq3=eqn.subs(x=x1,y=y1)
eq4=eqn.subs(x=z1,y=w1)
sol=solve([eq2,eq3,eq4],[k,a2,b2])
sol

Out[548]:

[[k == r97, a2 == 0, b2 == 0], [k == 1/2*(w1^2*x1^2 - y1^2*z1^2 - w1^2 + y1^2)/(w1*x1^2 - y1*z1^2 - w1 + y1), a2 == -1/4*(w1^4*x1^4 + y1^4*z1^4 - 2*w1^4*x1^2 + w1^4 - 4*w1*y1^3 + y1^4 - 2*(w1^2*x1^2 - 3*w1^2)*y1^2 + 2*(2*w1*y1^3 - y1^4 - (w1^2*x1^2 + w1^2)*y1^2)*z1^2 + 4*(w1^3*x1^2 - w1^3)*y1)/(w1*y1^3 + (w1^2*x1^2 - 2*w1^2)*y1^2 + (w1^2*y1^2 - w1*y1^3)*z1^2 - (w1^3*x1^2 - w1^3)*y1), b2 == 1/4*(w1^4*x1^4 + y1^4*z1^4 - 2*w1^4*x1^2 + w1^4 - 4*w1*y1^3 + y1^4 - 2*(w1^2*x1^2 - 3*w1^2)*y1^2 + 2*(2*w1*y1^3 - y1^4 - (w1^2*x1^2 + w1^2)*y1^2)*z1^2 + 4*(w1^3*x1^2 - w1^3)*y1)/(w1^2*x1^4 + y1^2*z1^4 - 2*w1^2*x1^2 - 2*((w1*x1^2 - w1)*y1 + y1^2)*z1^2 + w1^2 + 2*(w1*x1^2 - w1)*y1 + y1^2)]]

In [561]:

Eqn=eqn.subs(k=sol[1][0].rhs(),a2=sol[1][1].rhs(),b2=sol[1][2].rhs())
Eqn=Eqn.numerator().collect(x).collect(y)/4


Eqn

Out[561]:

-w1^3*x1^2*y1 + w1^2*x1^2*y1^2 + w1^2*y1^2*z1^2 - w1*y1^3*z1^2 + w1^3*y1 - 2*w1^2*y1^2 + w1*y1^3 + (w1^3*x1^2*y1 - w1^2*x1^2*y1^2 - w1^2*y1^2*z1^2 + w1*y1^3*z1^2 - w1^3*y1 + 2*w1^2*y1^2 - w1*y1^3)*x^2 + (w1^2*x1^4 - 2*w1*x1^2*y1*z1^2 + y1^2*z1^4 - 2*w1^2*x1^2 + 2*w1*x1^2*y1 + 2*w1*y1*z1^2 - 2*y1^2*z1^2 + w1^2 - 2*w1*y1 + y1^2)*y^2 - (w1^3*x1^4 - w1^2*x1^2*y1*z1^2 - w1*x1^2*y1^2*z1^2 + y1^3*z1^4 - 2*w1^3*x1^2 + w1^2*x1^2*y1 + w1*x1^2*y1^2 + w1^2*y1*z1^2 + w1*y1^2*z1^2 - 2*y1^3*z1^2 + w1^3 - w1^2*y1 - w1*y1^2 + y1^3)*y

In [565]:

Eqn.show()
Eqn.coefficient(y).factor()

Out[565]:

$\displaystyle -w_{1}^{3} x_{1}^{2} y_{1} + w_{1}^{2} x_{1}^{2} y_{1}^{2} + w_{1}^{2} y_{1}^{2} z_{1}^{2} - w_{1} y_{1}^{3} z_{1}^{2} + w_{1}^{3} y_{1} - 2 \, w_{1}^{2} y_{1}^{2} + w_{1} y_{1}^{3} + {\left(w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}\right)} x^{2} + {\left(w_{1}^{2} x_{1}^{4} - 2 \, w_{1} x_{1}^{2} y_{1} z_{1}^{2} + y_{1}^{2} z_{1}^{4} - 2 \, w_{1}^{2} x_{1}^{2} + 2 \, w_{1} x_{1}^{2} y_{1} + 2 \, w_{1} y_{1} z_{1}^{2} - 2 \, y_{1}^{2} z_{1}^{2} + w_{1}^{2} - 2 \, w_{1} y_{1} + y_{1}^{2}\right)} y^{2} - {\left(w_{1}^{3} x_{1}^{4} - w_{1}^{2} x_{1}^{2} y_{1} z_{1}^{2} - w_{1} x_{1}^{2} y_{1}^{2} z_{1}^{2} + y_{1}^{3} z_{1}^{4} - 2 \, w_{1}^{3} x_{1}^{2} + w_{1}^{2} x_{1}^{2} y_{1} + w_{1} x_{1}^{2} y_{1}^{2} + w_{1}^{2} y_{1} z_{1}^{2} + w_{1} y_{1}^{2} z_{1}^{2} - 2 \, y_{1}^{3} z_{1}^{2} + w_{1}^{3} - w_{1}^{2} y_{1} - w_{1} y_{1}^{2} + y_{1}^{3}\right)} y$

-(w1^2*x1^2 - y1^2*z1^2 - w1^2 + y1^2)*(w1*x1^2 - y1*z1^2 - w1 + y1)

In [11]:

plot(.5*t-ln(7-2),(1,4))

Out[11]:

In [9]:

var('t');find_root(.5*t-ln(7-2),1,4)

Out[9]:

3.2188758248682006

In [49]:

f=lambda t:t^t-49
f(3.2780315235960593)

Out[49]:

-8.437694987151190e-15

In [4]:

find_root(t^t-49,1,4)

Out[4]:

3.2780315235960593

In [7]:

plot(t^t-49,(1,3.5))

Out[7]:

In [0]:

In [476]:

k=sol[1][0].rhs()
a2=sol[1][1].rhs()
b2=sol[1][2].rhs()
#as2=a2.denominator().factor()/a2.numerator().factor()
#bs2=b2.denominator().factor()/b2.numerator().factor()
#ks=k.denominator().factor()/k.numerator().factor()
ks=2*(w1*x1^2 - y1*z1^2 - w1 + y1)/(w1^2*x1^2 - y1^2*z1^2 - w1^2 + y1^2)
as2=4*(w1*x1^2 - y1*z1^2 - w1 + y1)*(w1 - y1)*w1*y1/((w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))
bs2=4*(w1*x1^2 - y1*z1^2 - w1 + y1)^2/((w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))
show(ks,as2,bs2)

Out[476]:

$\displaystyle \frac{2 \, {\left(w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}\right)}}{w_{1}^{2} x_{1}^{2} - y_{1}^{2} z_{1}^{2} - w_{1}^{2} + y_{1}^{2}} \frac{4 \, {\left(w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}\right)} {\left(w_{1} - y_{1}\right)} w_{1} y_{1}}{{\left(w_{1} x_{1} + y_{1} z_{1} + w_{1} - y_{1}\right)} {\left(w_{1} x_{1} + y_{1} z_{1} - w_{1} + y_{1}\right)} {\left(w_{1} x_{1} - y_{1} z_{1} + w_{1} - y_{1}\right)} {\left(w_{1} x_{1} - y_{1} z_{1} - w_{1} + y_{1}\right)}} \frac{4 \, {\left(w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}\right)}^{2}}{{\left(w_{1} x_{1} + y_{1} z_{1} + w_{1} - y_{1}\right)} {\left(w_{1} x_{1} + y_{1} z_{1} - w_{1} + y_{1}\right)} {\left(w_{1} x_{1} - y_{1} z_{1} + w_{1} - y_{1}\right)} {\left(w_{1} x_{1} - y_{1} z_{1} - w_{1} + y_{1}\right)}}$

In [0]:

 $$k=\frac{2(w1 x1^2 - y1z1^2 - w1 + y1}{w1^2 x1^2 - y1^2*z1^2 - w1^2 + y1^2}$$
 $$a=\frac{ 4*(w1*x1^2 - y1*z1^2 - w1 + y1)*(w1 - y1)*w1*y1}{(w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1)}$$
$$b={ 4*(w1*x1^2 - y1*z1^2 - w1 + y1)^2/((w1*x1 + y1*z1 + w1 - y1}{w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))}$$

In [337]:

#b2.denominator().factor()/(a2.numerator().factor())
as2=4*(w1*x1^2 - y1*z1^2 - w1 + y1)*(w1 - y1)*w1*y1/((w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))

In [335]:

(4*(w1*x1^2 - y1*z1^2 - w1 + y1)*(w1 - y1)*w1*y1)/((w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))

Out[335]:

4*(w1*x1^2 - y1*z1^2 - w1 + y1)*(w1 - y1)*w1*y1/((w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))

In [378]:

help(Rot)

Out[378]:

Help on function Rot in module __main__:

Rot(Pol=[[0, 0], [1, 0], [1, 3]], cen=[0, 0], th=30, deg=True, dec=True)
    Rot(Pol=[[0,0],[1,0],[1,3]],cen=[0,0],th=30,deg=true) returns the rotation of Pol by th degrees about cen.

In [339]:

as2=(a2.denominator().factor())/(a2.numerator().factor())
bs=(b2.denominator().factor())/(b2.numerator().factor())
print(as2)
bs

Out[339]:

4*(w1*x1^2 - y1*z1^2 - w1 + y1)*(w1 - y1)*w1*y1/((w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))

4*(w1*x1^2 - y1*z1^2 - w1 + y1)^2/((w1*x1 + y1*z1 + w1 - y1)*(w1*x1 + y1*z1 - w1 + y1)*(w1*x1 - y1*z1 + w1 - y1)*(w1*x1 - y1*z1 - w1 + y1))

In [577]:

var('x1 y1 x y z1 w1')
EQ1=w1^2*x1^2*y - w1*x1^2*y^2 - w1^2*x^2*y1 + w1*x^2*y1^2 + y^2*y1*z1^2 - y*y1^2*z1^2 - w1^2*y + w1*y^2 + w1^2*y1 - y^2*y1 - w1*y1^2 + y*y1^2
EQ1=(EQ1.collect(y).collect(x)/(-w1^2*y1+w1*y1^2)).numerator().collect(y).collect(x)
EQ1.show()
A1=EQ1.coefficient(y^2);A1.show()
B1=-EQ1.coefficient(y);B1.show()
C1=EQ1.subs(x=0,y=0);C1.show()
D1=B1^2-4*A1*C1;D1.show()
Y1=(-B1+sqrt(D1))/(2*A1);show(Y1)

Out[577]:

$\displaystyle {\left(w_{1}^{2} y_{1} - w_{1} y_{1}^{2}\right)} x^{2} + {\left(w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}\right)} y^{2} - w_{1}^{2} y_{1} + w_{1} y_{1}^{2} - {\left(w_{1}^{2} x_{1}^{2} - y_{1}^{2} z_{1}^{2} - w_{1}^{2} + y_{1}^{2}\right)} y$

$\displaystyle w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}$

$\displaystyle w_{1}^{2} x_{1}^{2} - y_{1}^{2} z_{1}^{2} - w_{1}^{2} + y_{1}^{2}$

$\displaystyle -w_{1}^{2} y_{1} + w_{1} y_{1}^{2}$

$\displaystyle {\left(w_{1}^{2} x_{1}^{2} - y_{1}^{2} z_{1}^{2} - w_{1}^{2} + y_{1}^{2}\right)}^{2} + 4 \, {\left(w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}\right)} {\left(w_{1}^{2} y_{1} - w_{1} y_{1}^{2}\right)}$

$\displaystyle -\frac{w_{1}^{2} x_{1}^{2} - y_{1}^{2} z_{1}^{2} - w_{1}^{2} + y_{1}^{2} - \sqrt{{\left(w_{1}^{2} x_{1}^{2} - y_{1}^{2} z_{1}^{2} - w_{1}^{2} + y_{1}^{2}\right)}^{2} + 4 \, {\left(w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}\right)} {\left(w_{1}^{2} y_{1} - w_{1} y_{1}^{2}\right)}}}{2 \, {\left(w_{1} x_{1}^{2} - y_{1} z_{1}^{2} - w_{1} + y_{1}\right)}}$

In [0]:

In [40]:

var('a, b, c, d, e, k,f, g, x, y,x1,y1,z1,w1')
eqn=a*x^2+b*y^2+c*x+d*y+e*x*y+f
eqn=a*x^2+b*(y-k)^2-1
#eqn= x^2+b*(y-d)^2-1
#pol=poly4(.5,.9,1,.1)[0]
pol=[(1,0),(x1,y1),(-1,0),(z1,w1)]
eq1=eqn.subs(x=pol[0][0],y=pol[0][1])
eq2=eqn.subs(x=pol[1][0],y=pol[1][1])
eq3=eqn.subs(x=pol[2][0],y=pol[2][1])
eq4=eqn.subs(x=pol[3][0],y=pol[3][1])
sol=solve([eq1,eq2,eq3,eq4],[k,a,b]);print(len(sol));print(sol)
k=sol[0][0].rhs();a=sol[0][1].rhs();b=sol[0][2].rhs()
EQ1=eqn.subs(k=sol[0][0].rhs(),a=sol[0][1].rhs(),b=sol[0][2].rhs())


EQ1=w1^2*x1^2*y - w1*x1^2*y^2 - w1^2*x^2*y1 + w1*x^2*y1^2 + y^2*y1*z1^2 - y*y1^2*z1^2 - w1^2*y + w1*y^2 + w1^2*y1 - y^2*y1 - w1*y1^2 + y*y1^2
pol,pic=poly4(.8,.5,1,.3)
x1=pol[1][0];y1=pol[1][1];z1=pol[3][0];w1=pol[3][1]
EQ2=EQ1.subs(x1=x1,y1=y1,z1=z1,w1=w1)
sol=solve(EQ2.subs(x=0),y)
A1=vector([0,sol[0].rhs()])
A2=vector([0,sol[1].rhs()])
Cen=.5*(A1+A2)
sol=solve(EQ2.subs(y=Cen[1]),x)
B1=vector([sol[0].rhs(),Cen[1]])
B2=vector([sol[1].rhs(),Cen[1]])
wd=1.25;Pic2=pic+implicit_plot(EQ2,(-wd,wd),(-1.5,.8))+line([A1,A2],linestyle='dotted')+line([B1,B2],linestyle='dotted')+point([A1,B1,Cen,A2,B2])
 
show(Pic2,aspect_ratio=1,axes=false)

Out[40]:

1
[
[k == 1/2*(w1^2*x1^2 - y1^2*z1^2 - w1^2 + y1^2)/(w1*x1^2 - y1*z1^2 - w1 + y1), a == -4*(w1*y1^3 + (w1^2*x1^2 - 2*w1^2)*y1^2 + (w1^2*y1^2 - w1*y1^3)*z1^2 - (w1^3*x1^2 - w1^3)*y1)/(w1^4*x1^4 + y1^4*z1^4 - 2*w1^4*x1^2 + w1^4 - 4*w1*y1^3 + y1^4 - 2*(w1^2*x1^2 - 3*w1^2)*y1^2 + 2*(2*w1*y1^3 - y1^4 - (w1^2*x1^2 + w1^2)*y1^2)*z1^2 + 4*(w1^3*x1^2 - w1^3)*y1), b == 4*(w1^2*x1^4 + y1^2*z1^4 - 2*w1^2*x1^2 - 2*((w1*x1^2 - w1)*y1 + y1^2)*z1^2 + w1^2 + 2*(w1*x1^2 - w1)*y1 + y1^2)/(w1^4*x1^4 + y1^4*z1^4 - 2*w1^4*x1^2 + w1^4 - 4*w1*y1^3 + y1^4 - 2*(w1^2*x1^2 - 3*w1^2)*y1^2 + 2*(2*w1*y1^3 - y1^4 - (w1^2*x1^2 + w1^2)*y1^2)*z1^2 + 4*(w1^3*x1^2 - w1^3)*y1)]
]

In [51]:

def offset(Off=(.05,.05)):
    return vector(Off)

In [55]:

Pic3=Pic2+text("A",(1.05,0))+text("B",(pol[1]+offset((.05,.05))))+text("C",(-1.08,0))+text("D",pol[3]+offset((-.05,-.05)));#Pic3.show(axes=false)

save(Pic3, 'Pic3.png')

In [0]:

In [25]:

EQ1.collect(y ).collect(x)
soly=solve(EQ1.subs(x=0),y)
A1=vector([0,soly[0].rhs()])
A2=vector([0,soly[1].rhs()])
Cen=(A1+A2)/2
solx=solve(EQ1.subs(y=Cen[1]),x)
B1=vector([0,solx[0].rhs()])
B2=vector([0,solx[1].rhs()])
show(B2.simplify())

Out[25]:

$\displaystyle \left(0,\,\frac{1}{2} \, \sqrt{\frac{w_{1}^{4} x_{1}^{4}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} - \frac{2 \, w_{1}^{2} x_{1}^{2} y_{1}^{2} z_{1}^{2}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} + \frac{y_{1}^{4} z_{1}^{4}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} - \frac{2 \, w_{1}^{4} x_{1}^{2}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} + \frac{4 \, w_{1}^{3} x_{1}^{2} y_{1}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} - \frac{2 \, w_{1}^{2} x_{1}^{2} y_{1}^{2}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} - \frac{2 \, w_{1}^{2} y_{1}^{2} z_{1}^{2}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} + \frac{4 \, w_{1} y_{1}^{3} z_{1}^{2}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} - \frac{2 \, y_{1}^{4} z_{1}^{2}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} + \frac{w_{1}^{4}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} - \frac{4 \, w_{1}^{3} y_{1}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} + \frac{6 \, w_{1}^{2} y_{1}^{2}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} - \frac{4 \, w_{1} y_{1}^{3}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}} + \frac{y_{1}^{4}}{w_{1}^{3} x_{1}^{2} y_{1} - w_{1}^{2} x_{1}^{2} y_{1}^{2} - w_{1}^{2} y_{1}^{2} z_{1}^{2} + w_{1} y_{1}^{3} z_{1}^{2} - w_{1}^{3} y_{1} + 2 \, w_{1}^{2} y_{1}^{2} - w_{1} y_{1}^{3}}}\right)$

In [32]:

print(EQ1.collect(y).collect(x))

Out[32]:

-(w1^2*y1 - w1*y1^2)*x^2 - (w1*x1^2 - y1*z1^2 - w1 + y1)*y^2 + w1^2*y1 - w1*y1^2 + (w1^2*x1^2 - y1^2*z1^2 - w1^2 + y1^2)*y

In [507]:

var('a, b, c, d, e, k,f, g, x, y,x1,y1,z1,w1')
eqn=a*x^2+b*y^2+c*x+d*y+e*x*y+f
eqn=x^2+b*y^2+1*y+g*x*y-1

#eqn= x^2+b*(y-d)^2-1
#pol=poly4(.5,.9,1,.1)[0]
pol=[(1,0),(x1,y1),(-1,0),(z1,w1)]
eq1=eqn.subs(x=pol[0][0],y=pol[0][1])
eq2=eqn.subs(x=pol[1][0],y=pol[1][1])
eq3=eqn.subs(x=pol[2][0],y=pol[2][1])
eq4=eqn.subs(x=pol[3][0],y=pol[3][1])
sol=solve([eq1,eq2,eq3,eq4],[b,g ]);print(len(sol));print(sol)
EQ1=eqn.subs(b=sol[0][0].rhs(),g=sol[0][1].rhs())
show(EQ1)
b1 = (w1*x1^2 + (g*w1*x1 - g*w1*z1 - z1^2 + 1)*y1 - w1)/(w1^2*y1 - w1*y1^2); d1 = -(g*w1^2*x1*y1 + w1^2*x1^2 - (g*w1*z1 + z1^2 - 1)*y1^2 - w1^2)/(w1^2*y1 - w1*y1^2)
pol,pic=poly4(.6,.7,.8,.9)
wd=2.5;Pic2=pic+implicit_plot(EQ1.subs(x1=pol[1][0],y1=pol[1][1],z1=pol[3][0],w1=pol[3][1]),(-wd,wd),(-wd,wd))

Out[507]:

1
[
[b == (w1*x1^2*z1 - ((z1^2 + w1 - 1)*x1 - w1*z1)*y1 - w1*z1)/(w1^2*x1*y1 - w1*y1^2*z1), g == -(w1^2*x1^2 + w1^2*y1 - (z1^2 + w1 - 1)*y1^2 - w1^2)/(w1^2*x1*y1 - w1*y1^2*z1)]
]

$\displaystyle x^{2} - \frac{{\left(w_{1}^{2} x_{1}^{2} + w_{1}^{2} y_{1} - {\left(z_{1}^{2} + w_{1} - 1\right)} y_{1}^{2} - w_{1}^{2}\right)} x y}{w_{1}^{2} x_{1} y_{1} - w_{1} y_{1}^{2} z_{1}} + \frac{{\left(w_{1} x_{1}^{2} z_{1} - {\left({\left(z_{1}^{2} + w_{1} - 1\right)} x_{1} - w_{1} z_{1}\right)} y_{1} - w_{1} z_{1}\right)} y^{2}}{w_{1}^{2} x_{1} y_{1} - w_{1} y_{1}^{2} z_{1}} + y - 1$

In [0]:

b1

In [0]:

In [512]:

Pic2.show(axes=true)

Out[512]:

In [0]:

In [384]:

var('h,k,a2,b2,x,x1,y1,z1,w1')
Pol=poly4(.5,.1,1,.1)[0]
#Pol=Rot(pol,(0,0),20);print(Pol)
#x1,y1,z1,w1=pol[1][0],pol[1][1],pol[3][0],pol[3][1]
#pol=[(1,0),(x1,y1),(-1,0),(z1,w1)]
eqn=x^2/a2+(y-k)^2/b2-1
eq1=eqn.subs(x=Pol[0][0],y=Pol[0][1])
eq2=eqn.subs(x=Pol[1][0],y=Pol[1][1])
eq3=eqn.subs(x=Pol[2][0],y=Pol[2][1])
eq4=eqn.subs(x=Pol[3][0],y=Pol[3][1])
sol=solve([eq1,eq3,eq4],[k,a2,b2]);print(len(sol));show(sol)
j=1
ksol =sol[j][0].rhs().n() ; a2sol = sol[j][1].rhs().n(); b2sol = sol[j][2].rhs().n()
c=sqrt(abs(a2sol-b2sol))
if a2sol -b2sol >= 0:
    F1=vector([c,ksol]);F2=vector([-c,ksol])
else:
    F1=vector([0,c+ksol]);F2=vector([0,ksol-c])
cen=vector([0,ksol])
Eqn2=eqn.subs(k =sol[j][0].rhs().n(), a2 = sol[j][1].rhs().n(), b2 = sol[j][2].rhs().n())
show(Eqn2)
pic=polygon(pol,fill=false)+point(pol)+implicit_plot(Eqn2,(-2.1,2.1),(-1.1,2.6),figsize=5,color='red',aspect_ratio=1)+point([cen,F1,F2])
pic

Out[384]:

[(0.939692620785908, 0.342020143325669), (0.644770251820833, 0.267971995712102), (-0.939692620785908, -0.342020143325669), (-1.10446916245796, -0.701647244780572)]
2

$\displaystyle \left[\left[k = r_{66}, a_{2} = 0, b_{2} = 0\right], \left[k = 0, a_{2} = \left(\frac{895675310220197311064290413618058279128137323}{1151188606729446985558295821676365277620123890}\right), b_{2} = \left(-\frac{8061077791981775799578613722562524512153235907}{9297905218000288917006669051121867039586720000}\right)\right]\right]$

$\displaystyle 1.28527446675563 \, x^{2} - 1.15343201714897 \, y^{2} - 1$

In [434]:

var('h,k,a2,b2,x,x1,y1,z1,w1')
pol,pic=poly4(.5,.9,1,.6)
x1,y1,z1,w1=pol[1][0],pol[1][1],pol[3][0],pol[3][1]
#x1,y1,z1,w1=3/4,sqrt(1-4/16),-3/4,-sqrt(1-4/16)
#pol=[(1,0),(x1,y1),(-1,0),(z1,w1)]
eqn=x^2/a2+(y-k)^2/b2-1
eq1=eqn.subs(x=1,y=0)
eq2=eqn.subs(x=-1,y=0)
eq3=eqn.subs(x=x1,y=y1)
eq4=eqn.subs(x=z1,y=w1)
sol=solve([eq2,eq3,eq4],[k,a2,b2]);print(len(sol));show(sol)
j=1
ksol =sol[j][0].rhs().n() ; a2sol = sol[j][1].rhs().n(); b2sol = sol[j][2].rhs().n()
c=sqrt(abs(a2sol-b2sol))
if a2sol -b2sol >= 0:
    F1=vector([c,ksol]);F2=vector([-c,ksol])
else:
    F1=vector([0,c+ksol]);F2=vector([0,ksol-c])
cen=vector([0,ksol])
Eqn2=eqn.subs(k =sol[j][0].rhs().n(), a2 = sol[j][1].rhs().n(), b2 = sol[j][2].rhs().n())
show(Eqn2)
pic=polygon(pol,fill=false)+point(pol)+implicit_plot(Eqn2,(-2.1,2.1),(-1.1,2.6),figsize=5,color='red',aspect_ratio=1)+point([cen,F1,F2])
pic+Pic+Pic2

Out[434]:

2

$\displaystyle \left[\left[k = r_{87}, a_{2} = 0, b_{2} = 0\right], \left[k = \left(\frac{111142505458509796066655}{363200756838291748171074}\right), a_{2} = \left(\frac{1093837001075934519950566263497961}{1019758370227233636876860548687520}\right), b_{2} = \left(\frac{20266442424328830464818813449837148181196915225}{14657198863100881119036049257289261002663368164}\right)\right]\right]$

$\displaystyle 0.932276353080181 \, x^{2} + 0.723225051354137 \, {\left(y - 0.306008463269789\right)}^{2} - 1$

In [68]:

simple_nonconvex=poly4(-.5,1,1,.5,labels=true,model=false)[1]
save(simple_nonconvex,'images/simple_nonconvex.png')

In [67]:

poly4(-.5,1,1,.5,labels=true,model=false)[1]

Out[67]:

The code below was generated by GPT-3.5 using this prompt:

save a diagram in the png format to my directory

In [0]:

Source of the problem below

Product

Resources

Company