Как правильно расписать нахождение центра окружности на си в данной задаче?

Нахождение центра окружности и её радиуса по 3 точкам.

Даны 3 точки: A(x1,y1), B(x2,y2) и C(x3,y3);

Пусть точка D(x0,y0) - центр окружности, а точки M(xm,ym) и N(xn,yn) - середины хорд AB и BC.

1. В соответствие с теоремой о среднем перпендикуляре к хорде, линии DM и DN будут перпендикулярами.
2. Запишем уравннение прямой AB: (x-x1) / (x2-x1) = (y-y1) / (y2-y1)
3. Переходим от канонического уравнения прямой к уравнению с угловым коэффициентом (y = m_ab * x + b_ab): m_ab = (y2-y1) / (x2-x1); b_ab = y1 - x1(y2-y1)/(x2-x1)
4. Соответственно, угловой коэффициент перпендикуляра DM к прямой AB будет равен: _m_ab = -1/m_ab = (x1-x2) / (y2-y1)
5. Так же вычислим m_bc = (y3-y2) / (x3-x2) и _m_bc = (x2-x3) / (y3-y2)
6. Запишем m_dm = (ym-y0) / (xm-x0); m_dn = (yn-y0) / (xn-x0)
7. Составим систему { _m_ab = m_dm; _m_bc = m_dn }
8. Получим: { y0 = (x1-x2)/(y2-y1) * x0 + ym - xm * (x1-x2)/(y2-y1); y0 = (x2-x3)/(y3-y2) * x0 + yn - xn * (x2-x3)/(y3-y2) }
9. Пусть { y0 = am * x0 + dm; y0 = an * x0 + dn }
10. Получим: x0 = (dn - dm) / (am - an); А y0 = am * x0 + dm, соответсвенно.
11. Радиус равен дистанции между D и A; r = sqrt((x1-x0)^2 + (y1-y0)^2)

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

void swap(double* a, double* b) {
    double temp = *a;
    *a = *b;
    *b = temp;
}

bool has(double* beg, double* end, double x, double y) {
    for(double* it = beg; it < end; it += 2)
        if(*it == x && *(it + 1) == y)
            return true;
    return false;
}

double* scan(int* n) {
    double x, y;
    double *arr, *ptr;

    printf("Enter the number of points: ");
    scanf("%d", n);

    while(*n < 3) {
        printf("Wrong!!! Input must be greater or equal to 3!\nEnter the number of points: ");
        scanf("%d", n);
    }

    arr = (double*) malloc((*n << 1) * sizeof(double));
    ptr = arr;

    printf("Enter the coordinates of the points:\n");
    for(int i = 0; i < *n; ++i, ptr += 2) {
        // Считываем x и y
        printf("Point %d: ", i + 1);
        scanf("%lf %lf", &x, &y);

        // Проверка была ли точка уже введена
        if(has(arr, ptr + 1, x, y)) {
            printf("Point already exists! Input another one!\n");
            ptr -= 2;
            --i;
            continue;
        }

        // Записываем точку в массива
        *ptr = x;
        *(ptr + 1) = y;
    }
    *n = *n << 1;
    return arr;
}

// Проблема состоит в том, что иногда y2-y1 = 0 или y3-y2 = 0
// что приводит к x0 = NaN, y0 = NaN и r = NaN
// функция prep проходится по массиву и записывает в первые 6 ячеек
// массива значения которые не вызовут ошибок
void prep(double* arr, int n) {
    for(int i = 2; i < 6; i += 2) {
        if(arr[i - 1] != arr[i + 1]) { continue; }

        for(int j = 0; j < n; j += 2) {
            if(j == i || arr[i - 1] == arr[j + 1]) { continue; }
            
            (void) swap(&arr[i], &arr[j]);
            (void) swap(&arr[i + 1], &arr[j + 1]);
            i = 2;
        }
    }
}

// По теореме Пифагора расстояние между двумя точками:
// dist = sqrt((xn - x0)^2 + (yn - y0)^2)
double dist(double x0, double y0, double xn, double yn) {
    // Вычисляем (xn - x0)^2
    xn -= x0; xn *= xn;
    // Вычисляем (yn - y0)^2
    yn -= y0; yn *= yn;
    return sqrt(xn + yn);
}

void calc(double* arr, double* radius, double* x0, double* y0) {
    double xm = (arr[0] + arr[2]) / 2;
    double ym = (arr[1] + arr[3]) / 2;
    double am = (arr[0] - arr[2]) / (arr[3] - arr[1]);
    double dm = ym - xm * am;

    double xn = (arr[2] + arr[4]) / 2;
    double yn = (arr[3] + arr[5]) / 2;
    double an = (arr[2] - arr[4]) / (arr[5] - arr[3]);
    double dn = yn - xn * an;

    *x0 = (dn - dm) / (am - an);
    *y0 = am * (*x0) + dm;
    *radius = dist(*x0, *y0, arr[0], arr[1]);
}

int main() {
    int n;
    double *arr;
    double new_r, r, x0, y0;

    arr = scan(&n);
    (void) prep(arr, n);
    (void) calc(arr, &r, &x0, &y0);

    new_r = r;
    for(int i = 6; i < n; i += 2) {
        // Находим дистанцию от каждой новой точки до предполагаемого центра
        new_r = dist(x0, y0, arr[i], arr[i + 1]);

        // Если найденная дистанция не равна радиусу окружности
        // значит  данная точка не лежит на данной окружности
        // выходим из цикла
        if(new_r != r) { break; }
    }

    if(new_r == r) {
        printf("All points are on the same circle!\n");
        printf("Radius: %lf; Center: (%lf, %lf)\n", r, x0, y0);
    }
    else { printf("Not all points are on the same circle!\n"); }

    free(arr);
    return 0;
}

Решение

Вот так можно найти центр окружности по трём точкам:

...
    double a, b, c, d;
    buildCircle3(x1, y1, x2, y2, x3, y3, &a, &b, &c, &d);
    if (c != 0) {
        double r = sqrt(a * a + b * b - 4 * c * d) / (2 * c);
        double x = -a / (2 * c);
        double y = -b / (2 * c);
    }
...

void buildCircle3(
    double x1, double y1,
    double x2, double y2,
    double x3, double y3,
    double *a, double *b, double *c, double *d
) {
    const double r1 = x1 * x1 + y1 * y1;
    const double r2 = x2 * x2 + y2 * y2;
    const double r3 = x3 * x3 + y3 * y3;

    const double x2y3 = x2 * y3 - x3 * y2;
    const double x3y1 = x3 * y1 - x1 * y3;
    const double x1y2 = x1 * y2 - x2 * y1;

    *a = r1 * (y3 - y2) + r2 * (y1 - y3) + r3 * (y2 - y1);
    *b = r1 * (x2 - x3) + r2 * (x3 - x1) + r3 * (x1 - x2);

    *c =        x2y3 +      x3y1 +      x1y2 ;
    *d = -(r1 * x2y3 + r2 * x3y1 + r3 * x1y2);
}

Объяснение

Я буду изображать определители таблицами. Извините.

Утверждение: если точки (x₁, y₁), (x₂, y₂), (x₃, y₃) не лежат на одной прямой, то определитель ниже равен нулю тогда и только тогда, когда точка (x, y) лежит на окружности проходящей через (x₁, y₁), (x₂, y₂), (x₃, y₃).

Убедится в этом не сложно. Правый столбец из единиц можно использовать чтобы прибавлять одну и ту же константу ко всем иксам. Значение определителя не поменяется. В свою очередь столбец из иксов можно прибавлять к столбцу из квадратов. Следите за руками. Удвоенный столбец иксов добавляется к столбцу квадратов:

Cтолбец единиц добавляется к столбцу иксов и к столбцу квадратов:

Соберем квадраты в скобки:

Получилось, что все иксы в определителе увеличились на единицу, а его значение не изменилось. То есть, определитель сохраняет значение при трансляции по x на единицу. На самом деле определитель сохраняет значение при любой трансляции по x и по y.

Обозначим (x_c, y_c) - центр окружности проходящей через (x₁, y₁), (x₂, y₂), (x₃, y₃). Выполним трансляцию:

Нижние три значения в столбце квадратов - квадрат радиуса окружности. Обозначим r² = (x₁ - x_c)² + (y₁ - y_c)² = (x₂ - x_c)² + (y₂ - y_c)² = (x₃ - x_c)² + (y₃ - y_c)²:

С помощью столбца единиц вычтем r² из столбца квадратов:

Разложим определитель по верхней строке. Миноры при x - x_c , y - y_c и 1 равны нулю - в каждом есть столбец из нулей. Единственное не нулевое слагаемое это ((x - x_c)² + (y - y_c)² - r²) умноженный на минор:

Убираем из него трансляцию на (-x_c, -y_c):

Этот минор может быть равен нулю только если x₁, x₂, x₃ на одной прямой - это и есть уравнение прямой. Но у нас точки на окружности, не на прямой. Значит нулём должен стать другой множитель: ((x - x_c)² + (y - y_c)² - r²) = 0. А он может стать нулём только если точка (x, y) находится на расстоянии r от центра окружности (x_c, y_c). Что и требовалось доказать.

Теперь осталось понять что код в начале ответа разлагает исходный определитель по первой строке и получает уравнение ax + by + c(x² + y²) + d = 0, где a, b, c, d - соответствующие миноры.

У нас есть способ по трём точкам строить уравнение окружности. Если в этом уравнении c = 0, исходные точки лежат на одной прямой, решения нет. Если c ≠ 0, то из уравнения можно извлечь центр окружности и радиус:

r = √(a² + b² - 4cd) / (2c)

x_c = -a / (2c)

y_c = -b / (2c)

Прелесть этого решения в его численной устойчивости. До момента вычисления последних трёх выражений, использовались только сложение, вычитание и умножение. Проверка вырожденности тройки точек встроена в само уравнение.