📄 main.cpp
字号:
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include "mymatrix.h"
const long step1=75;
const long step2=90;////////////
const long step3=75;
const long step4=90;
const long max=step1*step2*step3*step4+1;/////
const float h1=5.05/step1;///////
const float h2=6.06/step2;///////
const float h3=5.05/step3;
const float h4=6.06/step4;
const float PI=3.14159265358979;
float x01=-2.525;
float x02=-3.03;
float x03=-2.525;
float x04=-3.03;
long gn; // 总组号
long ugn; ///总的不稳定组号
long g_attr;
struct dom
{
long total;
long n;
};
struct gcmdata
{
long Iz;
long Jz;
long gr;
long CpLoc;
long PreLoc;
};
struct cpdata
{
long Cz;
float Pz;
};
class lin
{
public:
lin *next;
long data;
lin();
};
lin::lin()
{
next=NULL;
}
class xlink
{
public:
float x;
int num;
xlink *next;
xlink();
};
xlink::xlink()
{
next=NULL;
}
void main()
{
gn=0;
ugn=0;
long map(long z,long j,long n);
void getbasedata();
void ab_per();
void determine();
void outputdata();
void fenzu();
//////////////////////////////////////////////////////////////////
///*
getbasedata();
ab_per();
/////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//gn=3;
//*
cout<<" start calculu"<<endl;
g_attr=gn;
ofstream of("per_groupn.txt");
of<<"persistent group numbers: "<<gn<<endl; /////////读取数据时要参考
////////////////////////////
if(gn<2)
{
cout<<"there 's no stable zone"<<endl;
exit(0);
}
////////////////////////////
//*/
outputdata();
//determine();
//fenzu();
//*/
}
long map(long m,long j)
{
//*
float x1,x2,x3,x4;
long z1,z2,z3,z4;
long j1,j2,j3,j4;
long number=2;
j4=(long)((j-1)/(number*number*number))+1;
j=j-(j4-1)*number*number*number;
j3=(long)((j-1)/(number*number))+1;
j=j-(j3-1)*(number*number);
j2=(long) ((j-1)/number)+1;
j1=j-(j2-1)*number;
z4=(long)((m-1)/(step1*step2*step3))+1;
m=m-(z4-1)*step1*step2*step3;
z3=(long)((m-1)/(step1*step2))+1;
m=m-(z3-1)*step1*step2;
z2=(long) ((m-1)/step1)+1;
z1=m-(z2-1)*step1;
x1=z1*h1-h1+x01+j1*h1/number-0.5*h1/number;
x2=z2*h2-h2+x02+j2*h2/number-0.5*h2/number;
x3=z3*h3-h3+x03+j3*h3/number-0.5*h3/number;
x4=z4*h4-h4+x04+j4*h4/number-0.5*h4/number;
//float x1,x2,x3,x4;
float h=0.08;
float mu=0.25;
float yt=0.04;
float v=0.1;
float t;
float k1,k2,k3,k4,l1,l2,l3,l4,s1,s2,s3,s4,d1,d2,d3,d4;
////////////////////////////
for(long i=0;i<30;i++)
{
k1=x2;
l1=mu*(1-x1*x1)*x2-(1+v)*x1+v*x3;
s1=x4;
d1=mu*(1-x3*x3)*x4+v*x1-(1+yt+v)*x3;
k2=x2+h*l1/2;
l2=mu*(1-(x1+h*k1/2)*(x1+h*k1/2))*(x2+h*l1/2)-(1+v)*(x1+h*k1/2)+v*(x3+h*s1/2);
s2=x4+d1*h/2;
d2=mu*(1-(x3+h*s1/2)*(x3+h*s1/2))*(x4+d1*h/2)+v*(x1+h*k1/2)-(1+yt+v)*(x3+h*s1/2);
k3=x2+h*l2/2;
l3=mu*(1-(x1+h*k2/2)*(x1+h*k2/2))*(x2+h*l2/2)-(1+v)*(x1+h*k2/2)+v*(x3+h*s2/2);
s3=x4+d2*h/2;
d3=mu*(1-(x3+h*s2/2)*(x3+h*s2/2))*(x4+d2*h/2)+v*(x1+h*k2/2)-(1+yt+v)*(x3+h*s2/2);
k4=x2+h*l3/2;
l4=mu*(1-(x1+h*k3)*(x1+h*k3))*(x2+h*l3)-(1+v)*(x1+h*k3)+v*(x3+h*s3);
s4=x4+d3*h/2;
d4=mu*(1-(x3+h*s3)*(x3+h*s3))*(x4+d3*h)+v*(x1+h*k3)-(1+yt+v)*(x3+h*s3);
x1=x1+h*(k1+2*k2+2*k3+k4)/6.0;
x2=x2+h*(l1+2*l2+2*l3+l4)/6.0;
x3=x3+h*(s1+2*s2+2*s3+s4)/6.0;
x4=x4+h*(d1+2*d2+2*d3+d4)/6.0;
}
z1=(long)((x1-x01)/h1)+1;
z2=(long)((x2-x02)/h2)+1;
z3=(long)((x3-x03)/h3)+1;
z4=(long)((x4-x04)/h4)+1;
if(z1>=step1 || z1<=0 || z2 >=step2 || z2<=0 || z3>=step3 || z3<=0 || m==0 || z4>=step4 || z4<=0)
return 0;
else
{
return (z4-1)*step1*step2*step3+(z3-1)*step1*step2+(z2-1)*step1+z1;
}
//*/
}
void getbasedata()
{
fstream iogcm,iocp,iopre;
iocp.open("cp.dat",ios::out|ios::binary|ios::in);
iogcm.open("gcm.dat",ios::in|ios::out|ios::binary);
//ifstream igcm("gcm.dat",ios::in);
iopre.open("pre.dat",ios::out|ios::binary|ios::in);
gcmdata gcm;
cpdata cp;
///*
long tz=0;
long snm=16;
long tempz[16][2]; ////////////记录
//////////////////////
for(long i=0;i<max;i++)
{
if(i%20==0)
cout<<"get Iz CpLoc cp.dat"<<i<<endl;
for(long k=0;k<snm;k++) /////////////每次都要置0 -1
{
for(long e=0;e<2;e++)
{
if(e==0)
{
tempz[k][e]=-1;
}
else
{
tempz[k][e]=0;
}
}
}
long r=0;
for(long j=0;j<snm;j++)
{
tz=map(i,j);
r=0;
while(tempz[r][0]!=tz && tempz[r][0]!=-1) ////////搜索到tz 或 tz为新值
{
r++;
}
tempz[r][0]=tz;
tempz[r][1]=tempz[r][1]+1;
}
r=0;
while(r<snm && tempz[r][0]!=-1)
{
r++;
}
//Iz[i]=r;
gcm.Iz=r;
long preIz;
if(i==0)
{
//CpLoc[0]=0;
gcm.CpLoc=0;
}
else
{
gcm.CpLoc=gcm.CpLoc+preIz;
}
///
preIz=r;
gcm.Jz=0;
gcm.PreLoc=0;
gcm.gr=0;
//outgcm.seekp(i*sizeof(gcmdata));
iogcm.write((char *)(&gcm),sizeof(gcmdata));
for(long t=0;t<r;t++)
{
cp.Cz=tempz[t][0];
cp.Pz=(float)(tempz[t][1])/snm;
///////////////////
iocp.write((char *)(&cp),sizeof(cpdata));
}
}
cout<<"The following part will get preimage data,please wait,first Jz"<<endl;
////////////////////// preimage ////////////////////////////////////
///iogcm.seekg(0+sizeof(long));
for(long b=0;b<max;b++)
{
if(b%2000==0)
cout<<"get Jz "<<b<<endl;
iogcm.seekg(b*sizeof(gcmdata));
iogcm.read((char *)(&gcm),sizeof(gcmdata));
long I=gcm.Iz;
for(long n=0;n<I;n++)
{
long jz;
long loc=gcm.CpLoc+n;
iocp.seekg(loc*sizeof(cpdata));
long cel;
iocp.read((char *)(&cel),sizeof(long));
iogcm.seekg(cel*sizeof(gcmdata)+sizeof(long),ios::beg);
iogcm.read((char *)(&jz),sizeof(long));///////先读取jz
jz++;
iogcm.seekp(cel*sizeof(gcmdata)+sizeof(long),ios::beg);
iogcm.write((char *)(&jz),sizeof(long));
}
}
//////// 计算 preloc /////
cout<<endl<<"get preLoc "<<endl;
long preloc;
for(b=0;b<max;b++)
{
if(b%2000==0)
cout<<"get preLoc "<<b<<endl;
if(b==0)
{
preloc=0;
iogcm.seekp(0+4*sizeof(long));//// preLoc的位置
iogcm.write((char *)(&b),sizeof(long));
}
else
{
iogcm.seekg((b-1)*sizeof(gcmdata)+sizeof(long));//////// prejz
long prejz;
iogcm.read((char *)(&prejz),sizeof(long));
preloc=preloc+prejz;
iogcm.seekp(b*sizeof(gcmdata)+4*sizeof(long));
iogcm.write((char *)(&preloc),sizeof(long));
}
}
///////////////////////////////////// 保存 preimage cell ///
cout<<" preimage cell "<<endl;
iogcm.seekg((max-1)*sizeof(gcmdata)+sizeof(long));
long jzmax;
iogcm.read((char *)(&jzmax),sizeof(long));
long totalpre=jzmax+preloc+1;///////////////////////////////////
/////////////////////////////////////////////////////////////////算法很糟糕/////////////////
cout<<"write preimage"<<endl;
long *pret;
pret=new long [totalpre];
for(b=0;b<totalpre;b++)
{
if(b%10000==0)
cout<<"write pre.dat "<<b<<endl;
pret[b]=-1;
}
cout<<" get preimage "<<endl;
//////////////////////////////////////////read prel data
long *ptt;
ptt=new long [max];
iogcm.seekg(0);
for(long tt=0;tt<max;tt++)
{
iogcm.read((char *)(&gcm),sizeof(gcmdata));
ptt[tt]=gcm.PreLoc;
}
iocp.seekg(0);
iogcm.seekg(0);
for(b=0;b<max;b++)
{
if(b%500==0)
cout<<"get preimage "<<b<<endl;
//cout<<b<<" "<<b<<endl;
iogcm.read((char *)(&gcm),sizeof(gcmdata));
long I=gcm.Iz;
for(long n=0;n<I;n++)
{
long cel;
iocp.read((char *)(&cp),sizeof(cpdata));
cel=cp.Cz;
long preL;
preL=ptt[cel];
long t=0;
long pce;
pce=pret[preL+t];
while(pce!=-1)
{
t=t+1;
pce=pret[preL+t];
}
pret[preL+t]=b;
}
}
////////////write data
for(tt=0;tt<totalpre;tt++)
{
if(tt%2000==0)
cout<<"save predata "<<tt<<endl;
iopre.write((char *)(&pret[tt]),sizeof(long));
}
delete [] ptt;
delete [] pret;
iogcm.close();
iocp.close();
iopre.close();
}
void ab_per()
{
////select simple map group
fstream iosgr,iogcm1,iocp1,ioscm,iopersistent;
iosgr.open("sgr.dat",ios::in|ios::out|ios::binary);
iogcm1.open("gcm.dat",ios::in|ios::out|ios::binary);
iocp1.open("cp.dat",ios::in|ios::out|ios::binary);
ioscm.open("tempscm.dat",ios::in|ios::out|ios::binary);////保存scm的结果
iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);
gcmdata gcm;
cpdata cp;
///*
cout<<"simple cell method"<<endl;
long sgr;
long sgn=0;
{
long temz[10000000]; //==========================????????????????高维时注意
for(long i=0;i<max;i++)
{
sgr=0;
iosgr.write((char *)(&sgr),sizeof(long));
//sgr[i]=0;
}
long m,b;
long g=0;
for(long j=0;j<max;j++)
{
if(j%1000==0)
cout<<j<<" simple cell method, get sgn "<<endl;
m=0;
iosgr.seekg(j*sizeof(long));
iosgr.read((char *)(&sgr),sizeof(long));
if(sgr!=0)
{
continue; ////////////////不是virgin cell,重新选取
}
else
{
sgr=-1;
iosgr.seekp(j*sizeof(long));
iosgr.write((char *)(&sgr),sizeof(long));
b=j;
temz[m]=b;
do
{
sgr=-1;
iosgr.seekp(b*sizeof(long));
iosgr.write((char *)(&sgr),sizeof(long));
m=m+1;
////////////////
iogcm1.seekg(b*sizeof(gcmdata));
iogcm1.read((char *)(&gcm),sizeof(gcmdata));
long loc=gcm.CpLoc;
iocp1.seekg(loc*sizeof(cpdata));
iocp1.read((char *)(&cp),sizeof(cpdata));
b=cp.Cz;
//b=Cz[CpLoc[b]];//////select the first image cell
temz[m]=b;
iosgr.seekg(b*sizeof(long));
iosgr.read((char *)(&sgr),sizeof(long));
}while(sgr==0); ///
if(sgr==-1) ///////////////新周期解
{
///////////寻找i/////////（刚好在周期轨道上）
g++;
sgn=g;
long u=0;
while(temz[u]!=b)
{
u++;
}
for(long q=0;q<m;q++)
{
if(q<u)
{
iosgr.seekg(temz[q]*sizeof(long));
long c=-2;
iosgr.write((char *)(&c),sizeof(long));
//sgr[temz[q]]=-2;
}
else
{
iosgr.seekp(temz[q]*sizeof(long));
iosgr.write((char *)(&g),sizeof(long));
long tt=temz[q];
ioscm.write((char *)(&g),sizeof(long));/////////先写group
ioscm.write((char *)(&tt),sizeof(long));
}
}
}
else //////吸引到已知一点上
{
for(long k=0;k<m;k++)
{
iosgr.seekg(temz[k]*sizeof(long));
long c=-2;
iosgr.write((char *)(&c),sizeof(long));
}
}
}
}
}
cout<<"search for persistent groups"<<endl;
////////////////////////////////////////////////////////////////////////////////////////
{
long temp[10000000]; ////======================??????????????save candidate cells
long n=0;
long flag=1;
//ioscm.clear();
for(long i=1;i<=sgn;i++)
{
cout<<sgn<<" "<<i<<endl;
///////////////////////////////先判断组号是否已经是absorbing group了
n=0;
ioscm.seekg(0);
long cel,gg;
ioscm.read((char *)(&gg),sizeof(long));
ioscm.read((char *)(&cel),sizeof(long));
long cur,pend;
cur=ioscm.tellg();
ioscm.seekg(0,ios::end);
pend=ioscm.tellg();
ioscm.seekg(cur);
while(pend!=ioscm.tellg())
{
if(gg==i)
{
temp[n]=cel;
n++;
}
ioscm.read((char *)(&gg),sizeof(long));
ioscm.read((char *)(&cel),sizeof(long));
}
ioscm.seekg(0);
cur=ioscm.tellg();
ioscm.seekg(0,ios::end);
pend=ioscm.tellg();
ioscm.seekg(cur);
n=0;
while(pend!=ioscm.tellg())
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -