源程序如下
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <conio.h>
#include <iostream.h>

double white()                                                
{
    static double    w0 = 2053.0;
    static double    s0 = 65536.0;
    static double    v0 = 13849.0;
    static double    r0 = 0.0;
    double            t0, EGx;
    int                m0, i;

    t0 = 0.0;
      
    for(i = 0; i < 12; i++)
    {
        r0 = w0 * r0 + v0;
        m0 = (int)(r0 / s0);
        r0 = r0 - m0 * s0;
        t0 = t0 + r0 /  s0;
    }
    EGx = t0 - 6.0;
    return(EGx);
}

void main()
{
    const int m=21,n=21;
    const double af=0.2;
    const int nn = (n-1)*(n-1)+1;
    const int mm = (m-1)*(m-1)+1;
    int i,j,i0,j0,xl,yl;
    int x[n],y[m];
    int cm,bz2,bm,dz2,dm;
    int ii,jj;
    double xv,yv;
    double cz,bz1,dz1;
    double dn[n-1][m-1];
    double z[n][m],xx[nn][mm],yy[nn][mm],zz[nn][mm];
    double a[n-1],b[n-1],c[n-1],d[n-1];
    double cc[n-1][m-1],dd[n-1][m-1],bb[n-1][m-1],kk[n-1][m-1];
    
    for(i = 0; i < n; i++)
        for(j = 0; j<m; j++)
        {
            z[i][j]=white();
        }
    for(i=0; i < n; i++)
    {
        x[i]=i;
        y[i]=i;
    }
    xl=x[n-1] - x[0];
    yl=y[n-1] - y[0];
    for(i=0;i<n-1;i++)
    {
        a[i]=(x[i+1]-x[i])*xl;
        b[i]=(x[n-1]*x[i]-x[0]*x[i+1])*xl;
        c[i]=(y[i+1]-y[i])*yl;
        d[i]=(y[n-1]*y[i]-y[0]*y[i+1])*yl;
    }
    cz=z[0][0]+z[m-1][n-1]-z[0][n-1]-z[m-1][0];
    cm=x[0]*y[0]+x[n-1]*y[m-1]-x[n-1]*y[0]-x[0]*y[m-1];
    bz1=z[0,0]-z[0,n-1];
    bz2=x[0]*y[0]-x[n-1]*y[0];
    bm=x[0]-x[n-1];
    dz1=z[0,0]-z[m-1,0];
    dz2=x[0]*y[0]-x[0]*y[m-1];
    dm=y[0]-y[m-1];
    for(i = 0; i < n-1; i++)
        for(j = 0; j < m-1; j++)
        {
            dn[i][j]= af;
            cc[i][j]=(z[i][j]-z[i][j+1]-z[i+1][j]+z[i+1][j+1]-dn[i][j]*cz) * (double)(1/cm);
            bb[i][j]=(z[i][j]-z[i][j+1]-dn[i][j]*bz1-cc[i][j]*bz2) * (double)(1/bm);
            dd[i][j]=(z[i][j]-z[i+1][j]-dn[i][j]*dz1-cc[i][j]*dz2) * (double)(1/dm);
            kk[i][j]=z[j+1][i+1]-bb[i][j]*x[n-1]-dd[i][j]*y[m-1]-dn[i][j]*z[m-1][n-1]-cc[i][j]*x[n-1]*y[m-1];
        }
    for(j=0;j<m-1;j++)
        for(j0=0;j0<m;j++)
        {
            yv=c[j]*y[j0]+d[j];
            jj=(j-1)*(m-1)+j0;
            for(i=0;i<n-1;i++)
                for(i0=0;i0<n;i++)
                {
                    ii=(i-1)*(n-1)+i0;
                    xv=a[i]*x[i0]+b[i];
                    //z方向压缩变换
                    zz[jj][ii]=bb[j][i]*x[i0]+dd[j][i]*y[j0]+cc[j][i]*x[i0]*y[j0]+dn[j][i]*z[j0][i0]+kk[j][i]+100;
                }
        }

        for(i=0;i<nn;i++)
            for(j=0;j<mm;j++)
            {
                xx[i][j] = 0.00125 * j + 38;
                yy[i][j] = 0.00125 * i + 120;
            }
            cout << endl << "Done" << endl;
}