牛顿迭代法、高精度牛顿迭代法开平方根（整数）

BB

此时，一个菜逼再次路过 QAQ

 

Introduction

牛顿迭代法 可以求在[a, b]上连续且单调的函数f(x)，f(x) = 0的近似解

 

The Roots - UVA 10428

Description
求解多项式f(x) = a0 + a1*x1 + a2*x2^2 + ... + an*xn^n = 0的所有根
Sample Input

2 1 -9.5750000000 -179.5585140000
4 1 53.3120000000 958.2510390000 6677.1763593480 15733.6254955064

Sample Output

Equation 1: -9.4420 19.0170
Equation 2: -22.8060 -18.1170 -6.7350 -5.6540

Solution 1 - 牛顿迭代法
从x=-inf开始使用牛顿迭代法求解该方程的根，假定求得的根为x0，这一根为最接近刚开始假定的初始解的根（不要问蒟蒻为什么，蒟蒻也不知道QAQ）。将多项式除以(x-x0)后，重复上述过程，直到多项式为1
这里除以(x-x0)是一定可以整除的，因为假定根为x0, x2, ..., xk，则原多项式实际上可以表述为(x-x0)(x-x1)...(x-xk)

// scanf printf push_back emplace_back memset sizeof
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const double eps = 1e-5;
const int N = (int)2e6 + 15;
const int inf = 0x3f3f3f3f;
const int MOD = 998244353;

inline int dcmp(double x) {
	return fabs(x) < eps ? 0 : (x > 0 ? 1 : -1);
}

inline double fun(const vector<double>& a, double x) {
	double res = 0;
	for(int i = 0; i < (int)a.size(); i++) {
		res += pow(x, i) * a[i];
	}
	return res;
}

inline double newton(const vector<double>& a, const vector<double>& aDiff) {
	double x = -30;
	for(int i = 0; i < 1000000; i++) {
		double nx = x - fun(a, x) / fun(aDiff, x);
		if(dcmp(x - nx) == 0) {
			return nx;
		}
		x = nx;
	}
	return x;
}

inline vector<double> getRoots(vector<double> a) {
	vector<double> res;
	while(a.size() > 1) {
		vector<double> aDiff, na;
		for(int i = 1; i < (int)a.size(); i++) {
			aDiff.push_back(i * a[i]);
		}
		double x0 = newton(a, aDiff);
		res.push_back(x0);
		for(int i = a.size() - 1; i >= 1; i--) {
			a[i - 1] = a[i - 1] + x0 * a[i];
		}

		a.erase(a.begin());
	}
	return res;
}


int main() {
	int n, csn = 1;
	while(~scanf("%d", &n) && n) {
		vector<double> a(n + 1);
		for(int i = n; i >= 0; i--) {
			scanf("%lf", &a[i]);
		}
		a = getRoots(a);
		printf("Equation %d:", csn++);
		for(double x : a) {
			printf(" %.4f", x);
		}
		puts("");
	}
}

Solution 2 - 递归 + 二分
求解该多项式的根，可以先确定极值所在的x坐标，再在相邻的x坐标之间二分找出零点（最左的x坐标还需要考虑左边，最右的同理）
问题转换为求解极值所在的x坐标，这一坐标实质上为导函数的零点，因此可以求出导函数后递归求解
最后考虑递归出口，当多项式为一次函数时，直接用除法就可以得到零点

// scanf printf push_back emplace_back memset sizeof
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const double eps = 1e-11;
const int N = (int)2e6 + 15;
const int inf = 0x3f3f3f3f;
const int MOD = 998244353;

inline int dcmp(double x) {
	return fabs(x) < eps ? 0 : (x > 0 ? 1 : -1);
}

inline double fun(const vector<double>& a, double x) {
	double res = 0;
	for(int i = 0; i < (int)a.size(); i++) {
		res += pow(x, i) * a[i];
	}
	return res;
}

inline double solve(const vector<double>& a, double l, double r) {
	double tmp = fun(a, l);
	if(dcmp(tmp) == 0) {
		return l;
	}
	while(r - l > eps) {
		double m = (l + r) / 2;
		double res = fun(a, m);
		if(dcmp(res) == 0) {
			return m;
		}
		if((dcmp(fun(a, m)) > 0) ^ (dcmp(tmp) < 0)) {
			l = m;
		} else {
			r = m;
		}
	}
	return l;
}

vector<double> getRoots(vector<double> a) {
	if(a.size() == 2) {
		return vector<double>(1, -a[0] / a[1]);
	}
	vector<double> b, aDiff, res;
	for(int i = 1; i < (int)a.size(); i++) {
		aDiff.push_back(i * a[i]);
	}
	b = getRoots(aDiff);
	b.insert(b.begin(), -inf);
	b.push_back(inf);
	for(int i = 1; i < (int)b.size(); i++) {
		double t1 = fun(a, b[i - 1]);
		double t2 = fun(a, b[i]);
		if((dcmp(t1) > 0 && dcmp(t2) > 0) || (dcmp(t1) < 0 && dcmp(t2) < 0)) {
			continue;
		}
		res.push_back(solve(a, b[i - 1], b[i] - eps));
	}
	return res;
}

int main() {
	int n, csn = 1;
	while(~scanf("%d", &n) && n) {
		vector<double> a(n + 1);
		for(int i = n; i >= 0; i--) {
			scanf("%lf", &a[i]);
		}
		a = getRoots(a);
		printf("Equation %d:", csn++);
		for(double x : a) {
			printf(" %.4f", x);
		}
		puts("");
	}
}

 

Square root - UVA 10023

Description
对x^2 <= N，给定整数N，求整数X，需要高精度
Sample Input

1
7206604678144

Sample Output

2684512

Solution
（套板 （逃
（不要问蒟蒻为什么，蒟蒻也不知道QAQ

import java.math.BigInteger;
import java.util.Scanner;

class Main {
    public static BigInteger newton(BigInteger n) {
        BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2);
        boolean isDec = false;
        while(true) {
            BigInteger b = n.divide(a).add(a).shiftRight(1);
            if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && isDec)
                break;
            isDec = a.compareTo(b) > 0;
            a = b;
        }
        return a;
    }

    public static void main(String[] args) {
        Scanner cin = new Scanner(System.in);
        int t = cin.nextInt();
        for(int i = 0; i < t; i++) {
            BigInteger a = cin.nextBigInteger();
            if(i != 0) {
                System.out.println("");
            }
            System.out.println(newton(a));
        }
    }
}