%! % PostScript prologue for pst-eucl.tex. % Version 1.04 2020/09/29 % For distribution, see pstricks.tex. % /tx@EcldDict 40 dict def tx@EcldDict begin %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Prcision /epsilon 1E-5 def /epsilonMin 1E-6 def /epsilonMax 1E-3 def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Pi /Pi 3.14159265359 def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% e /E 2.718281828459045 def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x -> true (if |x| < epsilonMin) /ZeroEq { abs epsilonMin lt } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x1 y1 x2 y2 -> a b c (ax-by+c=0 with a^2+b^2=1) /EqDr { 4 copy 3 -1 roll sub 7 1 roll exch sub 5 1 roll 4 -1 roll mul 3 1 roll mul exch sub 2 index dup mul 2 index dup mul add sqrt 4 -1 roll 1 index div exch 4 -1 roll 1 index div exch 4 -1 roll 1 index div exch pop } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% orthogonal projection of M1 onto (OM2) %% x1 y1 x2 y2 -> x3 y3 /Project { 2 copy dup mul exch dup mul add 5 1 roll 2 copy 5 -1 roll mul exch 5 -1 roll mul add 4 -1 roll div dup 4 -1 roll mul exch 3 -1 roll mul } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% a b c (ax2+bx+c=0) -> x1 y1 /SolvTrin { /c exch def /b exch def /a exch def b dup mul a c mul 4 mul sub dup 0 lt { pop 0 0 } %% no solutions {sqrt dup b neg add a 2 mul div exch b add neg 2 a mul div } ifelse } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x1 y1 x2 y2 -> Dist /ABDist { 3 -1 roll sub dup mul 3 1 roll sub dup mul add sqrt } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x1 y1 x2 y2 -> x2-x1 y2-y1 /ABVect { 3 -1 roll exch sub 3 1 roll sub exch } bind def %/ABVect { 3 -1 roll sub 3 1 roll exch sub exch } bind def %% wrong version %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x1 y1 x2 y2 x3 y3 x4 y4 -> x y /InterLines { EqDr /D1c exch def /D1b exch def /D1a exch def EqDr /D2c exch def /D2b exch def /D2a exch def D1a D2b mul D1b D2a mul sub dup ZeroEq % { pop pop pop 0 0 } %% parallel lines % --- hv 20110714 { pop 0 0 } %% parallel lines --- hv 20110714 { /Det exch def D1b D2c mul D1c D2b mul sub Det div D1a D2c mul D2a D1c mul sub Det div } ifelse } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% a b c R -> x1 y1 x2 y2 /InterLineCircle { /CR exch def /Dc exch def neg /Db exch def /Da exch def ABVect /Vy exch def /Vx exch def %% Dc==0 then O belong to the line %% First project O on the line -> M (-ca;-cb) %% l'abscisse de M sur (OM) divisee par R donne le cosinus %Dc neg dup Db mul exch Da mul 2 copy 0 0 %ABDist dup CR gt { pop pop pop 0 0 0 0 } %{ ZeroEq { pop pop Db Da } if Atan /alpha exch def Dc abs CR gt { 0 0 0 0 } { Db neg Da neg Atan /alpha exch def Dc CR div dup dup mul 1 exch sub sqrt exch Atan /beta exch def alpha beta add dup cos CR mul exch sin CR mul alpha beta sub dup cos CR mul exch sin CR mul 4 copy ABVect Vy mul 0 le exch Vx mul 0 le and { 4 2 roll } if } ifelse } def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% R R' OO' -> x1 y1 x2 y2 /InterCircles { /OOP exch def /CRP exch def /CR exch def OOP dup mul CRP dup mul sub CR dup mul add OOP div 2 div dup dup mul CR dup mul exch sub dup 0 lt { pop pop 0 0 0 0 } { sqrt 2 copy neg } ifelse } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x y theta -> x' y' (rotation of theta) /Rotate { dup sin /sintheta exch def cos /costheta exch def /y exch def /x exch def x costheta mul y sintheta mul sub y costheta mul x sintheta mul add } def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% N -> x y /GetNode { tx@NodeDict begin tx@NodeDict 1 index known { load GetCenter } { pop 0 0 } ifelse end } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x -> ch(x) /ch { dup Ex exch neg Ex add 2 div } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x -> sh(x) /sh { dup Ex exch neg Ex sub 2 div } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x -> e^(x) /Ex { E exch exp } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% x f g -> x y n /NewtonSolving { /g exch def /f exch def 0 { %%% STACK: x0 n 1 add exch %% one more loop dup ZeroEq { dup 0.0005 add fgeval 1 index 0.0005 sub fgeval sub .001 div } { dup 1.0005 mul fgeval 1 index 0.9995 mul fgeval sub .001 2 index mul div } ifelse %%% STACK: n x0 fg'(x0) %%% compute x1=x0-fg(x0)/fg'(x0) 1 index fgeval exch div dup 4 1 roll sub exch %% stack: dx x0 n 3 -1 roll ZeroEq %% exit if root found 1 index 100 eq or { exit } if %% or looping for more than 100 times } loop dup 100 lt { exch dup /x exch def f } { pop 0 0 } ifelse 3 -1 roll } def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /fgeval { /x exch def f g sub } bind def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% calculate the line coefficents Ax+By+C=0 %% x1 y1 x2 y2 -> A B C /LineCoefABC { 0 index 3 index sub % A=y2-y1 4 index 3 index sub % B=x1-x2 3 index 5 index mul 6 index 4 index mul sub % C=x2y1-x1y2 7 3 roll pop pop pop pop } def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% calculate the 2-order determinant %% |a11 a12| %% |a21 a22| %% a11 a12, a21 a22 -> X /DeterminantTwo { 4 1 roll mul 3 1 roll mul exch sub } def %% calculate the 3-order determinant %% |a11 a12 a13| %% |a21 a22 a23| %% |a31 a32 a33| %% 8 7 6 5 4 3 2 1 0 %% a11 a12 a13, a21 a22 a23, a31 a32 a33 -> X /DeterminantThree { % |a22 a23, a32 a33| * (-1)^(1+1)a11 8 index abs epsilon lt { %a11=0 0 } { 4 index 4 index 3 index 3 index DeterminantTwo 9 index mul } ifelse % |a12 a13, a32 a33| * (-1)^(1+2)a21 6 index abs epsilon lt { %a12=0 0 sub } { 8 index 8 index 4 index 4 index DeterminantTwo 7 index mul sub } ifelse % |a12 a13, a22 a23| * (-1)^(1+3)a31 3 index abs epsilon lt { %a13=0 0 add } { 8 index 8 index 7 index 7 index DeterminantTwo 4 index mul add } ifelse 10 1 roll pop pop pop pop pop pop pop pop pop } def %% calculate the 4-order determinant %% |a11 a12 a13 a14| %% |a21 a22 a23 a24| %% |a31 a32 a33 a34| %% |a41 a42 a43 a44| %% 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 %% a11 a12 a13 a14, a21 a22 a23 a24, a31 a32 a33 a34, a41 a42 a43 a44 -> X /DeterminantFour { % |a22 a23 a24, a32 a33 a34, a42 a43 a44| * (-1)^(1+1)a11 15 index abs epsilon lt { %a11=0 0 } { 10 index 10 index 10 index 9 index 9 index 9 index 8 index 8 index 8 index DeterminantThree 16 index mul } ifelse % |a12 a13 a14, a32 a33 a34, a42 a43 a44| * (-1)^(1+2)a21 12 index abs epsilon lt { %a21=0 0 sub } { 15 index 15 index 15 index 10 index 10 index 10 index 9 index 9 index 9 index DeterminantThree 13 index mul sub } ifelse % |a12 a13 a14, a22 a23 a24, a42 a43 a44| * (-1)^(1+3)a31 8 index abs epsilon lt { %a31=0 0 add } { 15 index 15 index 15 index 14 index 14 index 14 index 9 index 9 index 9 index DeterminantThree 9 index mul add } ifelse % |a12 a13 a14, a22 a23 a24, a32 a33 a34| * (-1)^(1+4)a41 4 index abs epsilon lt { %a41=0 0 sub } { 15 index 15 index 15 index 14 index 14 index 14 index 13 index 13 index 13 index DeterminantThree 5 index mul sub } ifelse 17 1 roll pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop } def %% calculate the 5-order determinant %% |a11 a12 a13 a14 a15| %% |a21 a22 a23 a24 a25| %% |a31 a32 a33 a34 a35| %% |a41 a42 a43 a44 a45| %% |a51 a52 a53 a54 a55| %% 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 %% a11 a12 a13 a14 a15, a21 a22 a23 a24 a25, a31 a32 a33 a34 a35, a41 a42 a43 a44 a45, a51 a52 a53 a54 a55-> X /DeterminantFive { % |a22 a23 a24 a25, a32 a33 a34 a35, a42 a43 a44 a45, a52 a53 a54 a55| * (-1)^(1+1)a11 24 index abs epsilon lt { %a11=0 0 } { 18 index 18 index 18 index 18 index 17 index 17 index 17 index 17 index 16 index 16 index 16 index 16 index 15 index 15 index 15 index 15 index DeterminantFour 25 index mul } ifelse % |a12 a13 a14 a15, a32 a33 a34 a35, a42 a43 a44 a45, a52 a53 a54 a55| * (-1)^(1+2)a21 20 index abs epsilon lt { %a21=0 0 sub } { 24 index 24 index 24 index 24 index 18 index 18 index 18 index 18 index 17 index 17 index 17 index 17 index 16 index 16 index 16 index 16 index DeterminantFour 21 index mul sub } ifelse % |a12 a13 a14 a15, a22 a23 a24 a25, a42 a43 a44 a45, a52 a53 a54 a55| * (-1)^(1+3)a31 15 index abs epsilon lt { %a31=0 0 add } { 24 index 24 index 24 index 24 index 23 index 23 index 23 index 23 index 17 index 17 index 17 index 17 index 16 index 16 index 16 index 16 index DeterminantFour 16 index mul add } ifelse % |a12 a13 a14 a15, a22 a23 a24 a25, a32 a33 a34 a35, a52 a53 a54 a55| * (-1)^(1+4)a41 10 index abs epsilon lt { %a41=0 0 sub } { 24 index 24 index 24 index 24 index 23 index 23 index 23 index 23 index 22 index 22 index 22 index 22 index 16 index 16 index 16 index 16 index DeterminantFour 11 index mul sub } ifelse % |a12 a13 a14 a15, a22 a23 a24 a25, a32 a33 a34 a35, a42 a43 a44 a45| * (-1)^(1+5)a51 5 index abs epsilon lt { %a51=0 0 add } { 24 index 24 index 24 index 24 index 23 index 23 index 23 index 23 index 22 index 22 index 22 index 22 index 21 index 21 index 21 index 21 index DeterminantFour 6 index mul add } ifelse 26 1 roll pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop pop } def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Conic Intersections % The following macros implements the conic intersectionpoints in Asymptote % module geometry.asy provided by Philippe IVALDI. % http://www.piprime.fr/files/asymptote/geometry/modules/geometry.asy.html %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % num len [array] numIsInArray -> bool /numIsInArray { /arr ED /len ED /num ED /res false def 0 1 len 1 sub { /idx ED /val arr idx get def num val sub abs epsilonMax lt { /res true def exit } if } for res } def % find the real roots of quadratic equation ax^2+bx+c=0 % a b c -> [roots] nroots /QuadraticRealRoots { 15 dict begin % all variables are local /Coefc ED /Coefb ED /Coefa ED /nroots 0 def /roots 2 array def %[ Coefa Coefb Coefc ] == Coefa abs epsilonMin lt { Coefb abs epsilonMin gt { roots 0 Coefc Coefb div neg put /nroots 1 def } }{ /delta Coefb Coefb mul 4 Coefc Coefa mul mul sub def delta abs epsilon lt { roots 0 Coefb 2 Coefa mul div neg put /nroots 1 def }{ delta 0.0 gt { /delta delta sqrt def roots 0 Coefb neg delta add 2 Coefa mul div put roots 1 Coefb neg delta sub 2 Coefa mul div put /nroots 2 def } if } ifelse } ifelse roots nroots % push the roots on stack. end} def % find the real roots of cubic equation a*x^3+b*x^2+c*x+d=0 % a b c d CubicRealRoots -> [roots] nroots /CubicRealRoots { 15 dict begin % all variables are local /Coefd ED /Coefc ED /Coefb ED /Coefa ED % [ Coefa Coefb Coefc Coefd ] == /nroots 0 def /roots 3 array def Coefa abs epsilonMin lt { % quadratic case Coefb Coefc Coefd QuadraticRealRoots /nroots ED roots copy pop }{% true cubic % normalize to x^3+bx^2+cx+d=0 /Coefb Coefb Coefa div def /Coefc Coefc Coefa div def /Coefd Coefd Coefa div def % let x=y-b/3, we have % y^3+py+q=0 % where p=(3c-b^2)/3,q=(27d-9bc+2b^3)/27 % let y=u+v, where uv=-p/3, then % u^3+v^3=-q, u^3v^3=-p^3/27 % then u^3=-q/2-\sqrt(q^2/4+p^3/27) % u^3=-q/2+\sqrt(q^2/4+p^3/27) % http://eqworld.ipmnet.ru/en/solutions/ae/ae0103.pdf Coefb abs epsilonMin lt { /bthird 0 def /p Coefc def /q Coefd def }{ /bthird Coefb 3 div def /p 3 Coefc mul Coefb Coefb mul sub 3 div def /q 27 Coefd mul 9 Coefb mul Coefc mul sub 2 Coefb mul Coefb mul Coefb mul add 27 div def } ifelse /ppp p dup dup mul mul def /qq q dup mul def /delta qq 4 div ppp 27 div add def % p = q = delta = delta abs epsilon lt { q abs epsilon lt { /nroots 1 def roots 0 bthird neg put }{ /nroots 2 def q 0.0 gt { /r1 q 2 div 1 3 div exp def }{ /r1 q 2 div neg 1 3 div exp neg def } ifelse /r2 r1 -2 mul def roots 0 r1 bthird sub put roots 1 r2 bthird sub put } ifelse }{ delta 0.0 gt { /nroots 1 def /qhalfneg q 2 div neg def /deltasqrt delta sqrt def /pthirdneg p 3 div neg def p abs epsilon lt { /tu qhalfneg deltasqrt add def /tv qhalfneg deltasqrt sub def tu abs epsilon lt { /u 0.0 def }{ tu 0.0 gt { /u tu 1 3 div exp def }{ /u tu neg 1 3 div exp neg def } ifelse } ifelse tv abs epsilon lt { /v 0.0 def }{ tv 0.0 gt { /v tv 1 3 div exp def }{ /v tv neg 1 3 div exp neg def } ifelse } ifelse }{ /tu qhalfneg deltasqrt add def tu 0.0 gt { /u tu 1 3 div exp def }{ /u tu neg 1 3 div exp neg def } ifelse /v pthirdneg u div def } ifelse /r1 u v add def roots 0 r1 bthird sub put }{ /nroots 3 def /qhalfneg q 2 div neg def /pthirdnegsqrt p 3 div neg sqrt def /argalpha qhalfneg pthirdnegsqrt dup dup mul mul div Acos 3 div def /r1 pthirdnegsqrt 2 mul argalpha cos mul def /r2 pthirdnegsqrt -2 mul argalpha 60 add cos mul def /r3 pthirdnegsqrt -2 mul argalpha 60 sub cos mul def roots 0 r1 bthird sub put roots 1 r2 bthird sub put roots 2 r3 bthird sub put } ifelse } ifelse } ifelse roots nroots % push the roots on stack. end} def % find the real roots of quartic equation a*x^4+b*x^3+c*x^2+dx+e=0 % a b c d e QuarticRealRoots -> [roots] nroots /QuarticRealRoots { 15 dict begin % all variables are local /Coefe ED /Coefd ED /Coefc ED /Coefb ED /Coefa ED /nroots 0 def /roots 4 array def Coefa abs epsilonMin lt { % cubic case Coefb Coefc Coefd Coefe CubicRealRoots /nroots ED roots copy pop }{ Coefe abs epsilonMin lt { % cubic case Coefa Coefb Coefc Coefd CubicRealRoots /nroots ED roots copy pop roots nroots 0.0 put /nroots nroots 1 add def }{ % true quartic % normalize to x^4+bx^3+cx^2+dx+e=0 /Coefb Coefb Coefa div def /Coefc Coefc Coefa div def /Coefd Coefd Coefa div def /Coefe Coefe Coefa div def /qeval { /vx ED 1 vx mul Coefb add vx mul Coefc add vx mul Coefd add vx mul Coefe add % x(x(x(ax+b)+c)+d)+e } def % [1 Coefb Coefc Coefd Coefe ] == % x^4+bx^3+cx^2+dx+e=0 % (x^2+b/2x)^2=(b^2/4-c)x^2-dx-e % add (x^2+b/2x)y+y^2/4 in each side, % (x^2+b/2x+y/2)^2=(b^2/4+y-c)x^2+(by/2-d)x+(y^2/4-e) % choose y such that % (by/2-d)^2-4(b^2/4+y-c)(y^2/4-e)=0 % i.e, % y^3-cy^2+(bd-4e)y-eb^2+4ec-d^2=0 % let t=sqrt(b^2/4+y-c), then % (x^2+b/2x+y/2)^2=(tx+(by/2-d)/(2t))^2 % we have % x^2+(b/2+t)x+y/2+(by/2-d)/(2t)=0 % x^2+(b/2-t)x+y/2-(by/2-d)/(2t)=0 1 Coefc neg Coefb Coefd mul 4 Coefe mul sub 4 Coefc Coefe mul mul Coefb Coefb Coefe mul mul sub Coefd Coefd mul sub CubicRealRoots /CubicNumRoots ED /CubicRoots ED %CubicNumRoots = CubicRoots == 0 1 CubicNumRoots 1 sub { /idx ED /y0 CubicRoots idx get def /delta y0 Coefc sub Coefb Coefb mul 4 div add def % [ idx y0 delta ] == delta abs epsilonMax lt { /squareval y0 y0 mul 4 div Coefe sub def squareval abs epsilon lt { % x^2+b/2x+y/2=0 %[squareval (squareval=0) ] == 1 Coefb 2 div y0 2 div QuadraticRealRoots /nroots1 ED /roots1 ED %nroots1 = roots1 == 0 1 nroots1 1 sub { /idx ED /rv roots1 idx get def rv qeval abs epsilonMax lt { rv nroots roots numIsInArray not { roots nroots rv put /nroots nroots 1 add def } if } if } for }{ squareval 0.0 gt { %[squareval (squareval>0) ] == % x^2+b/2x+y/2\pm\sqrt(y^2/4-e)=0 /squareval squareval sqrt def 1 Coefb 2 div y0 2 div squareval add QuadraticRealRoots /nroots1 ED /roots1 ED %nroots1 = roots1 == 1 Coefb 2 div y0 2 div squareval sub QuadraticRealRoots /nroots2 ED /roots2 ED %nroots2 = roots2 == nroots1 0 gt nroots2 0 gt or { 0 1 nroots1 1 sub { /idx ED /rv roots1 idx get def rv qeval abs epsilonMax lt { rv nroots roots numIsInArray not { roots nroots rv put /nroots nroots 1 add def } if } if } for 0 1 nroots2 1 sub { /idx ED /rv roots2 idx get def rv qeval abs epsilonMax lt { rv nroots roots numIsInArray not { roots nroots rv put /nroots nroots 1 add def } if } if } for } if } if } ifelse }{ delta 0.0 gt { %[delta (delta>0) ] == /sqrtdelta delta sqrt def /Coefp Coefb 2 div sqrtdelta add def /Coefq y0 2 div y0 Coefb mul 2 div Coefd sub sqrtdelta div 2 div add def 1 Coefp Coefq QuadraticRealRoots /nroots1 ED /roots1 ED %nroots1 = roots1 == /Coefp Coefb 2 div sqrtdelta sub def /Coefq y0 2 div y0 Coefb mul 2 div Coefd sub sqrtdelta div 2 div sub def 1 Coefp Coefq QuadraticRealRoots /nroots2 ED /roots2 ED %nroots2 = roots2 == nroots1 0 gt nroots2 0 gt or { 0 1 nroots1 1 sub { /idx ED /rv roots1 idx get def rv qeval abs epsilonMax lt { rv nroots roots numIsInArray not { roots nroots rv put /nroots nroots 1 add def } if } if } for 0 1 nroots2 1 sub { /idx ED /rv roots2 idx get def rv qeval abs epsilonMax lt { rv nroots roots numIsInArray not { roots nroots rv put /nroots nroots 1 add def } if } if } for } if } if } ifelse nroots 4 eq { exit } if } for } ifelse } ifelse roots nroots % push the roots on stack. end} def % SortInters % sort the intersections from smallest to largest abscissa. % [Inters] nInters SortInters -> [Inters] nInters /SortInters { /nInters ED /Inters ED nInters 1 gt { /PairCmp { /PairB ED /PairA ED PairA 0 get PairB 0 get sub } def /PairSwap { /PairB ED /PairA ED /IdxB ED /IdxA ED Inters IdxA PairB put Inters IdxB PairA put } def 0 1 nInters 2 sub { /Idx ED 0 1 nInters 2 sub Idx sub { /Bubidx ED /PairA Inters Bubidx get def /PairB Inters Bubidx 1 add get def PairA PairB PairCmp 0 gt { Bubidx Bubidx 1 add PairA PairB PairSwap } if } for } for } if Inters nInters % push the intersections on stack. } def % ConicCircleInter % find all intersections of conic $cc(x,y)=ax^2+bxy+cy^2+dx+ey+f=0$ and % circle $c(x,y)=(x-x0)^2+(y-y0)^2=r^2$. % % a b c d e f x0 y0 r ConicCircleInter -> [Inters] nInters /ConicCircleInter { 15 dict begin % all variables are local /cr ED /cy0 ED /cx0 ED /ccf ED /cce ED /ccd ED /ccc ED /ccb ED /cca ED %[ cca ccb ccc ccd cce ccf ] == % x y cceval -> cc(x,y) /cceval { /ty ED /tx ED cca tx mul tx mul ccb tx mul ty mul add ccc ty mul ty mul add ccd tx mul add cce ty mul add ccf add } def /m 2 cca mul cx0 mul ccb cy0 mul add ccd add def /n 2 ccc mul cy0 mul ccb cx0 mul add cce add def /s cx0 cy0 cceval def /p s cr div ccc cr mul add def /t cca cr mul ccc cr mul sub def % [ m n s p t ] == /Qqa ccb ccb mul cr mul cr mul t t mul add def /Qqb 2 t mul m mul 2 n mul ccb mul cr mul add def /Qqc m m mul n n mul add 2 t mul p mul add ccb ccb mul cr mul cr mul sub def /Qqd 2 m mul p mul 2 n mul ccb mul cr mul sub def /Qqe p p mul n n mul sub def %[ Qqa Qqb Qqc Qqd Qqe ] == Qqa Qqb Qqc Qqd Qqe QuarticRealRoots /nroots ED /roots ED %nroots = roots == /nInters 0 def /Inters 4 array def /SaveInter { /tty ED /ttx ED nInters 4 lt { Inters nInters [ttx tty] put /nInters nInters 1 add def } if } def 0 1 nroots 1 sub { /idx ED /cosx roots idx get def /sinx 1.0 cosx cosx mul sub sqrt def /x cr cosx mul cx0 add def /y cr sinx mul cy0 add def /ccxy x y cceval def % [idx 0 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if /y cr sinx mul neg cy0 add def /ccxy x y cceval def % [idx 1 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if } for Inters nInters SortInters% push the intersections on stack. end} def % ConicEllipseInter % find all intersections of conic $cc(x,y)=ax^2+bxy+cy^2+dx+ey+f=0$ and % ellipse $e(x,y): (x-x0)^2/m^2+(y-y0)^2/n^2=1$. % % a b c d e f x0 y0 m n ConicEllipseInter -> [Inters] nInters /ConicEllipseInter { 15 dict begin % all variables are local /en ED /em ED /ey0 ED /ex0 ED /ccf ED /cce ED /ccd ED /ccc ED /ccb ED /cca ED %[ cca ccb ccc ccd cce ccf ] == % x y cceval -> cc(x,y) /cceval { /ty ED /tx ED cca tx mul tx mul ccb tx mul ty mul add ccc ty mul ty mul add ccd tx mul add cce ty mul add ccf add } def /p 2 cca mul em mul ex0 mul ccb em mul ey0 mul add ccd em mul add def /q 2 ccc mul en mul ey0 mul ccb en mul ex0 mul add cce en mul add def /r ex0 ey0 cceval def /s r ccc en mul en mul add def /t cca em mul em mul ccc en mul en mul sub def % [ p q r s t ] == /Qqa t t mul ccb ccb mul em mul em mul en mul en mul add def /Qqb 2 t mul p mul 2 q mul ccb mul em mul en mul add def /Qqc p p mul q q mul add 2 t mul s mul add ccb ccb mul em mul em mul en mul en mul sub def /Qqd 2 p mul s mul 2 q mul ccb mul em mul en mul sub def /Qqe s s mul q q mul sub def %[ Qqa Qqb Qqc Qqd Qqe ] == Qqa Qqb Qqc Qqd Qqe QuarticRealRoots /nroots ED /roots ED %nroots = roots == /nInters 0 def /Inters 4 array def /SaveInter { /tty ED /ttx ED nInters 4 lt { Inters nInters [ttx tty] put /nInters nInters 1 add def } if } def 0 1 nroots 1 sub { /idx ED /cosx roots idx get def /sinx 1.0 cosx cosx mul sub sqrt def /x em cosx mul ex0 add def /y en sinx mul ey0 add def /ccxy x y cceval def % [idx 0 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if /y en sinx mul neg ey0 add def /ccxy x y cceval def % [idx 1 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if } for Inters nInters SortInters% push the intersections on stack. end} def % ConicHyperbolaInter % find all intersections of conic $cc(x,y)=ax^2+bxy+cy^2+dx+ey+f=0$ and % hyperbola $h(x,y): (x-x0)^2/m^2-(y-y0)^2/n^2=1$. % % a b c d e f x0 y0 m n ConicHyperbolaInter -> [Inters] nInters /ConicHyperbolaInter { 15 dict begin % all variables are local /hn ED /hm ED /hy0 ED /hx0 ED /ccf ED /cce ED /ccd ED /ccc ED /ccb ED /cca ED %[ cca ccb ccc ccd cce ccf ] == % x y cceval -> cc(x,y) /cceval { /ty ED /tx ED cca tx mul tx mul ccb tx mul ty mul add ccc ty mul ty mul add ccd tx mul add cce ty mul add ccf add } def /p 2 cca mul hm mul hx0 mul ccb hm mul hy0 mul add ccd hm mul add def /q 2 ccc mul hn mul hy0 mul ccb hn mul hx0 mul add cce hn mul add def /r hx0 hy0 cceval def /s r ccc hn mul hn mul sub def /t cca hm mul hm mul ccc hn mul hn mul add def % [ p q r s t ] == /Qqe t t mul ccb ccb mul hm mul hm mul hn mul hn mul sub def /Qqd 2 t mul p mul 2 q mul ccb mul hm mul hn mul sub def /Qqc p p mul q q mul sub 2 t mul s mul add ccb ccb mul hm mul hm mul hn mul hn mul add def /Qqb 2 p mul s mul 2 q mul ccb mul hm mul hn mul add def /Qqa s s mul q q mul add def %[ Qqa Qqb Qqc Qqd Qqe ] == Qqa Qqb Qqc Qqd Qqe QuarticRealRoots /nroots ED /roots ED % nroots = roots == /nInters 0 def /Inters 4 array def /SaveInter { /tty ED /ttx ED nInters 4 lt { Inters nInters [ttx tty] put /nInters nInters 1 add def } if } def 0 1 nroots 1 sub { /idx ED /cosx roots idx get def /sinx 1.0 cosx cosx mul sub sqrt def /x hm cosx div hx0 add def /y hn sinx mul cosx div hy0 add def /ccxy x y cceval def % [idx 0 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if /y hn sinx mul cosx div neg hy0 add def /ccxy x y cceval def % [idx 1 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if } for Inters nInters SortInters% push the intersections on stack. end} def % ConicIHyperbolaInter % find all intersections of conic $cc(x,y)=ax^2+bxy+cy^2+dx+ey+f=0$ and % conjugate hyperbola $h(x,y): (y-y0)^2/m^2-(x-x0)^2/n^2=1$. % % a b c d e f x0 y0 m n ConicIHyperbolaInter -> [Inters] nInters /ConicIHyperbolaInter { 15 dict begin % all variables are local /hn ED /hm ED /hy0 ED /hx0 ED /ccf ED /cce ED /ccd ED /ccc ED /ccb ED /cca ED % [ cca ccb ccc ccd cce ccf ] == % x y cceval -> cc(x,y) /cceval { /ty ED /tx ED cca tx mul tx mul ccb tx mul ty mul add ccc ty mul ty mul add ccd tx mul add cce ty mul add ccf add } def /p 2 ccc mul hm mul hy0 mul ccb hm mul hx0 mul add cce hm mul add def /q 2 cca mul hn mul hx0 mul ccb hn mul hy0 mul add ccd hn mul add def /r hx0 hy0 cceval def /s r cca hn mul hn mul sub def /t cca hn mul hn mul ccc hm mul hm mul add def % [ p q r s t ] == /Qqe t t mul ccb ccb mul hm mul hm mul hn mul hn mul sub def /Qqd 2 t mul p mul 2 q mul ccb mul hm mul hn mul sub def /Qqc p p mul q q mul sub 2 t mul s mul add ccb ccb mul hm mul hm mul hn mul hn mul add def /Qqb 2 p mul s mul 2 q mul ccb mul hm mul hn mul add def /Qqa s s mul q q mul add def %[ Qqa Qqb Qqc Qqd Qqe ] == Qqa Qqb Qqc Qqd Qqe QuarticRealRoots /nroots ED /roots ED %nroots = roots == /nInters 0 def /Inters 4 array def /SaveInter { /tty ED /ttx ED nInters 4 lt { Inters nInters [ttx tty] put /nInters nInters 1 add def } if } def 0 1 nroots 1 sub { /idx ED /cosx roots idx get def /sinx 1.0 cosx cosx mul sub sqrt def /x hn sinx mul cosx div hx0 add def /y hm cosx div hy0 add def /ccxy x y cceval def % [idx 0 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if /x hn sinx mul cosx div neg hx0 add def /ccxy x y cceval def % [idx 1 cosx sinx x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if } for Inters nInters SortInters% push the intersections on stack. end} def % ConicParabolaInter % find all intersections of conic $cc(x,y)=ax^2+bxy+cy^2+dx+ey+f=0$ and % parabola $p(x,y): (x-x0)^2=2p(y-y0)$. % % a b c d e f x0 y0 p ConicParabolaInter -> [Inters] nInters /ConicParabolaInter { 15 dict begin % all variables are local /pp ED /py0 ED /px0 ED /ccf ED /cce ED /ccd ED /ccc ED /ccb ED /cca ED % [ cca ccb ccc ccd cce ccf ] == % x y cceval -> cc(x,y) /cceval { /ty ED /tx ED cca tx mul tx mul ccb tx mul ty mul add ccc ty mul ty mul add ccd tx mul add cce ty mul add ccf add } def /Qqa ccc 4 pp pp mul mul div def /Qqb ccb 2 pp mul div def /Qqc cca ccb px0 mul 2 pp mul div add ccc py0 mul pp div add cce 2 pp mul div add def /Qqd 2 cca px0 mul mul ccb py0 mul add ccd add def /Qqe px0 py0 cceval def Qqa Qqb Qqc Qqd Qqe QuarticRealRoots /nroots ED /roots ED %nroots = roots == /nInters 0 def /Inters 4 array def /SaveInter { /tty ED /ttx ED nInters 4 lt { Inters nInters [ttx tty] put /nInters nInters 1 add def } if } def 0 1 nroots 1 sub { /idx ED /argt roots idx get def /x argt px0 add def /y argt argt mul 2 pp mul div py0 add def /ccxy x y cceval def % [idx 0 argt x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if } for Inters nInters SortInters% push the intersections on stack. end} def % ConicIParabolaInter % find all intersections of conic $cc(x,y)=ax^2+bxy+cy^2+dx+ey+f=0$ and % conjugate parabola $p(x,y): (y-y0)^2=2p(x-x0)$. % % a b c d e f x0 y0 p ConicIParabolaInter -> [Inters] nInters /ConicIParabolaInter { 15 dict begin % all variables are local /pp ED /py0 ED /px0 ED /ccf ED /cce ED /ccd ED /ccc ED /ccb ED /cca ED % [ cca ccb ccc ccd cce ccf ] == % x y cceval -> cc(x,y) /cceval { /ty ED /tx ED cca tx mul tx mul ccb tx mul ty mul add ccc ty mul ty mul add ccd tx mul add cce ty mul add ccf add } def /Qqa cca 4 pp pp mul mul div def /Qqb ccb 2 pp mul div def /Qqc ccc ccb py0 mul 2 pp mul div add cca px0 mul pp div add ccd 2 pp mul div add def /Qqd 2 ccc py0 mul mul ccb px0 mul add cce add def /Qqe px0 py0 cceval def Qqa Qqb Qqc Qqd Qqe QuarticRealRoots /nroots ED /roots ED % nroots = roots == /nInters 0 def /Inters 4 array def /SaveInter { /tty ED /ttx ED nInters 4 lt { Inters nInters [ttx tty] put /nInters nInters 1 add def } if } def 0 1 nroots 1 sub { /idx ED /argt roots idx get def /y argt py0 add def /x argt argt mul 2 pp mul div px0 add def /ccxy x y cceval def % [idx 0 argt x y ccxy ] == ccxy abs epsilonMax lt { x y SaveInter } if } for Inters nInters SortInters% push the intersections on stack. end} def % ConicInter % find all intersections of conic $cc1(x,y)=ax^2+bxy+cy^2+dx+ey+f=0$ and % $cc2(x,y)=a'x^2+b'xy+c'y^2+d'x+e'y+f'=0$. % % a b c d e f a' b' c' d' e' f' ConicInter -> [Inters] nInters /ConicInter { 15 dict begin % all variables are local /cxf ED /cxe ED /cxd ED /cxc ED /cxb ED /cxa ED /ccf ED /cce ED /ccd ED /ccc ED /ccb ED /cca ED % [ cxa cxb cxc cxd cxe cxf ] == % [ cca ccb ccc ccd cce ccf ] == % x y cceval -> cc(x,y) /cceval { /ty ED /tx ED cca tx mul tx mul ccb tx mul ty mul add ccc ty mul ty mul add ccd tx mul add cce ty mul add ccf add } def /cxeval { /ty ED /tx ED cxa tx mul tx mul cxb tx mul ty mul add cxc ty mul ty mul add cxd tx mul add cxe ty mul add cxf add } def /nInters 0 def /Inters 4 array def /SaveInter { /tty ED /ttx ED nInters 4 lt { Inters nInters [ttx tty] put /nInters nInters 1 add def } if } def cca cxa sub abs epsilon gt ccb cxb sub abs epsilon gt ccc cxc sub abs epsilon gt or or { /Qqa -2 cca ccc cxa cxc mul mul mul mul cca ccc cxb cxb mul mul mul add cca ccb cxc cxb mul mul mul sub ccb ccb cxa cxc mul mul mul add ccc ccb cxa cxb mul mul mul sub cca cca cxc cxc mul mul mul add ccc ccc cxa cxa mul mul mul add def /Qqb ccc ccb cxa cxe mul mul mul neg ccc cce cxa cxb mul mul mul sub ccb ccd cxc cxb mul mul mul sub 2 cca ccc cxb cxe mul mul mul mul add cca ccb cxc cxe mul mul mul sub ccb ccb cxc cxd mul mul mul add 2 ccc ccd cxa cxc mul mul mul mul sub 2 cca ccc cxc cxd mul mul mul mul sub ccc ccd cxb cxb mul mul mul add ccc ccb cxb cxd mul mul mul sub 2 ccb cce cxa cxc mul mul mul mul add 2 ccc ccc cxa cxd mul mul mul mul add cca cce cxc cxb mul mul mul sub 2 cca ccd cxc cxc mul mul mul mul add def /Qqc ccd cce cxc cxb mul mul mul neg ccc ccf cxb cxb mul mul mul add ccb ccf cxc cxb mul mul mul sub ccb ccd cxc cxe mul mul mul sub ccb ccb cxc cxf mul mul mul add 2 ccc ccd cxc cxd mul mul mul mul sub 2 ccc ccc cxa cxf mul mul mul mul add 2 cca ccf cxc cxc mul mul mul mul add ccd ccd cxc cxc mul mul mul add 2 ccc ccf cxa cxc mul mul mul mul sub 2 ccb cce cxc cxd mul mul mul mul add ccc cce cxb cxd mul mul mul sub 2 cca ccc cxc cxf mul mul mul mul sub ccc ccc cxd cxd mul mul mul add 2 ccc ccd cxb cxe mul mul mul mul add ccc cce cxa cxe mul mul mul sub cce cce cxa cxc mul mul mul add ccc ccb cxd cxe mul mul mul sub ccc ccb cxb cxf mul mul mul sub cca cce cxc cxe mul mul mul sub cca ccc cxe cxe mul mul mul add def /Qqd cce ccf cxc cxb mul mul mul neg ccc ccd cxe cxe mul mul mul add 2 ccd ccf cxc cxc mul mul mul mul add ccc ccb cxe cxf mul mul mul sub ccc cce cxd cxe mul mul mul sub 2 ccc ccc cxd cxf mul mul mul mul add 2 ccc ccd cxc cxf mul mul mul mul sub ccd cce cxc cxe mul mul mul sub 2 ccc ccf cxc cxd mul mul mul mul sub ccc cce cxb cxf mul mul mul sub 2 ccb cce cxc cxf mul mul mul mul add ccb ccf cxc cxe mul mul mul sub cce cce cxc cxd mul mul mul add 2 ccc ccf cxb cxe mul mul mul mul add def /Qqe -2 ccc ccf cxc cxf mul mul mul mul cce cce cxc cxf mul mul mul add ccf ccf cxc cxc mul mul mul add cce ccf cxc cxe mul mul mul sub ccc ccf cxe cxe mul mul mul add ccc ccc cxf cxf mul mul mul add ccc cce cxe cxf mul mul mul sub def % [ Qqa Qqb Qqc Qqd Qqe ] == Qqa Qqb Qqc Qqd Qqe QuarticRealRoots /nQRoots ED /QRoots ED } { cce cxe sub abs epsilon gt { /D cxe cce sub dup mul def /Qda cca cxe cxe mul mul ccb ccd mul ccb cxd mul sub 2 cca cce mul mul mul sub cxe mul add ccc cxd cxd mul mul add ccb cce mul 2 ccc ccd mul mul sub cxd mul add cca cce cce mul mul add ccb ccd cce mul mul sub ccc ccd ccd mul mul add D div def /Qdb ccb cxe mul 2 ccc cxd mul mul sub ccb cce mul sub 2 ccc ccd mul mul add cxf mul ccd cxe cxe mul mul sub cce cxd mul ccb ccf mul sub ccd cce mul add cxe mul add 2 ccc ccf mul mul cce cce mul sub cxd mul add ccb cce mul 2 ccc ccd mul mul sub ccf mul add D div neg def /Qdc ccc ccf cxf sub dup mul mul D div cce ccf cxf sub mul cxe cce sub div add ccf add def Qda Qdb Qdc QuadraticRealRoots /nQRoots ED /QRoots ED } { ccd cxd sub abs epsilon gt { /D cxd ccd sub def /Qda ccc def /Qdb ccb cxf mul neg cce cxd mul add ccb ccf mul add ccd cce mul sub D div def /Qdc cca ccf cxf sub dup mul mul D dup mul div ccd ccf cxf sub mul D div add ccf add def Qda Qdb Qdc QuadraticRealRoots /nDRoots ED /DRoots ED 0 1 nDRoots 1 sub { /DIdx ED /DVal DRoots DIdx get def /Qea cca def /Qeb ccb DVal mul ccd add def /Qec ccc DVal dup mul mul cce DVal mul add ccf add def Qea Qeb Qec QuadraticRealRoots /nERoots ED /ERoots ED 0 1 nERoots 1 sub { /EIdx ED /XVal ERoots EIdx get def /ccxy XVal DVal cceval def /cxxy XVal DVal cxeval def ccxy abs epsilonMax lt cxxy abs epsilonMax lt and { XVal DVal SaveInter } if } for } for Inters nInters SortInters% push the intersections on stack. exit } { ccf cxf sub abs epsilon lt { (cannot find intersection of identical conics) = exit } if } ifelse } ifelse } ifelse 0 1 nQRoots 1 sub { /QIdx ED /QVal QRoots QIdx get def /Qea ccc def /Qeb ccb QVal mul cce add def /Qec cca QVal dup mul mul ccd QVal mul add ccf add def Qea Qeb Qec QuadraticRealRoots /nERoots ED /ERoots ED 0 1 nERoots 1 sub { /EIdx ED /YVal ERoots EIdx get def /ccxy QVal YVal cceval def /cxxy QVal YVal cxeval def ccxy abs epsilonMax lt cxxy abs epsilonMax lt and { QVal YVal SaveInter } if } for } for Inters nInters SortInters% push the intersections on stack. end} def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end % END ps-euclide.pro