[uncompleted] preliminary study on numerical calculation in OI

Time:2022-5-29

catalogue

  • newton iteration
  • Lagrange multiplier method
  • Gauss elimination method
  • Lagrange interpolation
  • Adaptive Simpson method
  • Monte Carlo method
  • More randomization algorithms
  • More

preface

Don’t pollute oi with mathematics Jpg

There should be no such cancer problem… Probably?

Update: less than half of the time

The more I write, the more I find that there are many things about numerical calculation… At first, I only planned to write Newton iteration and Lagrange multiplier.Startling giant pit

newton iteration

Theoretical part

Newton iteration is usually used to solve the approximate root of a function. Abel rufini theorem points out that there is no general algebraic solution for the quintic and higher-order algebraic equations, that is, there is no form of exact solution such as root formula, and only the approximate solution can be considered.

Taylor expansion

Gradually deviate from the topic (?)

If you don’t want to see this part, skip to the fourth level title, “what if I don’t want to see Taylor unfold?” Division.

(Le Tai Le Taylor expansion.jpg) (there should be a Kennedy expression pack here, but I didn’t do it)

Let’s imitate a function, such as\(g(x) = e^x\). The premise of imitation is to describe this function. We can describe the function by derivative.

However, only one first derivative is very weak. It is easy to be attacked. There is little information, so consider finding out\(n\)The derivative of order if there is one\(f(x)\)And\(g(x)\)First order to\(n\)If the order derivatives are all the same, it can be considered that the two functions are similar enough.

And then we found out\(f(x) = e^x\)This function can be continuously and continuously (unambiguously).

Then we can find the multiple derivatives by imitating a polynomial that can be derived all the time.

\[\large

f(x) = \lim_{n \rightarrow \infty} \sum_{i = 0}^{n}\frac{f^{(n)}(x_0)}{n!} (x – x_0)^i + R_n(x)

\]

But obviously it is impossible to calculate\(\infty\)Yes, but the later payano remainder will be smaller and smaller to an acceptable error range, which can be regarded as a solution.

So we just need the first item of Taylor’s expansion:

\[\large\begin{aligned}

f(x) = f(x_0) + f'(x_0)(x – x_0)

\end{aligned}\]

Then we ask for the zero point of the function:

\[\large\begin{aligned}

&f(x) = f(x_0) + f'(x_0)(x – x_0) = 0\\
&f'(x_0)(x – x_0) = -f(x_0)\\
&x = x_0 – \frac{f(x_0)}{f'(x_0)}

\end{aligned}\]

What if the first term of Taylor expansion is not accurate enough?

Iteration, each time the last solution is substituted, the zero point can be approached.

Thus, the iterative formula is obtained:

\[\large

x_{n + 1} = x_n – \dfrac{f(x_n)}{f'(x_n)}

\]

What if I don’t want to see Taylor unfold?

Let’s look at the function image to understand:

For a function\(f(x)\), we find the derivative\(f'(x)\)

Generate an initial solution\(x_0\), find one pass at a time\((x_0,f(x_0))\), function\(y = f(x)\)The slope of this tangent is\(f'(x_0)\), and\(x\)Axis intersection is\(x_0 – \dfrac{f(x_0)}{f'(x_0)}\), you can find that repeating this process is approaching the function and\(x\)Intersection of axes.

So Newton’s iteration has the following formula:

\[\large

x_{n + 1} = x_n – \dfrac{f(x_n)}{f'(x_n)}

\]

When Newton iteration is not possible

  • The initial iteration value is selected as the stagnation point, so that\(f'(x) = 0\)
  • The function does not converge and runs farther and farther

Newton iteration of multivariate function

Change the derivative togradient

Abrupt linear algebra

Newton iteration process

First we need\(f(x)\)and\(f'(x)\)Then select an initial solution and start the iteration.

Newton iterative method has the property of square convergence, that is to say, it is\(\mathcal{O} (\log n)\)of

Example: calculate\(\sqrt{x}\), first\(\sqrt{a}\)Can be understood as equation\(x^2 – a = 0\)For the original function\(f(x) = x^2 – a\)Derived\(f'(x) = 2x\), and then directly iterate.

Reference code:

#include 
#include 
#include 
using std::abs;
constexpr double eps = 1e-9;

inline double NewtonSqrt(double x) {
	double x0 = x;
	while(abs(x0 * x0 - x) > eps)
		x0 -= (x0 * x0 - x) / (2 * x0);
	return x0;
}

int main() {
	double x;
	scanf("%lf",&x);
	printf("Newton's method %.8lf\n",NewtonSqrt(x));
	printf("std::sqrt %.8lf\n",std::sqrt(x));
	return 0;
}

It can be found that the initial value of iteration is not very important, but if there is a good initial solution, the number of iterations will be significantly reduced. The famous fast inverse square root algorithm uses “magic numbers”
0x5f3759dfThen only one iteration is performed and the speed is greatly accelerated.

In Chris lomont’s paper, an excellent initial value with slightly better efficiency than the quake III source code is given:0x5f375a86

Paper link:http://www.lomont.org/papers/2003/InvSqrt.pdf

#include 
#include 
#include 
using std::abs;
constexpr double eps = 1e-9;

inline float NewtonSqrt(float x) {
   float halfx = 0.5 * x;
   int i = *((int *)&x);
   i = 0x5f3759df - (i>>1);
   x = *((float *)&i);
   x = x * (1.5 - (halfx * x * x));
   return 1.0 / x;
} 

int main() {
   double x;
   scanf("%lf",&x);
   printf("Fast Sqrt %.8lf\n",NewtonSqrt(x));
   printf("std::sqrt %.8lf\n",std::sqrt(x));
   return 0;
}

Of course, this method also hasdoubleVersion, the initial solution used is0x5fe6eb50c7aa19f9

#include 
#include 
#include 
using std::abs;
constexpr double eps = 1e-9;

inline double NewtonSqrt(double x) {
    double halfx = 0.5 * x;
    long long i = *((long long *)&x);
    i = 0x5fe6eb50c7aa19f9 - (i >> 1);
    x = *((double *)&i);
    x = x * (1.5 - (halfx * x * x));
    return 1.0 / x;
} 

int main() {
	double x;
	scanf("%lf",&x);
	printf("Fast Sqrt %.8lf\n",NewtonSqrt(x));
	printf("std::sqrt %.8lf\n",std::sqrt(x));
	return 0;
}

However, the square root reciprocal efficiency has beenrsqrtssHang itThe victory of hardware is the victory of von Neumann!

Example

Equilibrium point

\(\texttt{Link.}\)

Find the weighted Fermat point on the plane.

This question seems to be calledFermat Weber problem(?), I’m not sure. I dare not talk nonsense.

First we firstWithout rightsTake a look, that is, the ordinary Fermat problem:

about\(n = 3\)That is, triangular Fermat points. Mathematician and physicist Torricelli gave an excellent solution to plane geometry:

  • When all three angles of a triangle are less than\(120°\)The point is the equilateral center of the triangle.
  • When a triangle has an internal angle greater than or equal to\(120°\)The point is the vertex of the maximum interior angle of the triangle.

But for\(n\)For larger problems, there is no better solution.

Now, what do we need to consider in detail:

about\(n\)individual\((x_i,y_i)\), find one\((x,y)\)bring

\[\large
\sum_{i = 1}^{n} \sqrt{(x – x_i)^2 + (y – y_i)^2}
\]

The value of is minimized.

(the weighted Fermat point is\(\sum_{i = 1}^{n} w_i\sqrt{(x – x_i)^2 + (y – y_i)^2}\)

So now the problem is to find the extreme value of a binary function:

\[\large
\min f(x,y) = \sum_{i = 1}^{n} \sqrt{(x – x_i)^2 + (y – y_i)^2}
\]

First of all, in order to get the change rate and concavity and convexity of this function, we need to find the derivative of it. But for a binary function, we need to fix one independent variable unchanged, and then find the relationship between the change of another independent variable and the dependent variable, which is equivalent to finding the derivative of two unary functions. At this time, the two derivatives are the derivatives of this binary functionpartial derivative

The partial derivative is obtained:

\[\large\begin{aligned}

&f_x(x,y) = \sum_{i = 1}^{n} \frac{x – x_i}{\sqrt{(x – x_i)^2 + (y – y_i)^2}} \\

&f_y(x,y) = \sum_{i = 1}^{n} \frac{y – y_i}{\sqrt{(x – x_i)^2 + (y – y_i)^2}}

\end{aligned}\]

Knowable when\(f_x(x,y) = 0\)And\(f_y(x,y) = 0\)And get the stagnation point of this function.

Then we find the second derivative:

\[\large\begin{aligned}

&\frac{\partial^2 z}{\partial x^2} = \sum_{i = 1}^{n} \frac{(y – y_i)^2}{\left(\sqrt{(x – x_i)^2 + (y – y_i)^2} \right)^3}\\

&\frac{\partial^2 z}{\partial y^2} = \sum_{i = 1}^{n} \frac{(x – x_i)^2}{\left(\sqrt{(x – x_i)^2 + (y – y_i)^2} \right)^3}

\end{aligned}\]

It can be seen that the second partial derivative is constant greater than\(0\)So this function is convex.

(so we can divide the three points directly\(x,y\)I have solved this problem, which is very reasonable.)

The original function is a convex function, so the first partial derivative is\(0\)The point of is its minimum value point. Consider how to solve the following equation:

\[\large\begin{aligned}

&\sum_{i = 1}^{n} \frac{x – x_i}{\sqrt{(x – x_i)^2 + (y – y_i)^2}} = 0 \\

&\sum_{i = 1}^{n} \frac{y – y_i}{\sqrt{(x – x_i)^2 + (y – y_i)^2}} = 0

\end{aligned}\]

Solve using Newton iteration:

express\(x\)And\(y\)As follows:

\[\Large\begin{aligned}

x = \frac{\sum_{i = 1}^{n} \frac{x_i}{\sqrt{(x – x_i)^2 + (y – y_i)^2}} } {\sum_{i = 1}^{n} \frac{1}{\sqrt{(x – x_i)^2 + (y – y_i)^2}}}\\

y = \frac{\sum_{i = 1}^{n} \frac{x_i}{\sqrt{(x – x_i)^2 + (y – y_i)^2}} } {\sum_{i = 1}^{n} \frac{1}{\sqrt{(x – x_i)^2 + (y – y_i)^2}}}

\end{aligned}\]

Then iterate vigorously and multiply it with weight.

How to solve faster? Set the initial solution to the centroid:

\[\large\begin{aligned}

x_0 = \frac{\sum_{i = 1}^{n}w_ix_i}{n}\\

y_0 = \frac{\sum_{i = 1}^{n}w_iy_i}{n}

\end{aligned}\]

Then remember to judge the case with only one point.

Run fast and get to the first page of the optimal solution.

// Fermat Point

#define sq(x) ((x) * (x))
std::pair p[N];
db w[N];

int main() {
    InitIO();
	int n = read();
	db x0 = 0,y0 = 0;
	rep(i,1,n) {
		p[i].first = readDB();
		p[i].second = readDB();
		w[i] = read();
		x0 += p[i].first * w[i];
		y0 += p[i].second * w[i];
	}
	x0 /= n,y0 /= n;
	if(n != 1) {
    	while(1) {
    		db x1 = 0,x2 = 0,y1 = 0,y2 = 0;
    		rep(i,1,n) {
    			db g = std::sqrt(sq(x0 - p[i].first) + sq(y0 - p[i].second)) / w[i];
    			// You can use your Newton's Method Square Root instead OuO
    			x1 += p[i].first / g;
    			x2 += 1 / g;
    			y1 += p[i].second / g;
    			y2 += 1 / g;
    		}
    		x1 /= x2,y1 /= y2;
    		if(std::abs(x1 - x0) < eps && std::abs(y1 - y0) < eps)
    			break;
    		else x0 = x1,y0 = y1;
    	}
    	writeDB(x0,3),space;
		writeDB(y0,3),enter;
	}
	else {
		writeDB(p[1].first,3),space;
		writeDB(p[1].second,3),enter;
	}
    EndIO();
    return 0;
}

Then the problem can be solved by other numerical methods, such as@damocrisBroydent iteration written by the boss:https://www.luogu.com.cn/record/67174683

Solution of one variable cubic equation

\(\texttt{Link.}\)

(it is true that this question is not divided into three points.)

Polynomial derivation:

\(\dfrac{x}{\mathrm{d}x} x^k = kx^{k – 1}\)

For this problem, the object function is\(f(x) = ax^3 +bx^2 +cx + d\)

Derivative is\(f'(x) = 3ax^2 + 2bx + c\)

#include 
#include 
#include 
#include 
#include 
#include 

typedef double db;
constexpr db eps = 1e-3;

db a,b,c,d;
inline db f(db x) {
	return a * x * x * x + b * x * x + c * x + d;
}

inline db fder(db x) {
	return 3.0 * a * x * x + 2.0 * b * x + c;
}

struct Node {
	db val;
	Node(){}
	Node(db v) :val(v) {}
	inline bool operator == (const Node &oth) const {
		return std::abs(val - oth.val) < eps;
	}
	inline bool operator < (const Node &oth) const {
		return -(val - oth.val) > eps;
	}
};
std::set res;

int main() {
	scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
	for(int i = -100;i <= 100;++i) {
		db x0 = i;
		while(std::abs(f(x0)) > eps)
			x0 -= f(x0) / fder(x0);
		res.emplace(x0);
	}
	for(auto i : res)
		printf("%.2lf ",i.val);
	return 0;
}

It is even possible to solve a simple and man high-order equation:

#include 
#include 
#include 
#include 
#include 
#include 
#include 

typedef double db;
constexpr db eps = 1e-3;
constexpr int N = 105;

int n;
db a[N];
db d[N];
inline db f(db x) {
	db res = 0,rx = 1;
	for(int i = 0;i <= n;++i) {
		res += rx * a[i];
		rx *= x;
	}
	return res;
}

inline db fder(db x) {
	db res = 0,rx = 1;
	for(int i = 0;i < n;++i) {
		res += rx * a[i];
		rx *= x;
	}
	return res;
}

struct Node {
	db val;
	Node(){}
	Node(db v) :val(v) {}
	inline bool operator < (const Node &oth) const {
		return -(val - oth.val) > eps;
	}
};
std::set res;

int main() {
	scanf("%d",&n);
	for(int i = n;~i;--i)
		scanf("%lf",&a[i]);
	for(int i = n;i;--i)
		d[i - 1] = a[i] * i;
	for(int i = -100;i <= 100;++i) {
		db x0 = i;
		while(std::abs(f(x0)) > eps)
			x0 -= f(x0) / fder(x0);
		res.emplace(x0);
	}
	for(auto i : res) if(std::isnormal(i.val))
		printf("%.2lf ",i.val);
	return 0;
}

Polynomial Newton iteration

No polynomial.

If you can’t, you have to write too. People are addicted to vegetables.

Pre knowledge\(n \log n\)Polynomial multiplication.

For a\(n\)Degree polynomial\(F(x)\)satisfy\(G(F(x)) \equiv 0 \bmod{x^n}\)

Polynomial inversion

Polynomial logarithmic function

Polynomial exponential function

Polynomial root

Polynomial trigonometric function

Newton iteration summary (?)

It can solve some function zero point problems or transform the function extreme value problem into derivative zero point to find its stationary point. It is also widely used in polynomials.

Lagrange multiplier method

There are so few Lagrange multiplier problems… They are basically found by searching other people’s blogs on the Internet.

Lagrange multiplier methodUnder constraintsThe method of finding the extreme value of multivariate function.

Theoretical part

Desirability function\(f(x,y)\)In constraints\(\phi (x,y) = 0\)First, the auxiliary function is constructed\(F(x,y) = f(x,y) + \lambda \phi(x,y)\)The stagnation point is obtained when the partial derivative is zero, and the simultaneous equations are as follows:

\[\large\left\{\begin{matrix}
&f_x(x,y) + \lambda\phi_x(x,y) = 0\\
&f_y(x,y) + \lambda\phi_y(x,y) = 0\\
&\phi (x,y) = 0
\end{matrix}\right.\]

Is it difficult to understand? Yes.

Curve in the following figure\(L\)Is a constraint\(\phi(x,y) = 0\), others are\(f(x,y) = C\)Isoline of.

Example

yja

Source unknown.

On the two-dimensional plane\(n\)We know the distance between each point and the origin, but we don’t know the specific location. Determine the location of these points, and find the maximum convex hull area that can be formed.

river

Source unknown.

CF 185 B

\(\texttt{Link.}\)

Noi2012 cycling in Sichuan and Tibet

\(\texttt{Link.}\)

Noi is afraid of this cancer

The degree of cancer is comparable to the beauty of circulation. Swimming pool, lemon tree under the moon, and big calculation in the wilderness

Gauss elimination method

Lagrange interpolation

This should also be a numerical calculation (?)

Adaptive Simpson method

Monte Carlo method

Pollard Rho discrete logarithm is also a kind of Monte Carlo. Its advantage is that the space overhead is less than bsgs and faster.

More randomization algorithms

Noip annealing saved my life. I had no choice but to write a blog.

The advantage of randomization algorithm is that it can be used as long as the data range is small enough or you find some properties. As long as you dare to write, you will have a chance to get points.

More

How does a computer find a logarithmic function\(\log\)of

First of all, we have the formula for changing the bottom\(\dfrac{\log_a b}{\log_a c} = \log_c b\), then as long as the logarithm of a base can be found, there will be a general solution for all bases.

Recommended Today

Cloud Native Virtualization: Building Edge Computing Instances with Kubevirt

With the popularity of Kubernetes, more and more businesses are running on containers, but there are still some business forms that are more suitable for running on virtual machines. How to control virtual machines and containers at the same time has gradually become a mainstream demand in the cloud-native era. Kubevirt gave came up with […]