在处理非常接近零的数字或其他边缘情况时，如何处理 C 中的浮点精度

Question

我用 C 编写了一个程序，根据给定的三个点计算三个点。所以我们给定点 A、B、C 并计算点 A'、B'、C'。现在这些点可以通过三种方式在 2D 平面上形成形状（ABA'C、ABCB' a AC'BC）。接下来，程序确定这些点形成的形状类型（正方形、菱形、矩形），或者这些点是否只输出平行四边形。这是运行程序时的示例：

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

这里有一张图片可以更好地理解：

现在，该程序在大多数情况下都运行良好。但考虑到数字接近零或大数字的边缘情况。如：

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

这是 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;
}

我应该如何努力确保它能够处理这些数字非常小或非常大的极端边缘情况？

为了在程序中处理浮点精度的边缘情况，我尝试了以下操作：

我定义了一个小的容差值（epsilon）来比较浮点数，而不是检查是否完全相等。例如，我将 epsilon 设置为 1e-10，允许程序考虑计算中的微小误差。但我仍然得到了错误的输出。

我应该尝试实施什么？你能引导我走向正确的方向吗？非常感谢您的帮助。

Answer 1

这些点可以通过三种方式在 2D 平面上形成形状

我应该尝试实施什么。

代码可能会使用浮点数学和类型。 然而浮点具有对数分布，而 OP 的几何形状确实需要线性（2D）分布。

定点

考虑仔细使用定点数学（例如比例整数）来解决此任务。

调查有问题的 FP 点并形成量表。 然后舍入（通过

llround(x * scale)

）到大约 32 位值。 然后，各种

is_a_certina_gemotric_shape()

测试继续进行内部整数数学计算。

区域很有趣，因为它可能涉及

sqrt()

。 许多 area 测试可以替换为 area_sqaured 测试，从而不再需要

sqrt()

。