int sgn(double pa)
{
	if (pa>0.0)
	{
		return 1;
	}
	else
	{
		if (pa<0.0)
		{
			return -1;
		}
	}
	return 0;
}

void qrdcmp(double a[], int m, int n, double q[])
{
	int i,j,k;
	double s,t,h,f;
    for (i = 1; i<=m; i++)
	{
        for (j = 1; j<=m; j++)
		{
            q[(i-1)*n+j] = 0.0;
        }
        q[(i-1)*n+i] = 1.0;
    }
    for (k = 1; k<=m-1; k++)
	{
        s = 0.0;
        for (i = k; i<=m;i++)
		{
            s = s + fabs(a[(i-1)*n+k]);
        }
        if (s != 0.0)
		{
            t = 0.0;
            for (i = k; i<=m; i++)
			{
                a[(i-1)*n+k] = a[(i-1)*n+k] / s;
                t = t + a[(i-1)*n+k] * a[(i-1)*n+k];
            }
            t = -sqrt(t) * sgn(a[(k-1)*n+k]);
            a[(k-1)*n+k] = a[(k-1)*n+k] - t;
            h = -t * a[(k-1)*n+k];
            for (j = k + 1; j<=n; j++)
			{
                f = 0.0;
                for (i = k; i<=m; i++)
				{
                    f = f + a[(i-1)*n+k] * a[(i-1)*n+j];
                }
                f = f / h;
                for (i = k; i<=m; i++)
				{
                    a[(i-1)*n+j] = a[(i-1)*n+j] - a[(i-1)*n+k] * f;
                }
            }
            for (j = 1; j<=m;j++)
			{
                f = 0.0;
                for (i = k; i<=m; i++)
				{
                    f = f + a[(i-1)*n+k] * q[(i-1)*n+j];
                }
                f = f / h;
                for (i = k; i<=m; i++)
				{
                    q[(i-1)*n+j] = q[(i-1)*n+j] - a[(i-1)*n+k] * f;
                }
            }
            a[(k-1)*n+k] = t * s;
            for (i = k + 1; i<=m; i++)
			{
                a[(i-1)*n+k] = 0.0;
            }
        }
    }
}