HDU-ITMO – 第 3 页 – Pancake's Personal Website

Timus #1650. Billionaires

题目：有n个富豪（数据包含：姓名，开始时的所在地，拥有财富），某地财富即为在此的所有富豪财富相加。在接下来的m天，会有k个土豪在某一天晚上离开所在地，在第二天早上到达所在地。求每天拥有财富最多的城市，按照字母表顺序输出某城市财富最多的天数。ps：若某一天有两个及以上的城市财富相同则该天不算在这几个城市的天数中，若最终有城市天数为0则不输出此城市。

解析：该题目逻辑很简单，建立Person和City结构体，Person储存姓名，所在地，财富；City储存城市名，拥有财富。循环查找/插入/删除/更新City和Person中的数据，并找到City中财富最大的城市。

代码分析：你必须建立map或set容器或写自己的Tree结构来完成此题，否则必会超时。比如，如果你用顺序的方式遍历City中最大的城市/城市名必定会超时。
1. 我的解决方法是用map<string, City*>cityptr;用指针快速修改和查找。
2. 更加快捷的方式是map<string, pair<City*, long long>>Person;
3. 傻一点的方法是建立map<string, int>index;这种方法缺点是每次都要更新index。

#include<iostream>
#include<map>
#include<set>
using namespace std;

struct City {
	string name;
	long long money;
	City(string nameIN, long long moneyIN):name(nameIN),money(moneyIN) {}
	bool operator< (const City& a)const
	{
		if (this->money == a.money) return this->name < a.name;
		return this->money > a.money;
	}
};

void add(set<City>& city, map<string, City*>& cityptr, string cityname, long long money) {
	map<string, City*>::iterator ptr;
	ptr = cityptr.find(cityname);
	City * cityNow = new City(cityname, money);
	if(ptr == cityptr.end()){
		city.insert(*cityNow);
		cityptr.insert({ cityNow->name, cityNow });
	}
	else {
		cityNow->money += ptr->second->money;
		city.erase(*ptr->second);
		ptr->second = cityNow;
		city.insert(*cityNow);
	}

}

void insertResult(set<City>& city, map<string, int> &result, int day) {
	string name;
	if (city.size() >= 2) {
		set<City>::iterator ptr = city.begin();
		if (ptr->money - (++ptr)->money != 0) name = city.begin()->name;
		else return;
	}
	else if (city.size() == 0) return;
	else if (city.size() == 1) name = city.begin()->name;
	map<string, int>::iterator ptr;
	ptr = result.find(name);
	if (ptr != result.end()) ptr->second += day;
	else result.insert({ name, day });
}


void del(set<City>& city, map<string, City*>& cityptr, string cityName, long long money) {
	map<string, City*>::iterator ptr;
	ptr = cityptr.find(cityName);
	City *tempCity = new City(cityName, ptr->second->money - money);
	city.erase(*ptr->second);
	delete ptr->second;
	ptr->second = tempCity;
	if (tempCity->money == 0) {
		cityptr.erase(ptr);
		return;
	}
	city.insert(*tempCity);
}

int main() {
	map<string, pair<string, long long>>Person;
	map<string, int>result;
	map<string, City*>cityptr;
	set<City>city;
	int n, m, k, day = 0, preday;
	cin >> n;
	string namePerson;
	string nameCity;
	string namePreCity;
	long long money;
	for (int i = 0; i < n; i++) {
		cin >> namePerson >> nameCity >> money;
		Person.insert({ namePerson,make_pair(nameCity,money) });
		add(city, cityptr, nameCity, money);
	}
	cin >> m >> k;
	for (int i = 0; i < k; i++) {
		preday = day;
		cin >> day >> namePerson >> nameCity;
		insertResult(city, result, day - preday);
		map<string, pair<string, long long>>::iterator ptr;
		ptr = Person.find(namePerson);
		money = ptr->second.second;
		add(city, cityptr, nameCity, money);
		namePreCity = ptr->second.first;
		ptr->second.first = nameCity;
		del(city, cityptr, namePreCity, money);
	}
	insertResult(city, result, m - day);
	for (auto i : result) {
		if (i.second == 0)continue;
		cout << i.first << " " << i.second << endl;
	}
}

#include<iostream>

#include<map>

#include<set>

using namespace std;

struct City {

string name;

long long money;

City(string nameIN, long long moneyIN):name(nameIN),money(moneyIN) {}

bool operator< (const City& a)const

{

if (this->money == a.money) return this->name < a.name;

return this->money > a.money;

}

};

void add(set<City>& city, map<string, City*>& cityptr, string cityname, long long money) {

map<string, City*>::iterator ptr;

ptr = cityptr.find(cityname);

City * cityNow = new City(cityname, money);

if(ptr == cityptr.end()){

city.insert(*cityNow);

cityptr.insert({ cityNow->name, cityNow });

}

else {

cityNow->money += ptr->second->money;

city.erase(*ptr->second);

ptr->second = cityNow;

city.insert(*cityNow);

}

void insertResult(set<City>& city, map<string, int> &result, int day) {

string name;

if (city.size() >= 2) {

set<City>::iterator ptr = city.begin();

if (ptr->money - (++ptr)->money != 0) name = city.begin()->name;

else return;

}

else if (city.size() == 0) return;

else if (city.size() == 1) name = city.begin()->name;

map<string, int>::iterator ptr;

ptr = result.find(name);

if (ptr != result.end()) ptr->second += day;

else result.insert({ name, day });

}

void del(set<City>& city, map<string, City*>& cityptr, string cityName, long long money) {

map<string, City*>::iterator ptr;

ptr = cityptr.find(cityName);

City *tempCity = new City(cityName, ptr->second->money - money);

city.erase(*ptr->second);

delete ptr->second;

ptr->second = tempCity;

if (tempCity->money == 0) {

cityptr.erase(ptr);

return;

}

city.insert(*tempCity);

}

int main() {

map<string, pair<string, long long>>Person;

map<string, int>result;

map<string, City*>cityptr;

set<City>city;

int n, m, k, day = 0, preday;

cin >> n;

string namePerson;

string nameCity;

string namePreCity;

long long money;

for (int i = 0; i < n; i++) {

cin >> namePerson >> nameCity >> money;

Person.insert({ namePerson,make_pair(nameCity,money) });

add(city, cityptr, nameCity, money);

}

cin >> m >> k;

for (int i = 0; i < k; i++) {

preday = day;

cin >> day >> namePerson >> nameCity;

insertResult(city, result, day - preday);

map<string, pair<string, long long>>::iterator ptr;

ptr = Person.find(namePerson);

money = ptr->second.second;

add(city, cityptr, nameCity, money);

namePreCity = ptr->second.first;

ptr->second.first = nameCity;

del(city, cityptr, namePreCity, money);

}

insertResult(city, result, m - day);

for (auto i : result) {

if (i.second == 0)continue;

cout << i.first << " " << i.second << endl;

}

Python 牛顿法求解多元非线性方程组

代码使用方式可见代码中示例代码。需安装numpy和sympy库。

from math import *
from sympy import *
import sympy
import numpy as np


class Newton:
    def __init__(self, equations):
        self.equations = equations
        self.n = len(equations)
        self.x = sympy.symbols(" ".join("x{}".format(i) for i in range(self.n)) + " x{}".format(self.n), real=True)
        self.equationsSymbol = [equations[i](self.x) for i in range(self.n)]
        # 初始化 Jacobian 矩阵
        self.J = np.zeros(self.n * self.n, dtype=sympy.core.add.Add).reshape(self.n, self.n)
        for i in range(self.n):
            for j in range(self.n):
                self.J[i][j] = sympy.diff(self.equationsSymbol[i], self.x[j])

    def cal_J(self, x):
        dict = {self.x[i]: x[i] for i in range(self.n)}
        J = np.zeros(self.n * self.n).reshape(self.n, self.n)
        for i in range(self.n):
            for j in range(self.n):
                J[i][j] = self.J[i][j].subs(dict)
        return J

    def cal_f(self, x):
        f = np.zeros(self.n)
        for i in range(self.n):
            f[i] = self.equations[i](x)
        f.reshape(self.n, 1)
        return f

    def newtonByStep(self, x0, step):
        x0 = np.array(x0)
        for i in range(step):
            x0 = x0 - np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)
            print("Step {}:".format(i + 1), ", ".join(["x{} = {}".format(j + 1, x0[j]) for j in range(self.n)]))
        return x0

    def newtonByEpsilon(self, x0, epsilon):
        error = float("inf")
        while error >= epsilon:
            cal = np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)
            error = max(abs(cal))
            x0 = x0 - cal
        print(x0)
        return x0


if __name__ == "__main__":
    # equations 为 equation = 0 的方程组
    # x0 为初始值 step 为迭代次数 epsilon 为精度

    # 多元非线性使用方法
    equations = [lambda x: cos(0.4 * x[1] + x[0] ** 2) + x[0] ** 2 + x[1] ** 2 - 1.6,
                 lambda x: 1.5 * x[0] ** 2 - x[1] ** 2 / 0.36 - 1,
                 lambda x: 3 * x[0] + 4 * x[1] + 5 * x[2]]

    newton = Newton(equations)
    newton.newtonByStep([1, 1, 1], 5)
    newton.newtonByEpsilon([1, 1, 1], 0.001)

    # 一元非线性使用方法
    equations = [lambda x: cos(x[0]) + sin(x[0])]

    newton = Newton(equations)
    newton.newtonByStep([1], 5)
    newton.newtonByEpsilon([1], 0.001)

    # 考试题
    equations = [lambda x: 2 * x[0] ** 3 - x[1] ** 2 - 1,
                 lambda x: x[0] * x[1] ** 3 - x[1] - 4]

    newton = Newton(equations)
    newton.newtonByStep([1.2, 1.7], 5)

from math import *

from sympy import *

import sympy

import numpy as np

class Newton:

def __init__(self, equations):

self.equations = equations

self.n = len(equations)

self.x = sympy.symbols(" ".join("x{}".format(i) for i in range(self.n)) + " x{}".format(self.n), real=True)

self.equationsSymbol = [equations[i](self.x) for i in range(self.n)]

# 初始化 Jacobian 矩阵

self.J = np.zeros(self.n * self.n, dtype=sympy.core.add.Add).reshape(self.n, self.n)

for i in range(self.n):

for j in range(self.n):

self.J[i][j] = sympy.diff(self.equationsSymbol[i], self.x[j])

def cal_J(self, x):

dict = {self.x[i]: x[i] for i in range(self.n)}

J = np.zeros(self.n * self.n).reshape(self.n, self.n)

for i in range(self.n):

for j in range(self.n):

J[i][j] = self.J[i][j].subs(dict)

return J

def cal_f(self, x):

f = np.zeros(self.n)

for i in range(self.n):

f[i] = self.equations[i](x)

f.reshape(self.n, 1)

return f

def newtonByStep(self, x0, step):

x0 = np.array(x0)

for i in range(step):

x0 = x0 - np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)

print("Step {}:".format(i + 1), ", ".join(["x{} = {}".format(j + 1, x0[j]) for j in range(self.n)]))

return x0

def newtonByEpsilon(self, x0, epsilon):

error = float("inf")

while error >= epsilon:

cal = np.linalg.pinv(self.cal_J(x0)) @ self.cal_f(x0)

error = max(abs(cal))

x0 = x0 - cal

print(x0)

return x0

if __name__ == "__main__":

# equations 为 equation = 0 的方程组

# x0 为初始值 step 为迭代次数 epsilon 为精度

# 多元非线性使用方法

equations = [lambda x: cos(0.4 * x[1] + x[0] ** 2) + x[0] ** 2 + x[1] ** 2 - 1.6,

lambda x: 1.5 * x[0] ** 2 - x[1] ** 2 / 0.36 - 1,

lambda x: 3 * x[0] + 4 * x[1] + 5 * x[2]]

newton = Newton(equations)

newton.newtonByStep([1, 1, 1], 5)

newton.newtonByEpsilon([1, 1, 1], 0.001)

# 一元非线性使用方法

equations = [lambda x: cos(x[0]) + sin(x[0])]

newton = Newton(equations)

newton.newtonByStep([1], 5)

newton.newtonByEpsilon([1], 0.001)

# 考试题

equations = [lambda x: 2 * x[0] ** 3 - x[1] ** 2 - 1,

lambda x: x[0] * x[1] ** 3 - x[1] - 4]

newton = Newton(equations)

newton.newtonByStep([1.2, 1.7], 5)

计算数学课程参考 Python代码

求积分以及方程的解需要自己修改f(x)函数里面的返回值。求非线性方程的解需要安装sympy包。

牛顿法解一元非线性方程

from math import *
from sympy import *


def f(x):
    return x ** 2 - 3 * x ** 3


def newtonByStep(x0, step):
    x = symbols("x", real=True)
    _f = f(x)
    f_x = diff(_f, x)
    for i in range(step):
        x1 = x0 - f(x0) / f_x.subs({x: x0})
        print("Step{}: x = {}, f(x) = {}".format(i + 1, x1, f(x1)))
        x0 = x1


def newtonByEpsilon(x0, epsilon):
    x = symbols("x", real=True)
    _f = f(x)
    f_x = diff(_f, x)
    error = float("inf")
    while error >= epsilon:
        x1 = x0 - f(x0) / f_x.subs({x: x0})
        error = abs(x1 - x0)
        x0 = x1
    return x0


if __name__ == "__main__":
    x0 = 1.04
    newtonByStep(x0, 5)
    print(newtonByEpsilon(x0, 0.001))

from math import *

from sympy import *

def f(x):

return x ** 2 - 3 * x ** 3

def newtonByStep(x0, step):

x = symbols("x", real=True)

_f = f(x)

f_x = diff(_f, x)

for i in range(step):

x1 = x0 - f(x0) / f_x.subs({x: x0})

print("Step{}: x = {}, f(x) = {}".format(i + 1, x1, f(x1)))

x0 = x1

def newtonByEpsilon(x0, epsilon):

x = symbols("x", real=True)

_f = f(x)

f_x = diff(_f, x)

error = float("inf")

while error >= epsilon:

x1 = x0 - f(x0) / f_x.subs({x: x0})

error = abs(x1 - x0)

x0 = x1

return x0

if __name__ == "__main__":

x0 = 1.04

newtonByStep(x0, 5)

print(newtonByEpsilon(x0, 0.001))

牛顿法解二元非线性方程

from math import *
from sympy import *


def f(x, y):
    return 2 * x ** 3 - y ** 2 - 1


def g(x, y):
    return x * y ** 3 - y - 4


def newtonByStep(x0, y0, step):
    x, y = symbols("x y", real=True)
    _f = f(x, y)
    _g = g(x, y)
    f_x = diff(_f, x)
    f_y = diff(_f, y)
    g_x = diff(_g, x)
    g_y = diff(_g, y)
    for i in range(step):
        x1 = x0 + (f(x0, y0) * g_y.subs({x: x0, y: y0}) - g(x0, y0) * f_y.subs({x: x0, y: y0})) / (
                g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))
        y1 = y0 + (g(x0, y0) * f_x.subs({x: x0, y: y0}) - f(x0, y0) * g_x.subs({x: x0, y: y0})) / (
                g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))
        print("Step{}: x = {}, y = {}, f(x, y) = {}, g(x, y) = {}".format(i + 1, x1, y1, f(x1, y1), g(x1, y1)))
        x0 = x1
        y0 = y1


def newtonByEpsilon(x0, y0, epsilon):
    x, y = symbols("x y", real=True)
    _f = f(x, y)
    _g = g(x, y)
    f_x = diff(_f, x)
    f_y = diff(_f, y)
    g_x = diff(_g, x)
    g_y = diff(_g, y)
    error = float("inf")
    while error >= epsilon:
        x1 = x0 + (f(x0, y0) * g_y.subs({x: x0, y: y0}) - g(x0, y0) * f_y.subs({x: x0, y: y0})) / (
                g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))
        y1 = y0 + (g(x0, y0) * f_x.subs({x: x0, y: y0}) - f(x0, y0) * g_x.subs({x: x0, y: y0})) / (
                g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))
        error = max(abs(x1 - x0), abs(y1 - y0))
        x0 = x1
        y0 = y1
    return x0, y0


if __name__ == "__main__":
    x0 = 1.2
    y0 = 1.7
    newtonByStep(x0, y0, 5)
    print(newtonByEpsilon(x0, y0, 0.001))

from math import *

from sympy import *

def f(x, y):

return 2 * x ** 3 - y ** 2 - 1

def g(x, y):

return x * y ** 3 - y - 4

def newtonByStep(x0, y0, step):

x, y = symbols("x y", real=True)

_f = f(x, y)

_g = g(x, y)

f_x = diff(_f, x)

f_y = diff(_f, y)

g_x = diff(_g, x)

g_y = diff(_g, y)

for i in range(step):

x1 = x0 + (f(x0, y0) * g_y.subs({x: x0, y: y0}) - g(x0, y0) * f_y.subs({x: x0, y: y0})) / (

g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))

y1 = y0 + (g(x0, y0) * f_x.subs({x: x0, y: y0}) - f(x0, y0) * g_x.subs({x: x0, y: y0})) / (

g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))

print("Step{}: x = {}, y = {}, f(x, y) = {}, g(x, y) = {}".format(i + 1, x1, y1, f(x1, y1), g(x1, y1)))

x0 = x1

y0 = y1

def newtonByEpsilon(x0, y0, epsilon):

x, y = symbols("x y", real=True)

_f = f(x, y)

_g = g(x, y)

f_x = diff(_f, x)

f_y = diff(_f, y)

g_x = diff(_g, x)

g_y = diff(_g, y)

error = float("inf")

while error >= epsilon:

x1 = x0 + (f(x0, y0) * g_y.subs({x: x0, y: y0}) - g(x0, y0) * f_y.subs({x: x0, y: y0})) / (

g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))

y1 = y0 + (g(x0, y0) * f_x.subs({x: x0, y: y0}) - f(x0, y0) * g_x.subs({x: x0, y: y0})) / (

g_x.subs({x: x0, y: y0}) * f_y.subs({x: x0, y: y0}) - f_x.subs({x: x0, y: y0}) * g_y.subs({x: x0, y: y0}))

error = max(abs(x1 - x0), abs(y1 - y0))

x0 = x1

y0 = y1

return x0, y0

if __name__ == "__main__":

x0 = 1.2

y0 = 1.7

newtonByStep(x0, y0, 5)

print(newtonByEpsilon(x0, y0, 0.001))

5种积分方法

from math import *


def sgn(x):
    if x > 0:
        return 1
    elif x == 0:
        return 0
    else:
        return -1


def func(x):
    # return sgn(sin(pi / sin(x)))
    return x ** 2 + 2


def simpson(begin, between, i):
    a = begin + (i - 1) * between
    b = begin + i * between
    return between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6


def rightRectangle(begin, between, i):
    return between * func(begin + between * i)


def leftRectangle(begin, between, i):
    return between * func(begin + between * (i - 1))


def midRectangle(begin, between, i):
    return between * func(begin + between / 2 * (2 * i - 1))


def trapezoidRectangle(begin, between, i):
    return between / 2 * (func(begin + between * (i - 1)) + func(begin + between * i))


def cal_IntegralByEpsilon(begin, end, epsilon, method):
    n = 1
    result = 0
    preResult = float("inf")
    while abs(preResult - result) >= epsilon:
        preResult = result
        result = 0
        n *= 2
        between = (end - begin) / n
        for i in range(n):
            try:
                result += method(begin, between, i + 1)
            except:
                return "Integrated function has discontinuity or does not defined in current interval"
    return result


def cal_IntegralByN(begin, end, n, method):
    result = 0
    between = (end - begin) / n
    for i in range(n):
        try:
            result += method(begin, between, i + 1)
        except:
            return "Integrated function has discontinuity or does not defined in current interval"


if __name__ == "__main__":
    begin = 0.2
    end = 0.5
    epsilon = 0.0001
    print(cal_IntegralByEpsilon(begin, end, epsilon, simpson))
    print(cal_IntegralByEpsilon(begin, end, epsilon, rightRectangle))
    print(cal_IntegralByEpsilon(begin, end, epsilon, leftRectangle))
    print(cal_IntegralByEpsilon(begin, end, epsilon, midRectangle))
    print(cal_IntegralByEpsilon(begin, end, epsilon, trapezoidRectangle))

    begin = 0.2
    end = 0.5
    n = 10
    print(cal_IntegralByN(begin, end, n, simpson))
    print(cal_IntegralByN(begin, end, n, rightRectangle))
    print(cal_IntegralByN(begin, end, n, leftRectangle))
    print(cal_IntegralByN(begin, end, n, midRectangle))
    print(cal_IntegralByN(begin, end, n, trapezoidRectangle))

from math import *

def sgn(x):

if x > 0:

return 1

elif x == 0:

return 0

else:

return -1

def func(x):

# return sgn(sin(pi / sin(x)))

return x ** 2 + 2

def simpson(begin, between, i):

a = begin + (i - 1) * between

b = begin + i * between

return between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6

def rightRectangle(begin, between, i):

return between * func(begin + between * i)

def leftRectangle(begin, between, i):

return between * func(begin + between * (i - 1))

def midRectangle(begin, between, i):

return between * func(begin + between / 2 * (2 * i - 1))

def trapezoidRectangle(begin, between, i):

return between / 2 * (func(begin + between * (i - 1)) + func(begin + between * i))

def cal_IntegralByEpsilon(begin, end, epsilon, method):

n = 1

result = 0

preResult = float("inf")

while abs(preResult - result) >= epsilon:

preResult = result

result = 0

n *= 2

between = (end - begin) / n

for i in range(n):

try:

result += method(begin, between, i + 1)

except:

return "Integrated function has discontinuity or does not defined in current interval"

return result

def cal_IntegralByN(begin, end, n, method):

result = 0

between = (end - begin) / n

for i in range(n):

try:

result += method(begin, between, i + 1)

except:

return "Integrated function has discontinuity or does not defined in current interval"

if __name__ == "__main__":

begin = 0.2

end = 0.5

epsilon = 0.0001

print(cal_IntegralByEpsilon(begin, end, epsilon, simpson))

print(cal_IntegralByEpsilon(begin, end, epsilon, rightRectangle))

print(cal_IntegralByEpsilon(begin, end, epsilon, leftRectangle))

print(cal_IntegralByEpsilon(begin, end, epsilon, midRectangle))

print(cal_IntegralByEpsilon(begin, end, epsilon, trapezoidRectangle))

begin = 0.2

end = 0.5

n = 10

print(cal_IntegralByN(begin, end, n, simpson))

print(cal_IntegralByN(begin, end, n, rightRectangle))

print(cal_IntegralByN(begin, end, n, leftRectangle))

print(cal_IntegralByN(begin, end, n, midRectangle))

print(cal_IntegralByN(begin, end, n, trapezoidRectangle))

绝对误差与相对误差

def errorByreal(real, standard):
    absolute = abs(real - standard)
    realative = absolute / real
    print("Absolute error = %e \nRelative error = %e" % (absolute, realative))


def errorByEpsilon(real, epsilon):
    accuracy = len(str(epsilon).split(".")[-1])
    standard = round(real, accuracy)
    absolute = abs(real - standard)
    realative = absolute / real
    print("Absolute error = %e \nRelative error = %e" % (absolute, realative))


if __name__ == "__main__":
    real = 0.63840532
    epsilon = 0.01
    errorByEpsilon(real, epsilon)

    real = 0.63840532
    standard = 0.6384
    errorByreal(real, standard)

def errorByreal(real, standard):

absolute = abs(real - standard)

realative = absolute / real

print("Absolute error = %e \nRelative error = %e" % (absolute, realative))

def errorByEpsilon(real, epsilon):

accuracy = len(str(epsilon).split(".")[-1])

standard = round(real, accuracy)

absolute = abs(real - standard)

realative = absolute / real

print("Absolute error = %e \nRelative error = %e" % (absolute, realative))

if __name__ == "__main__":

real = 0.63840532

epsilon = 0.01

errorByEpsilon(real, epsilon)

real = 0.63840532

standard = 0.6384

errorByreal(real, standard)

最小二乘法一元多项式拟合

def gauss(matrix):
    array = []
    for i in range(0, len(matrix) - 1):
        aii = matrix[i][i]
        if aii == 0:
            continue
        for j in range(i + 1, len(matrix)):
            aji = matrix[j][i]
            for k in range(i, len(matrix[j])):
                matrix[j][k] = matrix[j][k] - aji / aii * matrix[i][k]
    for i in range(len(matrix) - 1, -1, -1):
        if matrix[i][i] == 0:
            exit("The system has no roots of equations or has an infinite set of them.")
        array.append(matrix[i][-1] / matrix[i][i])
        for j in range(0, i):
            matrix[j][-1] -= matrix[j][i] * array[-1]
    array.reverse()
    return array


def leastSquareRegression(points, order):
    result = []
    for i in range(order + 1):
        sum = 0
        for j in range(len(points)):
            sum += points[j][0] ** i * points[j][1]
        result.append(sum)
    matrix = []
    for i in range(order + 1):
        matrix.append([])
        for j in range(2 * order):
            sum = 0
            for k in range(len(points)):
                sum += points[k][0] ** (i + j)
            matrix[i].append(sum)
        matrix[i].append(result[i])
    return gauss(matrix)


def R2_func(y_test, y):
    return 1 - ((y_test - y) ** 2).sum() / ((y.mean() - y) ** 2).sum()


def test(points, result):
    import numpy as np
    approach = []
    for i in range(len(points)):
        sum = 0
        for j in range(order + 1):
            sum += points[i][0] ** j * result[j]
        approach.append(sum)
    print("R^2 result =", R2_func(np.array(approach), np.array([i[1] for i in points])))


if __name__ == "__main__":
    points = [[1, 2], [3, 4], [5, 77], [8, -3]]
    order = 2
    result = leastSquareRegression(points, order)
    print("y =", " + ".join(["{}X^{}".format(result[i], i) for i in range(order + 1)]))
    test(points, result)

def gauss(matrix):

array = []

for i in range(0, len(matrix) - 1):

aii = matrix[i][i]

if aii == 0:

continue

for j in range(i + 1, len(matrix)):

aji = matrix[j][i]

for k in range(i, len(matrix[j])):

matrix[j][k] = matrix[j][k] - aji / aii * matrix[i][k]

for i in range(len(matrix) - 1, -1, -1):

if matrix[i][i] == 0:

exit("The system has no roots of equations or has an infinite set of them.")

array.append(matrix[i][-1] / matrix[i][i])

for j in range(0, i):

matrix[j][-1] -= matrix[j][i] * array[-1]

array.reverse()

return array

def leastSquareRegression(points, order):

result = []

for i in range(order + 1):

sum = 0

for j in range(len(points)):

sum += points[j][0] ** i * points[j][1]

result.append(sum)

matrix = []

for i in range(order + 1):

matrix.append([])

for j in range(2 * order):

sum = 0

for k in range(len(points)):

sum += points[k][0] ** (i + j)

matrix[i].append(sum)

matrix[i].append(result[i])

return gauss(matrix)

def R2_func(y_test, y):

return 1 - ((y_test - y) ** 2).sum() / ((y.mean() - y) ** 2).sum()

def test(points, result):

import numpy as np

approach = []

for i in range(len(points)):

sum = 0

for j in range(order + 1):

sum += points[i][0] ** j * result[j]

approach.append(sum)

print("R^2 result =", R2_func(np.array(approach), np.array([i[1] for i in points])))

if __name__ == "__main__":

points = [[1, 2], [3, 4], [5, 77], [8, -3]]

order = 2

result = leastSquareRegression(points, order)

print("y =", " + ".join(["{}X^{}".format(result[i], i) for i in range(order + 1)]))

test(points, result)

Python 求积分的几种方法

不同的方法最大区别在于公式的不同/矩形区域的不同，依据公式而定，需要什么代码你可以自己进行更改，很简单。
本文代码为不严格代码，最严格的代码为先判断连续性。当你需要对大量函数积分（懒得写函数连续性检测的时候）可以使用本文代码用try…except…不严格替代。
看不懂公式的可以去 CSDN 查看图片。

Simpson积分法:
\(\int _ { a } ^ { b } f ( x ) d x \approx \frac { ( b – a ) ( f ( a ) + f ( b ) + 4 f ( \frac { a + b } { 2 } ) ) } { 6 }\)

 def calculate_integral(begin, end, func, epsilon):
        n = 1
        result = 0
        preResult = float("inf")
        while abs(preResult - result) > epsilon:
            preResult = result
            result = 0
            n *= 2
            between = (end - begin) / n
            for i in range(n):
                try:
                    a = begin + i * between
                    b = begin + (i + 1) * between
                    result += between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6
                except:
                    return "Integrated function has discontinuity or does not defined in current interval"
        return result

def calculate_integral(begin, end, func, epsilon):

n = 1

result = 0

preResult = float("inf")

while abs(preResult - result) > epsilon:

preResult = result

result = 0

n *= 2

between = (end - begin) / n

for i in range(n):

try:

a = begin + i * between

b = begin + (i + 1) * between

result += between * (func(a) + func(b) + 4 * func((a + b) / 2)) / 6

except:

return "Integrated function has discontinuity or does not defined in current interval"

return result

右矩形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { n } f ( a + \frac { b – a } { n } i )\)

def calculate_integral(begin, end, func, epsilon):
    n = 1
    result = 0
    preResult = float("inf")
    while abs(preResult - result) > epsilon:
        preResult = result
        result = 0
        n *= 2
        for i in range(n):
            try:
                result += (end - begin) / n * func(begin + (end - begin) / n * (i + 1))
            except:
                return "Integrated function has discontinuity or does not defined in current interval"
    return result

def calculate_integral(begin, end, func, epsilon):

n = 1

result = 0

preResult = float("inf")

while abs(preResult - result) > epsilon:

preResult = result

result = 0

n *= 2

for i in range(n):

try:

result += (end - begin) / n * func(begin + (end - begin) / n * (i + 1))

except:

return "Integrated function has discontinuity or does not defined in current interval"

return result

左矩形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { n } f ( a + \frac { b – a } { n } (i – 1) )\)

中矩形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { n } f ( a + \frac { (b – a)} { 2 n }(2i – 1 ) )\)

梯形积分公式
\(\int _ { a } ^ { b } f ( x ) d x = \lim _ { n \rightarrow \infty } \sum _ { i = 1 } ^ { n } \frac { b – a } { 2n } (f ( a + \frac { b – a } { n }( i – 1) ) + f ( a + \frac { b – a } { n } i ))\)

Timus #1628. White Streaks

题目：有m*n的矩形，每个方块默认为白色，现给k个方块图上黑色，问有多少连续的num*1或1*num的长方形。如果一个1*1的方块周围没有白色方块，那么这个方块也包含在答案里面。

解析：行列遍历。对行和列分别进行查找，对行和列查找过程中单独成一个方块的坐标分别记录在两个集合中，最终只需要在结果上加上两个集合的交集中元素的个数。

#include<iostream>
#include<algorithm>
#include<iterator>
#include<set>
using namespace std;

struct Point
{
	int x, y;
	bool operator< (const Point & a)const
	{
		if (this->x == a.x)return this->y < a.y;
		return this->x < a.x;
	}
}point[60000];

int fx(Point& a, Point& b) { 
	if (a.x == b.x)return a.y < b.y;
	return a.x < b.x; 
}

int fy(Point& a, Point& b) {
	if (a.y == b.y)return a.x < b.x;
	return a.y < b.y;
}

void insert(set<Point>& a, int x, int y) {
	Point tempPoint;
	tempPoint.x = x;
	tempPoint.y = y;
	a.insert(tempPoint);

}

int main() {
	int m, n, k, ans = 0;
	cin >> n >> m >> k;
	set<Point> A;
	set<Point> B;
	set<Point> C;
	for (int i = 0; i < k; i++)cin >> point[i].x >> point[i].y;
	int now = 0, row = 0, col = 0;

	if (n == 1) {
		sort(point, point + k, fx);
		if (point[0].y == 1) ans++; 
		for (int i = 1; i < k; i++) {
			if (point[i].y - point[i - 1].y == 1) ans++;
		}
		if (point[k - 1].y == m) ans++;
		cout << k + 1 - ans;
		return 0;
	}

	if (m == 1) {
		sort(point, point + k, fy);
		if (point[0].x == 1) ans++;
		for (int i = 1; i < k; i++) {
			if (point[i].x - point[i - 1].x == 1) ans++;
		}
		if (point[k - 1].x == n) ans++;
		cout << k + 1 - ans;
		return 0;
	}

	sort(point, point + k, fx);
	while (row < n) {
		if (point[now].x != row + 1) {
			ans++;
			row++;
			continue;
		}
		if (point[now].y == 2) insert(A, point[now].x, 1);
		else if (point[now].y > 2) ans++;
		now++;
		while (now < k)
		{
			if (point[now].x == point[now - 1].x) {
				if (point[now].y - point[now - 1].y == 2) insert(A, point[now].x, point[now].y - 1);
				else if (point[now].y - point[now - 1].y > 2) ans++;
			}
			else {
				break;
			}
			now++;
		}
		if (point[now - 1].y == m - 1) insert(A, point[now - 1].x, m);
		else if (point[now - 1].y < m - 1) ans++;
		row++;
	}

	sort(point, point + k, fy);
	now = 0;
	while (col < m) {
		if (point[now].y != col + 1) {
			ans++;
			col++;
			continue;
		}
		if (point[now].x == 2) insert(B, 1, point[now].y);
		else if (point[now].x > 2) ans++;
		now++;
		while (now < k)
		{
			if (point[now].y == point[now - 1].y) {
				if (point[now].x - point[now - 1].x == 2) insert(B, point[now].x - 1, point[now].y);
				else if (point[now].x - point[now - 1].x > 2) ans++;
			}
			else {
				break;
			}
			now++;
		}
		if (point[now - 1].x == n - 1) insert(B, n, point[now - 1].y);
		else if (point[now - 1].x < n - 1) ans++;
		col++;
	}

	set_intersection(A.begin(), A.end(), B.begin(), B.end(), inserter(C, C.begin()));
	ans += C.size();
	cout << ans << endl;

}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

#include<iostream>

#include<algorithm>

#include<iterator>

#include<set>

using namespace std;

struct Point

{

int x, y;

bool operator< (const Point & a)const

{

if (this->x == a.x)return this->y < a.y;

return this->x < a.x;

}

}point[60000];

int fx(Point& a, Point& b) {

if (a.x == b.x)return a.y < b.y;

return a.x < b.x;

}

int fy(Point& a, Point& b) {

if (a.y == b.y)return a.x < b.x;

return a.y < b.y;

}

void insert(set<Point>& a, int x, int y) {

Point tempPoint;

tempPoint.x = x;

tempPoint.y = y;

a.insert(tempPoint);

}

int main() {

int m, n, k, ans = 0;

cin >> n >> m >> k;

set<Point> A;

set<Point> B;

set<Point> C;

for (int i = 0; i < k; i++)cin >> point[i].x >> point[i].y;

int now = 0, row = 0, col = 0;

if (n == 1) {

sort(point, point + k, fx);

if (point[0].y == 1) ans++;

for (int i = 1; i < k; i++) {

if (point[i].y - point[i - 1].y == 1) ans++;

}

if (point[k - 1].y == m) ans++;

cout << k + 1 - ans;

return 0;

}

if (m == 1) {

sort(point, point + k, fy);

if (point[0].x == 1) ans++;

for (int i = 1; i < k; i++) {

if (point[i].x - point[i - 1].x == 1) ans++;

}

if (point[k - 1].x == n) ans++;

cout << k + 1 - ans;

return 0;

}

sort(point, point + k, fx);

while (row < n) {

if (point[now].x != row + 1) {

ans++;

row++;

continue;

}

if (point[now].y == 2) insert(A, point[now].x, 1);

else if (point[now].y > 2) ans++;

now++;

while (now < k)

{

if (point[now].x == point[now - 1].x) {

if (point[now].y - point[now - 1].y == 2) insert(A, point[now].x, point[now].y - 1);

else if (point[now].y - point[now - 1].y > 2) ans++;

}

else {

break;

}

now++;

}

if (point[now - 1].y == m - 1) insert(A, point[now - 1].x, m);

else if (point[now - 1].y < m - 1) ans++;

row++;

}

sort(point, point + k, fy);

now = 0;

while (col < m) {

if (point[now].y != col + 1) {

ans++;

col++;

continue;

}

if (point[now].x == 2) insert(B, 1, point[now].y);

else if (point[now].x > 2) ans++;

now++;

while (now < k)

{

if (point[now].y == point[now - 1].y) {

if (point[now].x - point[now - 1].x == 2) insert(B, point[now].x - 1, point[now].y);

else if (point[now].x - point[now - 1].x > 2) ans++;

}

else {

break;

}

now++;

}

if (point[now - 1].x == n - 1) insert(B, n, point[now - 1].y);

else if (point[now - 1].x < n - 1) ans++;

col++;

}

set_intersection(A.begin(), A.end(), B.begin(), B.end(), inserter(C, C.begin()));

ans += C.size();

cout << ans << endl;

}