Maybe error when delta is close to 0?

[asyvn]

I write code to calculate radius of circle thought 2 point and tangents to a line.

`import geometry; unitsize(1cm); defaultpen(fontsize(9pt));

point B=(0.0, 0); label("$B$", B, SW); point A=(1.3, 4.2); label("$A$", A, NW); point C=(6.2, 0); label("$C$", C, SE);

point P=C, Q=A;

line lm=line(B, (2,0)); circle [] res; int [] arr= {-1,1}; real xp=P.x, yp=P.y, xq=Q.x, yq=Q.y; real a1=lm.a, b1=lm.b, c1=lm.c;

int k1=1; real A1=a1, B1=b1, C1=-c1, D1=-k1*sqrt(a1^2+b1^2); real A2=2(xp-xq), B2=2(yp-yq), C2=(xp^2-xq^2)+(yp^2-yq^2), D2=0; real PP, QQ, RR, SS, AA, BB, CC;

if (abs(B1A2-B2A1) > EPS) { PP = (C1B2 - C2B1)/(-B1A2+B2A1); QQ = (D1B2 - D2B1)/(-B1A2+B2A1); RR = (C1A2 - C2A1)/(B1A2-B2A1); SS = (D1A2 - D2A1)/(B1A2-B2A1); AA = QQ^2+SS^2-1; BB = 2(PP-xp)QQ + 2(RR-yp)SS; CC = (PP-xp)^2 + (RR-yp)^2;

write(BB^2-4AA*CC); real [] rx = quadraticroots(AA, BB, CC); //radius of PPL

for (int i=0; i<rx.length; ++i) { if (rx[i] >0) { point cx=(PP+QQrx[i], RR+SSrx[i]); res.push(circle(cx, rx[i])); } } } draw(res[0], blue); draw(A--B--C--cycle); dot(A^^B^^C, Fill(white));

shipout(bbox(1mm, invisible));`

I try with point A, point C and line lm. if P=C, Q=A, quadraticroots(AA, BB, CC) have delta = BB^2-4AACC=0 and has roots. But change P=A, Q=C, delta = -7.105427357601e-15 < 0 then quadratic equation has no real roots!!!

How to fix it?

[johncbowman]

Can you repost your code using code mode so we can actually run it? There are many missing *s in what you posted.

In any case, you might need to use a true Circle (rather than the Bezier approximation called circle).

[asyvn]

CODE HERE: https://gist.github.com/asyvn/c57d12e8047845f472394e528f9a53ad

[johncbowman]

That link doesn't actually point to the correct location. It is better to include code inline by applying the <> icon to the code and then checking with Preview:

import geometry;
unitsize(1cm);
defaultpen(fontsize(9pt));

point B=(0.0, 0); label("$B$", B, SW);
point A=(1.3, 4.2); label("$A$", A, NW);
point C=(6.2, 0); label("$C$", C, SE);

point P=C, Q=A;

line lm=line(B, (2,0));
circle [] res;
int [] arr= {-1,1};
real xp=P.x, yp=P.y, xq=Q.x, yq=Q.y;
real a1=lm.a, b1=lm.b, c1=lm.c;

int k1=1;
real A1=a1, B1=b1, C1=-c1, D1=-k1*sqrt(a1^2+b1^2);
real A2=2(xp-xq), B2=2(yp-yq), C2=(xp^2-xq^2)+(yp^2-yq^2), D2=0;
real PP, QQ, RR, SS, AA, BB, CC;

if (abs(B1*A2-B2*A1) > EPS) {
  PP = (C1*B2 - C2*B1)/(-B1*A2+B2*A1);
  QQ = (D1*B2 - D2*B1)/(-B1*A2+B2*A1);
  RR = (C1*A2 - C2*A1)/(B1*A2-B2*A1);
  SS = (D1*A2 - D2*A1)/(B1*A2-B2*A1);
  AA = QQ^2+SS^2-1; BB = 2(PP-xp)*QQ + 2(RR-yp)*SS; CC = (PP-xp)^2 + (RR-yp)^2;

  write("DELTA = ", BB^2-4AA*CC);
  real [] rx = quadraticroots(AA, BB, CC); //radius of PPL

  for (int i=0; i<rx.length; ++i) {
    if (rx[i] >0) {
      point cx=(PP+QQ*rx[i], RR+SS*rx[i]);
      res.push(circle(cx, rx[i]));
    }
  }
}

draw(res[0], blue);

draw(A--B--C--cycle);
dot(A^^B^^C, Fill(white));

shipout(bbox(1mm, invisible));

[johncbowman]

This is just a simple numerical precision issue.

I see now that you are using the circle type defined in the geometry module, not the built-in Bezier approximation given by circle(pair c, real r).

A simple workaround is to adjust BB a bit immediately after defining it. For example, try: BB += BB*realEpsilon;