跪求能用c语言生成蒙特卡罗方法中的sobol序列的代码

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

   请问各位大佬有没有能够用用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);
}

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：蒙特卡罗方法 蒙特卡罗 蒙特卡 OBO C语言

你好，您的C语言版本的Sobol序列代码改好了吗，可以给我一份吗？谢谢