原创 单片机模拟笔算实现超大整数相乘

完整文档和Keil演示程序：
https://static.assets-stash.eet-china.com/album/old-resources/2009/7/14/670e726a-4f6c-4793-ab08-9e11712ab07d.zip

单片机模拟笔算实现超大整数相乘

目录
第1章 引言 2
第2章 详细设计过程 3
2.1 简化计算公式 3
2.2 模拟笔算的思路 3
2.3 C语言代码具体设计 3
2.4 仿真调试祥解 8
2.5 后记 9
第1章 引言
从单片机汇编时代走过来的老鸟们，也许至今回想起数值计算的时候，仍然会忍不住的发抖……^_^
然而，即使到了今天，可以方面的利用C51丰富的数学库函数来实现复杂的计算，我们仍然不时感到不爽。毕竟C51里只支持float型，7个有效位，引入的误差有时候实在令人难以容忍。什么，你大声反对AMO，说你知道可以在C51里定义double型浮点数变量？呵呵，你肯定被编译器骗了，不信你看看它的手册或者仿真一下，你会发现，它把double通通当作float了^_^
在以前的一段时间里，www.18ic.com出现了一阵子DDS热潮，很多网友在做毕业设计的时候选了DDS，然后用百度搜索到了18ic……也有几个网友干脆发email来问AMO，然而身为斑竹的 AMO多少显得有点不太称职……在最近的一个工程项目里，AMO终于忍无可忍了，决定彻底解决这个舍入误差问题。这是一个DDS（直接数字频率合成器）的应用模块，DDS芯片是AD9851。AD9851的频率字为32位，2的32次方等于4294967296，外接晶振频率为50MHz，那么频率分辨率为：50000000/(4294967296-1)=0.012Hz。AMO希望可以直接输入一个频率值F，然后通过公式F* (4294967296-1)/50000000计算出频率字，然后将频率字写入AD9851，从而输出期望频率的信号（这就是那几个网友发email来问的问题……不要朝AMO扔鸡蛋西红柿，呵呵）。愿望是美好的，实际上在C51里却是不可实现的，因为计算过程带来的舍入误差太大了。所以AMO绞尽脑汁想出了一种不会带来误差的算法：模拟笔算。

下面以51系列单片机＋Keil来阐述具体设计。

关键点：模拟笔算
第2章 详细设计过程
2.1 简化计算公式
实际上，大部分的计算过程都可以简化为加法和乘法。比如前言中计算AD9851频率字的公式F*(4294967296-1)/50000000就可以简化为：
F*4294967295/50000000
并且进一步简化为：F* 85.8993459
这样我们只需要做一次乘法运算就可以了。

2.2 模拟笔算的思路
模拟笔算，顾名思义，就是用程序来模拟我们手工运算的过程，可以想象，这个过程是不会带来额外的误差的（除非你的程序有问题^_^）。

两整数相乘，一般手工笔算过程如下：
         1234
      ×  121
      -----------
         1234
        2468
       1234
      -----------
       149314
用C语言模拟手工笔算实现大整数相乘，相对来说比较简单。做法是：
1． 先用sprintf()函数将乘数格式化成ASCII字符串；
2． 将ASCII字符串转换为十进制数，这样，数组中每个元素存储1位十进制数据；以上两步也可以改为直接构造十进制数组；
3． 通过循环结构来模拟手工笔算过程；
4． 根据实际需要对运算结果进行处理，如转换为变量，转换为显示码输出到终端上等等。
2.3 C语言代码具体设计
/******************************************************************************
* 函数原型： void BigNumProduct(char x[],char y[],char z[],int xLen,int yLen)
* 功能描述： 超大整数相乘：z=x*y
* 入口参数： 字符型数组x,y,z以及x,y的长度；存放格式为10进制，按日常习惯方向存放，
             如：x存放"2197"，则x[0]=2,x[1]=1,x[2]=9,x[3]=7。
             注意：存放的必须是十进制数0－9，而不是ascii字符'0'-'9'！
* 出口参数： 无，直接将输出结果写到z中。
* 作    者： AMO
* 电子邮件： amo73@126.com
* 备    注： 要认真检查数组的长度，确保不越界；为了便于理解，没有对代码进行优化
******************************************************************************/
void BigNumProduct(char x[],char y[],char z[],int xLen,int yLen)
{
    int i,j,k; //AMO提问：为什么是int，而不是unsigned int？好好想想^_^
    char temp1,temp2;
    for(i=0;i<xLen+yLen;i++)//清空用来存放结果的数组
    {
        z=0;
    }
    for(i=xLen-1;i>=0;i--)
    {
        for(j=yLen-1;j>=0;j--)
        {
            temp1 = x*y[j];      //相乘
            for(k=i+j+1;k>=0;k--)   //相加并移位
            {
                temp2 = z[k]+temp1; //相乘结果或者进位，与本位相加
                z[k] = temp2%10;    //取余，写入本位
                temp1 = temp2/10;   //取模，作为进位
                if(temp1 == 0)      //如果进位为0，则跳过
                {
                    break;
                }
            }
        }
    }
}

以上就是核心的算法（大整数相乘）了，应该能看明白吧？AMO觉得注释已经写得非常明白了^_^
接下来看看完整的示例代码吧：

＃i nclude <reg52.h>
＃i nclude <stdio.h>
＃i nclude <stdlib.h>
＃i nclude <string.h>

/******************************************************************************
用C语言模拟手工笔算实现大整数相乘
两整数相乘，一般手工笔算过程如下：
         1234
      ×  121
      ----------
         1234
        2468
       1234
      ----------
       149314
模拟该过程的C语言代码如下：
（整数在数组中正序存放，结果是正序输出的，数组中每个元素用来存放1个十进制位）
******************************************************************************/

/******************************************************************************
* 函数原型： void BigNumProduct(char x[],char y[],char z[],int xLen,int yLen)
* 功能描述： 超大整数相乘：z=x*y
* 入口参数： 字符型数组x,y,z以及x,y的长度；存放格式为10进制，按日常习惯方向存放，
             如：x存放"2197"，则x[0]=2,x[1]=1,x[2]=9,x[3]=7。
             注意：存放的必须是十进制数0－9，而不是ascii字符'0'-'9'！
* 出口参数： 无，直接将输出结果写到z中。
* 作    者： AMO
* 电子邮件： amo73@126.com
* 备    注： 要认真检查数组的长度，确保不越界；为了便于理解，没有对代码进行优化
******************************************************************************/
void BigNumProduct(char x[],char y[],char z[],int xLen,int yLen)
{
    int i,j,k; //AMO提问：为什么是int，而不是unsigned int？好好想想^_^
    char temp1,temp2;
    for(i=0;i<xLen+yLen;i++)//清空用来存放结果的数组
    {
        z=0;
    }
    for(i=xLen-1;i>=0;i--)
    {
        for(j=yLen-1;j>=0;j--)
        {
            temp1 = x*y[j];      //相乘
            for(k=i+j+1;k>=0;k--)   //相加并移位
            {
                temp2 = z[k]+temp1; //相乘结果或者进位，与本位相加
                z[k] = temp2%10;    //取余，写入本位
                temp1 = temp2/10;   //取模，作为进位
                if(temp1 == 0)      //如果进位为0，则跳过
                {
                    break;
                }
            }
        }
    }
}

/******************************************************************************
* 函数原型： unsigned long CalculateAD9851(unsigned long Frequency)
* 功能描述： 根据频率值（正整数）计算AD9851频率字（32位频率字）
* 入口参数： ulong型频率值
* 出口参数： ulong型32位频率字
* 作    者： AMO
* 电子邮件： amo73@126.com
* 备    注： 频率分辨率为1Hz，要更精确的话，请用CalculateAD9851Plus函数。
             为了便于理解，用了几个库函数，占用了较多资源，实际使用可以优化一下
******************************************************************************/
unsigned long CalculateAD9851(unsigned long Frequency)
{//晶振频率50MHz，不倍频；AD9851频率字为32位(2的32次方-1=4294967295)，则：
 //频率字=频率*4294967295/50000000=频率*85.8993459；
 //在此把85.8993459的小数点向右移7位，得到858993459（计算完以后要把小数点移回来）
    char x[10];//unsigned long型最大为4294967295（10位）
    char y[] = "858993459"; //9位
    char z[20] = {0};       //10+9=19,取整20
    int xLen,yLen;          //x,y的位数
    int i;
    xLen = sprintf(x,"%Ld",Frequency);//将频率转为ascii字符串，并返回长度
    yLen = strlen(y);       //求另一个乘数的长度
    i=0;
    while(x)
    {
        x = x-'0';    //ascii字符串转为10进制字符串
        i++;
    }
    i=0;
    while(y)
    {
        y = y-'0';    //ascii字符串转为10进制字符串
        i++;
    }
    BigNumProduct(x,y,z,xLen,yLen);//z=x*y，结果保存在数组z里存放格式与x,y相同
    for(i=0;i<xLen+yLen;i++)
    {
        z = z + '0';  //十进制字符串转为ascii字符串
    }
    z[xLen+yLen-7] = '\0';  //设置字符串结尾，实际上是屏蔽掉小数点后面的数据^_^
    return( atol(z) );      //提取出计算结果
}


/******************************************************************************
* 函数原型： void Init_Uart(void)
* 功能描述： 配置串行口为9600波特率
* 入口参数： 无
* 出口参数： 无
* 备注：     晶振为11.0592MHz；其余见代码注释
******************************************************************************/
void Uart_Init(void)
{
    SCON  = 0x40; //串行方式1，不允许串行接收 
    TMOD  = 0x20; //定时器T1,方式2
//**** 根据单片机时钟模式修改下面的代码 ****   
  //  TH1   = 0xFA; //6时钟模式
    TH1   = 0xFD; //12时钟模式
//******************************************
    TR1 = 1;  //启动定时器
    TI  = 1;  //发送第一个字符
}

void main(void)
{
    unsigned int i;
    unsigned long Frequency;//频率
    unsigned long Result,last;

    Uart_Init();//初始化串行口
    Frequency = 0;
 printf("\n清输入一个初始频率值：\n");//提示输入初始频率值
 REN=1;//允许串行接收
    scanf("%Ld",&Frequency);//等待从串口窗口输入起始频率值
 REN=0;//不允许串行接收
    while(1)
    {
        Result = CalculateAD9851(Frequency);//根据频率值，计算AD9851频率字
        //打印频率、对应的频率字以及频率字增量
        printf("\n%LdHz --->%LU;  +%Ld\n",Frequency,Result,Result-last);
        Frequency ++;
        last = Result;
        i=2;
        while(i--);//延时
    }
}

还是那句话，应该能看明白吧？AMO觉得注释已经写得非常明白了^_^
为什么里面没有加法？因为AMO觉得既然乘法你都可以看懂了，加法肯定不会有问题了^_^

2.4 仿真调试祥解
也许你手头没有硬件仿真器，没有烧录器，没有单片机开发板，没有直流电源，甚至连一只电阻都找不到……没关系，AMO已经想到这些问题了，所以AMO做好了一个可以直接使用的Keil工程（在压缩包里）。
1． 编译工程；
2． 启动调试（软件仿真模式）；
3． 打开串行窗口，如下图所示，点击一下红圈内的图标即可：
 
出来的界面如下：


4． 全速运行，程序将在scanf()函数那里等待你从串行窗口输入一个起始频率值（整数），此时切换到串行窗口，随便输入一个整数（比如说1），并敲一下回车键：
 
可以从串行窗口看出，频率每升高1Hz，频率字就会大约增加86。


2.5 后记
由于超大整数算法是用C语言编写，可以轻松应用到别的系统上（最低级的单片机都能实现了，更何况别的平台了^_^），所需要做的改动就是改变缓冲区长短了。当然了，AMO最主要的目的还是介绍一下这种模拟笔算的方法，大家完全可以编写出更实用更有效率的算法来^_^