#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define real double
#define RND ((real)rand()/RAND_MAX)

int qrdec(real**,real**,int,int);
int qrback(real**,real**,int,int,real*,real*);

real funs(int i, real x){
	if(i==0) return 1;
	if(i==1) return x;
	else {fprintf(stderr,"funs: wrong parameter"); return 0;}
}

int main(int argc, char** argv){
	int n=10,m=2;
	if(argc>1) {n=atoi(argv[1]); m=n;}
	if(argc>2) m=atoi(argv[2]);
	if(m>n) m=n;

	real* x = (real *)calloc(n,sizeof(real)); 
	real* y = (real *)calloc(n,sizeof(real)); 
	real* s = (real *)calloc(n,sizeof(real)); 
	for(int i=0;i<n;i++) x[i]=i+0.5;
	for(int i=0;i<n;i++) y[i]=i+cos(2.0*i);
	for(int i=0;i<n;i++) s[i]=1+pow(sin(2.0*i),2);

	real** a = (real **)calloc(m,sizeof(real *)); 
	for (int i=0; i<m ; i++) a[i] = (real *)calloc(n,sizeof(real)); 

	for(int row=0;row<n;row++) for(int col=0;col<m;col++)
		a[col][row]=funs(col,x[row])/s[row];

	real* b = (real *)calloc(n,sizeof(real)); 
	for(int i=0;i<n;i++) b[i]=y[i]/s[i];

	real** r = (real **)malloc(m*sizeof(real *)); 
	for (int i=0; i<m ; i++) r[i] = (real *)malloc(n*sizeof(real)); 

	qrdec(a,r,n,m);

	real* sol = (real *)calloc(n,sizeof(real)); 
	qrback(a,r,n,m,b,sol);

	for(int i=0;i<n;i++){
		printf("# m=1,S=28\n");
		printf("%g %g\n",x[i],y[i]-s[i]);
		printf("%g %g\n",x[i],y[i]);
		printf("%g %g\n",x[i],y[i]+s[i]);
		}

	printf("# m=1,S=0\n");
	for(int i=0;i<n;i++){
		printf("%g %g\n",x[i],sol[0]+x[i]*sol[1]);
		}

}