求救： f(x,y)=f(x-1,y)+f(x+1,y)+f(x,y-1)+f(x,y+1)，求f(x,y).

[goodenfeet]

其中，
f(x,0)=0,f(500,y)=50,f(x,500)=75,f(0,y)=100.
然后，
0<=x,y<=500,x,y是正整数。
要求是，
输入x,y，返回f(x,y).


向各位大侠求救了，小弟想了好几天，实在是没有办法了，就贴到这里来了。这个用简单的函数递归实现不了，
会出现重复递归的问题。听同学说可能要用带记录的递归，但是挺迷惑的。还望大侠们指教，不胜感激！

[jsjacky1101]

郁闷!

[Rivest]

to mathe:
你的算法有点问题

[mathe]

呵呵，我已经计算出来了

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
typedef double ftype;
typedef ftype ** MATRIX;
typedef ftype * VECTOR;
typedef const ftype ** 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){
    MATRIX x = (MATRIX)malloc(sizeof(ftype *)*n+sizeof(ftype)*n*n);
    int i;
    x[0]=(ftype *)(x+n);
    for(i=1;i<n;i++)x[i]=x[i-1]+n;
    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);
}

void matrix_sum(MATRIX A, CONST_MATRIX B, int n){
    int i,j;
    for(i=0;i<n;i++){
	for(j=0;j<n;j++){
	    A[i][j] += B[i][j];
	}
    }
}

void matrix_diff(MATRIX A, CONST_MATRIX B, int n){
    int i,j;
    for(i=0;i<n;i++){
	for(j=0;j<n;j++){
	    A[i][j] -= B[i][j];
	}
    }
}

void vector_sum(VECTOR a, CONST_VECTOR b, int n){
    int i;
    for(i=0;i<n;i++)a[i]+=b[i];
}

void vector_diff(VECTOR a, CONST_VECTOR b, int n){
    int i;
    for(i=0;i<n;i++)a[i]-=b[i];
}

void vector_rdiff(VECTOR a, CONST_VECTOR b, int n){
    int i;
    for(i=0;i<n;i++)a[i]=b[i]-a[i];
}

void matrix_mul(MATRIX out, CONST_MATRIX in1, CONST_MATRIX in2, int n){
    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;
    }
}
//-H*in
void vector_mul_mH(VECTOR out, CONST_VECTOR in, int n){
    int i;
    for(i=0;i<n;i++){
        ftype sum = in[i];
        if(i>0)sum-=in[i-1];
        if(i<n-1)sum-=in[i+1];
        out[i]=sum;
    }
}
//-H*in
void matrix_mul_mH(MATRIX out, CONST_MATRIX in, int n){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<n;j++){
	ftype sum=in[i][j];
	if(i>0)sum-=in[i-1][j];
	if(i<n-1)sum-=in[i+1][j];
	out[i][j]=sum;
    }
}
//H*in
void matrix_mul_H(MATRIX out, CONST_MATRIX in, int n){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<n;j++){
	ftype sum=-in[i][j];
	if(i>0)sum+=in[i-1][j];
	if(i<n-1)sum+=in[i+1][j];
	out[i][j]=sum;
    }
}

void matrix_mul_vector(VECTOR out, CONST_MATRIX m, CONST_VECTOR in, int n){
    int i,j;
    for(i=0;i<n;i++){
	ftype sum=0;
	for(j=0;j<n;j++){
	    sum+=m[i][j]*in[j];
	}
	out[i]=sum;
    }
}

void H_mul_vector(VECTOR out, CONST_VECTOR in, int n){
    int i;
    for(i=0;i<n;i++){
	ftype sum=-in[i];
	if(i>0)sum+=in[i-1];
	if(i<n-1)sum+=in[i+1];
        out[i]=sum;
    }
}

void vector_init_const(VECTOR v, ftype c, int n){
    int i;
    for(i=0;i<n;i++)v[i]=c;
}

void matrix_init_O(MATRIX o, int n){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<n;j++)o[i][j]=0.0;
}

void matrix_init_const(MATRIX m, ftype c, int n){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<n;j++)m[i][j]=c;
}

void matrix_init_E(MATRIX e, int n){
   int j;
   matrix_init_O(e,n);
   for(j=0;j<n;j++)e[j][j]=1.0;
}

void matrix_init_H(MATRIX h, int n){
   int i;
   matrix_init_O(h,n);
   for(i=0;i<n;i++){
       if(i>=1)h[i][i-1]=1.0;
       h[i][i]=-1.0;
       if(i<n-1)h[i][i+1]=1.0;
   }
}

void matrix_copy(MATRIX m, CONST_MATRIX k, int n){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<n;j++)m[i][j]=k[i][j];
}

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 i,j;
    for(i=0;i<n;i++){
	for(j=0;j<n;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;
VECTOR cn, sn;
MATRIX T,F;
void fini()
{
    vector_free(cn);
    vector_free(sn);
    matrix_free(T);
    matrix_free(F);
}

void init()
{
    int i,j;
    ftype mul=sqrt(2.0/(N+1));
    cn=vector_alloc(N);
    sn=vector_alloc(N);
    T=matrix_alloc(N);
    F=matrix_alloc(N+1);
    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<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){
            T[i][j]=0;
        }else{
            T[i][j]=sn[index-1]*mul*is_minus;
        }
    }
    for(i=0;i<N;i++){///Calculate F[j,i+1]
        double x=-2.0*cn[i]+1.0;
        F[0][i]=1.0;
        F[1][i]=x;
        for(j=2;j<=N;j++){
            F[j][i]=x*F[j-1][i]-F[j-2][i];
        }
    }
}

void DST(VECTOR out, const VECTOR in){
    matrix_mul_vector(out, T, in, N);///It could be improved to O(nlogn) by using Fast Fourier Transform
}

void Solve(MATRIX x,const VECTOR init){///Init holding initial value of r(n)
    int i,j;
    int failed=0;
    VECTOR R;
    MATRIX y;
    R=vector_alloc(N);
    y=matrix_alloc(N);
    DST(R,init);
    for(i=0;i<N;i++)R[i]=-R[i];
    for(i=0;i<N;i++){
        double v=F[N][i];
        if(fabs(v)<ERROR){
            if(fabs(R[i])>=ERROR){
                failed|=1;
            }
            failed|=2;
            y[0][i]=0.0;
        }else{
            y[0][i]=-R[i]/v;
        }
        for(j=1;j<N;j++){
            y[j][i]=y[0][i]*F[j][i];
        }
    }
    for(i=0;i<N;i++){
        DST(x[i],y[i]);
    }
    if(failed){
        if(failed&1){
            fprintf(stderr,"WARNING: No Valid Solution\n");
        }else{
            fprintf(stderr,"WARNING: Solution not unique\n");
        }
    }
    matrix_free(y);
    vector_free(R);
}

void Usage(char *name){
    printf("Usage:\n\t%s n\n\t\tWhere n is positive integer\n"
              "Or\n\t%s d1 d2 ... dn\n\t\tWhere di is a floating pointer number\n",
              name, name);
    exit(-1);
}

void dump_diff(const MATRIX x, const VECTOR r)
{
    int i;
    VECTOR d=vector_alloc(N);
    for(i=0;i<N;i++){
        ftype sum=0.0;
        sum=x[N-1][i]-x[N-2][i];
        if(i>=1)sum-=x[N-1][i-1];
        if(i<N-1)sum-=x[N-1][i+1];
        d[i]=sum;
    }
    printf("Required Boundary:\n");
    vector_output(r,N);
    printf("Result Boundary:\n");
    vector_output(d,N);
    vector_diff(d,r,N);
    printf("Diff:\n");
    vector_output(d,N);
    vector_free(d);
}

int main(int argc, char *argv[]){
    int i;
    MATRIX x;
    VECTOR r;
    if(argc<=1){
        Usage(argv[0]);
    }else if(argc==2){
        N=atoi(argv[1]);
        if(N<=0){
            Usage(argv[0]);
        }
        if(N==1){
            printf("1\n");
            return 0;
        }
    }else{
        N=argc-1;
    }
    
    init();
    x=matrix_alloc(N);
    r=vector_alloc(N);
    if(argc==2){
            for(i=0;i<N;i++)r[i]=1.0;
    }else{
        for(i=0;i<N;i++){
            r[i]=atof(argv[i+1]);
        }
    }
    Solve(x,r);
    matrix_output(x,N);
    dump_diff(x,r);
    vector_free(r);
    matrix_free(x);
    fini();
    return 0;
}

[mathe]

上面程序由于没有使用快速离散正弦变换，时间复杂度为O(n^3)

[mathe]

上面程序都假设初始数据三个边界全为0
有两种不同输入数据
i)一个整数作为命令行参数
计算出另外一条边界全部是1的情况。
比如参数为499就对应本题。
ii)一串浮点数作为参数表示另外一条边界的情况。
比如 1.0 1.0 2.0作为参数，就是一个3*3情况，3条边界全部为0，另外一条边界为1.0 1.0 2.0
输出:
首先是一个n*n的矩阵，表示各个格子应该填入的数据
然后输出用户要求的边界
在输出计算结果要求的边界和误差
程序还会判断是否有解或者解是否唯一。如果无解或解不唯一，会给出一条警告信息。

[mathe]

例子：
$ ./bulb3 1 1 -1 -1
WARNING: Solution not unique
2 1.25 -1.25 -2
0.75 0.5 -0.5 -0.75
-1.75 -1 1 1.75
-1.5 -0.75 0.75 1.5
Required Boundary:
1 1 -1 -1
Result Boundary:
1 1 -1 -1
Diff:
0 -1.11022e-16 1.11022e-16 1.11022e-16

可以看出，输出的4*4矩阵中，除了最后一行以外，其他各个数据都是边上数据之和。而最后一行，则还需要加上边界条件中的对应数据。

[goodenfeet]

信息量太大了，我慢慢看，待会打分。我现在就想说一句话：mathe,偶像！

[dchilli]

mark

[mathe]

昨天晚上分析了一下，这个题目竟然可以用O(n^2log(n))的算法达到，
而其中离散正弦变换竟然可以得到应用。也许，这个题目一开始就应该用傅立叶变换来分析。
还是采用我前面的记号，分别用a,b,c,d代替四条边界上的取值（其实每条边界上不需要要求所有点的初始值相同）。
首先我们可以发现，这个题目满足可加性。
我们可以先计算出(a,b,c,d)=(a,0,0,0)的解，然后计算(a,b,c,d)=(0,b,0,0)的解，然后再计算另外两个只有一条边界为0的解。算出来后，将4个解对应位置相加，就是所求的解了。所以我们可以先只分析(0,0,0,d)这组情况。
记函数f(1,x)=1,f(2,x)=x, f(n+1,x)=x*f(n,x)-f(n-1,x)
那么f(n,2x)就是第二类切皮雪夫多项式.
记矩阵T为一个n阶对称阵，其中
T[i,j]=sqrt(2/n)*sin((i*j*PI)/(n+1)) (1<=i<=n,1<=j<=n)
可以知道T[i,j]是正交阵
而且对于任意一个长度为n的向量X, T*X其实就是X的离散正弦变换，而计算离散正弦变换可以有O(n*log(n))的快速算法。

现在回到题目,n=499.
由于a=b=c=0,所以题目中r1,r2,...,r(n-1)都是0，只有r(n)=(-d,....,-d)不等于0. (其实对于这个边界条件，可以直接求解不同的位置不同的情况，只要修改r(n)的值就可以了）
所以得到上面递推式中u1=u2=...=u(n)=0, u(n+1)=r(n),而V(n+1)=f(n,-B),所以我们得到
f(n,-B)x1=-r(n)
而
TBT=diag{2cos(1*Pi/(n+1))-1,2cos(2*Pi/(n+1))-1,...,2cos(n*Pi/(n+1))-1}
所以
如果假设Txk=yk, Tr(n)=-R(n)
那么我们得到
f(n,(-2cos(h*Pi/(n+1))+1))*y1(h)=-R(n)(h)
其中y1(h)和R(n)(h)分别表示这两个向量的第h项。
所以我们可以
i)首先计算r(n)的离散正弦变换R(n). (如果所有d相同，我们可以发现，R(n)的所有奇数项都是0）
ii)计算f(n, -2cos(h*Pi/(n+1))+1) (h=1,2,...,n)
这个可以采用f(n,.)的迭代式计算，计算每个时间复杂度为O(n)，总共O(n^2)
iii)解出y1,如果对于某个h,对应f(n,-2cos(h*Pi/(n+1))+1)=0, 而R(h)!=0,那么方程无解。
而如果对于f(n,-2cos(h*Pi/(n+1))+1)=0的h都有，R(h)=0,那么方程有无穷个解，这是y(1,h)可以任意取.
而对于f(n,-2cos(h*Pi/(n+1))+1)=0的情况，不论那种情况，取y(1,h)=0都是x1的最小二乘解。
而在我们的题目中，d这条边界对应所有数据相同，所以对于h是偶数的情况R(n)才为0
而对于n=499,我们发现h=100,200,300,400的情况f(n,-2cos(h*Pi/(n+1))+1)=0. 所以我们知道，对于这个题目，有无穷的，但是我们可以取y(1,100)=y(1,200)=y(1,300)=y(1,400)=0,求出最小二乘解。
iv)通过再次离散正弦变换（逆变换)计算出x1
上面部分计算x1总共花费了O(n^2)的时间。
此后，我们就可以通过方程
x(k)=f(k-1,-B)x1来计算x(k)
同样做离散正弦变换，上面等式变成
y(k,h)=f(k-1,-2cos(h*Pi/(n+1))+1) y(1,h)
而其中f(k-1,-2cos(h*Pi/(n+1))+1)在ii)中已经产生，直接使用就可以解出yk
所以最后我们只要在通过n-1次离散正弦变换(逆变换）可以计算出x2,x3,...,xn
这一步需要花费O(n*nlog(n))的时间
所以总共时间复杂度O(n^2log(n)).

[mathe]

上面分析计算出边界条件是(0,0,0,d)的情况
而对于另外三种边界条件，我们不需要重新计算，只要将计算结果旋转一下乘上一个不同的系数就可以了。
最后将计算结果累加就可以了。

>>>>>>>>>>>>加分要选择“管理帖子按钮”
发表于：2007-11-08 17:36:3720楼 得分:0
怎么评分啊，我怎么似乎没有权限给大家加分呢。。。

[jason1019]

mathe是数学牛人啊，佩服。

[mathe]

改用离散正弦变换除了速度快以外，还有一个好处应该就是解很稳定，我们总可以求得一个最小二乘解，不会象前面解方程一样，计算误差会引起解的很大变换，而且还有象这个题目，方程有无穷组解，解方程的方法也处理不了。
而计算精度为什么会引起这么大的问题，其实很好理解
比如现在我们计算两个相差不大的数的差值，比如只保留4位有效数字。
比如计算(2/5-3/7)*7的差
2/5=0.4000
3/7=0.4286
2/5-3/7=-0.0286
(2/5-3/7)*7=-0.2002
所以这里同理论值－0.2差了0.1%,
也就是精度是万分之一的计算经过了一步两个比较接近的数的减法后误差达到千分之一了。
而我们的题目中，会经常有数百个绝对值相差不大的数做加减法，误差就非常大了。
比如上面题目中，理论上T^2是单位阵TBT是对角阵
对于这个结论，我写过代码来检验，结果发现T^2总是非常接近单位阵
但是TBT同理论值的误差比较大，在N=10时就达到0.178
而在N=500时越为0.004.虽然这个数值看起来不大，但是考虑到T中每个元素平均数值也就在1/sqrt(500)=0.044左右，所以这个相对误差也在10％左右，是非常大的误差了。

[kojie_chen]

牛人啊,我都没信心了

[feixue_XXXX]

这个题目有问题啊：
从f(x,0)=0;
可知 当y=0 时 x取任何值此函数为 0
但 f(500,y)=50 ！=0

在从f(0,y)=100
可知 当x=0 时 不管y 取任何值都为100
但 f(x,500)75 !=100
所以这个题目本身就是错的！

[flyingscv]

满足f(x,y)=f(x-1,y)+f(x+1,y)+f(x,y-1)+f(x,y+1)的N*N方阵
可以列出N*N条N*N元1次方程式，因为每个都不能用替他带入,所以只有一组解，就是0解
也就是说方阵中任意f(x,y)<>0，则方阵无解

[flyingscv]

忽律了f(x,y)可以是负数
这样到是可能有解的

解501*501+4个501*501元一次方程式，也只有用计算机来做了

大概需要用链表实现方程式替换，复杂度<o(n2)[n=501*501],等周末试试有没有解

[kendan]

study

[jixingzhong]

f(x,y)=f(x-1,y)+f(x+1,y)+f(x,y-1)+f(x,y+1)

题目不是给定了么，
然后给定几个映射：
f(x,0)=0,f(500,y)=50,f(x,500)=75,f(0,y)=100.

那么就递归好了么 ...

信号分析学的不好，
从程序设计的角度来看，
似乎没那么复杂

[xkw365]

楼上的，五颗星，羡慕啊

[mathe]

其实这道题目应该放在算法版讨论。
再算法版里面还有一道关于给数扫雷的算法，算法和这个很相似，只是在那里，我分析的没有这么彻底。
此外，过去我还分析过一个叫关灯游戏的问题，也很类似，只是那个题目是在二阶域上的。本来一位关灯游戏是最简单的，因为二阶域的计算比较简单，没有想到反而是实数域中问题更加简单，可以转化为离散正弦变化问题，简化到O(n^2log(n))的算法。而二阶域就比较难了。