Constructible

def is_constructible(a): '''This function takes in a number a and tells the user whether or not a is constructible. Not complete. Author: bleveque''' if a in QQ: print a, 'is constructible' else: is_a_algebraic = 1 try: minimal_polynomial(a) except TypeError: print 'a is not constructible -- a is transcendental' is_a_algebraic = 0 if is_a_algebraic == 1: K.<t>=NumberField(minimal_polynomial(a)) d=K.degree() if d==2: print a, 'is constructible' elif log(d,2) in ZZ: print a, 'is possibly constructible. Check that the Number Field it lies in a series of quadratic extensions of QQ.' else: print a, 'is not constructible' 

is_constructible(cos(pi/9)) 

cos(1/9*pi) is not constructible

cos(1/9*pi) is not constructible

is_constructible(sqrt(3)) 

sqrt(3) is constructible

sqrt(3) is constructible

is_constructible(3^(1/4)) 

3^(1/4) is possibly constructible. Check that the Number Field it lies
in a series of quadratic extensions of QQ.

3^(1/4) is possibly constructible. Check that the Number Field it lies in a series of quadratic extensions of QQ.

is_constructible(3) 

3 is constructible

3 is constructible

is_constructible(pi) 

a is not constructible -- a is transcendental

a is not constructible -- a is transcendental

############################################################### 

############################################################### 

def FareyGraph(xrange,yrange,fraction_limit=3): p=line([]) for i in range(xrange[0],xrange[1]): p+=line([(i,0),(i,yrange[1])]) p+=farey_fill_in(i,i+1,fraction_limit) p+=line([(xrange[1],0),(xrange[1],yrange[1])]) return p ################################################### def farey_fill_in(x0,x1,flimit): if flimit <= 0: return line([]) elif flimit < 0: 'the input flimit must be non-negative' return line([]) else: q=line([]) for j in range(1,flimit+1): if j==1: q+=plot_semi_circle(x0,x1) else: j=float(j) q+=plot_semi_circle(x0,x0+(x1-x0)/j) q+=plot_semi_circle(x1-(x1-x0)/j,x1) if j>2: for k in range(1,j+1): k=float(k) if gcd(j,k)==1: nextx=x1 for h in range(int(j)-1,0,-1): for g in range(1,h+1): g,h=float(g),float(h) if (x0+k/j)<(x0+g/h)<nextx: nextx=x0+g/h q+=plot_semi_circle(x0+(x1-x0)*k/j,nextx) q+=plot_semi_circle(x1-(nextx-x0),x1-(x1-x0)*k/j) return q ############################################################ def plot_semi_circle(x0,x1): if x0==x1: return line([]) else: x0,x1=min(x0,x1),max(x0,x1) t=var('t') r=(x1-x0)/2 xmid=x0+r p=parametric_plot((r*cos(t)+xmid,r*sin(t)),(t,0,pi)) return p 

def slaction(P,a11,a12,a21,a22,plotcolor='green'): z=P[0]+P[1]*i Q=point((real((a11*z+a12)/(a21*z+a22)),imag((a11*z+a12)/(a21*z+a22))),color=plotcolor,pointsize=5) return Q ########################################### @interact def fgraph_interact(slide=slider(0,10,1,0,label='Detail of Farey Graph: ',display_value=True), x_min_box=input_box(default=-2,type=int,label='x-min: ',width=5), x_max_box=input_box(default=2,type=int,label='x-max: ',width=5), A=input_grid(2,2,default=[[1,0],[0,1]],label='A: '), ExtraPoints=input_box(default=[],type=list,label='Any other points? ',width=50), pointcheck=checkbox(default=False,label='Show points? '), check=checkbox(default=False,label='Show action on points by matrix A? ')): print 'The matrix A is: ' print matrix(ZZ,[[A[0][0], A[0][1]],[A[1][0], A[1][1]]]) if A[0][0]*A[1][1]-A[0][1]*A[1][0] == 1: print 'A is in SL(2,Z)' else: print 'A is not in SL(2,Z)' p=FareyGraph((x_min_box,x_max_box),(-1,1),slide) if slide>0 and pointcheck==True: a=(0,.2) b=(0,.4) c=(0,.6) d=(0,.8) e=(0,1) f=(1,.2) g=(1,.4) h=(1,.6) i=(1,.8) j=(1,1) k=(.5,.5) l=(.5,1) m=(.25,.75) n=(.75,.75) o=(.5,.75) examplepoints=[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o] for expoi in examplepoints: p+=point(expoi,color='red') for extrapoint in ExtraPoints: p+=point(extrapoint,color='red') if A[0][0]*A[1][1]-A[0][1]*A[1][0] == 1: if check == True: for slpoints in examplepoints: p+=slaction(slpoints,A[0][0],A[0][1],A[1][0],A[1][1]) for extraslpoints in ExtraPoints: p+= slaction(extraslpoints,A[0][0],A[0][1],A[1][0],A[1][1],'orange') show(p,figsize=(10,6)) 

############################################################### 

############################################################### 

def pointfind(b,d,e,xval,sign): if xval^3+b*xval^2+d*xval+e >= 0: z = n(sqrt(xval^3+b*xval^2+d*xval+e)) if sign>=0: return (xval,z) else: return (xval,-z) ######################################## def compute_intersection(b,d,e,P1,P2): if P1==P2: m=(3*P1[0]^2+2*b*P1[0]+d)/float(2*P1[1]) v=P1[1]-m*P1[0] x=(P1[0]^4-2*d*P1[0]^2-8*e*P1[0]+d^2-4*b*e)/float(4*P1[0]^3+4*b*P1[0]^2+4*d*P1[0]+4*e) return(x,m*x+v) else: m=(P2[1]-P1[1])/float(P2[0]-P1[0]) v=P1[1]-m*P1[0] x=m^2-b-P1[0]-P2[0] return (x,m*x+v) ######################################### def point_text(P,name,rgbcolor): return text(name,(P[0]-.3,P[1]+1),rgbcolor=rgbcolor) ######################################### def plot_grouplaw(b,d,e,xmin,xmax,ymin,ymax,P1,P2): PA=linedraw(P1,P2,xmin,xmax,color='green') PA+=point(P1,rgbcolor=(0,0,0),pointsize=40) PA+=point(P2,rgbcolor=(0,0,0),pointsize=40) C=compute_intersection(b,d,e,P1,P2) PA+=point(C,rgbcolor=(0,1,0),pointsize=40) #PA+=linedraw(P1,C,xmin,xmax,color='green') PA+=line([C,(C[0],-C[1])],color='green',linestyle='--') PA+=line([(C[0],-C[1]),(C[0],ymin)],color='green',linestyle='--') PA+=line([(C[0],C[1]),(C[0],ymax)],color='green',linestyle='--') PA+=point((C[0],-C[1]),rgbcolor=(0,1,0),pointsize=40) PA+=point_text(C,"$P*Q$",rgbcolor=(1,0,0)) PA+=point_text((C[0],-C[1]),"$P+Q$",rgbcolor=(1,0,0)) return PA ########################################## def linedraw(P1,P2,xmin,xmax,**keywords): var('x') if P1==P2: m=(3*P1[0]^2+2*b*P1[0]+d)/float(2*P1[1]) v=P1[1]-m*P1[0] else: m=(P2[1]-P1[1])/float(P2[0]-P1[0]) v=P1[1]-m*P1[0] f(x)=m*x+v return plot(f(x),(xmin,xmax),**keywords) ########################################## def is_on_curve(b,d,e,x): return x^3+b*x^2+d*x+e >= 0 ########################################## @interact def ellplot(xmin=slider(-20,-1,1,-3,label='$x$-min on plot: '), xmax=slider(0,30,1,5,label='$x$-max on plot: '), b=slider(-10,10,1,0,label='coefficient of $x^2$: '), d=slider(-10,10,1,0,label='coefficient of $x$: '), e=slider(-10,10,1,1,label='constant coefficient: '), xP=input_box(default=4,type=float,label='$x$-coord of P: '), xPsign=checkbox(default=True,label='desired sign of $y$-coord of P is positive? '), xQ=input_box(default=0,type=float,label='$x$-coord of Q: '), xQsign=checkbox(default=True,label='desired sign of $y$-coord of Q is positive? ')): if (4*d^3+27*e^2)==0: print 'please enter valid coefficients for a non-singular curve' else: if (not is_on_curve(b,d,e,xP)) or (not is_on_curve(b,d,e,xQ)): print 'at least one of the x-coordinates entered does not define a point on the curve over RR' else: if xPsign == True: xpsign=1 else: xpsign=-1 if xQsign == True: xqsign=1 else: xqsign=-1 E=EllipticCurve([0,b,0,d,e]) P=pointfind(b,d,e,xP,xpsign) Q=pointfind(b,d,e,xQ,xqsign) p=plot(E,xmin,xmax) p+=point_text(P,"$P$",rgbcolor=(1,0,0)) p+=point_text(Q,"$Q$",rgbcolor=(1,0,0)) ymin=p.ymin() ymax=p.ymax() p+=plot_grouplaw(b,d,e,xmin,xmax,ymin,ymax,P,Q) print E show(p) print P,' + ',Q,' = ',(compute_intersection(b,d,e,P,Q)[0],-compute_intersection(b,d,e,P,Q)[1])