// convolve a gaussian with the non-first columns of data

#include <fstream>
#include <sstream>
#include <iostream>
#include <vector>
#include <cmath>

using namespace std;


int main (int argc, char *argv[])
{
    if (argc != 3) {
        cerr << "usage: convolve <sigma> <filename>\n";
        return 1;
    }
    double sigma = atof(argv[1]);
    const char* fname = argv[2];

    ifstream fstr(fname);
    if (!fstr) {
        cerr << "Failed to open " << fname << endl;
        return 2;
    }
    
    vector<vector<double> > column;

    char buf[1024];
    string s;

    int line;
    for (line=0; ! fstr.eof(); ++line) {
        fstr.get(buf,1024,'\n');
        char dummy;
        if (! fstr.get(dummy)) break; // get newline or eof

        if (!line) {
            istringstream s(buf);
            double x;
            while (s>>x) { column.push_back(vector<double>()); }
        }

        istringstream iss(buf);
        double tmp;
        for (int col = 0; ! iss.eof(); ++col) {
            iss>>tmp;
            column.reserve(col+1);
            column[col].push_back(tmp);
        }
    }

    int ncols = column.size();
    int nrows = column[0].size();

    const double sqrt2pi = sqrt(2.0*3.14159);

    vector<vector<double> > newcol;
    for (int col=1; col<ncols; ++col) {
        vector<double> c;
        c.push_back(0.0);       // j=0 not allowed
        for (int j=1; j<nrows; ++j) {
            double c_j=0;

            for (int i=0; i<nrows; ++i) {
                double e = (1.0-(1.0*i)/(1.0*j))/sigma;
                c_j += column[col][i]*exp(-e*e/2.0);
            }
            c_j /= sqrt2pi*sigma*j;
            
            c.push_back(c_j);
        }
        newcol.push_back(c);
    }

    for (int row=0; row<nrows;++row) {
        cout << column[0][row];
        for (int col=0; col<ncols-1; ++col) 
            cout << " " << newcol[col][row];
        cout << endl;
    }

    return 0;
}