#include <cstdio>   // For printf, scanf, fopen, fclose, perror, getchar
#include <cstdlib>  // For exit
#include <cstring>  // For strcpy, strcmp
#include <cmath>    // For cos, sin, M_PI (from cmath)
#include <cfloat>   // For DBL_MAX
#include <cctype>   // For toupper

// Define MMAX as a constant for array sizing
#define MMAX_DEFINE 300 // Renamed to avoid conflict with potential variable mmax

// Global variables (consider encapsulating these in a class/struct for larger programs)
// C++ arrays are 0-indexed. If the original BASIC used 1-based indexing
// (e.g., X(1) to X(1024)), we need 1025 elements for X[0] to X[1024].
// For Y, in gs1420, it accesses y[i-1] which implies y[0] is the first usable index.
// If ND is the number of data points, and they are stored from x[1] to x[ND],
// then x and y need to be at least of size (ND + 1) for 1-based indexing.
// The MMAX for fpe, r, rr, rfpe also suggests 1-based indexing for these,
// so they need MMAX_DEFINE + 1 size.
double x[1025]; // Max ND is assumed to be 1024 based on original BASIC DIM X(1024)
double y[1025]; // Max ND is assumed to be 1024. Accessed as y[i-1], so needs to handle i=1 (y[0]) up to y[ND-1]
double fpe[MMAX_DEFINE + 1];
double r[MMAX_DEFINE + 1];
double rr[MMAX_DEFINE + 1];
double rfpe[MMAX_DEFINE + 1];

double rm, wnmin, wnmax, sum1, sum2, sum, sumn, sumd, wnint;
int isw_global; // Renamed to avoid conflict with local 'isw' in main
int nd, mmax_global, ns, minm; // mmax_global renamed for clarity
double dt, pmm, pm, fpemin, av;
double ci_const = 2 * M_PI; // Use M_PI from cmath for better precision

// Function prototypes
void gs1880();
void gs1800(int m);
void gs1340();
void gs1420();
void gs1960();
void gs2110();
int mem(); // mem returns int to indicate success/failure

//---
//## サブルーチン: gs1880
//`minm`に`mmax`を、`pmm`に`pm`を設定し、`rfpe`配列に`r`配列の値をコピーします。
void gs1880() {
    minm = mmax_global; // Use the global mmax
    pmm = pm;
    for (int i = 1; i <= mmax_global; i++) { // Loop up to mmax_global
        rfpe[i] = r[i];
    }
}

//---
//## サブルーチン: gs1800
//モデル次数`m`におけるFPE (Final Prediction Error) を計算し、
//それが最小値であれば`fpemin`、`minm`、`pmm`、`rfpe`を更新します。
void gs1800(int m) {
    fpe[m] = static_cast<double>(nd + m + 1) / static_cast<double>(nd - m - 1) * pm;
    if (fpe[m] > fpemin) return;
    fpemin = fpe[m];
    minm = m;
    pmm = pm;
    for (int i = 1; i <= m; i++) {
        rfpe[i] = r[i];
    }
}

//---
//## サブルーチン: gs1340 (データ入力)
//指定されたファイルからデータを読み込み、`x`配列に格納します。
void gs1340() {
    char nm[256];
    FILE *fp;
    double bbb, ccc; // These variables are read but not used, matching original BASIC.

    printf(" *** データファイルメイ    =");
    if (scanf("%255s", nm) != 1) { // Safely read string, limit length
        fprintf(stderr, "ファイル名の読み込みに失敗しました。\n");
        exit(1);
    }
    
    fp = fopen(nm, "r");
    if (fp == nullptr) { // Check for nullptr instead of NULL for C++ style
        perror(nm); // Prints system error message
        exit(1); // Exit with an error code
    }

    for (int i = 1; i <= nd; i++) {
        if (fscanf(fp, "%lf %lf %lf", &x[i], &bbb, &ccc) != 3) { // Expect 3 doubles
            fprintf(stderr, "警告: データファイルの行 %d の読み込みに失敗しました。不足している可能性があります。\n", i);
            // Optionally, handle error more strictly or fill with default
            x[i] = 0.0; // Set to a default value if read fails
        }
        printf("I=%d  X=%g\n", i, x[i]);
    }
    fclose(fp);
}

//---
//## サブルーチン: gs1420 (ARモデル推定 - Burg法)
//Burg法を用いて自己回帰係数を推定します。
void gs1420() {
    double z;
    
    // Calculate average and subtract it from data
    sum = 0.0;
    for (int i = 1; i <= nd; i++) { // Assuming x[1] to x[nd] contain data
        sum += x[i];
    }
    av = sum / static_cast<double>(nd);

    sum = 0.0;
    for (int i = 1; i <= nd; i++) { // Assuming x[1] to x[nd] contain data
        z = x[i] - av;
        x[i] = z; // x[i] now holds detrended data
        y[i-1] = z; // y[0] to y[nd-1] will hold detrended data, used in subsequent steps
    }
    
    // Correctly calculate sum of squares after detrending
    // The original BASIC had sum += z*z after Y(I-1)=Z
    // Here we need to re-sum for the correct pm
    double detrended_sum_sq = 0.0;
    for (int i = 1; i <= nd; i++) {
        detrended_sum_sq += x[i] * x[i];
    }
    pm = detrended_sum_sq / static_cast<double>(nd);

    fpemin = static_cast<double>(nd + 1) / static_cast<double>(nd - 1) * pm;
    fpe[0] = fpemin;

    for (int m = 1; m <= mmax_global; m++) {
        sumn = 0.0;
        sumd = 0.0;
        for (int i = 1; i <= nd - m; i++) {
            sumn += x[i] * y[i]; // Original C code used y[i] here, likely intending current values of y
            sumd += x[i] * x[i] + y[i] * y[i];
        }
        
        if (sumd == 0.0) { // Prevent division by zero
            rm = 0.0;
        } else {
            rm = -2.0 * sumn / sumd;
        }
        r[m] = rm;
        pm *= (1.0 - rm * rm);

        if (m > 1) { // Replaced `goto label1680` with a conditional
            for (int i = 1; i <= m - 1; i++) {
                r[i] = rr[i] + rm * rr[m - i];
            }
        }
        for (int i = 1; i <= m; i++) {
            rr[i] = r[i];
        }

        // Update forward and backward prediction errors (Burg's algorithm core)
        // This part is a common source of off-by-one errors or misinterpretation
        // in Burg's algorithm implementations translated from different languages.
        // The core idea is that X and Y are updated to become the forward and backward
        // prediction errors for the *next* order.
        double temp_xi; 
        for (int i = 1; i <= nd - m - 1; i++) {
            temp_xi = x[i]; // Store old x[i] before it's updated
            x[i] = x[i] + rm * y[i];
            y[i] = y[i+1] + rm * temp_xi; // Corrected: y[i+1] (old backward error) + rm * old x[i] (old forward error)
                                          // This aligns with standard Burg's method where b_i = b_{i+1} + k * f_i
        }

        if (isw_global == 0) gs1800(m);
    }
    if (isw_global == 1) gs1880();
}

//---
//## サブルーチン: gs1960 (パワースペクトル推定)
//推定されたAR係数と予測誤差分散を用いてパワースペクトルを計算します。
void gs1960() {
    static bool ci_initialized = false; // Flag to ensure ci_const * dt is done only once per `mem` call for the initial `ci_const`
    if (!ci_initialized) {
        ci_const *= dt; // Multiply ci_const by dt only once
        ci_initialized = true;
    }
    
    int i = 0;
    // Iterate from wnmin to wnmax with step wnint
    // Using a for loop with a small epsilon for float comparison to avoid precision issues
    for (double f = wnmin; f <= wnmax + (wnint / 2.0); f += wnint) { // Add a small epsilon
        sum1 = 1.0;
        sum2 = 0.0;
        for (int j = 1; j <= minm; j++) {
            sum1 += rfpe[j] * cos(ci_const * f * j);
            sum2 += rfpe[j] * sin(ci_const * f * j);
        }
        
        double denominator = (sum1 * sum1 + sum2 * sum2);
        if (denominator == 0.0) { // Prevent division by zero
            x[i] = DBL_MAX; // Assign DBL_MAX (from <cfloat>)
        } else {
            x[i] = pmm * dt / denominator;
        }
        i++;
        if (i >= 1025) { // Prevent out-of-bounds access for x array (0 to 1024)
            fprintf(stderr, "警告: スペクトル点数が配列サイズを超過しました。処理を中断します。\n");
            break;
        }
    }
}

//---
//## サブルーチン: gs2110 (データ出力)
//結果のスペクトルデータをファイルに保存します。
void gs2110() {
    char nm[256];
    FILE *fp;

    printf("   Out Data File Name: ");
    if (scanf("%255s", nm) != 1) {
        fprintf(stderr, "ファイル名の読み込みに失敗しました。\n");
        // Clear input buffer on error
        while (getchar() != '\n' && getchar() != EOF); 
        return; 
    }

    fp = fopen(nm, "w");
    if (fp == nullptr) {
        perror(nm);
        return; 
    }

    for (int i = 0; i < ns; i++) { // Assuming x[0] to x[ns-1] for spectrum data
        fprintf(fp, "%g\n", x[i]);
    }
    fclose(fp);
    printf("データがファイル '%s' に保存されました。\n", nm);
}

//---
//## メインルーチン: mem
//プログラムの主要な処理を制御します。
int mem() {
    char an_char; 
    char p_str[256], n_str[256]; 

    printf("      ***** MEM *****\n\n");

    printf(" *** サンプル テンスウ    =");
    if (scanf("%d", &nd) != 1 || nd <= 0) {
        fprintf(stderr, "無効なサンプル点数です。正の整数を入力してください。\n");
        return 0;
    }
    printf("nd:%d\n", nd);

    printf(" *** サンプル カンカク    =");
    if (scanf("%lf", &dt) != 1 || dt <= 0) {
        fprintf(stderr, "無効なサンプル間隔です。正の数を入力してください。\n");
        return 0;
    }
    printf("dt:%g\n", dt);

    while (true) { 
        printf(" モデル ジスウ ヲ アタエマスカ (Y/N)?\n");
        // Clear input buffer first to avoid reading leftover newline from previous scanf.
        // Then read the actual character.
        while (getchar() != '\n' && getchar() != EOF); 
        an_char = getchar(); 
        an_char = toupper(an_char); // Convert to uppercase for consistent comparison

        if (an_char == 'Y') {
            isw_global = 1;
            strcpy(p_str, " *** モデル ジスウ          =");
            break;
        } else if (an_char == 'N') {
            isw_global = 0;
            strcpy(p_str, " *** サイダイ モデル ジスウ    =");
            break;
        } else {
            printf("YまたはNを入力してください。\n");
        }
    }
    printf("%s", p_str);
    if (scanf("%d", &mmax_global) != 1 || mmax_global <= 0) {
        fprintf(stderr, "無効なモデル次数です。正の整数を入力してください。\n");
        return 0;
    }
    printf("mmax:%d\n", mmax_global);

    gs1340(); // データ入力

    ns = nd;
    wnmin = 0.0;
    wnmax = static_cast<double>(nd); 
    strcpy(n_str, " Data ");

    gs1420(); // Burg法による自己回帰係数と予測誤差の推定

    while (true) { 
        printf("   サイショウ シュウハスウ =");
        if (scanf("%lf", &wnmin) != 1) {
            fprintf(stderr, "無効な入力です。数値を入力してください。\n");
            // Clear input buffer on error
            while (getchar() != '\n' && getchar() != EOF); 
            continue;
        }
        printf("   サイダイ シュウハスウ =");
        if (scanf("%lf", &wnmax) != 1) {
            fprintf(stderr, "無効な入力です。数値を入力してください。\n");
            while (getchar() != '\n' && getchar() != EOF);
            continue;
        }
        if (wnmax <= wnmin) {
            fprintf(stderr, "最大周波数は最小周波数より大きい必要があります。\n");
            continue;
        }
        break; 
    }
    printf("wnmin,max:%g , %g\n", wnmin, wnmax);

    printf("   スペクトル テンスウ  =");
    if (scanf("%d", &ns) != 1 || ns <= 0) {
        fprintf(stderr, "無効なスペクトル点数です。正の整数を入力してください。\n");
        return 0;
    }
    
    // Check if ns exceeds the pre-allocated x array size (1025 for 0-1024).
    if (ns > 1025) {
        fprintf(stderr, "警告: スペクトル点数 (%d) が許容範囲 (1025) を超えています。プログラムを終了します。\n", ns);
        return 0;
    }

    wnint = (wnmax - wnmin) / static_cast<double>(ns);
    printf("ns:%d    wnint:%g\n", ns, wnint);
    
    gs1960(); // スペクトルの計算

    strcpy(n_str, " MEM Spectrum ");
    gs2110(); // スペクトル推定結果のグラフ (ファイル出力のみ)

    while (true) { 
        printf("   モウイチド シマスカ (Y/N)? ");
        while (getchar() != '\n' && getchar() != EOF); 
        an_char = getchar(); 
        an_char = toupper(an_char);

        if (an_char == 'Y') {
            printf("\n"); 
            
            while (true) {
                printf("   サイショウ シュウハスウ =");
                if (scanf("%lf", &wnmin) != 1) {
                    fprintf(stderr, "無効な入力です。数値を入力してください。\n");
                    while (getchar() != '\n' && getchar() != EOF);
                    continue;
                }
                printf("   サイダイ シュウハスウ =");
                if (scanf("%lf", &wnmax) != 1) {
                    fprintf(stderr, "無効な入力です。数値を入力してください。\n");
                    while (getchar() != '\n' && getchar() != EOF);
                    continue;
                }
                if (wnmax <= wnmin) {
                    fprintf(stderr, "最大周波数は最小周波数より大きい必要があります。\n");
                    continue;
                }
                break;
            }
            printf("wnmin,max:%g , %g\n", wnmin, wnmax);

            printf("   スペクトル テンスウ  =");
            if (scanf("%d", &ns) != 1 || ns <= 0) {
                fprintf(stderr, "無効なスペクトル点数です。正の整数を入力してください。\n");
                return 0;
            }
            if (ns > 1025) {
                fprintf(stderr, "警告: スペクトル点数 (%d) が許容範囲 (1025) を超えています。プログラムを終了します。\n", ns);
                return 0;
            }

            wnint = (wnmax - wnmin) / static_cast<double>(ns);
            printf("ns:%d    wnint:%g\n", ns, wnint);
            gs1960();
            strcpy(n_str, " MEM Spectrum ");
            gs2110();
            continue; 
        } else if (an_char == 'N') {
            break; 
        } else {
            printf("YまたはNを入力してください。\n");
        }
    }
    return 1; 
}

//---
//## メイン関数
//プログラムのエントリポイントです。
int main() { 
    if (mem() == 1) {
        printf("プログラムが正常に終了しました。\n");
        return 0; 
    } else {
        printf("プログラムがエラーで終了しました。\n");
        return 1; 
    }
}