一、大数的运算算法
使用原子的STM32F4开发板制作一个计算器，类似电脑上的标准型计算器。
在制作界面之前，首先要解决大数运算问题。因为STM32为32位MCU，直接计算最大不过32位，按十进制也就是10位数，显然不够用。
为了方便计算，使用数组储存各位数据，并且低位在前。下面数据结构中，使用指针代替数组，这是为了方便以后与界面函数部分对接。
#define THIS_MAXLEN 50 //数组长度
#define SHOW_MAXLEN 30 //自我调整中使用的最大长度
#define DIV_PREC 20 //除法精度相关
#define RESULT_TYPE signed char
#define RESULT_OK 0
#define RESULT_ERR_MEM 1
#define RESULT_ERR_DIVZERO 2
#define RESULT_ERR_MINUS 3
#define RESULT_ERR_OVERLEN 4
typedef struct
{
signed char len;//数据长度
signed char *num;
} bignumBASE_st;//正整数结构
typedef struct
{
bignumBASE_st baseT;
signed char symbol; //0正，1负
signed char exponent;//指数，10^exp次方
} bignumESB_st;//包括浮点，正负
二、加法、减法
四则运算中，加法减法是最容易，甚至都不需要额外的临时大数组进行辅助。基本方法是从低位开始计算，有进位借位则记录，下一个循环计上。唯一需要注意的就是，类似a=a+1这种源参数与目的参数其实相同的情况，要避免影响。再一个，设计函数时，因为怕MCU栈空间不够，所以源参数也使用指针传递，而非值传递。其实这样做反而有些麻烦，需要手工做更多的保护处理，千万不要以为加了cons修饰t就没事了。
由于有符号加法与减法是相互调用转化的，为避免出现递归调用（虽然实际上应该不会），我对函数进行多层封装。
/*
说明：正整数加法 ,加法最内层
参数：pdst 存放结果，pa pb乘数
返回：无
*/
static void _bigADD(bignumBASE_st *const pdst,const bignumBASE_st *pa,const bignumBASE_st *pb)
{
unsigned char i=0,len = (pa->len)>(pb->len)?pa->len:pb->len;
signed char snum;
signed char cc=0; //进位标志
const bignumBASE_st *pl,*ps;
if( pa->len > pb->len )
{
pl=pa,ps=pb;
}else
{
pl=pb,ps=pa;
}
for(i=0; i< len; i++)
{ //避免数组溢出
snum=( i < ps->len ) ? ps->num[i] : 0 ;
pdst->num [i] = cc+ pl->num[i] + snum ; //加法
if( pdst->num[i] > 9 )
{
cc=1;
pdst->num [i]-=10;
}
else cc=0;
}
pdst->num [i] =cc;
pdst->len = len + cc;
}
//无符号大数加法，处理浮点
#define _BIG_ABS(a) (a<0?-a:a)
RESULT_TYPE nsbigADD(bignumESB_st *const pdst,const bignumESB_st *pa,const bignumESB_st *pb)
{
signed char ec=0;
ec=pa->exponent-pb->exponent ;
if( ec != 0)
{
const bignumESB_st *pt,*pn;//临时变量，pt需要调整指数
bignumESB_st ts;
ts.baseT.num = (signed char*)malloc(THIS_MAXLEN);//申请内存
if(ts.baseT.num == NULL)return RESULT_ERR_MEM;
ec>0 ?( pt=pa,pn=pb):(pt=pb,pn=pa);
bigNumTenfold( &ts.baseT,&pt->baseT, _BIG_ABS(ec) ); //移动数组与长度
ts.exponent = pn->exponent ; //对齐指数
_bigADD( & pdst->baseT , &ts.baseT, &pn->baseT ) ;
pdst->exponent = ts.exponent ;
}
else
{
_bigADD( &(pdst->baseT), &(pa->baseT), &(pb->baseT)) ;
pdst->exponent=pa->exponent;
}
return RESULT_OK ;
}
//无符号大整数减法 ，大减去小, 减法最内层
static void _bigSUB(bignumBASE_st *pdst,const bignumBASE_st *pa,const bignumBASE_st *pb)
{
unsigned char i=0 ,cflag=0 ,len;
signed char cc=0 ; //借位标志
len= pa->len ;
for(i=0; i< len; i++)
{
signed char bnum=0,danum=0;
bnum = (i < pb->len )?pb->num[i]:0;//避免地址溢出
danum = pa->num[i] - cc;
cc=0;
if( danum > bnum )
{
pdst->num[i] = danum - bnum ;
cflag=i;
}
else if( danum == bnum )
{
pdst->num[i] = 0;
}
else
{
cc=1; //借位
pdst->num[i] = 10 + danum - bnum;//借10
if( pdst->num[i] != 0)cflag=i;
}
}
pdst->len = cflag+1 ;
}
//无符号大数减法,处理浮点
RESULT_TYPE nsbigSUB(bignumESB_st *const pdst,const bignumESB_st *pa,const bignumESB_st *pb)
{
signed char ec=0;
signed char res;
ec = pa->exponent-pb->exponent ;//比较浮点长度
if( ec != 0)
{
const bignumESB_st *pt,*pn; //pt为需要调整指数的一方
const bignumESB_st *pxa,*pxb;
bignumESB_st ts;//用于调整指数
ts.baseT.num =(signed char*) malloc(THIS_MAXLEN);//申请内存
if(ts.baseT.num == NULL)return RESULT_ERR_MEM;
ec>0 ?( pt=pa,pn=pb):(pt=pb,pn=pa); //浮点长的一方为pt
bigNumTenfold( &(ts.baseT),&(pt->baseT), _BIG_ABS(ec) ); //移动数组与长度
ts.exponent = pn->exponent ; //对齐指数
if( ec>0 )
{
pxa= &ts,pxb=pb;
}else
{
pxa=pa,pxb=&ts;
}
res = arrcmp( &pxa->baseT , &pxb->baseT );//计较大小
if( res == -1)
{
_bigSUB(&pdst->baseT,&pxb->baseT,&pxa->baseT);
pdst->symbol=1;
pdst->exponent = ts.exponent ;//指数与多多的一方相同
}
else if( res == 0)
{
pdst->baseT.num[0]=0;
pdst->baseT.len=1 ;
pdst->symbol=0;
pdst->exponent = 0;
}
else if( res == 1)
{
_bigSUB(&pdst->baseT,&pxa->baseT ,&pxb->baseT );
pdst->symbol=0;
pdst->exponent = ts.exponent ;//指数与多的一方相同
}
}
else
{
res = arrcmp( & (pa->baseT) , & (pb->baseT) );//比较大小
if( res == -1)
{
_bigSUB(&pdst->baseT,&pb->baseT ,&pa->baseT );
pdst->symbol=1;
pdst->exponent = pa->exponent ;
}
else if( res == 0)
{
pdst->baseT.num[0]=0;
pdst->baseT.len=1 ;
pdst->symbol=0;
pdst->exponent = 0 ;
}
else if( res == 1)
{
_bigSUB(& pdst->baseT,& pa->baseT,& pb->baseT);
pdst->symbol=0;
pdst->exponent = pa->exponent ;
}
}
return RESULT_OK;
}
上面几个内部加法减法函数都是为下面两个最外层的函数准备
//有符号大数加法
RESULT_TYPE bigNumADD(bignumESB_st *const pdst,const bignumESB_st *pa,const bignumESB_st *pb)
{
RESULT_TYPE err;
if( pa->symbol == pb->symbol)
{
err=nsbigADD( pdst, pa, pb);
pdst->symbol = pa->symbol ;
}else
{
if ( pa->symbol ==0 ) // 正
err=nsbigSUB( pdst, pa, pb);
else
err=nsbigSUB( pdst, pb, pa);
}
removeTailZero( pdst->baseT.num , &pdst->baseT.len , 0 ); //去掉多余的零
bignumst_adjust( pdst, SHOW_MAXLEN ); //对数据进行标准化调整
return err;
}//有符号大数减法
RESULT_TYPE bigNumSUB(bignumESB_st *const pdst,const bignumESB_st *pa,const bignumESB_st *pb)
{
RESULT_TYPE err;
if( pa->symbol == pb->symbol)
{
err=nsbigSUB( pdst, pa, pb);
if( pa->symbol == 1) pdst->symbol = ! (pdst->symbol);
}
else
{
err=nsbigADD( pdst, pa, pb);
pdst->symbol = pa->symbol;
}
removeTailZero( pdst->baseT.num , &pdst->baseT .len , 0 ); //去掉多余的零
bignumst_adjust( pdst, SHOW_MAXLEN ); //对数据进行标准化调整
return err ;
}
三、乘法
乘法相对加减法稍微复杂那么一点点，方法依然是从低位到高位算，并且会用到加法。
总的来说，就是先单乘，然后逐步将结果相加
/*
说明：单位整数与大数相乘 ,乘法最内层
参数：pdst 存放结果，pa 乘数所在 ，singledigit 单位乘数0-9，pa的长度
返回：整数位数，长度
*/
static signed char _bigMUL_single(signed char *pdst,const signed char *pa,const unsigned char singledigit,const unsigned char len)
{
//即使源参数与目的参数一样，也不会影响
unsigned char i=0 ,tlen=len;
signed char cc=0;//进位结果
if( singledigit == 0 )
{
pdst[0]=0;
return 1;
}
for(i=0; i<len; i++)
{
pdst[i] = pa[i]*singledigit + cc;
if( pdst[i] > 9)
{
cc = pdst[i]/10; //进位
pdst[i] = pdst[i]%10; //取模
}
else cc=0;
}
pdst[i]=cc;
cc=(cc==0)?0:1;
tlen+= cc ;
return tlen;
}
/*
说明：无符号大数乘法，处理指数浮点
参数：pdst 存放结果，pa pb 乘数所在
返回：错误类型
*/
static RESULT_TYPE nsbigMUL( bignumESB_st *const pdst, const bignumESB_st *pa,const bignumESB_st *pb)
{
unsigned char i=0,lenlong,lenshort;
bignumBASE_st temp,*ptemp=&temp;//保存单次结果
bignumESB_st tdst,*ptdst=&tdst;//开辟内存,避免因源参数与目的参数一样导致的BUG
const bignumESB_st *plong,*pshort;
//申请内存
ptemp->num =(signed char*) malloc(THIS_MAXLEN);
ptdst->baseT.num =(signed char*) malloc(THIS_MAXLEN);
if( ptemp->num==NULL || ptdst->baseT.num == NULL )return RESULT_ERR_MEM;
//指数处理
pdst->exponent = pa->exponent+pb->exponent ;
//直接相乘
if( pa->baseT.len + pb->baseT.len <=9 )
{
unsigned int val=0;
val =arr2int(pa->baseT.num, pa->baseT.len);
val *=arr2int(pb->baseT.num, pb->baseT.len);
pdst->baseT.len =int2arr( val, pdst->baseT.num );
return RESULT_OK;
}
//清零
bignumst_MEMinit( ptdst) ;
//调整，长乘以短
if( pa->baseT.len > pb->baseT.len ) plong=pa,pshort=pb;
else plong=pb,pshort=pa;
lenlong=plong->baseT.len ;
lenshort=pshort->baseT.len;
//第一次,省掉复制
ptdst->baseT.len = _bigMUL_single( ptdst->baseT.num, plong->baseT.num , pshort->baseT.num[0], lenlong) ;
for(i=1; i < lenshort ; i++)
{
ptemp->len = _bigMUL_single( ptemp->num, plong->baseT.num , pshort->baseT.num[i], lenlong) ;
bigNumTenfold( ptemp,ptemp,i) ;//10的i次方 ,即数组右移i位
_bigADD( &ptdst->baseT, &ptdst->baseT, ptemp) ;
}
//复制
bignumst_cpyto(pdst,ptdst);
return RESULT_OK ;
}
/* 说明：无符号大数乘法，处理指数浮点 参数：pdst 存放结果，pa pb 乘数所在 返回：错误类型 */RESULT_TYPE bigNumMUL(bignumESB_st *const pdst,const bignumESB_st *pa,const bignumESB_st *pb)
{
RESULT_TYPE err;
err=nsbigMUL( pdst, pa, pb); //真正的乘法
pdst->symbol = pa->symbol ^ pb->symbol ; //处理正负号
removeTailZero( pdst->baseT.num , &pdst->baseT.len , 0 ); //去掉多余的零
bignumst_adjust( pdst, SHOW_MAXLEN ) ; //数据标准化
return err ;
}
四、除法
上面的加减乘法，都不难，我花了一天。但除法我整整花了4天，实在可恶，并且方法十分粗暴，特别是浮点精度处理，由于我缺乏相关的数学理论知识，实在不知道除数与商的之间的精度关系，或许有空可以研究研究。但目前我只能尽量申请尽量多的内存，以保证精度。这对于一向对内存吝啬的我简直是不可忍受。不过我还是忍了，这就是成长。
除法过程模拟手工算式，从高位算起，这里要用到单乘函数，然后使用正整数比较函数，正整数减法函数；过程不复杂，就是要特别注意精度处理。
/*
说明：无符号大数除法，处理指数浮点
参数：pdst 存放结果，pa 被除数所在 ， pb 除数所在
返回：错误类型
*/
static RESULT_TYPE nsbigDIV(bignumESB_st *const pdst,const bignumESB_st *pa,const bignumESB_st *pb)
{
signed char div_preci= DIV_PREC < THIS_MAXLEN ? DIV_PREC : THIS_MAXLEN ;
const signed char tblen=pb->baseT.len ;
const signed char talen=pa->baseT.len ;
bignumESB_st tdst,*ptdst=&tdst;//开辟内存,避免因源参数与目的参数一样导致的BUG
//错误情况 除数=0
if( IsZERO( &pb->baseT ) == 0 ) return RESULT_ERR_DIVZERO;
//特殊情况=1，直接复制
if( pb->baseT.len == 1 && pb->exponent == 0 )
{
if( pb->baseT.num[0] == 1 )
{
bignumst_cpyto( pdst, pa );
return RESULT_OK;
}
}
//申请内存
ptdst->baseT.num =(signed char*) malloc(THIS_MAXLEN);
if( ptdst->baseT.num == NULL )return RESULT_ERR_MEM ;
bignumst_MEMinit(ptdst);//深度初始化结构体，赋代数值0
{
bignumESB_st tm,txa; //两个中间变量，一个单乘的结果，一个每次减法的余数
bignumESB_st *ptm=&tm ,*ptxa=&txa;
signed char i,j;
//申请内存
ptm->baseT.num =(signed char*) malloc(THIS_MAXLEN);
ptxa->baseT.num =(signed char*) malloc(THIS_MAXLEN);
if( ptm->baseT.num == NULL || ptxa->baseT.num == NULL )return RESULT_ERR_MEM ;
{//用被除数赋值给临时变量
signed char len;
len = (talen < tblen) ? talen:tblen ; //取其小者
for( i=0;i < len ; i++ )
ptxa->baseT.num[i] = pa->baseT.num[ talen - len + i ]; //除法是从高开始的
ptxa->baseT.len= len ;
}
ptm->baseT.len=tblen; //赋初值
for(i=0; i <= (signed char)talen- tblen + div_preci ; i++)
{
signed char res ;
for(j=1;j<=10;j++)
{
ptm->baseT.len = _bigMUL_single( ptm->baseT.num, pb->baseT.num, j, tblen ) ; // m = b * j
res = arrcmp( &ptm->baseT , &ptxa->baseT ) ; // m > a?
if( res > 0 ) //若大于
{
ptdst->baseT.num[ talen - tblen-i + div_preci ] = j-1 ; //结果，商
//退一步
//为下一轮准备被除数 ptxa
ptm->baseT.len = _bigMUL_single( ptm->baseT.num, pb->baseT.num, j-1, tblen ) ; // m = b * j-1
_bigSUB( &ptxa->baseT ,& ptxa->baseT, &ptm->baseT ) ; // a = a - m
bigNumTenfold( &ptxa->baseT, &ptxa->baseT, 1) ;// 递推进一位
if( talen - tblen -1 >= i)
ptxa->baseT.num[0] = pa->baseT.num[ talen - tblen - i-1] ;// 头填pa ,next
else
{
ptxa->baseT.num[0] = 0;
}
break;
}
else if( res == 0 )//若相等
{
ptdst->baseT.num[ talen - tblen - i + div_preci ] = j; //结果，商
//为下一轮准备被除数 ptxa
if( talen - tblen -1 >= i)
ptxa->baseT.num[0] = pa->baseT.num[ talen - tblen - i-1] ;// 头填pa ,next
else ptxa->baseT.num[0] = 0;
ptxa->baseT.len =1;
break;
}
} ///end for(j;;)
if( i > talen - tblen )
{
if( IsZERO( &(ptxa->baseT) ) == 0 )
{
break;//除尽了，跳出
}
}
} ///end for(i;;)
{
signed char dstlen;
dstlen = talen - tblen + 1 + div_preci;
ptdst->baseT.len = dstlen ;
dstlen = (talen < tblen)?tblen - talen:0;
ptdst->exponent = (pa->exponent) - (pb->exponent) - div_preci + dstlen ;//
}
}
//复制
bignumst_cpyto(pdst,ptdst);
return RESULT_OK ;
}
五、开方
开方我使用的算法是牛顿逼近法，网上很多，总结过来就是一条公式：xn+1 = ( xn + c/xn )/2 ；其中C是需要开方的数；
首先针对上面公式中出现的除以2，对之前的除法进行精简
//除以个位数，single :1-9
static RESULT_TYPE bigNumDIV_single(bignumESB_st *const pdst,const bignumESB_st *pa,const unsigned char single)
{
signed char div_preci= DIV_PREC < THIS_MAXLEN ? DIV_PREC : THIS_MAXLEN ;
const signed char talen=pa->baseT.len ;
bignumESB_st tdst,*ptdst=&tdst;//开辟内存,避免因源参数与目的参数一样导致的BUG
//特殊情况=1，直接复制
if( single == 1 )
{
bignumst_cpyto( pdst, pa );
return RESULT_OK;
}
//符号
pdst->symbol=pa->symbol;
//申请内存
ptdst->baseT.num =(signed char*) malloc(THIS_MAXLEN);
if( ptdst->baseT.num == NULL )return RESULT_ERR_MEM ;
bignumst_MEMinit(ptdst);//深度初始化结构体，赋代数值0
{
signed char i,j;
signed char tm=0,txa=pa->baseT.num[ talen -1 ];//两个中间变量，一个单乘的结果，一个每次减法的余数
for(i=0; i <= (signed char)talen- 1 + div_preci ; i++)
{
signed char res ;
for(j=1;j<=10;j++)
{
tm = single * j ; // m = b * j
res = tm-txa ; // m > a?
if( res > 0 ) //大于
{
ptdst->baseT.num[ talen - 1-i + div_preci ] = j-1 ; //结果，商
//退一步
//为下一轮准备被除数 ptxa
tm = single * (j-1) ; // m = b * j-1
txa= txa - tm ; // a = a - m
txa=txa*10 ;// x10
if( talen - 1 -1 >= i)
txa += pa->baseT.num[ talen - 1 - i-1] ;// 头填pa ,next
else txa+=0;
break;
}
else if( res == 0 )//相等
{
ptdst->baseT.num[ talen - 1 - i + div_preci ] = j; //结果，商
//为下一轮准备被除数 ptxa
if( talen - 1 - 1 >= i)
txa = pa->baseT.num[ talen - 1 - i-1] ;// 头填pa ,next
else txa=0;
break;
}
} ///end for(j;;)
if( i > talen -1 )
{
if( txa == 0 )
{
//ptdst->baseT.len = i+1;
//ptdst->exponent = (pa->exponent) + i+1 - talen ;//
//goto lab_return;//除尽了，跳出
break;
}
}
} ///end for(i;;)
{
ptdst->baseT.len = talen + div_preci ;
ptdst->exponent = (pa->exponent) - div_preci ;//
}
}
//去掉无效的0
removeTailZero( ptdst->baseT.num , &ptdst->baseT.len , 0 ); //去掉多余的零
//调整指数
bignumst_adjust( ptdst ,SHOW_MAXLEN) ;
//复制
bignumst_cpyto(pdst,ptdst);
return RESULT_OK ;
}
这里固定迭代次数，实际中应该通过比较精度，然后决定是否继续迭代。//根号2
static RESULT_TYPE nsbigNumSQR(bignumESB_st *const pdst,const bignumESB_st *psrc)
{
//特殊情况=0,1,4,9
if( psrc->baseT.len == 1 && psrc->exponent == 0)
{
if( psrc->baseT.num[0] == 0 || psrc->baseT.num[0] == 1)
{
bignumst_setval( pdst, psrc->baseT.num[0] );
return RESULT_OK;
}
else
if( psrc->baseT.num[0] == 4 )
{
bignumst_setval( pdst, 2 );
return RESULT_OK;
}else
if( psrc->baseT.num[0] == 9 )
{
bignumst_setval( pdst, 3 );
return RESULT_OK;
}
}
{
bignumESB_st tsrc,*ptsrc=&tsrc;//开辟内存,避免因源参数与目的参数一样导致的BUG
bignumESB_st tm,*ptm=&tm;
//申请内存
ptsrc->baseT.num =(signed char*) malloc(THIS_MAXLEN);
ptm->baseT.num =(signed char*) malloc(THIS_MAXLEN);
if( ptsrc->baseT.num == NULL || ptm->baseT.num == NULL )return RESULT_ERR_MEM ;
//深度初始化结构体，赋代数值0
bignumst_MEMinit(pdst);
bignumst_MEMinit(ptm);
bignumst_cpyto(ptsrc,psrc);//复制源至临时变量
//给牛顿法的初始值
{
signed char zero = ptsrc->exponent + ptsrc->baseT.len;
bignumst_cpyto(pdst,ptsrc);
// zero 区间[0.1~10)不变
if( zero > 1 ) // zero是整数部分的长度
{
pdst->exponent -= (zero-1)/2 + (zero+1)%2;
}
else if( zero < 0 )//0.0nnnnn,-zero就是小数点后连续0的个数
{
pdst->exponent -= (zero-1)/2 ;
}
bignumst_valprec(pdst,7);//取精度
}
//牛顿迭代
{
signed char i=0;
for(;i<12;i++)
{
nsbigDIV( ptm , ptsrc, pdst );
bignumst_valprec(ptm,7);//取精度,过长反而影响除法结果
nsbigADD( pdst, ptm , pdst ) ;
bignumst_valprec(pdst,8);//取精度
bigNumDIV_single(pdst,pdst,2) ;
bignumst_valprec(pdst,7);//取精度
bignumst_adjust( pdst, SHOW_MAXLEN );//
//ShowBignum( pdst ) ;printf("x/2");
}
}
}
return RESULT_OK ;
}
//根号2
RESULT_TYPE bigNumSQR(bignumESB_st *const pdst,const bignumESB_st *psrc)
{
RESULT_TYPE err;
if( psrc->symbol == 1)return RESULT_ERR_MINUS;
err=nsbigNumSQR( pdst, psrc);
return err;
}
六、测试
VS2010测试
代码下载：http://download.csdn.net/download/wangzibigan/9997298
扫一扫
441
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
