任意积分函数（采用龙贝格算法）

函数接口为 double romberg_integer(char string[],double up,double down,double e)
   up 为上限 down为下限  e误差
   函数支持 ax^n 三角函数和其反函数  指数函数 对数函数
   注意：使用时诸如 x+1.2sinx/x  要写为 1x+1.2sinx/1x
        用反函数时用 sinnx cosnx tannx 表示
        对数函数用log(b)x表示

        函数除对数函数有括号，其余不可有括号

这里给出龙贝格迭代式：

T（i）=T（0）=（b-a)/2*(f(a)+f(b))

T（j）=T（i）/2+(b-a)/（2^k） *  ∑（j=0->2^(k-1)-1)  f(a+(2i+1)*(b-a)/2^k)                                           

由于要实现任意积分，而f(x)的组合是千变万化的，所以必须将一个大函数分解成一一个基本函数（如 ax^n , bx , c, asinx 等）然后进行模式匹配转变为相应代码，进行求值。

这部分代码包含在deal_function 函数中

#include<string.h>            
#define max_string_size 20
int failure[max_string_size];
char xn[5]="ax^n",
    ax[3]="ax",
    lnx[5]="alnx",
    logx[10]="alog(b)x",
    SIN[6]="asinx",
    COS[6]="acosx",
    TAN[6]="atanx",
    COSN[6]="acosnx",
    SINN[6]="asinnx",
    TANN[6]="atannx";           /* 要匹配的字符串 */
static void fail(char *pat)     /* fail pmatch 函数为KMP 字符匹配算法 */
{
    int n=strlen(pat),i,j;
    failure[0]=-1;
    for(j=1;j<n;j++)
    {
        i=failure[j-1];
        while((*(pat+j)!=*(pat+i+1))&&(i>=0)){
            i=failure[i];
        }
        if(pat[j]==pat[i+1]){
             failure[j]=i+1;
        }
        else {
            failure[j]=-1;
        }
    }
}
static int pmatch(char *string,char *pat)
{
    void fail(char *pat);
    int i=0,j=0;
    int lens=strlen(string);
    int lenp=strlen(pat);
    fail(pat);
    while(i<lens&&j<lenp)
    {
        if(string[i]==pat[j]){
            i++;j++;
        }
        else if(j==0) i++;
        else j=failure[j-1]+1;
    }
    return ((j==lenp)?1:0);
}
static double deal_function(char string[],double x) /*本函数对任意函数求值 返回f(x) f(x)可以为ax^n 三角函数和其反函数  指数函数 对数函数
                                                      将诸如1.7x^2+2.5x+1.2 这样的字符串为1.7x^2  2.5x  1.2这样的小字符串

                                                      然后将1.7x^2  2.5x  1.2 分别转化为 ax^n bx c   然后进行字符串匹配 求值

例：string=1.2x^2+2x 

string 会按 + - * / 切分，并放在tmp_string 然后将tmp_string 中的数字转化为变量 a , b ,n中 再形成新的字符串newtmp_string

newtmp_string 进行模式匹配

例： string:   1.2x^2+5x+3 

           tmp_string    1.2x^2   5x   3(string 被切分为3部分,分成三步求值）

            newtmp_string   ax^2   ax    a*/(a=1.2 a=5  a=3)*/{
    int pmatch(char *string,char *pat);
    char tmp_string[10],newtmp_string[10],form_number[10];
    char lastoptr='+',nextoptr;
    double tmpresult=0,finresult=0,a=1,b=1,n=1;
    int i=0,j,k,t,len;
    short buttom=1,flags=1;
    if(*string=='-')  {
        lastoptr='-';
        i=1;
    }
    len=strlen(string);
    *(string+len)='+';
    *(string+len+1)='\0';
    for(j=0;*(string+i)!='\0';i++){
        switch(*(string+i)){
        case '+':
        case '-':
        case '*':

        case '/':buttom=0;nextoptr=*(string+i);break; /* 当buttom为0时 即字符串已被切分开  */

        default:*(tmp_string+j++)=*(string+i);break;
        }
        if(buttom==0)     {
            *(tmp_string+j)='\0';
            for(j=0,k=0,t=0;;){ /* form_number 的作用是将1.75x^2 中的数字字符转化为数字 */
                for(t=0;(*(tmp_string+j)>='0')&&(*(tmp_string+j)<='9')||(*(tmp_string+j)=='.');t++,j++){
                    *(form_number+t)=*(tmp_string+j);
                }
                *(form_number+t)='\0';
                if((*(tmp_string+j)>='a')&&(*(tmp_string+j)<='z')) {
                    if(*(form_number)!=0){
                        a=atof(form_number);
                        *(newtmp_string+k++)='a';
                    }
                }
                else if(*(tmp_string+j)=='^'){ /* ax^n的情况  */
                    flags=0;                
                }
                else if(*(tmp_string+j)==')'){ /* log(b)x的情况 */
                    b=atof(form_number);
                    *(newtmp_string+k++)='b';
                }
                else if(*(tmp_string+j)=='\0'&&(flags==0)){
                    n=atof(form_number);
                    *(newtmp_string+k++)='n';
                    break;
                }
                else if(flags!=0&&*(tmp_string+j)=='\0'){ /* tmp_string为常数情况 */
                    tmpresult=atof(form_number);
                    break;
                }
                *(newtmp_string+k++)=*(tmp_string+j++);   
            }
            *(newtmp_string+k)='\0';
            if(pmatch(newtmp_string,xn)){
                tmpresult=pow(x,n)*a;  
            }
            else if(pmatch(newtmp_string,ax)){
                tmpresult=a*x;
            }
            else if(pmatch(newtmp_string,lnx)){
                tmpresult=a*log(x);
            }
            else if(pmatch(newtmp_string,logx)){
                tmpresult=log10(x)/log10(b)*a;
            }
            else if(pmatch(newtmp_string,SIN)){
                tmpresult=a*sin(x);
            }
            else if(pmatch(newtmp_string,COS)){
                tmpresult=a*cos(x);
            }
            else if(pmatch(newtmp_string,TAN)){
                tmpresult=a*tan(x);
            }
            else if(pmatch(newtmp_string,SINN)){
                tmpresult=a*asin(x);
            }
            else if(pmatch(newtmp_string,COSN)){
                tmpresult=a*acos(x);
            }
            else if(pmatch(newtmp_string,TANN)){
                tmpresult=a*atan(x);
            }
            switch(lastoptr){
            case '+': finresult+=tmpresult;break;
            case '-': finresult-=tmpresult;break;
            case '*': finresult*=tmpresult;break;
            case '/': finresult/=tmpresult;break;
            default: printf ("you got a trouble in deal_function ,line 149\n");
            }
            a=b=n=buttom=flags=1; /* 回复初始状态 */
            lastoptr=nextoptr;
            j=0;
        }
    }
    return finresult;
}
double romberg_integer(char string[],double up,double down,double e) /*龙贝格迭代公式
                                                                       T_i=T（0）=（b-a)/2*(f(a)+f(b))
                                                                      T_j=T_i/2+(b-a)/2^k +F
                                                                      F=f(a+(2i+1)*(b-a)/2^k[])  F从0累加到2^(k-1)
                                                                      |T_i-T_j|<=e
                                                                     */
{
    double deal_function(char string[],double x);
    double T_i,T_j,tmpresult,middle_num,x;
    int i,limit=1,k=1;
    if(down>up){
        printf ("error,uplimit can not smaller than lowlimit\n");
        return 0;
    }
    middle_num=(up-down)/2;
    T_i=middle_num*(deal_function(string,up)+deal_function(string,down));
    do{
        for(i=0,tmpresult=0;i<=(limit-1);i++){
            x=down+(2*i+1)*middle_num;
            tmpresult+=deal_function(string,x);
        }
        T_j=T_i/2+middle_num*tmpresult;
        k++;
        middle_num/=2;
        limit*=2;
        if(abs(T_i-T_j)<=e){
            break;
        }
        else{
            T_i=T_j;
        }
    }while(1);
    return T_j;
}