请问各位大佬有没有能够用用c语言生成蒙特卡罗方法中的sobol序列的代码?我现在在写相关的论文,急需这样的代码,之前找到了相关的代码,但是输出的结果全部都是0,我不明白是哪里出错了,还有想请问一下sobol序列中的简单多项式的次数怎么来确定,系数如何确定? 我怀疑是我所用的sobol_para.text文本文件里的数据不对,麻烦各位大佬帮忙看看,感激不尽。
#include"math.h"
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#define MAXBIT 30
#define MAXDIM 160
#define S 2
int IMIN(int a,int b)
{
return a<b?a:b;
}
void sobseq(int *n,double x[])
{ FILE *frl;
  int j,k,l;
  unsigned long i,im,ipp;
  static double fac;
  static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
  static unsigned long mdeg[MAXDIM+1];
  static unsigned long ip[MAXDIM+1];
  static unsigned long iv[MAXDIM*MAXBIT+1];
  long temp;
if(*n<0)
{ //read mdeg[]and ip[]
  if((frl=fopen("sobol_para.txt","r"))==NULL)
 {
     printf("\nFile_not_found!");
     exit(0);
  }
  for(i=1;i<=MAXDIM;i++)
    fscanf(frl,"%d",&mdeg);//read mdeg[]
  for(i=1;i<=MAXDIM;i++)
    fscanf(frl,"%d",&ip);//read ip[]
  fclose(frl);
  //set the values of iv[]//
  srand((unsigned)time(NULL));//give the rand seed
  for(i=1;i<=mdeg[MAXDIM];i++)
  {
for(j=1;j<=MAXDIM;j++)
{temp=2*(int((1L<<(i-1))*float(rand())/RAND_MAX)+1)-1;
iv[(i-1)*MAXDIM+j]=temp;
 }
}
/*Initialize; don't return a vector*/
for(k=1;k<=MAXDIM;k++)
     ix[k]=0;
 in=0;
 if(iv[1]!=1)
         return;
 fac=1.0/(1L<<MAXBIT);
 for(j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM)
     iu[j]=&iv[k];
 for(k=1;k<=MAXDIM;k++)
 {
 for(j=1;j<=mdeg[k];j++)
      iu[j][k]<<=(MAXBIT-j);
    /*Stored values only require normalization*/
for(j=mdeg[k]+1;j<=MAXBIT;j++)
{  /*Use the recurrence to get other values.*/
   ipp=ip[k];
   i=iu[j-mdeg[k]][k];
   i ^=(i>>mdeg[k]);
   for(l=mdeg[k]-1;l>=1;l--)
  {
      if(ipp&1)
         i ^=iu[j-1][k];
      ipp>>=1;
   }
     iu[j][k]=i;
 }
}
}
else
{
/*Calculate the next vector in the sequence*/
im=in++;
for(j=1;j<=MAXBIT;j++)
{
   /*Find the rightmost=ero bit*/
   if(!(im&1))
      break;
   im>>=1;
 }
 if(j>MAXBIT)
    printf("\nMAXBIT_too_small_in_sobseq\n");
 im=(j-1)*MAXDIM;
for(k=1;k<=IMIN(*n,MAXDIM);k++)
{ 
    ix[k] ^=iv[im+k];
    x[k]=ix[k]*fac;
 }
}
}
int main()
{FILE *fwl;
 static double xsob[S+1];
 int i,j;
 int ini,run;
 ini=-1;
 sobseq(&ini,xsob);
 if(( fwl=fopen("sobol_points.txt","w"))==NULL)
 { printf( "The_file_'sobol_points'_was_not_opened\n");
  exit(0);
  }
 for(i=1;i<=1024;i++)
 { run=S;
  sobseq(&run,xsob);
  for(j=1;j<=S;j++)
     fprintf(fwl,"%f\t",xsob[j]);
  fprintf(fwl,"\n");
 }//end of i
  fclose(fwl);
}