How can I handle floating-point precision in C when dealing with numbers very close to zero or other edge cases

I wrote a program in C that calculates three points based on three points that are given. So we are given points A, B, C and calculate points A’, B’, C’. Now there are three ways that these points form a shape on a 2D plane (ABA’C, ABCB’ a AC’BC). Next, the program determines what type of shape was formed with these points either (square, rhombus, rectangle) or if none of those just output parallelogram. Here’s an example when running the program:

Point A:
[0,0]
Point B:
[-3,4]
Point C:
[4,3]
A': [1,7], square
B': [7,-1], parallelogram
C': [-7,1], parallelogram

Here’s an image for better understanding:

Now, the program works well for most cases. But given an edge case with numbers close to zero or big numbers. Such as:

Point A:
[78890.469e-63,0.669e-63]
Point B:
[2876.019e-63,78838.689e-63]
Point C:
[0.008e-63,0.669e-63]

My program outputs this:
Point A:
Point B:
Point C:
A': [-7.601444e-59,7.883869e-59], parallelogram
B': [7.601446e-59,-7.883735e-59], parallelogram
C': [8.176648e-59,7.883869e-59], square

//But it should output

Point A:
Point B:
Point C:
A': [-7.6014442e-59,7.8838689e-59], parallelogram
B': [7.6014458e-59,-7.8837351e-59], parallelogram
C': [8.176648e-59,7.8838689e-59], rhombus

Here is the whole program in C:

#include <stdio.h>
#include <math.h>
#include <float.h>

// Struct for Point
typedef struct {
    double x;
    double y;
} Point;

// Struct for Vector
typedef struct {
    double u1;
    double u2;
} Vector;

// Function to compute area
double area(Vector v) {
    return sqrt((v.u1 * v.u1) + (v.u2 * v.u2));
}

// Function to check for a right angle
double rightAngle(Vector u, Vector v) {
    if (fabs((u.u1 * v.u1) + (u.u2 * v.u2)) < 100 * DBL_EPSILON * fabs(u.u1 * v.u1) + 0.000000001) {
        return 1;
    }
    return 0;
}

// Function to read a coordinate
int readCoord(const char* prompt, Point* p) {
    char z1, z2, z3;
    printf("%s", prompt);
    if (scanf(" %c %lf %c %lf %c", &z1, &p->x, &z2, &p->y, &z3) != 5 || z1 != 91 || z2 != 44 || z3 != 93) {
        printf("Incorrect input.\n");
        return 0;
    }
    return 1;
}

// Function to check if three points are parallel
int areTheyParallel(Point a, Point b, Point c) {
    if ((a.x == b.x && a.y == b.y) || (b.x == c.x && b.y == c.y) || (a.x == c.x && a.y == c.y)) {
        printf("Cannot create shapes.\n");
        return 1;
    }

    double s1 = (b.y - a.y) / (b.x - a.x);
    double s2 = (c.y - b.y) / (c.x - b.x);

    if (fabs(s1 - s2) < 100 * DBL_EPSILON * s2) {
        printf("Cannot create shapes.\n");
        return 1;
    }

    return 0;
}

// Function to determine and print the type of quadrilateral
void determineType(const char* point, Vector u, Vector v, Point p) {
    double size_u = area(u);
    double size_v = area(v);
    double isRightAngle = rightAngle(u, v);

    if (fabs(size_u - size_v) <= 100 * DBL_EPSILON * size_v && isRightAngle == 1) {

        printf("%s: [%.13g,%.13g], square\n", point, p.x, p.y);

    } else if (fabs(size_u - size_v) <= 100 * DBL_EPSILON * size_v && isRightAngle == 0) {

        printf("%s: [%.13g,%.13g], rhombus\n", point, p.x, p.y);

    } else if (fabs(size_u - size_v) >= 100 * DBL_EPSILON * size_v && isRightAngle == 1) {

        printf("%s: [%.13g,%.13g], rectangle\n", point, p.x, p.y);

    } else if (fabs(size_u - size_v) >= 100 * DBL_EPSILON * size_v && isRightAngle == 0) {

        printf("%s: [%.13g,%.13g], parallelogram\n", point, p.x, p.y);
    }
}

Point computeNewPoint(Point base, Vector translation) {
    Point newPoint;
    newPoint.x = base.x + translation.u1;
    newPoint.y = base.y + translation.u2;
    return newPoint;
}

Vector computeVector(Point from, Point to) {
    Vector v;
    v.u1 = to.x - from.x;
    v.u2 = to.y - from.y;
    return v;
}

int main() {
    Point A, B, C, newA, newB, newC;
    Vector u, v;

    if (!readCoord("Point A:\n", &A)) return 1;
    if (!readCoord("Point B:\n", &B)) return 1;
    if (!readCoord("Point C:\n", &C)) return 1;

    if (areTheyParallel(A, B, C)) return 1;

    // A'
    u = computeVector(A, C);
    newA = computeNewPoint(B, u);

    u = computeVector(A, B);
    v = computeVector(A, C);
    determineType("A'", u, v, newA);

    // B'
    u = computeVector(B, A);
    newB = computeNewPoint(C, u);

    u = computeVector(C, B);
    v = computeVector(A, B);
    determineType("B'", u, v, newB);

    // C'
    u = computeVector(C, A);
    newC = computeNewPoint(B, u);

    u = computeVector(A, C);
    v = computeVector(B, C);
    determineType("C'", u, v, newC);

    return 0;
}

What should I try to ensure it will work with these extreme edge cases where there are very small numbers or very big numbers?

To handle the edge case with floating-point precision in my program, I tried the following:

I defined a small tolerance value (epsilon) to compare floating-point numbers instead of checking for exact equality. For instance, I set epsilon to 1e-10, allowing the program to account for minor inaccuracies in calculations. But still I got the wrong output.

What should I try to implement. Could you guide me in the right direction. Help will be much appreciated.

You need to sign in to view this answers