Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <bits/stdc++.h>
- #include <experimental/array>
- #define mar experimental::make_array
- using namespace std;
- array<array<long double, 100>, 50> table {};
- array<long double, 50> val {};
- int n{};
- int mfd=4;
- long double h=0.1;
- long double y {1.683635};
- array<long double, 3> arX{0.151316, 0.815359, 0.281586};
- void finDiff(int k);
- long double intBegTable(long double x);
- long double intEndTable(long double x);
- long double intCenterTable(long double x);
- long double revIntTable(long double y);
- void crTable();
- void solveForX();
- void solveForY();
- long double findT(long double y, int ind);
- long double fact(int a);
- int main(){
- freopen("in","r",stdin);
- while(cin >> val[n] >> table[n<<1][0] && ++n){}
- crTable();
- solveForX();
- solveForY();
- return 0;
- }
- void solveForX(){
- for(int i=0; i<3; i++){
- double res=-1337;
- for(int j=0; j<n; j++){
- if (arX[i]<val[j]){
- if (j<n/3)
- res=intBegTable(arX[i]);
- if (j>=n/3 && j<2*n/3)
- res=intCenterTable(arX[i]);
- if (j>=2*n/3)
- res=intEndTable(arX[i]);
- break;
- }
- }
- cout << "\np("<< arX[i] << ") ~=~ " << setprecision(7)<< res << "\n\n";
- if(res!=-1337)
- continue;
- res=intEndTable(arX[i]);
- }
- double res{};
- res=intBegTable(arX[2]);
- cout << "p("<< arX[2] << ") ~=~ " << setprecision(7)<< res << "\n\n";
- res=intEndTable(arX[2]);
- cout << "\np("<< arX[2] << ") ~=~ " << setprecision(7) << res << "\n\n";
- }
- void solveForY(){
- double res{},
- resX{};
- res=revIntTable(y);
- cout << "\n\np^-1("<< y << ") ~=~ " << setprecision(7) << res << "\n\n";
- cout << "Check: \n";
- for(int j=0; j<n; j++){
- if (res<val[j]){
- if (j<n/3)
- resX=intBegTable(res);
- if (j>=n/3 && j<2*n/3)
- resX=intCenterTable(res);
- if (j>=2*n/3)
- resX=intEndTable(res);
- break;
- }
- }
- cout << "\np("<< res << ") ~=~ " << resX << "\n\n";
- cout << y << " == " << resX << "?\n\n\n\n";
- }
- void finDiff(int k){
- for(int i=k; i<(n<<1)-k; i++)
- table[i][k]=table[i+1][k-1]-table[i-1][k-1];
- }
- void crTable(){
- mfd=min(mfd, n-1);
- cout << "x\t\t px\t ";
- for(int i=1; i<=mfd; i++)
- finDiff(i),
- cout << "d^{" << i << "}" << mar("\t", "\n\n")[i==mfd];
- for(int i=0; i<(n<<1)-1; i++){
- cout << mar("\t", to_string(val[i>>1]))[i%2==0] << " |\t";
- for(int j=0; j<=mfd; j++)
- cout << " \a"[table[i][j]<0] << mar("\t", to_string(table[i][j]))[table[i][j]!=0] << "\a\n"[j==mfd];
- }
- }
- long double intBegTable(long double x){
- int ind{};
- for(ind=0; ind<n; ind++)
- if (val[ind]>x){
- ind--;
- break;
- }
- long double t=(x-val[ind])/h,
- res=0,
- tt=1;
- cout << "\n\nInterp BEGIN table start, x=" << x << ", x_0="<<val[ind] << ", t=" << t << endl;
- for(int i=0; i<5; i++){
- //cout << table[(ind<<1)+i][i] << ' ' << tt << endl;
- res+=tt*table[(ind<<1)+i][i]/fact(i), tt=tt*(t-i);
- }
- return res;
- }
- long double intEndTable(long double x){
- int ind{};
- for(ind=n-1; ind>=0; ind--)
- if (val[ind]<x){
- ind++;
- break;
- }
- long double t=(x-val[ind])/h,
- res=0,
- tt=1;
- cout << "Interp END table start, x=" << x << ", x_0="<<val[ind] << ", t=" << t << endl;
- for(int i=0; i<5; i++){
- //cout << table[(ind<<1)-i][i] << ' ' << tt << endl;
- res+=tt*table[(ind<<1)-i][i]/fact(i), tt=tt*(t+i);
- }
- return res;
- }
- long double intCenterTable(long double x){
- int ind{};
- for(ind=0; ind<n; ind++)
- if (val[ind]>x){
- if (fabs(val[ind]-x)>fabs(val[ind-1]-x))
- ind--;
- break;
- }
- long double t=(x-val[ind])/h,
- res=0,
- tt=1;
- cout << "Interp CENTER table start, x=" << x << ", x_0="<<val[ind] << ", t=" << t << endl;
- for(int i=0; i<5; i++){
- //cout << table[(ind<<1)+1*(i%2==1)][i] << ' ' << tt << endl;
- res+=tt*table[(ind<<1)+1*(i%2==1)][i]/fact(i), tt=tt*(t+(((i+1)>>1))*pow(-1,i));
- }
- return res;
- }
- long double revIntTable(long double y){
- int ind{};
- for(ind=0; ind<(n<<1)-1; ind++)
- if (table[ind<<1][0]>y){
- ind--;
- break;
- }
- long double res=y,
- tt=1,
- t{};
- t=findT(y, ind);
- res=t*h+val[ind];
- return res;
- }
- long double findT(long double y, int ind){
- long double res{-1},
- t{1},
- tt;
- while(fabs(res-t)>1e-4){
- t=res;
- cout << "t=" << t << endl;
- res=0;
- tt=1;
- res=y-table[ind<<1][0]-t*(t-1)/2.*table[(ind<<1)+2][2]-t*(t-1)*(t-2)/6.*table[(ind<<1)+3][3];
- res/=table[(ind<<1)+1][1];
- }
- return res;
- }
- long double fact(int a){
- if (a<2)
- return 1.;
- return a*fact(a-1);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement