restart;
with(geometry);
with(plots);
_EnvHorizontalName = 'x';
_EnvVerticalName = 'y';
point(A, -3, 9);
point(B, -5, 0);
point(C, 6, 0);
xA := -3;
yA := 9;
xB := -5;
yB := 0;
xC := 6;
yC := 0;
triangle(ABC, [A, B, C]);
midpoint(M2, A, C);
midpoint(M1, B, C);
midpoint(M3, A, B);
coordinates(M1);
coordinates(M2);
coordinates(M3);
segment(sA, [A, M1]);
segment(sB, [B, M2]);
segment(sC, [C, M3]);
centroid(G, ABC);
midpoint(J1, A, G);
midpoint(J2, B, G);
midpoint(J3, C, G);
coordinates(J1);
coordinates(J2);
coordinates(J3);
zA := xA + yA*I;
zB := xB + yB*I;
zC := xC + yC*I:
diff((z - zA)*(z - zB)*(z - zC), z)*I;);
allvalues(%, z)
;I am trying to find focus of Steiner ellipse with Marden's theoreme. Thank you.

I would like to apply the pappus theorem to an ellipse. Pascale's theorem 

I am trying to decompose an isprime into sum of 2 squares.
Can you tell me why yhse procedure are not goog.
                       

Sumof2Squares:= proc(p::And(prime, satisfies(p-> irem(p,4)=1)))
local x, y:= 1;
   x:= mods(Roots(x^2+y^2), p)[2,1];
   while x^2+y^2 > p do
      (x,y):= FermatDescent(x,y,p)
   end do;
   (x,y)
end proc:

FermatDescent:= proc(x::posint, y::posint, p::posint)
local 
   m:= (x^2+y^2)/p,
   a:= mods(x,m),  
   b:= mods(y,m);

   (abs((a*x+b*y)/m), abs((a*y-b*x)/m))
end proc:
   
trace(FermatDescent);

Sumof2Squares(1973);
Thank you.

I should like to find the list of Friday 13 of a year

restart;  
with(geometry):  
with(plots):  
_EnvHorizontalName = 'x':  _EnvVerticalName = 'y':
point(A, -1, 9):                                                                                                       
point(B, -5, 0):
point(C, 6, 0):
triangle(ABC,[A,B,C]):
midpoint(M1,A,C): midpoint(M2,B,C):midpoint(M3,A,B):
rotation(J, C, Pi/2, 'counterclockwise', M1):triangle(AJC,[A,J,C]):
rotation(Ii, C, Pi/2, 'counterclockwise', M2):triangle(BIC,[B,Ii,C]):
rotation(K, A, Pi/2, 'counterclockwise', M3):triangle(AKB,[A,K,B]):
midpoint(O1,K,J): coordinates(O1):
midpoint(O2,A,Ii): coordinates(O2):  
poly:=[coordinates(A),coordinates(J),coordinates(Ii),coordinates(K)]:   

display(draw([A(color = black, symbol = solidcircle, symbolsize = 12), 
B(color = black, symbol = solidcircle, symbolsize = 12), 
C(color = black, symbol = solidcircle, symbolsize = 12), 
J(color = black, symbol = solidcircle, symbolsize = 12), 
polygonplot(poly,color = "DarkGreen", transparency = 0.5),
ABC(color = red ),
BIC(color = green),
AKB(color = grey),
AJC(color =blue)]),
textplot([[coordinates(A)[], "A"],[coordinates(J)[], "J"],[coordinates(Ii)[], "I"],   
[coordinates(B)[], "B"], [coordinates(K)[], "K"], 
[coordinates(C)[], "C"]], 
align = [above, right]),  axes = none);
Error, (in geometry:-draw) the option must be of type equation or name. I don't see how to correct this error/