基于8位加法的初等函数计算器设计
2004年初，我在一家台资小公司上班，公司老板想在电子行业有一个持久的，产品来维系公司的资金来源，开始着手于计算器的研制。最后把这个任务分派给了我，由我写出初等函数，多项表达式由老板自己写。于是我开始一项痛苦的摸索过来，有两上问题我必须解决：
1． 成本问题（即要达到高精度，也要考虑成本问题，主芯片选用6502,下Mask 2-3元人发币，工程可实现性，计算一个初等函数最多允许1/20秒）.
2． 计算方法，在上大学的时候用的是泰勒级数，但是展开式太多了，是不是别他国外的公司也是用这个级数展开了，这是一个问题。
我看了许多的国外的网站，讲到了卡西欧 CASIO的设计,但是也没有说到具体的底层应该怎样构建。我查阅了一些关于数据计算方法的书籍，其中有一个叫密比雪夫多项式的数学方法吸引了我，直觉告诉我，这就是我要找到的。后来计算器的设计，我用的是密比雪夫多项式的系数，而不是用泰勒级数展开得的系数。
整个程序的过程是基于8位的加减法进行的，浮点数是BCD码。我把它做成了Class,用C++编写，其实也同struct 一样的。对于浮点数的乘除法，我采用了建倍数表的方法。
设计的整过程序是基于{AddSub8. ADC8();AddSub8.SBC8()}8位的加减法去模拟6502,
用它搭建浮点数的加减乘除: FPAddSub(),FPMul();FPDiv();其中Cchebft 用于计算多项式的系数.
笔者把这几年收集的关于国人对计算器的研究**”计算器的计算误差与计算精确度的研讨”,”怎样用没有开立方根按键的计算器开立方根和高次方根”等等文件及计算器源程序的工程文件放置网站: http://www.sokutek.com/document.asp,如果您感兴趣，请去下载。如有疑问也可以来电，我们一起讨论:chenshiyangyi@163.com;QQ:3637323.
还有一个问题，就是对数函数，并不是一直都用切比雪夫系数计算，在一定泛围，用了帕德函数计算。具体的方法，可能要做适当的改动.
源程序如下：
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <math.h>
#define CHM 7
#define ENTER printf("\n");
#define PI 3.141592653589792384626433832795
//程序有一些问题：可能在运算前先判断末尾连续为零的个数可提高效率
class AddSub8
{
public:
bool c;//进位、借位标志
short ADC8(short,short);//8bits相加
short SBC8(short,short);//8bits相减
void CLC(){c=0;};
void SEC(){c=1;};
};
short AddSub8::ADC8(short a1,short a2)
{
short result;
if(c==true) result=a1+a2+1;
else result=a1+a2;
if(result>99) { c=1;return result-100;}
else {c=0;return result;}
}
short AddSub8::SBC8(short a1,short a2)
{
short result;
if(c==true) result=a1-a2;
else result=a1-a2-1;
if(result<0) {c=0;return result+100;}
else {c=1;return result;}
};
class FP
{
public:
bool Sign;//符号 (1)为负
short Exponenet;//以二进制形式表示的指数
short Mantissa[CHM];//CHM BCD码形式的尾数
void Equ(short *data) {for(int i=0;i<CHM;i++) Mantissa=data;}
// bool ExponenetSign() {if((Exponenet&128)>0) return true; return false;}
FP& operator =(FP & fp)
{
Sign=fp.Sign;
Exponenet=fp.Exponenet;
for(int i=0;i<CHM;i++) Mantissa=fp.Mantissa;
return *this;}
friend FP operator +(FP,FP);
friend FP operator -(FP,FP);
friend FP operator *(FP,FP);
friend FP operator /(FP,FP);
};
FP operator +(FP num1,FP num2)
{
FP num3;
void FPAddSub(FP,FP,FP&);
FPAddSub(num1,num2,num3);
return num3;
}
FP operator -(FP num1,FP num2)
{
FP num3;
void FPAddSub(FP,FP,FP &);
num2.Sign=!num2.Sign;
FPAddSub(num1,num2,num3);
return num3;
}
FP operator *(FP num1,FP num2)
{
FP num3;
void FPMul(FP,FP,FP&);
FPMul(num1,num2,num3);
return num3;
}
FP operator /(FP num1,FP num2)
{
FP num3;
bool FPDiv(FP,FP,FP&);
FPDiv(num1,num2,num3);
return num3;
}
class Vector
{
public:
void PrintfShort(short *data,int n,bool flag=true)
{if(flag==true) printf("\n");
for(int i=0;i<n;i++)
{ if(data==0) {printf("00");continue;}
if(data<10) printf("0");
printf("%d",data);
}
}
void BCDLeftShift(short *,short *,int);
void BCDLeftShift(short *,int);
void PrintfFp(FP);
void Equ(short *data1,short *data2,int n){for(int i=0;i<n;i++) data1=data2;}
void PrintZero(int n){for(int i=0;i<n;i++) printf("0");}
void ClrFP(FP &data) {data.Exponenet=0;data.Sign=false;for(int i=0;i<CHM;i++) data.Mantissa=0;}
void ClrShort(short *data,int n){for(int i=0;i<n;i++) data=0;}
void BCD8ToBCD4(short *,int,short *);
void BCDRightShift(short *,int);//左移一位，高BCD移出，末低BCD添零
double FP2Double(FP);//把FP数据类型转化为double型
void Double2FP(double,FP&);
FP Double2FP(double);
void Norm(FP &);
void PrintfDouble(double *,int);
};
FP Vector::Double2FP(double a)
{
FP num1;
Double2FP(a,num1);
return num1;
}
void Vector::PrintfDouble(double *data,int num)
{
for(int i=0;i<num;i++) printf("%.15e\n",data);
}
void Vector::Norm(FP &x)
{
for(int i=0;i<2*CHM;i++)
if(x.Mantissa[0]/10==0) {x.Exponenet--;BCDRightShift(x.Mantissa,CHM);}
else
break;
if(i==2*CHM) x.Exponenet=0;
}
void Vector::Double2FP(double a,FP &fp)
{
if(a>0) fp.Sign=false;
else
{
fp.Sign=true;
a=fabs(a);
}
fp.Exponenet=0;
if(a<1e-110) {ClrFP(fp);return;}
int exp=0;
for(int i=0;i<110;i++)
{
if(floor(a)==0) {a*=10;fp.Exponenet--;}
else
{
a/=10;fp.Exponenet++;
}
if(floor(a)>=1 && floor(a)<=9) break;
}
a=a*10;
for(i=0;i<CHM;i++)
{
fp.Mantissa=(short)floor(a);
a=(a-floor(a))*100;
}
}
double Vector::FP2Double(FP num)
{
double a=0,b=10;
for(int i=0;i<CHM;i++)
{
a+=num.Mantissa/b;
b*=100;
}
a*=pow(10,num.Exponenet);
if(num.Sign==false) return a;
else
return a*(-1);
}
void Vector::BCDRightShift(short *num,int n)//左移一位BCD码
{
for(int i=0;i<n-1;i++)
num=(num%10)*10+num[i+1]/10;
num[n-1]=(num[n-1]%10)*10;
}
void Vector::BCDLeftShift(short *num,int n)
{
for(int i=n-1;i>=1;i--)
num=num/10+(num[i-1]%10)*10;
num[0]=num[0]/10;
}
void Vector::BCD8ToBCD4(short *num1,int n,short *num)
{
for(int i=0;i<n;i++)
{
num[2*i]=num1/10;//在汇编程序的实现: &0xF0
num[2*i+1]=num1%10;//&0x0F
}
}
void Vector::BCDLeftShift(short *num1,short *num2,int n)
{
short a=0,b=0;
for(int i=0;i<n;i++)
{
a=num1/10;
num2=b*10+a;
b=num1%10;
}
num2[n]=b*10;
};
void Vector::PrintfFp(FP data)
{
if(data.Sign==true) printf("-");
printf("%d.%d",data.Mantissa[0]/10,data.Mantissa[0]%10);
for(int i=1;i<CHM;i++)
{
if(data.Mantissa<=9) printf("0%d",data.Mantissa);
else
printf("%d",data.Mantissa);
}
printf("E%d",data.Exponenet);
}
class FileFP
{
public:
void SaveDouble(double);
void SaveFP(FP);
void SaveCoAsm(FP);
FP ReadFP();
FileFP(char *,bool);
~FileFP(){fclose(fp);}
private:
FILE *fp;
};
FileFP::FileFP(char *str,bool flag=true)
{
if(flag==true) fp=fopen(str,"w");
else
fp=fopen(str,"r");
if(fp==NULL) exit(0);
}
FP FileFP::ReadFP()
{
char c1,c2;
short n=0,n1,n2=0;
FP x;
Vector vect; vect.ClrFP(x);
c1=getc(fp);
while(c1<'0' || c1>'9')
{
if(c1=='-') x.Sign=true;
c1=getc(fp);
}
c2=getc(fp);//小数点
while(c1!='e' && c1!='E')
{
c2=getc(fp);
x.Mantissa[n]=(c1-'0')*10+c2-'0';
n++;
c1=getc(fp);
}
fscanf(fp,"%d",&n1);
x.Exponenet=n1;
return x;
}
void FileFP::SaveDouble(double x)
{
fprintf(fp,"%.13e\n",x);
}
void FileFP::SaveFP(FP x)
{
if(x.Sign==true) fprintf(fp,"-");
fprintf(fp,"%d.%d",x.Mantissa[0]/10,x.Mantissa[0]%10);
for(int i=1;i<CHM;i++)
{
if(x.Mantissa<=9) fprintf(fp,"0%d",x.Mantissa);
else
fprintf(fp,"%d",x.Mantissa);
}
fprintf(fp,"e%d\n",x.Exponenet);
}
void FileFP::SaveCoAsm(FP x)
{
fprintf(fp,"\n .DB %d",x.Exponenet);
if(x.Sign==true) fprintf(fp,",1");
else
fprintf(fp,",0");
for(int i=0;i<CHM;i++)
if(x.Mantissa<=9) fprintf(fp,",$0%d",x.Mantissa);
else
fprintf(fp,",$%d",x.Mantissa);
}
void FPAddSub(FP num1,FP num2,FP & num3)
{
short shift=0,i;
bool flag=false;//false 表示num1为主运算区,num2进行移位操作 ,c表示有无进位
short num[CHM+2];//进行移位对齐时应该注意任零两加数都不应该为零
if(num1.Mantissa[0]==0) {num3=num2;return;}
if(num2.Mantissa[0]==0) {num3=num1;return;}//不是很规范；判零
num[CHM]=0;
if(num1.Exponenet==num2.Exponenet) //指数相等，比较尾数
{
for(i=0;i<CHM;i++)
{ if(num1.Mantissa==num2.Mantissa) continue;
if(num1.Mantissa>num2.Mantissa) flag=false;
else
flag=true;
break;
}
}
else//指数不相等，比较指数
{
shift=num1.Exponenet-num2.Exponenet;
if(shift>0) flag=false;
else
{
flag=true;shift=shift*(-1);
}
}
//确定主运算区(num3),加数(num)
Vector vect;
if(flag==false)
{num3=num1;
if(shift%2==0) vect.Equ(num,num2.Mantissa,CHM);
else
vect.BCDLeftShift(num2.Mantissa,num,CHM);//BCD右移一位
}
else
{num3=num2;
if(shift%2==0)
vect.Equ(num,num1.Mantissa,CHM);
else
vect.BCDLeftShift(num1.Mantissa,num,CHM);//BCD码右移一位
}
AddSub8 add;
shift=shift/2;
if(shift>CHM+1) return;
short *result;
result=new short[4*CHM];//可能太浪费空间
vect.ClrShort(result,3*CHM);
if((num1.Sign^num2.Sign)==false) //加法
{
add.CLC();
for(i=shift+CHM;i>=0;i--)
if(i>=CHM) result[i+1]=add.ADC8(0,num[i-shift]); //有几步是多余的
else if(i<shift) result[i+1]=add.ADC8(num3.Mantissa,0);
else
result[i+1]=add.ADC8(num3.Mantissa,num[i-shift]);
result[0]=add.ADC8(0,0);//以上完成了数值的加法
if(result[0]==0)
vect.Equ(num3.Mantissa,result+1,CHM);//没有进位
else//否则右移一位
{
num3.Exponenet++;
//vect.BCDLeftShift(result+1,num3.Mantissa,CHM-1);//一次左移更干净
//num3.Mantissa[0]+=10*result[0];//???
vect.BCDRightShift(result,CHM+1);
vect.Equ(num3.Mantissa,result,CHM);
}
delete []result;
return;
}
else//减法
{
add.SEC();
for(i=shift+CHM;i>=0;i--)
if(i>=CHM) result=add.SBC8(0,num[i-shift]); //有几步是多余的:
//至少不应该保存起来(没有多余，有可能前面几位为零)
//即使这种情况没有，但可以为后续的开发利用
else if(i<shift) result=add.SBC8(num3.Mantissa,0);
else
result=add.SBC8(num3.Mantissa,num[i-shift]);
for(i=0;i<2*CHM;i++)//因为有可能被减数与减数的前几位是相等的
if(result==0) {num3.Exponenet-=2;continue;}
else
break;
if(i>2*CHM-1) {vect.ClrFP(num3);return;}//如果结果等于零，返回(这一步可以在函数开始处进行判断)
if(result>9) vect.Equ(num3.Mantissa,result+i,CHM);
else
{
// vect.BCDLeftShift(result+i+1,num3.Mantissa,CHM-1);//???????????
// num3.Mantissa[0]+=result*10;
num3.Exponenet--;
vect.BCDRightShift(result+i,CHM+1);
vect.Equ(num3.Mantissa,result+i,CHM);
}
}
delete []result;
};
void FPMul(FP num1,FP num2,FP& num3)
{
//首先应文盲判断是否为零
num3.Exponenet=num1.Exponenet+num2.Exponenet;//指数相加
num3.Sign=num1.Sign^num2.Sign;//符号异或
//尾数移位相加
int i,j;
short Rnum1[10][CHM+1];//一张乘法表:Rnum1[0]:进位标志
for(i=1;i<10;i++) Rnum1[0]=0;
for(i=0;i<=CHM;i++) Rnum1[0]=0;
AddSub8 add;//8位的BCD加法器
Vector vect;
for(i=1;i<=9;i++)//计算被乘数的1-9倍，以便于查表计算结果
{
add.CLC();
for(j=CHM-1;j>=0;j--)
Rnum1[j+1]=add.ADC8(num1.Mantissa[j],Rnum1[i-1][j+1]);
Rnum1[0]=add.ADC8(0,Rnum1[i-1][0]);
}
//为了减少移位，我们分为二个结果，只需移一次即可:移8位运算(0,2,4,6,8;1,3,5,7,9)
short Result1[2*CHM],Result2[2*CHM];//这里定义的数过于大
short Address[2*CHM];
vect.BCD8ToBCD4(num2.Mantissa,CHM,Address);//为了便于查表(可能在机器上实现起来更节约时间)：
vect.ClrShort(Result1,2*CHM);
vect.ClrShort(Result2,2*CHM);
for(i=CHM-1;i>=0;i--)
{
add.CLC();
if(Address[i*2]!=0)
for(j=CHM;j>=0;j--)
{
Result1[j+i]=add.ADC8(Result1[j+i],Rnum1[Address][i*2]][j]);
}
add.CLC();
if(Address[i*2+1]!=0)
for(j=CHM;j>=0;j--)
{
Result2[j+i]=add.ADC8(Result2[j+i],Rnum1[Address][i*2+1]][j]);
}
}
//将Result2右移一位，然后两数相加
short Result3[2*CHM+1];
vect.ClrShort(Result3,2*CHM+1);
vect.BCDLeftShift(Result2,Result3,2*CHM);
add.CLC();
for(i=2*CHM-1;i>=0;i--)
Result2=add.ADC8(Result1,Result3);
if(Result2[0]==0) vect.Equ(num3.Mantissa,Result2+1,CHM);//进行归一化时可以编一个子程序
else
{
// vect.BCDLeftShift(Result2+1,num3.Mantissa,CHM-1);//???????????:倒不如一次左移来得干净
// num3.Mantissa[0]+=Result2[0]*10;
num3.Exponenet++;
vect.BCDRightShift(Result2,CHM+1);
vect.Equ(num3.Mantissa,Result2,CHM);
}
}
bool FPDiv(FP num1,FP num2,FP & num3)
{
/* 除法运算：符号异或，阶数相减：尾数
(1.先检查除数的所谓的有效长度；
2.两数相减，判断进位标志
3.保存商
)采取的运算法则是基于人工的算法，精度可以任意长，我们作了一张除数的倍数表*/
num3.Exponenet=num1.Exponenet-num2.Exponenet;
num3.Sign=num1.Sign^num2.Sign;
short L=CHM;//表示除数的长度：
int i,j,k;
for(i=CHM-1;i>=0;i--)
if(num2.Mantissa==0) {L--;continue;}
else
break;
if(L==0) {num3.Exponenet=120;return false;}//除数为零
short *Dividend,*Divisor[10];//被除数，除数
Dividend=new short[L+1];
for(i=0;i<10;i++) Divisor=new short[L+1];
Vector vect;
for(i=0;i<L+1;i++) Divisor[0]=0;
AddSub8 add;//8位的BCD加法器
for(i=1;i<=9;i++)//计算除数的1-9倍，以便于查表计算结果
{
add.CLC();
for(j=L-1;j>=0;j--)
Divisor[j+1]=add.ADC8(num2.Mantissa[j],Divisor[i-1][j+1]);
Divisor[0]=add.ADC8(0,Divisor[i-1][0]);
}
//被除数初始化
Dividend[0]=0;
vect.Equ(Dividend+1,num1.Mantissa,L);
short Result[CHM+1],m;
for(i=0;i<CHM+1;i++) //修改i的上限值可以动态调节精度
{
for(m=0;m<2;m++)//BCD的高低字节
{
for(j=9;j>0;j--)
{
add.SEC();//进位标志置1
for(k=L;k>=0;k--)
add.SBC8(Dividend[k],Divisor[j][k]);
if(add.c==true) break;
}
if(j!=0)
for(k=L;k>=0;k--)
Dividend[k]=add.SBC8(Dividend[k],Divisor[j][k]);//更新被除数
if(m==0) Result=j*10;//保存结果
else
Result+=j;
vect.BCDRightShift(Dividend,L+1);//被除数左移一位
if(L+i<CHM)
{
if(m==0) Dividend[L]+=num1.Mantissa[L+i]/10;
else
Dividend[L]+=num1.Mantissa[L+i]%10;
}
}
}
if(Result[0]>9) vect.Equ(num3.Mantissa,Result,CHM);
else
{
num3.Exponenet--;
// vect.BCDLeftShift(Result+1,num3.Mantissa,CHM);//左移一次可以节约时间:有一个问题：Mantissa溢出
// num3.Mantissa[0]+=Result[0]*10;
vect.BCDRightShift(Result,CHM+1);
vect.Equ(num3.Mantissa,Result,CHM);
}
delete []Dividend;
for(i=0;i<10;i++) delete []Divisor;
return true;
}
//虽然底层的程序都做好了，但是有一点：我做了两个模块：
//左移、右移；其实如果我考虑得周密一点的话，左移右移是一回事
void Cchebft(double a,double b,double *c,int N,int math)
{
int k,j;
double fac,bpa,bma,*f,y,sum;
f=new double[N];
for(j=0;j<N;j++) f[j]=0;
bma=0.5*(b-a);
bpa=0.5*(b+a);
for(k=0;k<N;k++)
{
y=cos(PI*(k+0.5)/N);
switch(math)
{
case 1:
f[k]=pow(10,y*bma+bpa); break;
case 2:
f[k]=atan(y*bma+bpa)/(y*bma+bpa);break;
case 3:
f[k]=sin(y*bma+bpa)/(y*bma+bpa);break;
case 4:
f[k]=log10(y*bma+bpa+1);break;
case 5:
f[k]=acos(y*bma+bpa) ;break;
case 6:
f[k]=log(y*bma+bpa);break;
case 7:
f[k]=log10(y*bma+bpa)/(y*bma+bpa-1);break;
case 8:
f[k]=asin(y*bma+bpa)/(y*bma+bpa);break;
case 9:
f[k]=acos(y*bma+bpa);break;
}
}
fac=2.0/N;
for(j=0;j<N;j++)
{
sum=0;
for(k=0;k<N;k++)
sum+=f[k]*cos(PI*j*(k+0.5)/N);
c[j]=fac*sum;
}
}
void Cchebpc(double *c,double *q,int N)
{
int k,j;
double sv,*Q1;
Q1=new double [N];
for(j=0;j<N;j++)
{q[j]=0;Q1[j]=0;}
q[0]=c[N-1];
for(j=N-2;j>=1;j--)
{
for(k=N-j;k>=1;k--)
{
sv=q[k];
q[k]=2.0*q[k-1]-Q1[k];
Q1[k]=sv;
}
sv=q[0];
q[0]=-Q1[0]+c[j];
Q1[0]=sv;
}
for(j=N-1;j>=1;j--) q[j]=q[j-1]-Q1[j];
q[0]=-Q1[0]+0.5*c[0];
delete []Q1;
}
void Cpcshift(double a,double b,double *coeff,int N)
{
int k,j;
double fac,cnst;
cnst=2.0/(b-a);
fac=cnst;
//for(j=0;j<N;j++) coeff[j]=1;
for(j=1;j<N;j++)
{
coeff[j]*=fac;
fac *=cnst;
}
cnst=0.5*(a+b);
for(j=0;j<=N-2;j++)
for(k=N-2;k>=j;k--)
coeff[k]-=cnst*coeff[k+1];
}
void PrintfVect(double *c,int n)
{
for(int i=0;i<n;i++) printf("\n%.12e",c);
}
double cddpoly(double x,double *c,int n)
{
double p=0;
p=c[n-1];
for(int j=n-2;j>=0;j--)
p=p*x+c[j];
return p;
}
FP poly(FP x,FP *c,int n)
{
FP p;
Vector vect;
vect.ClrFP(p);
p=c[n-1];
for(int j=n-2;j>=0;j--)
{ p=p*x+c[j];
}
// printf("\n%.11e\n",0.11111111111111111111);
return p;
}
FP chebyshev(double a,double b,int N,int n,FP x)
{
double c[40],q[40];
Cchebft(a,b,c,N,n);
Cchebpc(c,q,N);
Cpcshift(a,b,q,N);
Vector vect;
FP fp[40];
for(int i=0;i<N;i++)
vect.Double2FP(q,fp);
return poly(x,fp,N);
}
FP pow10(FP x)
{
FP p;
Vector vect;
vect.ClrFP(p);
short k=0;
//判断是否溢出:1,0,<0(指数存在的情况)
if(x.Exponenet>1 && x.Sign==false) {p.Mantissa[0]=10;p.Exponenet=100;return p;}
if(x.Exponenet>1 && x.Sign==true) {p.Mantissa[0]=10;p.Exponenet=-100;return p;}
if(x.Exponenet==1)
{
k=x.Mantissa[0];
x.Mantissa[0]=0;
}
if(x.Exponenet==0)
{
k=x.Mantissa[0]/10;
x.Mantissa[0]=x.Mantissa[0]%10;
}
vect.Norm(x);//进行归一化处理
if(x.Exponenet==0) //是否是整数
{
p.Mantissa[0]=10;
if(x.Sign==true) p.Exponenet=k*(-1);
else
p.Exponenet=k;
return p;
}
if(x.Sign==true) { p.Mantissa[0]=10;x=p+x;k=k*(-1)-1;}
double a,b;
int N=12;
if(x.Mantissa[0]/10>5 && x.Exponenet==-1) {a=0.5;b=1;} else {a=0;b=0.5;}
// a=0;
// if(x.Exponenet==-1) a=x.Mantissa[0]/10;
// a=a*0.1;
// b=a+0.1;
p=chebyshev(a,b,N,1,x);
p.Exponenet=k;
return p;
}
FP arctan(FP x)
{
if(x.Mantissa[0]==0) return x;
FP p,pi,one;
bool flag1=false,flag2=false,flag3=false;;
Vector vect;
vect.ClrFP(one);
one.Mantissa[0]=10;//置一
vect.Double2FP(PI*0.5,pi);
if(x.Exponenet>=0) {x=one/x;flag2=true;}
if(x.Sign==true) {x.Sign=false;flag1=true;}
double a,b;
int N=12,u;
if(x.Exponenet<-1) u=0;
else
u=x.Mantissa[0]/25;//不对
if(x.Exponenet==0) u=3;//x=1:可以直接赋值退出子程序
a=0.25*u;
b=a+0.25;
p=chebyshev(a,b,N,2,x);//x.exponenet<-2进行修正
if(x.Exponenet<-2)
flag3=true;
p=p*x;
if(flag2==true && flag1==true) p=p-pi;
else
if(flag2==true && flag1==false) p=pi-p;
else if(flag1==true) p.Sign=true;
FP w;
if(flag3==true)//进行修正:对末位有效数字修正
{
vect.ClrFP(w);
w.Mantissa[0]=10;
w.Sign=p.Sign;
w.Exponenet=p.Exponenet-CHM*2+1;
p=p+w;
}
return p;
}
FP Sin(FP x)//-pi/2----pi/2
{
FP p,pi4;
bool flag=false;
if(x.Sign==true) {x.Sign=false;flag=true;}//绝对值
double a,b;
int N=12;
a=0;b=PI/4;
Vector vect;
vect.Double2FP(PI/4.00,pi4);
pi4=x-pi4;
if(pi4.Sign==false) {a=PI/4;b=PI/2;}
p=chebyshev(a,b,N,3,x);
if(flag==true) p.Sign=true;
return p*x;
}
FP Cos(FP x)//0---pi
{
FP pi2;
Vector vect;
vect.Double2FP(PI/2,pi2);
x=pi2-x;
return Sin(x);
}
FP Tan(FP x)
{
return Sin(x)/Cos(x);//Cos!=0
}
FP Sqrt(FP x)
{
FP p,p1;
double k=0.5;
Vector vect;
vect.Double2FP(k,p1);
vect.ClrFP(p);
if(x.Mantissa[0]==0) return p;
p.Mantissa[0]=10;p.Exponenet=x.Exponenet/2;
for(int i=0;i<15;i++)//迭代法:这里应该不断判断两次的结果是否已经收敛
{
p=p+x/p;
p=p*p1;
}
return p;
}
FP arcsin(FP x)
{
FP p,p1;
Vector vect;
bool flag=false,flag1=false;
if(x.Sign==true) {flag1=true;x.Sign=false;}
p=x;
if(x.Exponenet>=0) {vect.ClrFP(p);p.Exponenet=110;return p;}
if(x.Exponenet==-1 && x.Mantissa[0]>50)
{flag=true;
vect.ClrFP(p1);
vect.Double2FP(0.5,p1);
x=Sqrt((p1-p1*x));
}
double a=0,b=0.25;
int N=12;
if(x.Exponenet==-1 && x.Mantissa[0]>25) {a=0.25;a=0.5;}
p=chebyshev(a,b,N,8,x)*x;
if(flag==true)
{
vect.ClrFP(p1);p1.Mantissa[0]=20;
p=p*p1;//p=p*2
vect.Double2FP(PI/2,p1);
p=p1-p;
}
if(flag1==true) p.Sign=true;
return p;
}
FP arccos(FP x)
{
FP p,p1;
Vector vect;
bool flag=false,flag1=false;
if(x.Sign==true) {x.Sign=false;flag1=true;}
p=x;
if(x.Exponenet>=0 && x.Mantissa[0]!=0) {vect.ClrFP(p);p.Exponenet=110;return p;}//这里判断不是很准确：x=1;x=-1;x>1;x<1
if(x.Exponenet==-1 && x.Mantissa[0]>50)
{flag=true;
vect.ClrFP(p1);
vect.Double2FP(0.5,p1);
x=Sqrt((p1-p1*x));
}
double a=0,b=0.25;
int N=12;
if(x.Exponenet==-1 && x.Mantissa[0]>25) {a=0.25;a=0.5;}
p=chebyshev(a,b,N,5,x);
if(flag==true)
{
vect.ClrFP(p1);p1.Mantissa[0]=20;
p=p*p1;//p=p*2
vect.Double2FP(PI,p1);
p=p1-p;
}
if(flag1==true)
{vect.Double2FP(PI,p1);p=p1-p;}
return p;
}
FP exp(FP x)
{
FP p;
double log10E=0.4342944819032518;//x=0.5有点问题
FP q;
Vector vect;
vect.Double2FP(log10E,q);
x=x*q;
p=pow10(x);
return p;
}
FP Log10(FP x)//x必须为正：此程序没有判断
{
//log10(a*10^b)=log10(a)+b=log10(a*c)-log10(c)
Vector vect;
FP One,x1;
vect.ClrFP(One);One.Mantissa[0]=10;
x1=x-One;
//x1趋于零时，C语言计算的结果不一定是准确的
if(x1.Exponenet<-5)
return x1*(vect.Double2FP(6)+x1)/(vect.Double2FP(6.0)+vect.Double2FP(4.0)*x1)*vect.Double2FP(0.4342944819032518);
if(x1.Exponenet<-2)
return (((x1+vect.Double2FP(21))*x1+vect.Double2FP(30))*x1)/((vect.Double2FP(9)*x1+vect.Double2FP(36))*x1+vect.Double2FP(30))*vect.Double2FP(0.4342944819032518);
if(x.Exponenet==0 && x.Mantissa[0]<15) //1<=x<=1.5
return chebyshev(1.0,1.50,12,7,x)*(x-One);
if(x.Exponenet==-1 && x.Mantissa[0]>=80)// 0.8<=x<1
return chebyshev(0.8,1,12,7,x)*(x-One);
FP p,p1;
double c;
//vect.PrintfFp(x);
if(x.Mantissa[0]/10>=5) c=1.0;
else
// c=6.0-x.Mantissa[0]/10;
switch(x.Mantissa[0]/10)
{
case 1:c=5.0;break;
case 2:c=3.0;break;
case 3:c=2.0;break;
case 4: c=2.0;break;//???????
}
vect.Double2FP((double)x.Exponenet,p1);
vect.Double2FP(log10(c),p);
p=p1-p;
vect.Double2FP(c,p1);
x.Exponenet=0;
x=x*p1;
double a,b;
int N=8;
// if(x.Mantissa[0]/10==5) N=10;
//if(x.Mantissa[0]>75) {a=7.5;b=10;} else {a=5;b=7.5;}
// a=x.Mantissa[0]/10;b=a+1;
a=x.Mantissa[0]/10+((x.Mantissa[0]%10*10+x.Mantissa[1]/10)/25)*0.25;
b=a+0.25;
N=7;
p=chebyshev(a,b,N,4,x)+p;
return p;
}
FP Ln(FP x)
{
FP p;
double p4=0.4342944819032518;
Vector vect;
vect.Double2FP(1/p4,p);
return Log10(x)*p;
}
FP Pow(FP x,FP y)
{
FP p;
p=Log10(x);
//p.Mantissa[CHM-1]=0;
return pow10(y*p);//x的y次方
}
FP ain(FP in)
{
FileFP s1("c:\\COEF.txt");
double a,b,c[100],e=2.71828182845904;
double q[20]={
3.979400086720374e-001,
1.737177927611043e-001,
-3.474355855214964e-002,
9.264948969143740e-003,
-2.779484694009623e-003,
8.894344021844690e-004,
-2.964780808962932e-004,
1.016592169579166e-004,
-3.558111680579244e-005,
1.258718493447536e-005,
-4.530049685422902e-006,
1.851579095212842e-006,
-6.806462801319921e-007,0},q1[20]={0};
int N=13;
a=0;b=0.3;
FP nine;
Vector vect;
vect.ClrFP(nine);
nine.Exponenet=0;
nine.Mantissa[0]=25;
nine.Mantissa[1]=00;
//Cchebft(a,b,c,N,7);
//Cchebpc(c,q,N);
//Cpcshift(a,b,q,N);
double s=6;
FP input,fp[20];
vect.Double2FP(s,input);
for(int i=0;i<N;i++) { vect.Double2FP(q,fp);}
//for(i=N-1;i>=0;i--) s1.SaveCoAsm(fp);
s1.SaveCoAsm(vect.Double2FP(log(10)*6));
s1.SaveCoAsm(vect.Double2FP(log(10)*4));
exit(0);
//ENTER;
//vect.PrintfFp(input);
FP one;vect.ClrFP(one);one.Mantissa[0]=10;
//ENTER;
//vect.PrintfFp(poly(in,fp,N));
//printf("\n\n%.15e %.15e %.15e \n",s,cddpoly(s,q,N),log10(s));
in=in-nine;
return poly(in,fp,N);
}
void CheckDesign()
{
int i;
double x1,y1,e=2.71828182845904;
FP x2,y2;
Vector vect;
while(1)
{
system("cls");
printf("a:sin(x) b:cos(x) c:tan(x) d:log(x) e:ln(x)\n");
printf("f:asin(x) g:acos(x) h:atan(x) i:10^x j:e^x\n");
printf("k:x^y l:sqrt(x) x:Exit\n\n\n\n\n");
printf("\nChoice:");
i=getchar();
if(i=='x') return;
x1=0;y1=0 ;
printf("x=");scanf("%lf",&x1);vect.Double2FP(x1,x2);
x1=vect.FP2Double(x2);
if(i=='k')
{
printf("y=");scanf("%lf",&y1);vect.Double2FP(y1,y2);
printf("x=%.15e ",x1);vect.PrintfFp(x2);
printf("\ny=%.15e ",y1);vect.PrintfFp(y2);
}
else
{
printf("x=%.15e ",x1);vect.PrintfFp(x2);
}
switch(i)
{
case 'a':printf("\n%.15e\n",sin(x1));vect.PrintfFp(Sin(x2));break;
case 'b':printf("\n%.15e\n",cos(x1));vect.PrintfFp(Cos(x2));break;
case 'c':printf("\n%.15e\n",tan(x1));vect.PrintfFp(Tan(x2));break;
case 'd':printf("\n%.15e\n",log10(x1));vect.PrintfFp(Log10(x2));break;
case 'e':printf("\n%.15e\n",log(x1));vect.PrintfFp(Ln(x2));break;
case 'f':printf("\n%.15e\n",asin(x1));vect.PrintfFp(arcsin(x2));break;
case 'g':printf("\n%.15e\n",acos(x1));vect.PrintfFp(arccos(x2));break;
case 'h':printf("\n%.15e\n",atan(x1));vect.PrintfFp(arctan(x2));break;
case 'i':printf("\n%.15e\n",pow(10,x1));vect.PrintfFp(pow10(x2));break;
case 'j':printf("\n%.15e\n",pow(e,x1));vect.PrintfFp(exp(x2));break;
case 'k':printf("\n%.15e\n",pow(x1,y1));vect.PrintfFp(Pow(x2,y2));break;
case 'l':printf("\n%.15e\n",sqrt(x1));vect.PrintfFp(Sqrt(x2));break;
}
getch();getchar();getchar();getch();
}
}
void main()
{
Vector vect;
CheckDesign();
//return;
int num=9999;
double step=0,begin=2,end=3,x1;
FileFP s1("c:\\d.txt");
FileFP s2("c:\\f.txt");
FileFP s("c:\\input.txt");
FP x2;
step=(end-begin)/num;
for(int i=0;i<num;i++)
{
x1=begin+i*step; x2=vect.Double2FP(x1);
x1=vect.FP2Double(x2);
s.SaveFP(x2);
// s1.SaveDouble(pow(10,x1));
s1.SaveDouble(log10(x1));
s2.SaveFP(ain(x2));
}
}
|
楼主
标题置顶
标题高亮
