URAL 1828. Approximation by a Progression 数学

1828. Approximation by a Progression

Time limit: 0.5 second
Memory limit: 64 MB

Your are given a sequence of integers  a 1, …,  an. Find an arithmetic progression  b 1, …,  bn for which the value ∑( ai −  bi) 2 is minimal. The elements of the progression can be non-integral.

Input

The first line contains the number  n of elements in the sequence (2 ≤  n ≤ 10 4). In the second line you are given the integers  a 1, …,  an; their absolute values do not exceed 10 4.

Output

Output two numbers separated with a space: the first term of the required arithmetic progression and its difference, with an absolute or relative error of at most 10 −6. It is guaranteed that the answer is unique for all input data.

Samples

input	output
4 0 6 10 15	0.400 4.900
4 -2 -2 -2 -2	-2 0

题意 ： 给定整数数组 ai。求一个等差数列bi。 等差数列bi要能使 ∑( ai  −  bi ) 2   最小的。输出 b1，d。两者可以是小数。

做法：推公式比较麻烦。步骤如下：

为方便计算 设 第零为b0，方差为d。  那么 bi=b0+i*d 。b1=b0+d。

然后带入 ∑ ( ai  −  bi ) 2

带入 bi=b0+i*d 。得到 ∑（b0^2+(i*d)^2+2*b0*i*d-2*b0*ai-2*i*d*ai+ai^2）

把b0提出来  得到 ∑ （ b0^2+b0*(2*i*d-2*ai)+(i*d)^2-2*i*d*ai+ai^2）

把d提出来   得到∑ （ i^2* d^2+d*(2*b0*i-2*i*ai)+ b0^2-2*b0*ai+ai^2）

可以发现 两个 分别关于 b0 和d 的二元一次函数都是开口向上的。所以极值落在 -b/2a上。

把∑化到括号里面。

设常数 c1 等于∑bi   ，c2 =∑（i*bi）  c3= ∑i   c4=∑i^2   ∑中i都是从1到n的累加

最后可以化简得到 b0的极值点 b0=-b/2a=(-d*c3+c1)/n，d的极值点  d= =-b/2a=(-b0*c3+c2)/c4。

这两个二元一次方程。 两个式子，可以解出 d =(（c1*c3）/n-c2)/(c3*c3/n-c4)   b0=(-d*c3+c1)/n

方程里a1就是等于上面的b0.

#include <stdio.h>
#include <stdlib.h>
#include <string.h> 
#include <string>
#include <iostream>
#include <algorithm>
using namespace std; 
#include <vector> 


int b[11000];
int main()
{
	double n;
	while(scanf("%lf",&n)!=EOF)
	{
		for(int i=1;i<=n;i++)
		{
			cin>>b[i];
		}


		//c1 bi he
		double c1=0;
		for(int i=1;i<=n;i++)
		{
			c1+=b[i];
		}
		//c2 ibi he
		double c2=0;
		for(int i=1;i<=n;i++)
		{
			c2+=b[i]*i;
		}

		//c3 i he
		double c3=0;
		for(int i=1;i<=n;i++)
		{
			c3+=i;
		}

		//c4 i^2 he
		double c4=0;
		for(int i=1;i<=n;i++)
		{
			c4+=i*i;
		}

		double d=1.0*((c1*c3)/n-c2)/((c3*c3)/n-c4);
		double a=1.0*(-d*c3+c1)/n;
		a+=d;

		if(d<0.000001&&d>-0.000001)
			d=0;
		printf("%lf %lf\n",a,d);

	}
	return 0;
}