#include <math.h>
#include <stdio.h>
#include <stdlib.h>
typedef double ftype;
typedef ftype ** MATRIX;
typedef ftype * VECTOR;
typedef const ftype *const*const CONST_MATRIX;
typedef const ftype * CONST_VECTOR;
#define INIT_a 0.0
#define INIT_b 100.0
#define INIT_c 75.0
#define INIT_d 50.0
#define ERROR 0.000001
MATRIX matrix_alloc(int n, int m){
MATRIX x = (MATRIX)malloc(sizeof(ftype *)*n+sizeof(ftype)*n*m);
int i;
x[0]=(ftype *)(x+n);
for(i=1;i<n;i++)x[i]=x[i-1]+m;
return x;
}
void matrix_free(MATRIX x){
free(x);
}
VECTOR vector_alloc(int n){
VECTOR x = (VECTOR)malloc(sizeof(ftype)*n);
return x;
}
void vector_free(VECTOR x){
free(x);
}
//A+=B; n rows m columns
void matrix_sum(MATRIX A, CONST_MATRIX B, int n, int m){
int i,j;
for(i=0;i<n;i++){
for(j=0;j<m;j++){
A[i][j] += B[i][j];
}
}
}
//A-=B; n rows m columns
void matrix_diff(MATRIX A, CONST_MATRIX B, int n, int m){
int i,j;
for(i=0;i<n;i++){
for(j=0;j<m;j++){
A[i][j] -= B[i][j];
}
}
}
//a+=b;
void vector_sum(VECTOR a, CONST_VECTOR b, int n){
int i;
for(i=0;i<n;i++)a[i]+=b[i];
}
//a-=b;
void vector_diff(VECTOR a, CONST_VECTOR b, int n){
int i;
for(i=0;i<n;i++)a[i]-=b[i];
}
//a = a-b;
void vector_rdiff(VECTOR a, CONST_VECTOR b, int n){
int i;
for(i=0;i<n;i++)a[i]=b[i]-a[i];
}
//out = in1*in2; the size of in1 is n*c, the size of in2 is c*m
void matrix_mul(MATRIX out, CONST_MATRIX in1, CONST_MATRIX in2, int n, int c, int m){
int i,j,k;
for(i=0;i<n;i++)for(j=0;j<n;j++){
ftype sum=0.0;
for(k=0;k<n;k++){
sum+=in1[i][k]*in2[k][j];
}
out[i][j]=sum;
}
}
//out = M*in, where the size of 'M' is n*m and the size of 'in' is m
void matrix_mul_vector(VECTOR out, CONST_MATRIX M, CONST_VECTOR in, int n, int m){
int i,j;
for(i=0;i<n;i++){
ftype sum=0;
for(j=0;j<m;j++){
sum+=M[i][j]*in[j];
}
out[i]=sum;
}
}
//v[i]=c for all i
void vector_init_const(VECTOR v, ftype c, int n){
int i;
for(i=0;i<n;i++)v[i]=c;
}
//O[i][j]=0 for all (i,j)
void matrix_init_O(MATRIX O, int n, int m){
int i,j;
for(i=0;i<n;i++)for(j=0;j<m;j++)O[i][j]=0.0;
}
//extract the k-th column of matrix M (n rows) to vector out
void extract_matrix_column(VECTOR out, CONST_MATRIX M, int n, int k)
{
int i;
for(i=0;i<n;i++){
out[i]=M[i][k];
}
}
//Change the k-th column of matrix M (n rows) to vector in
void update_matrix_column(MATRIX M, const VECTOR in, int n, int k)
{
int i;
for(i=0;i<n;i++){
M[i][k]=in[i];
}
}
//M[i][j]=c for all (i,j)
void matrix_init_const(MATRIX M, ftype c, int n,int m){
int i,j;
for(i=0;i<n;i++)for(j=0;j<m;j++)M[i][j]=c;
}
//Square matrix E =diag{1,1,...,1}
void matrix_init_E(MATRIX E, int n){
int j;
matrix_init_O(E,n,n);
for(j=0;j<n;j++)E[j][j]=1.0;
}
//M=K; the size of matrix is n*m
void matrix_copy(MATRIX M, CONST_MATRIX K, int n,int m){
int i,j;
for(i=0;i<n;i++)for(j=0;j<m;j++)M[i][j]=K[i][j];
}
//c=w; for size of vector is n
void vector_copy(VECTOR v, CONST_VECTOR w, int n){
int i;
for(i=0;i<n;i++)v[i]=w[i];
}
void matrix_output(CONST_MATRIX M, int n,int m){
int i,j;
for(i=0;i<n;i++){
for(j=0;j<m;j++){
printf("%g\t",M[i][j]);
}
printf("\n");
}
}
void vector_output(VECTOR v, int n){
int i;
for(i=0;i<n;i++)printf("%g\t",v[i]);
printf("\n");
}
int N, M;
VECTOR cn, sn;///Array to hold cos(k*Pi/(N+1)), sin(k*Pi/(N+1)), 1<=k<=N
VECTOR cm, sm;///Array to hold cos(k*Pi/(M+1)), sin(k*Pi/(M+1)), 1<=k<=M
MATRIX TN,TM;
void fini()
{
vector_free(cn);
vector_free(sn);
vector_free(cm);
vector_free(sm);
matrix_free(TN);
matrix_free(TM);
}
void init()
{
int i,j;
ftype mul=sqrt(2.0/(N+1));
ftype mul2=sqrt(2.0/(M+1));
cn=vector_alloc(N);
sn=vector_alloc(N);
cm=vector_alloc(M);
sm=vector_alloc(M);
TN=matrix_alloc(N,N);
TM=matrix_alloc(M,M);
for(i=0;i<N;i++){
cn[i]=cos((i+1)*M_PI/(ftype)(N+1));
sn[i]=sin((i+1)*M_PI/(ftype)(N+1));
}
for(i=0;i<M;i++){
cm[i]=cos((i+1)*M_PI/(ftype)(M+1));
sm[i]=sin((i+1)*M_PI/(ftype)(M+1));
}
for(i=0;i<N;i++)for(j=0;j<N;j++){
int index=(i+1)*(j+1);
int is_minus=1;
index%=2*(N+1);
if(index>N){
is_minus=-1;
index-=N+1;
}
if(index==0){
TN[i][j]=0;
}else{
TN[i][j]=sn[index-1]*mul*is_minus;
}
}
for(i=0;i<M;i++)for(j=0;j<M;j++){
int index=(i+1)*(j+1);
int is_minus=1;
index%=2*(M+1);
if(index>M){
is_minus=-1;
index-=M+1;
}
if(index==0){
TM[i][j]=0;
}else{
TM[i][j]=sm[index-1]*mul2*is_minus;
}
}
}
void TWODDST(MATRIX out, const MATRIX in)
{
int i;
VECTOR tmp1=vector_alloc(N);
VECTOR tmp2=vector_alloc(N);
for(i=0;i<N;i++){
matrix_mul_vector(out[i], TM, in[i], M, M);
}
for(i=0;i<M;i++){
extract_matrix_column(tmp1, out, N, i);
matrix_mul_vector(tmp2, TN, tmp1, N, N);
update_matrix_column(out, tmp2, N, i);
}
vector_free(tmp2);
vector_free(tmp1);
}
int zero_candidate(int x, int y){
if(x*2==N+1&&y*2==M+1)
return 1;
if(x*5==N+1&&y*5==3*(M+1))
return 1;
if(x*5==3*(N+1)&&y*5==M+1)
return 1;
return 0;
}
//Return nonzero when there's solution for equation.
int Solve(MATRIX x, const MATRIX init)
{
MATRIX tmp = matrix_alloc(N,M);
int i,j;
int failed=0;
int zero_count=0;
TWODDST(tmp, init);
for(i=0;i<N;i++)for(j=0;j<M;j++){
if(zero_candidate(i+1,j+1)){
if(fabs(tmp[i][j])>=ERROR){//No solution, this element must be 0.
failed=1;
}
tmp[i][j]=0.0;
zero_count++;
}else{
tmp[i][j]/=(2.0*cn[i]+1.0)*(2.0*cm[j]+1.0)-1.0;
}
}
if(failed)return 0;//No solution available;
TWODDST(x, tmp);
return 1;
}
int IsValidSolution(MATRIX x)
{
int i,j;
for(i=0;i<N;i++)for(j=0;j<M;j++){
if(fabs(x[i][j])>ERROR&&fabs(x[i][j]-1)>ERROR){
return 0;
}
}
return 1;
}
void GaussSolve(const MATRIX kmatrix, VECTOR r,int order)
{
int i,j,k;
MATRIX tmp=matrix_alloc(order,order);
for(i=0;i<order;i++)for(j=0;j<order;j++){
tmp[j][i]=kmatrix[i][j];
}
for(i=0;i<order;i++){
ftype f=tmp[i][i];
for(j=i;j<order;j++)
tmp[i][j]/=f;
r[i]/=f;
for(j=0;j<order;j++){
ftype div=tmp[j][i];
if(j==i)continue;
for(k=i;k<order;k++){
tmp[j][k]-=div*tmp[i][k];
}
r[j]-=div*r[i];
}
}
matrix_free(tmp);
}
void SolutionOutput(MATRIX x)
{
int i,j;
for(i=0;i<N;i++){
for(j=0;j<M;j++){
int v=(int)(x[i][j]+0.5);
printf("%d",v);
}
printf("\n");
}
}
