C/C++求解两圆交点源代码

前段时间用到了圆定位算法，通过求其交点算出定位坐标，网上找了许久也没找到一个心仪的算法，索性自己写，下面是源代码，亲测无bug,喜欢的小伙伴可以参考一下。需要源文件的可以q我 442499086
#include "iostream"

#define square(val) ((val) * (val))
#define equal(c1, c2) (((c1.axis.x) == (c2.axis.x)) && ((c1.axis.y) == (c2.axis.y)) && ((c1.r) == (c2.r)))
#define yequal(c1, c2) ((c1.axis.y) == (c2.axis.y))
#define xequal(c1, c2) ((c1.axis.x) == (c2.axis.x))

typedef struct
{
    double x;
    double y;
}CircleAxis;

typedef struct
{
    CircleAxis axis;
    double r;
}Circle;

typedef struct
{
    CircleAxis axis[2];
}*PCircularIntersection, CircularIntersection;

int GetIntersectionPointCoordinate(PCircularIntersection ci, Circle c1, Circle c2);

int main()
{
    using namespace std;
    CircularIntersection ci;
    Circle c1, c2;
    c1.r = 4;
    c1.axis.x = 0;
    c1.axis.y = 0;

    c2.r = 4;
    c2.axis.x = 2;
    c2.axis.y = 0;

    GetIntersectionPointCoordinate(&ci, c1, c2);
    
    cout << "IntersectionPointCoordinate: " << endl;
    cout << "x1:" <<  ci.axis[0].x << " " << "y1:" << ci.axis[0].y << endl;
    cout << "x2:" <<  ci.axis[1].x << " " << "y2:" << ci.axis[1].y << endl;
    return 0;
}

int GetIntersectionPointCoordinate(PCircularIntersection ci, Circle c1, Circle c2)
{
    double error, A, B, C, K, S, x1, y1, x2, y2, x, y;
    if (equal(c1, c2) || ci == NULL)
        return 0;
    else
    {
        if (yequal(c1, c2))
        {
            x = (square((c1.axis.x)) + square((c1.axis.y)) + square((c2.r)) - square((c2.axis.x)) - square((c2.axis.y)) - square((c1.r))) / (2 * (c1.axis.x) - 2 * (c2.axis.x));
            A = 1;
            B = (-2) *  (c1.axis.y);
            C = square(x) + square((c1.axis.x)) + square((c1.axis.y)) - square((c1.r)) - 2 * (c1.axis.x)*x;
            error = square(B) - 4 * A*C;

            if (error < 0)
                return -1;
            else
            {
                x1 = x2 = x;
                y1 = ((-B) + sqrt(error)) / (2 * A);
                y2 = ((-B) - sqrt(error)) / (2 * A);
            }
        }
        else
            if (xequal(c1, c2))
            {
                y = (square((c1.axis.x)) + square((c1.axis.y)) + square((c2.r)) - square((c2.axis.x)) - square((c2.axis.y)) - square((c1.r))) / (2 * (c1.axis.y) - 2 * (c2.axis.y));
                A = 1;
                B = (-2) * (c1.axis.x);
                C = square((c1.axis.x)) + square((c1.axis.y)) + square(y) - square((c1.r)) - 2 * y* (c1.axis.y);
                error = square(B) - 4 * A*C;

                if (error < 0)
                    return -1;
                else
                {
                    y1 = y2 = y;
                    x1 = ((-B) + sqrt(error)) / (2 * A);
                    x2 = ((-B) - sqrt(error)) / (2 * A);
                }
            }
            else
            {
                K = (2 * (c2.axis.x) - 2 * (c1.axis.x)) / (2 * (c1.axis.y) - 2 * (c2.axis.y));
                S = (square((c1.axis.x)) + square((c1.axis.y)) + square((c2.r)) - square((c2.axis.x)) - square((c2.axis.y)) - square((c1.r))) / (2 * (c1.axis.y) - 2 * (c2.axis.y));

                //解二次方程
                A = (1 + square(K));
                B = 2 * K*S - 2 * (c1.axis.x) - 2 * (c1.axis.y)*K;
                C = square(S) + square((c1.axis.x)) + square((c1.axis.y)) - square((c1.r)) - 2 * (c1.axis.y)*S;
                error = square(B) - 4 * A*C;

                if (error < 0)
                    return -1;
                else
                {
                    x1 = ((-B) + sqrt(error)) / (2 * A);
                    x2 = ((-B) - sqrt(error)) / (2 * A);
                    y1 = K * x1 + S;
                    y2 = K * x2 + S;
                }
            }
        ci->axis[0].x = x1;
        ci->axis[0].y = y1;

        ci->axis[1].x = x2;
        ci->axis[1].y = y2;
        return 1;
    }
}