C++ converted to Rust, compiled ok but got a bad result. Which step is wrong?

My C++ code is here:

//最小二乘曲面拟合.cpp
  #include <iostream>
  #include <cmath>
  using namespace std;
//x[n]       存放给定数据点的X坐标。
//y[n]       存放给定数据点的Y坐标。
//z[n][m]    存放给定n*m个网点上函数值。
//n          X坐标个数。
//m          Y坐标个数。
//a[p][q]    返回二元拟合多项式的系数。
//p          拟合多项式中x的最高次为p-1。要求p<=min(n,20)。
//q          拟合多项式中y的最高次为q-1。要求q<=min(n,20)。
//dt[3]      分别返回误差平方和,误差绝对值之和与误差绝对值的最大值。
  void pir2(double x[], double y[], double z[], int n, int m,
              double a[], int p, int q, double dt[])
  {
      int i,j,k,l,kk;
      double apx[20],apy[20],bx[20],by[20],u[20][20];
      double t[20],t1[20],t2[20],d1,d2,g,g1,g2;
      double x2,dd,y1,x1,*v;
      v = new double[20*m];
      for (i=0; i<=p-1; i++)
      {
          l=i*q;
          for (j=0; j<=q-1; j++) a[l+j]=0.0;
      }
      if (p>n) p=n;
      if (p>20) p=20;
      if (q>m) q=m;
      if (q>20) q=20;
      d1=1.0*n; apx[0]=0.0;
      for (i=0; i<=n-1; i++)  apx[0]=apx[0]+x[i];
      apx[0]=apx[0]/d1;
      for (j=0; j<=m-1; j++)
      {
          v[j]=0.0;
          for (i=0; i<=n-1; i++)  v[j]=v[j]+z[i*m+j];
          v[j]=v[j]/d1;
      }
      if (p>1)
      {
          d2=0.0; apx[1]=0.0;
          for (i=0; i<=n-1; i++)
          {
              g=x[i]-apx[0];
              d2=d2+g*g;
              apx[1]=apx[1]+x[i]*g*g;
          }
          apx[1]=apx[1]/d2;
          bx[1]=d2/d1;
          for (j=0; j<=m-1; j++)
          {
              v[m+j]=0.0;
              for (i=0; i<=n-1; i++)
              {
                  g=x[i]-apx[0];
                  v[m+j]=v[m+j]+z[i*m+j]*g;
              }
              v[m+j]=v[m+j]/d2;
          }
          d1=d2;
      }
      for (k=2; k<=p-1; k++)
      {
          d2=0.0; apx[k]=0.0;
          for (j=0; j<=m-1; j++) v[k*m+j]=0.0;
          for (i=0; i<=n-1; i++)
          {
              g1=1.0; g2=x[i]-apx[0];
              for (j=2; j<=k; j++)
              {
                  g=(x[i]-apx[j-1])*g2-bx[j-1]*g1;
                  g1=g2; g2=g;
              }
              d2=d2+g*g;
              apx[k]=apx[k]+x[i]*g*g;
              for (j=0; j<=m-1; j++) v[k*m+j]=v[k*m+j]+z[i*m+j]*g;
          }
          for (j=0; j<=m-1; j++) v[k*m+j]=v[k*m+j]/d2;
          apx[k]=apx[k]/d2;
          bx[k]=d2/d1;
          d1=d2;
      }
      d1=m; apy[0]=0.0;
      for (i=0; i<=m-1; i++) apy[0]=apy[0]+y[i];
      apy[0]=apy[0]/d1;
      for (j=0; j<=p-1; j++)
      {
          u[j][0]=0.0;
          for (i=0; i<=m-1; i++)  u[j][0]=u[j][0]+v[j*m+i];
          u[j][0]=u[j][0]/d1;
      }
      if (q>1)
      {
          d2=0.0; apy[1]=0.0;
          for (i=0; i<=m-1; i++)
          {
              g=y[i]-apy[0];
              d2=d2+g*g;
              apy[1]=apy[1]+y[i]*g*g;
          }
          apy[1]=apy[1]/d2;
          by[1]=d2/d1;
          for (j=0; j<=p-1; j++)
          {
              u[j][1]=0.0;
              for (i=0; i<=m-1; i++)
              {
                  g=y[i]-apy[0];
                  u[j][1]=u[j][1]+v[j*m+i]*g;
              }
              u[j][1]=u[j][1]/d2;
          }
          d1=d2;
      }
      for (k=2; k<=q-1; k++)
      {
          d2=0.0; apy[k]=0.0;
          for (j=0; j<=p-1; j++) u[j][k]=0.0;
          for (i=0; i<=m-1; i++)
          {
              g1=1.0;
              g2=y[i]-apy[0];
              for (j=2; j<=k; j++)
              {
                  g=(y[i]-apy[j-1])*g2-by[j-1]*g1;
                  g1=g2; g2=g;
              }
              d2=d2+g*g;
              apy[k]=apy[k]+y[i]*g*g;
              for (j=0; j<=p-1; j++) u[j][k]=u[j][k]+v[j*m+i]*g;
          }
          for (j=0; j<=p-1; j++) u[j][k]=u[j][k]/d2;
          apy[k]=apy[k]/d2;
          by[k]=d2/d1;
          d1=d2;
      }
      v[0]=1.0; v[m]=-apy[0]; v[m+1]=1.0;
      for (i=0; i<=p-1; i++)
          for (j=0; j<=q-1; j++) a[i*q+j]=0.0;
      for (i=2; i<=q-1; i++)
      {
          v[i*m+i]=v[(i-1)*m+(i-1)];
          v[i*m+i-1]=-apy[i-1]*v[(i-1)*m+i-1]+v[(i-1)*m+i-2];
          if (i>=3)
              for (k=i-2; k>=1; k--)
                  v[i*m+k]=-apy[i-1]*v[(i-1)*m+k]+
                      v[(i-1)*m+k-1]-by[i-1]*v[(i-2)*m+k];
          v[i*m]=-apy[i-1]*v[(i-1)*m]-by[i-1]*v[(i-2)*m];
      }
      for (i=0; i<=p-1; i++)
      {
          if (i==0) { t[0]=1.0; t1[0]=1.0;}
          else
          {
              if (i==1)
              {
                  t[0]=-apx[0]; t[1]=1.0;
                  t2[0]=t[0]; t2[1]=t[1];
              }
              else
              {
                  t[i]=t2[i-1];
                  t[i-1]=-apx[i-1]*t2[i-1]+t2[i-2];
                  if (i>=3)
                      for (k=i-2; k>=1; k--)
                          t[k]=-apx[i-1]*t2[k]+t2[k-1]
                              -bx[i-1]*t1[k];
                  t[0]=-apx[i-1]*t2[0]-bx[i-1]*t1[0];
                  t2[i]=t[i];
                  for (k=i-1; k>=0; k--)
                  { t1[k]=t2[k]; t2[k]=t[k];}
              }
          }
          for (j=0; j<=q-1; j++)
              for (k=i; k>=0; k--)
                  for (l=j; l>=0; l--)
                      a[k*q+l]=a[k*q+l]+u[i][j]*t[k]*v[j*m+l];
      }
      dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
      for (i=0; i<=n-1; i++)
      {
          x1=x[i];
          for (j=0; j<=m-1; j++)
          {
              y1=y[j];
              x2=1.0; dd=0.0;
              for (k=0; k<=p-1; k++)
              {
                  g=a[k*q+q-1];
                  for (kk=q-2; kk>=0; kk--) g=g*y1+a[k*q+kk];
                  g=g*x2; dd=dd+g; x2=x2*x1;
              }
              dd=dd-z[i*m+j];
              if (fabs(dd)>dt[2]) dt[2]=fabs(dd);
              dt[0]=dt[0]+dd*dd;
              dt[1]=dt[1]+fabs(dd);
          }
      }
      delete[] v;
      return;
  }

//最小二乘曲面拟合例
#include <iostream>
#include <iomanip>
#include <cmath>
//  #include "最小二乘曲面拟合.cpp"
using namespace std;
int main()
{
    int i,j;
    double x[11],y[21],z[11][21],a[6][5],dt[3];
    for (i=0; i<=10; i++) x[i]=0.2*i;
    for (i=0; i<=20; i++) y[i]=0.1*i;
    for (i=0; i<=10; i++)
        for (j=0; j<=20; j++)
            z[i][j]=exp(x[i]*x[i]-y[j]*y[j]);
    pir2(x,y,&z[0][0],11,21,&a[0][0],6,5,dt);
    cout <<"二元拟合多项式系数矩阵 : " <<endl;
    for (i=0; i<=5; i++)
    {
        for (j=0; j<=4; j++)
            cout <<setw(14) <<a[i][j];
        cout <<endl;
    }
    cout <<"误差平方和 = " <<dt[0] <<endl;
    cout <<"误差绝对值和 = " <<dt[1] <<endl;
    cout <<"误差绝对值最大值 = " <<dt[2] <<endl;
    return 0;
}

The result is :

二元拟合多项式系数矩阵 :
0.888696 0.096095 -1.38677 0.903053 -0.171411
8.59716 0.929614 -13.4154 8.73604 -1.65821
-45.7529 -4.94728 71.3951 -46.492 8.82477
85.1474 9.20702 -132.868 86.5229 -16.4231
-62.8532 -6.79634 98.0792 -63.8686 12.1231
16.99 1.83714 -26.5121 17.2645 -3.27702
误差平方和 = 5.80139
误差绝对值和 = 25.3761
误差绝对值最大值 = 0.555972

This result should be right.

The Rust code is here:

#[no_mangle]
extern "C" fn fit_surface_least_squares_m( n:usize, m:usize,mut p:usize,mut q:usize,x:&[f64;100],y:&[f64;100],z:&[[f64;100];100],a:&mut [[f64;20];20],dt:&mut [f64;3])->bool
{
    //let mut l:usize;
    //let mut kk:usize;
    let mut apx=[0.0;20];
    let mut apy=[0.0;20];
    let mut bx=[0.0;20];
    let mut by=[0.0;20];
    let mut t=[0.0;20];
    let mut t1=[0.0;20];
    let mut t2=[0.0;20];
    //let mut d2:f64;
    let mut d=[0.0;2];
    let mut g=[0.0;3];
    //let mut g1:f64;
    //let mut g2:f64;
    let mut x2:f64;
    let mut dd:f64;
    let mut y1:f64;
    let mut x1:f64;
    let mut v=[[0.0;100];20];
    let mut u=[[0.0;20];20];
    //matrix<_Ty> v(20,m), u(20,20);

    //int p = a.GetRowNum();       //拟合多项式中变量x的最高次数加1
    //并要求p<=n且p<=20; 否则在本函数中自动取p=min{n,20}处理

    //int q = a.GetColNum();        //拟合多项式中变量y的最高次数加1
    //并要求q<=m且q<=20; 否则在本函数中自动取q=min{m,20}处理

    for i in 0..p
    {
        for j in 0..q
        {
            a[i][j]=0.0;
        }
    }   

    if p>n {p=n;}
    if p>20 {p=20;}
    if q>m {q=m;}
    if q>20 {q=20;}

    let mut xx:f64=0.0;


    for i in 0..n
    {
        xx+=  x[i] / (n as f64);
    }
    let mut yy:f64=0.0;
    for i in 0..m
    {
        yy+= y[i] / (m as f64);
    }
    d[0]=n as f64;
    apx[0]=0.0;
    for i in 0..n
    {
        apx[0]+=x[i] - xx;
    }
    apx[0] = apx[0] / d[0];
    for j in 0..m
    {
        v[0][j]=0.0;
        for i in 0..n
        {
            v[0][j]=v[0][j]+z[i][j];
        }  
        v[0][j]=v[0][j]/d[0];
    }
    if p>1
    {
        d[1]=0.0;
        apx[1]=0.0;
        for i in 0..n
        {
            g[0]=x[i]-xx-apx[0];
            d[1]+=g[0]*g[0];
            apx[1]+=(x[i]-xx)*g[0]*g[0];
        }
        apx[1]=apx[1]/d[1];
        bx[1]=d[1]/d[0];
        for j in 0..m
        {
            v[1][j]=0.0;
            for i in 0..n
            {
                g[0]=x[i]-xx-apx[0];
                v[1][j]=v[1][j]+z[i][j]*g[0];
            }
            v[1][j]/=d[1];
        }
        d[0]=d[1];
    }
    for k in 2..p
    {
        d[1]=0.0;
        apx[k]=0.0;
        for j in 0..m
        {
            v[k][j]=0.0;
        }
        for i in 0..n
        {
            g[1]=1.0;
            g[2]=x[i]-xx-apx[0];

            for j in 2..=k
            {
                g[0]=(x[i]-xx-apx[j-1])*g[2]-bx[j-1]*g[1];
                g[1]=g[2];
                g[2]=g[0];
            }
            d[1]+=g[0]*g[0];
            apx[k]=apx[k]+(x[i]-xx)*g[0]*g[0];
            for j in 0..m    
            {
                v[k][j]+=z[i][j]*g[0];
            }
        }
        for j in 0..m
        {
            v[k][j]/=d[1];
        }  
        apx[k]=apx[k]/d[1];
        bx[k]=d[1]/d[0];
        d[0]=d[1];
    }
    d[0]=m as f64;
    apy[0]=0.0;
    for i in 0..m
    {
        apy[0]+=y[i]-yy;
    }
    apy[0]=apy[0]/d[0];
    for j in 0..p
    {
        u[j][0]=0.0;
        for i in 0..m
        {
            u[j][0]+=v[j][i];
        }  
        u[j][0]/=d[0];
    }
    if q>1
    {
        d[1]=0.0;
        apy[1]=0.0;
        for i in 0..m
        {
            g[0]=y[i]-yy-apy[0];
            d[1]+=g[0]*g[0];
            apy[1]+=(y[i]-yy)*g[0]*g[0];
        }
        apy[1]=apy[1]/d[1];
        by[1]=d[1]/d[0];
        for j in 0..p
        {
            u[j][1]=0.0;
            for i in 0..m
            {
                g[0]=y[i]-yy-apy[0];
                u[j][1]+=v[j][i]*g[0];
            }
            u[j][1]/=d[1];
        }
        d[0]=d[1];
    }
    for k in 2..q
    {
        d[1]=0.0;
        apy[k]=0.0;
        for j in 0..p
        {
            u[j][k]=0.0;
        }
        for i in 0..m
        {
            g[1]=1.0;
            g[2]=y[i]-yy-apy[0];
            for j in 2..=k
            {
                g[0]=(y[i]-yy-apy[j-1])*g[2]-by[j-1]*g[1];
                g[1]=g[2];
                g[2]=g[0];
            }
            d[1]+=g[0]*g[0];
            apy[k]=apy[k]+(y[i]-yy)*g[0]*g[0];
            for j in 0..p
            {
                u[j][k]+=v[j][i]*g[0];
            }    
        }
        for j in 0..p
        {
            u[j][k]/=d[1];
        }  
        apy[k]/=d[1];
        by[k]=d[1]/d[0];
        d[0]=d[1];
    }
    v[0][0]=1.0;
    v[1][0]=-apy[0];
    v[1][1]=1.0;
    for i in 0..p
    {
        for j in 0..q
        {
            a[i][j]=0.0;
        }    
    }

    for i in 2..q
    {
        v[i][i]=v[i-1][i-1];
        v[i][i-1]=-apy[i-1]*v[i-1][i-1]+v[i-1][i-2];
        if i>=3
        {
            for k in (1..=(i-2)).rev()
            {
                v[i][k]=-apy[i-1]*v[i-1][k]+v[i-1][k-1]-by[i-1]*v[i-2][k];
            }

        }

        v[i][0]=-apy[i-1]*v[i-1][0]-by[i-1]*v[i-2][0];
    }
    for i in 0..p
    {
        if i==0
        {
            t[0]=1.0;
            t1[0]=1.0;
        }
        else
        {
            if i==1
            {
                t[0]=-apx[0];
                t[1]=1.0;
                t2[0]=t[0];
                t2[1]=t[1];
            }
            else
            {
                t[i]=t2[i-1];
                t[i-1]=-apx[i-1]*t2[i-1]+t2[i-2];
                if i>=3
                {
                    for k in (1..=(i-2)).rev()
                    {
                        t[k]=-apx[i-1]*t2[k]+t2[k-1]-bx[i-1]*t1[k];

                    }

                }

                t[0]=-apx[i-1]*t2[0]-bx[i-1]*t1[0];
                t2[i]=t[i];
                for k in (0..=(i-1)).rev()
                {
                    t1[k]=t2[k];
                    t2[k]=t[k];
                }
            }
        }
        for j in 0..q
        {
            for k in (0..=i).rev()
            {
                for l in (0..=j).rev()
                {
                    a[k][l]+=u[i][j]*t[k]*v[j][l];
                }
            }

        }

    }
    for i in 0..3
    {
        dt[i]=0.0;
    }
    //计算dt
    //dt[0]存放拟合多项式与数据点误差的平方和
    //dt[1]存放拟合多项式与数据点误差的绝对值之和
    //dt[2]存放拟合多项式与数据点误差绝对值的最大值
    for i in 0..n
    {
        x1=x[i]-xx;
        for j in 0..m
        {
            y1=y[j]-yy;
            x2=1.0;
            dd=0.0;
            for k in 0..p
            {
                g[0]=a[k][q-1];
                for kk in (0..=(q-2)).rev()
                {
                    g[0]=g[0]*y1+a[k][kk];
                }
                g[0]=g[0]*x2;
                dd+=g[0];
                x2*=x1;
            }
            dd-=z[i][j];
            if dd.abs()>dt[2]
            {
                dt[2]=dd.abs();
            }
            dt[0]=dt[0]+dd*dd;
            dt[1]=dt[1]+dd.abs();
        }
    }
    true
}

When I compile it to C dll and call it by Ruby FFI I got the following result:

系数矩阵a=[[12.610828495549994, -6.4031174330927225, -48.421984673330485, 7.374037165505316, 38.79011866159077], [-6.304552141689332, 3.2011209842825186, 24.207682095178615, -3.68651447611775, -19.39240754646785], [-1.5801195738843747, 0.8023018625207973, 6.06720849590951, -0.9239567779291769, -4.860348849570485], [1.3068055530349345, -0.6635273345643454, -5.01776060806868, 0.7641395424233058, 4.019652038542775], [0.043835634614448844, -0.022257436637848578, -0.16831633450534478, 0.025632384020965664, 0.13483566673677283], [-0.12004056416126747, 0.06095030388610434, 0.46092152946540504, -0.07019234159007005, -0.3692372574619461]]

误差平方和 = 14546.87281940942
误差绝对值和 =1105.0604796924251
误差绝对值最大值 =38.4058183240485

The Ruby FFI function is as follows:

module Qin
extend FFI::Library
….
attach_function :fit_surface_least_squares_m,[:int,:int,:int,:int,:pointer,:pointer,:pointer,:pointer,:pointer],:bool
def self.fit_surface_least_squares(p,q,x_coord,y_coord,z_coord)
    #extern "C" fn fit_surface_least_squares_m( n:usize, m:usize,mut p:usize,mut q:usize,x:&[f64;100],y:&[f64;100],z:&[[f64;100];100],a:&mut [[f64;100];100],dt:&mut [f64;3])->bool
    #int p = a.GetRowNum();        //拟合多项式中变量x的最高次数加1
    #int q = a.GetColNum();          //拟合多项式中变量y的最高次数加1

    n=x_coord.length
    m=y_coord.length

    pointer_x=FFI::MemoryPointer.new :double,100
    pointer_y=FFI::MemoryPointer.new :double,100
    pointer_x.put_array_of_double(0,x_coord)
    pointer_y.put_array_of_double(0,y_coord)


    pointer_z=FFI::MemoryPointer.new :double,100*100
    pointer_a=FFI::MemoryPointer.new :double,20*20
    pointer_dt=FFI::MemoryPointer.new :double,3
    z_coord.each_with_index do |e,i|
        pointer_z.put_array_of_double(i*100*8,e)  ##8-->f64--float,n=100
    end

    if n<=100 and m<=100
        if self.fit_surface_least_squares_m(n,m,p,q,pointer_x,pointer_y,pointer_z,pointer_a,pointer_dt)
            a=[]
            (0..(p-1)).each do |i|
                a[i]=pointer_a.get_array_of_double(i*20*8,q)
            end
            dt=pointer_dt.get_array_of_double(0,3)
            return [a,dt]
        end
    else
        nil
    end

end
…
end #end module
end

x_coord=(0..10).to_a.map{|e| e*0.2}
y_coord=(0..20).to_a.map{|e| e*0.1}
z_coord=Array.new(21) { Array.new(11) { |i|  0.0}  }
for j in 0..20
    for i in 0..10
        z_coord[j][i]=Math.exp(x_coord[i]**2-y_coord[j]**2)
    end
end
#puts z_coord.inspect

a=Qin.fit_surface_least_squares(p=6,q=5,x_coord,y_coord,z_coord)
puts a[0].inspect
puts a[1]

The conversion seems to be right but the Ruby result seems to be wrong. Why?

Please format your code according to this topic, it's mostly unreadable now.

(upd: somehow I inserted the link to this topic itself, thanks eminence for fixing)

Fixed URL for formatting instructions: Forum Code Formatting and Syntax Highlighting

You should use testing to help figure out why you get a different result than expected.

I found it. The Ruby code is wrong; test panics because different language's system. The result both are allowable.

This topic was automatically closed 90 days after the last reply. We invite you to open a new topic if you have further questions or comments.