## 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: