Numerical rounding errors in intersection code

Question

I hope this question is in the right place, as it is as much about programming as it is about math. I'm trying to find the intersection between a circle and a line. My implementation of the algorithm works fine, until the line segment becomes very short. Take for example the following circle: center (-0.0177,0.0098) radius 0.00985 and line: position (-0.0275498730956292,0.00985) length 1.269043708E-7 direction (-1,0)

My algorithm (based on https://github.com/jpkinney/meshmorph/blob/master/doc/IntersectionLine2Circle2.pdf):

PGVector DELTA=PGVectorSubtract(line.position, arc.center);
CGFloat DDELTA=PGVectorDotProd(line.direction, DELTA);
CGFloat lengthDELTA2=PGVectorDotProd(DELTA, DELTA);
CGFloat delta=DDELTA*DDELTA-(lengthDELTA2-arc.radius*arc.radius);

if (CGFloatIsSmallerThanZero(delta)) return PGIntersectionNone;
else if (CGFloatIsBiggerThanZero(delta)) {
    CGFloat c,d;
    c=(-DDELTA-sqrt(delta))/line.length;
    d=(-DDELTA+sqrt(delta))/line.length;
    return PGIntersecionMultiple;
} else {
    *s1=-DDELTA/line.length;
    return PGIntersectionSingle;
}

It calculates that d=1,3669532933245911E-11. It should be zero... My float compare functions use 1000*DBL_EPSILON as its thresholds.

I've already solved some of these problems by storing the line direction separate from its length., so the delta calculation is more reliable. But it seems the problem has only shifted.

How do I solve problems like these?

Thanks in advance.

Kind regards,

Remco Poelstra