一、求方程的根

牛顿法

1、求数字的算术平方根

// 牛顿法求 num 的平方根
// 令f(x) = x^2 - num , 求 f(x) 的零点

// initial_num  (x0)     last_num: x(1)

#include<iostream>
using namespace std;
#include<cmath>

double getroot(const double & num){
    if(num==0||num==1)return num;

    double error = 1e-17;
    double initial_num = num/2.0;
    double last_num = (num/initial_num+initial_num)/2.0;

    while ( fabs(last_num - initial_num) > error)
    {
        initial_num = last_num;
        last_num = (num/initial_num+initial_num)/2.0;
    }
    
    return last_num;
}

int main(){
    double enter_num;
    cin>>enter_num;
    cout<< "the root is :"<<getroot(enter_num)<<endl;
    return 0;

}

2、求x 2 − ln ⁡ x = 0的根

原文链接：https://blog.csdn.net/weixin_42301220/article/details/126816206

#include<iostream>
using namespace std;

constexpr double esp = 1e-6;


int main() {
        double initial_num = 2;
        double last_num = initial_num - (initial_num * initial_num + log(initial_num)) / (2.0 * initial_num + (1 / initial_num));

        while (abs(initial_num-last_num)>=esp)
        {
                initial_num = last_num;
                last_num = initial_num - (initial_num * initial_num + log(initial_num)) / (2.0 * initial_num + (1 / initial_num));

        }
        cout << last_num << endl;
        return 0;

}

二、计算几何学

Geometry

Point：

class Point {
public:
        Point() = default;
        Point(double x_in, double y_in) :x(x_in), y(y_in) {};  //拷贝构造函数

        Point(const Point& p) :x(p.x), y(p.y) {};
        Point& operator=(const Point& p) {  //赋值运算符
                x = p.x;
                y = p.y;
                return *this;
        }

        Point operator+(const Point& p) {
                return { x + p.x, y + p.y };
        }
        Point operator-(const Point& p) const {
                return { x - p.x, y - p.y };
        }
        double operator*(const Point& p) const {
                return x * p.x + y * p.y;
        }
        Point operator*(double k)const {
                return { x * k, y * k };
        }
        friend Point operator*(double k, const Point& p) {
                return { p.x * k, p.y * k };
        }
        bool operator==(const Point& p)const {
                return p.x == x && p.y == y;
        }
        bool operator!=(const Point& p)const {
                return !(p.x == x && p.y == y);
        }

        double modulus()const {
                return sqrt(x * x + y * y);
        }

        double DistanceTo(const Point& other)const {
                double dx = x - other.x;
                double dy = y - other.y;
                return sqrt(dx * dx + dy * dy);
        }

        friend std::ostream& operator<<(std::ostream& out, const Point& p) {
                out << "(" << p.x << ", " << p.y << ")";
                return out;
        }

public:
        double x;
        double y;
};

直线line：

//Define line segment.
class Line {
public:
        Line() = default;
        Line(Point p1_in, Point p2_in) : p1(p1_in), p2(p2_in), direction(p2_in - p1_in) {
        };

        friend std::ostream& operator<<(std::ostream& out, const Line& line) {
                out << "Line: " << line.p1 << " ---> " << line.p2;
                return out;
        }

public:
        Point p1;
        Point p2;
        Point direction;
};

线段segment：

class Segment {
public:
        Segment() = default;
        Segment(Point start_in, Point end_in) : start(start_in), end(end_in), direction(end - start) {
        }

        Segment(const Segment& s) : start(s.start), end(s.end), direction(end - start) {}

        Segment& operator=(const Segment& s) {
                start = s.start;
                end = s.end;
                return *this;
        }

        //Segment operator+(const Segment& rhs) const {
        //        return { start + rhs.start, end + rhs.end };
        //}
        Segment operator-(const Segment& rhs)const {
                return { start - rhs.start, end - rhs.end };
        }

        double length()const {
                return direction.modulus();
        }

        Point unit_direction()const {
                int len = length();
                if (len != 0) {
                        return { len / direction.x , len / direction.y };
                }
                else throw std::runtime_error("Cannot calculate unit direction for a segment with zero length.");
        }
        friend std::ostream& operator<<(std::ostream& out, const Segment& s) {
                out << "Segment: " << s.start << " ---> " << s.end;
                return out;
        }

public:
        Point start;
        Point end;
        Point direction;
};

折线段polyline：

class Polyline {
        Polyline() = default;

        Polyline(const std::vector<Point>& pts) :points(pts) {
                for (int i = 1; i < points.size(); ++i) {
                        segs.push_back(Segment(points[i - 1], points[i]));
                }
        }

        Polyline(const std::vector<Segment>& segs_) : segs(segs_) {
                for (int i = 0; i < segs.size(); ++i) {
                        points.push_back(segs[i].start);
                }
                points.push_back(segs[segs.size() - 1].end);
        }

        void append(const Segment& seg) {
                if (!segs.empty() && segs.back().end != seg.start) {
                        throw std::invalid_argument("Disconnected Segment");
                }
                segs.push_back(seg);
                points.push_back(seg.end);
        }
        void append(const Point& p) {
                const auto seg = Segment(points.back(), p);
                points.push_back(p);
                segs.push_back(seg);
        }

        Polyline operator+(const Polyline& other) const {
                Polyline result;
                result.segs = this->segs;
                result.points = this->points;
                for (auto& seg : other.segs) {
                        result.append(seg);
                }

                return result;
        }

        Segment GetSegmentByIndex(int index) {
                if (index < 0 || index >= segs.size()) {
                        throw std::out_of_range("Index out of range");
                }
                return segs[index];
        }
        std::vector<Point> Points() const {
                return points;
        }
        std::vector<Segment> Segments()const {
                return segs;
        }

private:
        vector<Point> points;
        vector<Segment> segs;
};

Algorithm

基本运算

点积：

constexpr double Epsilon = 1e-6;


using Vec = Point;   //类型别名

// 点积运算
double DotProduct(const Vec& v1, const Vec& v2) {
        return v1.x * v2.x + v1.y + v2.y;
}

叉乘：

// 叉乘运算
double CrossProduct(const Vec& v1, const Vec& v2) {
        return v1.x * v2.y - v2.x * v1.y;
}

int Sign(double x) {
        return x < 0 ? -1 : x>0 ? 1 : 0;
}

Projection

点到线段的距离：

// 计算点到线段的投影长度 x   (p1p*cos(theta))
double ComputeProjectionLength(const Point& p, const Segment& seg) {
        const auto& p1p = p - seg.start;
        return DotProduct(p1p, seg.unit_direction());
}

点到线段的投影点：

// Compute projection point of point p.
Point ComputeProjection(const Point& p, const Segment& segment) {
    double projection_length = ComputeProjectionLength(p, segment);
    Point ans;
    ans.x = segment.start.x + segment.unit_direction().x * projection_length;
    ans.y = segment.start.y + segment.unit_direction().y * projection_length;
    return ans;
 
}

Distance-求距离

point to point：

// Get distance between point p1 and point p2.
double GetDistance(const Point& p1, const Point& p2) {
        return p1.DistanceTo(p2);
}

point to line：

// Get distance between point p and a straight line.
double GetDistance(const Point & p, const Line& line) {
        Segment p1p2(line.p1, line.p2);
        Segment p1p(line.p1, p);

        /*p1p_.x = p1p.end.x - p1p.start.x;
        p1p_.y = p1p.end.y - p1p.start.y;*/
        return abs(CrossProduct(p1p.direction, p1p2.unit_direction()));
}

point to segment：

// Get distance between point p and segment(p1,p2).
double GetDistance(const Point& p, const Segment& seg) {
        Segment p1p(seg.start, p);
        Segment p2p(seg.end, p);
        const auto c1 = DotProduct(p1p.direction,seg.direction);
        const auto c2 = DotProduct(p2p.direction,seg.direction);
        if (c1 <= 0)return GetDistance(p, seg.start);
        if (c2 >= 0)return GetDistance(p, seg.end);
        return abs(CrossProduct(p1p.direction, seg.unit_direction()));// 投影点在线段上
}

segment to segment：

四个端点到线段的最小值

double GetDistance(const Segment& seg1, const Segment& seg2) {
        const double d1 = GetDistance(seg1.start, seg2);
        const double d2 = GetDistance(seg1.end, seg2);
        const double d3 = GetDistance(seg2.start, seg1);
        const double d4 = GetDistance(seg1.end, seg1);
        return min(min(d1, d2), min(d3, d4));

}

Side-求相对位置关系

enum class Side {
    // When the segment's length is zero, it's useless to determine the side, so we use DEFAULT_NO_SIDE to show.
    DEFAULT_NO_SIDE = 0,
    LEFT,
    RIGHT,
    // The three below states mean that the point is in line.
    ON_AFTER,
    ON_BEFORE,
    WITHIN
};

也就是说点与线段的相对位置关系包括以下5种：

• 点在线段的左边
• 点在线段的右边
• 点在线段所在的直线上点在线段前面点在线段后面点在线段内部

判断点跟一条线段的相对位置关系：

//判断点和线段的位置关系（5种）
// Determine which side of the segment the point is.
Side GetSide(const Point& p, const Segment& s) {
    const auto& p0 = s.start;
    const auto& p2 = s.end;
    const auto& a = p - p0;
    const auto& b = p2 - p0;
    const auto cross_product = CrossProduct(a, b);
    if (cross_product != 0) {
        // Returns LEFT if p0,p,p2 are clockwise position, returns RIGHT means p0,p,p2 are counter-clockwise position.