函数接口为 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表示
#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这样的小字符串
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 '*':
}
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;
}
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部分,分成三步求值)
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时 即字符串已被切分开 */
}
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;
}