ACM内部函数3

8.组合序列

语法：

m_of_n(int m, int n1, int m1, int* a, int head)

参数：

m：组合数C的上参数

n1：组合数C的下参数

m1：组合数C的上参数，递归之用

*a：1～n的整数序列数组

head：头指针

返回值：

null

注意：

*a需要自行产生

初始调用时，m=m1、head=0

调用例子：求C(m,n)序列：m_of_n(m,n,m,a,0);

源程序：

void m_of_n(int m, int n1, int m1, int* a, int head)

{

    int i,t;

    if(m1<0 || m1>n1) return;

 

    if(m1==n1)

        {

        for(i=0;i<m;i++) cout<<a[i]<<' '; // 输出序列

        cout<<'/n';

        return;

        }

    m_of_n(m,n1-1,m1,a,head); // 递归调用

    t=a[head];a[head]=a[n1-1+head];a[n1-1+head]=t;

    m_of_n(m,n1-1,m1-1,a,head+1); // 再次递归调用

    t=a[head];a[head]=a[n1-1+head];a[n1-1+head]=t;

}

9.快速傅立叶变换（FFT）

语法：

kkfft(double pr[],double pi[],int n,int k,double fr[],double fi[],int l,int il);

参数：

pr[n]：输入的实部

pi[n]：数入的虚部

n，k：满足n=2^k

fr[n]：输出的实部

fi[n]：输出的虚部

l：逻辑开关，0 FFT，1 ifFT

il：逻辑开关，0 输出按实部/虚部；1 输出按模/幅角

返回值：

null

注意：

需要 math.h

源程序：

void kkfft(pr,pi,n,k,fr,fi,l,il)

int n,k,l,il;

double pr[],pi[],fr[],fi[];

{

    int it,m,is,i,j,nv,l0;

    double p,q,s,vr,vi,poddr,poddi;

    for (it=0; it<=n-1; it++)

        {

         m=it; is=0;

        for (i=0; i<=k-1; i++)

            {j=m/2; is=2*is+(m-2*j); m=j;}

        fr[it]=pr[is]; fi[it]=pi[is];

        }

    pr[0]=1.0; pi[0]=0.0;

    p=6.283185306/(1.0*n);

    pr[1]=cos(p); pi[1]=-sin(p);

    if (l!=0) pi[1]=-pi[1];

    for (i=2; i<=n-1; i++)

        {

       p=pr[i-1]*pr[1];

       q=pi[i-1]*pi[1];

        s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);

        pr[i]=p-q; pi[i]=s-p-q;

        }

    for (it=0; it<=n-2; it=it+2)

        {

       vr=fr[it]; vi=fi[it];

        fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];

        fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];

        }

    m=n/2; nv=2;

    for (l0=k-2; l0>=0; l0--)

        {

        m=m/2; nv=2*nv;

        for (it=0; it<=(m-1)*nv; it=it+nv)

            for (j=0; j<=(nv/2)-1; j++)

                {

               p=pr[m*j]*fr[it+j+nv/2];

                q=pi[m*j]*fi[it+j+nv/2];

                s=pr[m*j]+pi[m*j];

                s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);

                poddr=p-q; poddi=s-p-q;

                fr[it+j+nv/2]=fr[it+j]-poddr;

                fi[it+j+nv/2]=fi[it+j]-poddi;

                fr[it+j]=fr[it+j]+poddr;

                fi[it+j]=fi[it+j]+poddi;

                }

        }

    if (l!=0)

        for (i=0; i<=n-1; i++)

            {

           fr[i]=fr[i]/(1.0*n);

            fi[i]=fi[i]/(1.0*n);

            }

    if (il!=0)

            for (i=0; i<=n-1; i++)

            {

           pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);

            if (fabs(fr[i])<0.000001*fabs(fi[i]))

                {

               if ((fi[i]*fr[i])>0) pi[i]=90.0;

                else pi[i]=-90.0;

                }

            else

                pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;

            }

    return;

}

10.Ronberg算法计算积分

语法：

result=integral(double a,double b);

参数：

a：积分上限

b：积分下限

function f：

积分函数

返回值：

f在（a,b）之间的积分值

注意：

function f(x)需要自行修改，程序中用的是sina(x)/x

需要 math.h

默认精度要求是1e-5

源程序：

double f(double x)

{

    return sin(x)/x; //在这里插入被积函数

}

 

double integral(double a,double b)

{

    double h=b-a;

    double t1=(1+f(b))*h/2.0;

    int k=1;

    double r1,r2,s1,s2,c1,c2,t2;

loop:

    double s=0.0;

    double x=a+h/2.0;

    while(x<b)

        {

        s+=f(x);

        x+=h;

        }

    t2=(t1+h*s)/2.0;

    s2=t2+(t2-t1)/3.0;

    if(k==1)

      {

        k++;h/=2.0;t1=t2;s1=s2;

        goto loop;

        }

    c2=s2+(s2-s1)/15.0;

    if(k==2){

        c1=c2;k++;h/=2.0;

        t1=t2;s1=s2;

        goto loop;

        }

    r2=c2+(c2-c1)/63.0;

    if(k==3){

        r1=r2; c1=c2;k++;

        h/=2.0;

        t1=t2;s1=s2;

        goto loop;

        }

    while(fabs(1-r1/r2)>1e-5){

        r1=r2;c1=c2;k++;

        h/=2.0;

        t1=t2;s1=s2;

        goto loop;

        }

    return r2;

}