## The procedure that does all the job: read CP4; ############################################################ ## Some testing code: ## print(`The procedure C4 takes a 4th degree polynomial as argument and returns a list. The first element is Argmax: the point where the maximum is attained. The second element is the conjugate of P. If polynomial P is not convex there are remaining elements: u, v, x1, x2, P'(x1), and P'(x2)`); ##The unitary case... lprint(`The unitary case...`): P:=x^4; L:=C4(P); assume(s,real);lprint(`The argmax is `):L[1](s); lprint(` simplifying gives:`):simplify("): convert(",abs):convert(",piecewise,s):simplify("); lprint(`The conjugate is `):L[2](s); ## other examples... lprint(`Other examples...`): P:='P': for P in [5*(x^2-1)^2,3*(x+4)*(x+2)*(x-3)*(x-10),sqrt(3)*(x+4)^2*(x-3)*(x-10)] do print(P): L:=C4(P): if nops(L)>2 then d:=(L[6]-L[5])/8: print(plot({P,L[3]*x-L[4]},x=L[5]-d..L[6]+d,title='P')); d:=(L[8]-L[7])/4: print(plot(L[1](s),s=L[7]-d..L[8]+d,title=`X(s):=Argmax`)); p1:=plot(L[2](s),s=L[7]-d..L[8]+d,title=`P*(s)`,color=red,thickness=3): t:='t':dP:=diff(P,x):d:=(L[6]-L[5])/8: p2:=plot([dP,x*dP-P,x=L[5]-d..L[6]+d],color=green): print(plots[display]([p2,p1])); else print(plot(P,x=-100..100,title='P')); print(plot(L[1](s),s=-100..100,title=`X(s):=Argmax`)); p1:=plot(L[2](s),s=-100..100,title=`P*(s)`,color=red,thickness=3): t:='t':dP:=diff(P,x): p2:=plot([dP,x*dP-P,x=L[1](-100)..L[1](100)],color=green): print(plots[display]([p1,p2])); fi: od: ## Convex case lprint(`Random convex polynomials`): i:='i':for i to 3 do a:=rand(-100..100):a0:=a(): b:=rand(floor(3*a0^4/8)..2*abs(a0)^4):b0:=b(): P:=a()*(x^4+a0^2*x^3+b0*x^2+a()*x+a());print(`P=`,eval(P)): print(`convex =`,evalb(3*a0^4-8*b0<0)); print(P): L:=C4(P): if nops(L)>2 then d:=(L[6]-L[5])/8: print(plot({P,L[3]*x-L[4]},x=L[5]-d..L[6]+d,title='P')); d:=(L[8]-L[7])/4: # print(plot(L[1](s),s=L[7]-d..L[8]+d,title=`X(s):=Argmax`)); p1:=plot(L[2](s),s=L[7]-d..L[8]+d,title=`P*(s)`,color=red,thickness=3): t:='t':dP:=diff(P,x):d:=(L[6]-L[5])/8: p2:=plot([dP,x*dP-P,x=L[5]-d..L[6]+d],color=green): print(plots[display]([p2,p1])); else print(plot(P,x=-100..100,title='P')); print(plot(L[1](s),s=-100..100,title=`X(s):=Argmax`)); p1:=plot(L[2](s),s=-100..100,title=`P*(s)`,color=red,thickness=3): t:='t':dP:=diff(P,x): p2:=plot([dP,x*dP-P,x=L[1](-100)..L[1](100)],color=green): print(plots[display]([p1,p2])); fi: od: ## Non-convex case lprint(`Random polynomials`): i:='i':for i to 3 do P:=sort(randpoly(x,degree=4)); print(P): L:=C4(P): if nops(L)>2 then d:=(L[6]-L[5])/8: print(plot({P,L[3]*x-L[4]},x=L[5]-d..L[6]+d,title='P')); d:=(L[8]-L[7])/4: print(plot(L[1](s),s=L[7]-d..L[8]+d,title=`X(s):=Argmax`)); p1:=plot(L[2](s),s=L[7]-d..L[8]+d,title=`P*(s)`,color=red,thickness=3): t:='t':dP:=diff(P,x):d:=(L[6]-L[5])/8: p2:=plot([dP,x*dP-P,x=L[5]-d..L[6]+d],color=green): print(plots[display]([p2,p1])); else print(plot(P,x=-100..100,title='P')); # print(plot(L[1](s),s=-100..100,title=`X(s):=Argmax`)); p1:=plot(L[2](s),s=-100..100,title=`P*(s)`,color=red,thickness=3): t:='t':dP:=diff(P,x): p2:=plot([dP,x*dP-P,x=L[1](-100)..L[1](100)],color=green): print(plots[display]([p1,p2])); fi: od: