Untitled

#include<bits/stdc++.h>
#define ll long long
const double dt = 0.5;
const double dx = 0.125;
const double du = 1;
const double dv = du*dt;
const double tmax = 40; 
const double xmax = 100;
const double umin = -7;
const double umax = 0;
const double vmax = 15;
using namespace std;
double J[90][35][807];
pair<int,int> back[90][35][561];
double phi(double v1, double u)
{
    return u*dt*(v1+0.5*u*dt);
}
double beta(int k)
{
    return 1-dt/(tmax-k*dt);
}
int main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    vector<double>v,x,u;
    for(double i=0;i<=vmax;i+=dv) v.push_back(i);
    for(double i=umin;i<=umax;i+=du) u.push_back(i);
    for(double i=0;i<=xmax;i+=dx) x.push_back(i);
    int n=tmax/dt-1;
    for(int i=0;i<v.size();i++)
    {
        for(int j=0;j<x.size();j++)
        {
            if(j==0) J[n][i][j]=0;
            else J[n][i][j]=INT_MAX;
        }
    }
    double newj=0;
    for(int k=n-1;k>=0;k--)
    {
        for(int j=0;j<x.size();j++)
        {
            for(int i=0;i<v.size();i++)
            {
                J[k][i][j]=INT_MAX;
                for(double a=umin;a<=umax;a+=du)
                {
                    if(x[j]-0.5*dt*dt*a+v[i]*dt < 0 || v[i]+a*dt < dv || v[i]+a*dt>vmax || x[j]>71)continue;
                    newj=J[k+1][(int)((v[i]+a*dt)/dv)][(int)((x[j]-0.5*dt*dt*a-v[i]*dt)/dx)]*beta(k)+phi(v[i],a);
                    if(J[k][i][j]>newj) 
                    {
                        J[k][i][j]=newj;
                        back[k][i][j]=make_pair((v[i]+a*dt)/dv,(x[j]-0.5*dt*dt*a-v[i]*dt)/dx);
                    }
                }
            }
        }
    }
    cout<<J[0][18][560]<<endl;
    pair<int,int>p=make_pair(18,560);
    for(int i=0;i<n;i++)
    {
        cout<<v[p.first]<<"\n"<<x[p.second]<<"\n---\n";
        p=back[i][p.first][p.second];
    }
}