#include <iostream>
#include <math.h>

using namespace std;

const double gamma = 0;
const double w0 = 10;
const double A = 3.1;
const double h = 0.01;
const int N = 1025;

struct point {
	double x;
	double y;
};

point F(point& v) {
	point ans;

    ans.x = v.y;
    ans.y = -2 * gamma * v.y - w0 * w0 * sin(v.x);

    return ans;
}

double linear(double t) {
	return A * exp(-gamma * t) * cos(sqrt(w0 * w0 - gamma * gamma) * t);
}

void FFTAnalysis(double AVal[], double FTvl[]) {
	int n = (N - 1) * 2, m = 0;

	double Tmvl[2 * N];

	for (int i = 0; i < n; i += 2) {
		Tmvl[i] = 0;
		Tmvl[i + 1] = AVal[i / 2];
	}

	for (int i = 1, j = 1; i < n; i += 2, j += m) {
		if (j > i) {
			swap(Tmvl[i]    , Tmvl[j]);
            swap(Tmvl[i + 1], Tmvl[j + 1]);
        }

        m = (N - 1);
		while ((m >= 2) && (j > m)) {
			j -= m;
             m /= 2;
		}
	}

	for (int Mmax = 2; Mmax < n; Mmax *= 2) {
		double Theta = -M_PI * 2 / Mmax;
		double Wpr = sin(Theta / 2) * sin(Theta / 2) * 2;
		double Wr = 1, Wi = 0; m = 1;

		while (m < Mmax) {
			m = m + 2;
            double Tmpr = Wr, Tmpi = Wi;
			Wr -= Tmpr * Wpr + Tmpi * sin(Theta);
			Wi += Tmpr * sin(Theta) - Tmpi * Wpr;

			for (int i = m - 2; i < n; i += Mmax * 2) {
				int temp = i + Mmax;
				Tmpr = Wr * Tmvl[temp] - Wi * Tmvl[temp - 1];
				Tmpi = Wi * Tmvl[temp] + Wr * Tmvl[temp - 1];

				Tmvl[temp] = Tmvl[i] - Tmpr; Tmvl[temp - 1] = Tmvl[i - 1] - Tmpi;
				Tmvl[i] = Tmvl[i] + Tmpr; Tmvl[i - 1] = Tmvl[i - 1] + Tmpi;
			}
		}
	}

	for (int i = 0; i < (N - 1); i++) {
        FTvl[i] = 2 * sqrt(pow(Tmvl[i * 2], 2) + pow(Tmvl[i * 2 + 1], 2)) / (N - 1);
	}
}

int main() {
	point V[N];

	V[0].x = A;
	V[0].y = 0;

    point k1, k2, k3, k4;
	for (int i = 0; i < N - 1; i++) {
        point temp;

		k1 = F(V[i]);

        temp.x = V[i].x + (h / 2) * k1.x;
        temp.y = V[i].y + (h / 2) * k1.y;

        k2 = F(temp);

		temp.x = V[i].x + (h / 2) * k2.x;
		temp.y = V[i].y + (h / 2) * k2.y;

        k3 = F(temp);

        temp.x = V[i].x + h * k3.x;
		temp.y = V[i].y + h * k3.y;

        k4 = F(temp);

		V[i + 1].x = V[i].x + h / 6 * (k1.x + 2 * k2.x + 2 * k3.x + k4.x);
		V[i + 1].y = V[i].y + h / 6 * (k1.y + 2 * k2.y + 2 * k3.y + k4.y);
	}

	double coords[N - 1], FFTans[N - 1];

	for (int i = 0; i < N - 1; i++) {
		coords[i] = V[i].x;
		FFTans[i] = 0;
	}

	FFTAnalysis(coords, FFTans);

	for (int i = 0; i < N - 1; i++) {
		cout << V[i].x << " " << linear(h * i) << " " << FFTans[i] << " " << h * i << '\n';
    }

	return 0;
}