序列求和

Description

求A的第1项到第N项的序列之和，精心设计以避免损失精度。

Input Description

第一行是包含一个整数M，表示输入有M个样例，接下来的M行，每行表示一个样例，包含一个整数，表示求和的项数N（N<=1000000000）。

Output Description

输出为对应每组输入的序列求和结果，使用双精度浮点数（double）表示，要求相对误差小于10^(-12)，即为正确答案。请统一使用”%.12e\n”格式输出结果。

Input Sample

3
1
2
3

Output Sample

2.121320343560e+000
4.609396906025e+000
7.938599833904e+000

Idea

上述求和公式的计算结果精度主要受到两方面影响：

当求和公式中的两个加数异号并且绝对值接近时，计算结果的精度会降低；

当求和过程中出现大数加小数的情况时，如果两个操作数的比值低于机器精度，计算结果的精度也会降低

当n很大时，可以用调和序列的欧拉公式代替原序列。

Code

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define Euler 0.577215664901533

int main(){
    int M;
    cin>>M;
    int input;
    double sum,sum_big,sum_small,sum_har;
    for(int i=0;i<M;i++){   
        cin>>input;
        if(input<=100000){
            sum_big=sum_small=0;
            for(int j=1;j<=input;j++){
                double jd=static_cast<double>(j);
                sum_big+=sqrt(jd*jd+1);     
                sum_small+=sqrt(jd*jd-1);
            }
            sum=(3*sum_big-sum_small)/2;
        }
        else if(input<=1e7){
            sum_big=sum_small=0;
            for(int j=1;j<=100000;j++){
                double jd=static_cast<double>(j);
                sum_big+=sqrt(jd*jd+1)+sqrt(jd*jd-1);
            }
            for(int j=input;j>=1;j--){
                double jd=static_cast<double>(j);
                double one=0.000000000000001e15;
                sum_small+=2/(sqrt(jd*jd+one)+sqrt(jd*jd-one));
            }
            sum=sum_big/2+sum_small;
        }
        else{
            sum_big=sum_small=sum_har=0;
            for(int j=1;j<=100000;j++){
                double jd=static_cast<double>(j);
                sum_big+=sqrt(jd*jd+1)+sqrt(jd*jd-1);
            }
            for(int j=1e6;j>=1;j--){
                double jd=static_cast<double>(j);
                double one=0.000000000000001e15;
                sum_small+=2/(sqrt(jd*jd+one)+sqrt(jd*jd-one));
                sum_har+=1/jd;
            }
            sum=sum_small-sum_har+log(input)+Euler+sum_big/2;
        }
        printf("%.12e\n",sum);
    }

    return 0;
}