发新帖我要提问
返回列表
打印

神奇的一篇算法**-大数阶乘

[复制链接]
2958|7
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
gaoyang9992006|  楼主 | 2016-6-14 13:50 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
计算器, 代码, 精度, 高精度, fft
大数阶乘
序
大数阶乘的计算是一个有趣的话题，从中学生到大学教授，许多人都投入到这个问题的探索和研究之中，并发表了他们自己的研究成果。如果你用阶乘作关键字在google上搜索，会找到许多此类文章，另外，如果你使用google学术搜索，也能找到一些计算大数阶乘的学术论文。但这些文章和论文的深度有限，并没有给出一个高速的算法和程序。
我和许多对大数阶乘感兴趣的人一样，很早就开始编制大数阶乘的程序。从2000年开始写第一个大数阶乘程序算起，到现在大约己有6-7年的时光，期间我写了多个版本的阶乘计算器，在阶乘计算器的算法探讨和程序的编写和优化上，我花费了很大的时间和精力，品尝了这一过程中的种种甘苦，我曾因一个新算法的实现而带来速度的提升而兴奋，也曾因费了九牛二虎之力但速度反而不及旧的版本而懊恼，也有过因解一个bug而通宵敖夜的情形。我写的大数阶乘的一些代码片段散见于互联网络，而算法和构想则常常萦绕在我的脑海。自以为，我对大数阶乘计算器的算法探索在深度和广度上均居于先进水平。我常常想，应该写一个关于大数阶乘计算的系列文章，一为整理自己的劳动成果，更重要的是可以让同行分享我的知识和经验，也算为IT界做一点儿贡献吧。
　
　 我的第一个大数阶乘计算器始于2000年，那年夏天，我买了一台电脑，开始在家专心学习VC，同时写了我的第一个VC程序，一个仿制windows界面的计算器。该计算器的特点是高精度和高速度，它可以将四则运算的结果精确到6万位以内，将三角、对数和指数函数的结果精确到300位以内，也可以计算开方和阶乘等。当时，我碰巧看到一个叫做实用计器的软件。值得称颂的是，该计算器的作者姜边竟是一个高中生，他的这个计算器功能强大，获得了各方很高的评价。该计算的功能之一是可以计算9000以内阶乘的精确值，且速度很快。在佩服之余，也激起我写一个更好更快的程序的决心，经过几次改进，终于使我的计算器在做阶乘的精确计算时 (以9000!为例)，可比实用计算器快10倍。而当精确到30多位有效数字时，可比windows自带的计算器快上7500倍。其后的2001年1月，我在csdn上看到一个贴子，题目是“有谁可以用四行代码求出1000000的阶乘”，我回复了这个贴子，给出一个相对简洁的代码，这是我在网上公布的第一个大数阶乘的程序。这一时期，可以看作是我写阶乘计算器的第一个时期。

　　我写阶乘计算器的第二个时期始于2003年9月，在那时我写了一组专门计算阶乘的程序，按运行速度来分，分为三个级别的版本，初级版、中级版和高级版。初级版本的算法许多人都能想到，中级版则采用大数乘以大数的硬乘法，高级版本在计算大数乘法时引入分治法。期间在csdn社区发了两个贴子，“擂台赛：计算n!(阶乘)的精确值，速度最快者2000分送上”和“擂台赛：计算n!(阶乘)的精确值，速度最快者2000分送上（续）”。其高级算的版本完成于2003年11月。此后，郭先强于2004年5月10日也发表了系列贴子，“擂台：超大整数高精度快速算法”、“擂台：超大整数高精度快速算法(续)”和“擂台：超大整数高精度快速算法(续2)”, 该贴重点展示了大数阶乘计算器的速度。这个贴子一经发表即引起了热列的讨论，除了我和郭先强先生外，郭雄辉也写了同样功能的程序，那时，大家都在持续改进自己的程序，看谁的程序更快。初期，郭先强的稍稍领先，中途郭子将apfloat的代码应用到阶乘计算器中，使得他的程序胜出，后期（2004年8月后）在我将程序作了进一步的改进后，其速度又稍胜于他们。在这个贴子中，arya提到一个开放源码的程序，它的大数乘法采用FNT＋CRT(快速数论变换＋中国剩余定理)。郭雄辉率先采用apflot来计算大数阶乘，后来郭先强和我也参于到apfloat的学习和改进过程中。在这点上，郭先强做得非常好，他在apfloat的基础上，成功地优化和改时算法，并应用到大数阶乘计算器上，同时他也将FNT算法应用到他的＜超大整数高精度快速算法库＞中，并在2004.10.18正式推出V3.0.2.1版。此后，我在2004年9月9日也完成一个采用FNT算法的版本，但却不及郭先强的阶乘计算器快。那时，郭先强的程序是我们所知的运算速度最快的阶乘计算器，其速度超过久负盛名的数学软件Mathematica和Maple,以及通用高精度算法库GMP。
　　
　　我写阶乘计算器的第三个时间约开始于2006年，在2005年8月收到北大刘楚雄老师的一封e-mail，他提到了他写的一个程序在计算阶乘时比我们的更快。这使我非常吃惊，在询问后得知，其核心部分使用的是ooura写的FFT函数。ooura的FFT代码完全公开，是世界上运行的最快的FFT程序之一，从这点上，再次看到了我们和世界先进水平的差距。佩服之余，我决定深入学习FFT算法，看看能否写出和ooura速度相当或者更快的程序，同时一个更大计划开始形成，即写一组包括更多算法的阶乘计算器，包括使用FFT算法的终极版和使用无穷级数的stirling公式来计算部分精度的极速版，除此之外，我将重写和优化以前的版本，力争使速度更快，代码更优。这一计划的进展并不快，曾一度停止。
　　
　　目前，csdn上blog数量正在迅速地增加，我也萌生了写blog的计划，借此机会，对大数阶乘之计算作一个整理，用文字和代码详述我的各个版本的算法和实现，同时也可能分析一些我在网上看到的别人写的程序，当然在这一过程中，我会继续编写未完成的版本或改写以前己经实现的版本，争取使我公开的第一份代码都是精品，这一过程可能是漫长的，但是我会尽力做下去。
菜鸟篇
程序1，一个最直接的计算阶乘的程序

#include "stdio.h"
#include "stdlib.h"

int main(int argc, char* argv[])
{
         long i,n,p;
         printf("n=?");
         scanf("%d",&n);
         p=1;
         for (i=1;i<=n;i++)
                  p*=i;
         printf("%d!=%d/n",n,p);
         return 0;
}

程序2，稍微复杂了一些，使用了递归，一个c++初学者写的程序

#include   <iostream.h>   
  long   int   fac(int   n);   
  void   main()   
  {   
          int   n;   
          cout<<"input   a   positive   integer:";   
          cin>>n;   
          long   fa=fac(n);   
          cout<<n<<"!   ="<<fa<<endl;   
  }   
  long   int   fac(int   n)   
  {   
          long   int   p;   
          if(n==0)   p=1;   
          else   
              p=n*fac(n-1);   
          return   p;   
  }   
  
程序点评，这两个程序在计算12以内的数是正确，但当n>12,程序的计算结果就完全错误了，单从算法上讲，程序并没有错，可是这个程序到底错在什么地方呢？看来程序作者并没有意识到，一个long型整数能够表示的范围是很有限的。当n>=13时，计算结果溢出，在C语言，整数相乘时发生溢出时不会产生任何异常，也不会给出任何警告。既然整数的范围有限，那么能否用范围更大的数据类型来做运算呢？这个主意是不错，那么到底选择那种数据类型呢？有人想到了double类型，将程序1中long型换成double类型，结果如下:
#include "stdio.h"
#include "stdlib.h"
int main(int argc, char* argv[])
{
   double i,n,p;
   printf("n=?");
   scanf("%lf",&n);
   p=1.0;
   for (i=1;i<=n;i++)
            p*=i;
   printf("%lf!=%.16g/n",n,p);
   return 0;
}
运行这个程序，将运算结果并和windows计算器对比后发现，当于在170以内时，结果在误差范围内是正确。但当N>=171,结果就不能正确显示了。这是为什么呢？和程序1类似，数据发生了溢出，即运算结果超出的数据类型能够表示的范围。看来C语言提供的数据类型不能满足计算大数阶乘的需要，为此只有两个办法。1.找一个能表示和处理大数的运算的类库。2.自己实现大数的存储和运算问题。方法1不在本文的讨论的范围内。本系列的后续文章将围绕方法2来展开。


使用特权

评论回复 赏

相关下载

相关帖子
沙发
gaoyang9992006|  楼主 | 2016-6-14 13:51 | 只看该作者
程序3.
#include "stdafx.h"
#include "math.h"
#define MAX_N 10000000.00          //能够计算的最大的n值，如果你想计算更大的数对数，可将其改为更大的值
#define MAX_MANTISSA   (1e308/MAX_N)    //最大尾数
struct bigNum
{
    double n1;     //表示尾数部分
     int n2;   //表示指数部分
  };
void calcFac(struct bigNum *p,int n)
{
     int i;
     double 　MAX_POW10_LOG=(floor(log10(1e308/MAX_N))); //最大尾数的常用对数的整数部分,
     double MAX_POW10=    (pow(10.00, MAX_POW10_LOG)); // 10 ^ MAX_POW10_LOG
        p->n1=1;
     p->n2=0;
     for (i=1;i<=n;i++)
     {
          if (p->n1>=MAX_MANTISSA)
          {
               p->n1 /= MAX_POW10;
               p->n2 += MAX_POW10_LOG;
          }
          p->n1 *=(double)i;
     }
}
void printfResult(struct bigNum *p,char buff[])
{
     while (p->n1 >=10.00 )
     {p->n1/=10.00; p->n2++;}
     sprintf(buff,"%.14fe%d",p->n1,p->n2);
}
int main(int argc, char* argv[])
{
     struct bigNum r;
     char buff[32];
     int n;
   
     printf("n=?");
     scanf("%d",&n);
     calcFac(&r,n);           //计算n的阶乘
     printfResult(&r,buff);   //将结果转化一个字符串
        printf("%d!=%s/n",n,buff);
     return 0;
}
以上代码中的数的表示方式为：数的值等于尾数乘以 10 ^ 指数部分，当然我们也可以表示为：尾数 乘以 2 ^ 指数部分，这们将会带来这样的好处，在调整尾数部分和指数部分时，不用除法，可以依据浮点数的格式直读取阶码和修改阶码(上文提到的指数部分的标准称呼)，同时也可在一定程序上减少误差。为了更好的理解下面的代码，有必要了解一下浮点数的格式。浮点数主要分为32bit单精度和64bit双精度两种。本方只讨论64bit双精度(double型)浮点数的格式，一个double型浮点数包括8个字节(64bit)，我们把最低位记作bit0,最高位记作bit63,则一个浮点数各个部分定义为:第一部分尾数：bit0至bit51，共计52bit，第二部分阶码:bit52-bit62,共计11bit，第三部分符号位:bit63，0:表示正数，1表示负数。如一个数为0.xxxx * 2^ exp,则exp表示指数部分，范围为－1023到1024，实际存储时采用移码的表示法，即将exp的值加上0x3ff，使其变为一个0到2047范围内的一个值。函数GetExpBase2 中各语句含义如下：1.“(*pWord & 0x7fff)”，得到一个bit48-bit63这个16bit数，最高位清0。2.“>>4”,将其右移4位以清除最低位的4bit尾数，变成一个11bit的数（最高位5位为零）。3(rank-0x3ff)”, 减去0x3ff还原成真实的指数部分。以下为完整的代码。
程序4：
#include "stdafx.h"
#include "math.h"
#define MAX_N 10000000.00      //能够计算的最大的n值，如果你想计算更大的数对数，可将其改为更大的值
#define MAX_MANTISSA   (1e308/MAX_N) //最大尾数
typedef unsigned short WORD;
struct bigNum
{
double n1;     //表示尾数部分
int n2;   //表示阶码部分
};
short GetExpBase2(double a) // 获得 a 的阶码
  {
// 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu
        WORD *pWord=(WORD *)(&a)+3;
    WORD rank = ( (*pWord & 0x7fff) >>4 );
    return (short)(rank-0x3ff);
}
double GetMantissa(double a) // 获得 a 的 尾数
{
// 按照IEEE 754浮点数格式，取得尾数，仅仅适用于Intel 系列 cpu
       WORD *pWord=(WORD *)(&a)+3;
*pWord &= 0x800f; //清除阶码
*pWord |= 0x3ff0; //重置阶码
return a;
}
void calcFac(struct bigNum *p,int n)
{
int i;
p->n1=1;
p->n2=0;
for (i=1;i<=n;i++)
{
if (p->n1>=MAX_MANTISSA)　//继续相乘可能溢出，调整之
      {
           p->n2 += GetExpBase2(p->n1);
           p->n1 = GetMantissa(p->n1);
      }
      p->n1 *=(double)i;
}
}
void printfResult(struct bigNum *p,char buff[])
{
double logx=log10(p->n1)+ p->n2 * log10(2);//求计算结果的常用对数
int logxN=(int)(floor(logx)); //logx的整数部分
sprintf(buff,"%.14fe%d",pow(10,logx-logxN),logxN);//转化为科学计算法形式的字符串
}
int main(int argc, char* argv[])
{
struct bigNum r;
char buff[32];
int n;
printf("n=?");
scanf("%d",&n);
calcFac(&r,n);           //计算n的阶乘
printfResult(&r,buff);   //将结果转化一个字符串
printf("%d!=%s/n",n,buff);
return 0;
}
程序优化的威力：
程序4是一个很好的程序，它可以计算到1千万（当n更大时，p->n2可能溢出）的阶乘，但是从运行速度上讲，它仍有潜力可挖，在采用了两种优化技术后，（见程序5），速度竟提高5倍，甚至超出笔者的预计。
第一种优化技术，将频繁调用的函数定义成inline函数，inline是C++关键字，如果使用纯C的编译器，可采用MACRO来代替。
第二种优化技术，将循环体内的代码展开，由一个循环步只做一次乘法，改为一个循环步做32次乘法。
实际速度：计算10000000！，程序4需要0.11667秒，程序5只需要0.02145秒。测试环境为迅驰1.7G,256M RAM。关于程序优化方面的内容，将在后续**专门讲述。下面是被优化后的部分代码，其余代码不变。
程序5的部分代码：
inline short GetExpBase2(double a) // 获得 a 的阶码
{
// 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu
    WORD *pWord=(WORD *)(&a)+3;
    WORD rank = ( (*pWord & 0x7fff) >>4 );
    return (short)(rank-0x3ff);
}
inline double GetMantissa(double a) // 获得 a 的 尾数
{
// 按照IEEE 754浮点数格式，取得尾数，仅仅适用于Intel 系列 cpu
    WORD *pWord=(WORD *)(&a)+3;
   *pWord &= 0x800f; //清除阶码
    *pWord |= 0x3ff0; //重置阶码

    return a;
}
void calcFac(struct bigNum *p,int n)
{
   int i,t;
   double f,max_mantissa;
   p->n1=1;p->n2=0;t=n-32;
   for (i=1;i<=t;i+=32)
   {
        p->n2 += GetExpBase2(p->n1); p->n1 = GetMantissa(p->n1);
        f=(double)i;
        p->n1 *=(double)(f+0.0);      p->n1 *=(double)(f+1.0);
        p->n1 *=(double)(f+2.0);      p->n1 *=(double)(f+3.0);
        p->n1 *=(double)(f+4.0);      p->n1 *=(double)(f+5.0);
        p->n1 *=(double)(f+6.0);      p->n1 *=(double)(f+7.0);
        p->n1 *=(double)(f+8.0);      p->n1 *=(double)(f+9.0);
        p->n1 *=(double)(f+10.0);          p->n1 *=(double)(f+11.0);
        p->n1 *=(double)(f+12.0);          p->n1 *=(double)(f+13.0);
        p->n1 *=(double)(f+14.0);          p->n1 *=(double)(f+15.0);
        p->n1 *=(double)(f+16.0);          p->n1 *=(double)(f+17.0);
        p->n1 *=(double)(f+18.0);          p->n1 *=(double)(f+19.0);
        p->n1 *=(double)(f+20.0);          p->n1 *=(double)(f+21.0);
        p->n1 *=(double)(f+22.0);          p->n1 *=(double)(f+23.0);
        p->n1 *=(double)(f+24.0);          p->n1 *=(double)(f+25.0);
        p->n1 *=(double)(f+26.0);          p->n1 *=(double)(f+27.0);
        p->n1 *=(double)(f+28.0);          p->n1 *=(double)(f+29.0);
        p->n1 *=(double)(f+30.0);          p->n1 *=(double)(f+31.0);
   }
  
   for (;i<=n;i++)
   {
        p->n2 += GetExpBase2(p->n1);
        p->n1 = GetMantissa(p->n1);
        p->n1 *=(double)(i);
   }
}
注1：10＾0，表示10的0次方


使用特权

评论回复 赏
板凳
gaoyang9992006|  楼主 | 2016-6-14 13:51 | 只看该作者
近似计算之二
在《阶乘之计算从入门到精通－近似计算之一》中，我们采用两个数来表示中间结果，使得计算的范围扩大到1千万，并可0.02秒内完成10000000！的计算。在保证接近16位有效数字的前提下，有没有更快的算法呢。很幸运，有一个叫做Stirling的公式，它可以快速计算出一个较大的数的阶乘，而且数越大，精度越高。有http://mathworld.wolfram.com查找Stirling’s Series可找到2个利用斯特林公式求阶乘的公式，一个是普通形式，一个是对数形式。前一个公式包含一个无穷级数和的形式，包含级数前5项的公式为n!=sqrt(2*PI*n)*(n/e)n*(1+ 1/12/n+ 1/288/n2 –139/51840/n3-571/2488320/n4+…),这里的PI为圆周率，e为自然对数的底。后一个公式也为一个无穷级数和的形式，一般称这个级数为斯特林(Stirling)级数,包括级数前3项的n!的对数公式为：ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n3 + 1/1260/n5)-…,下面我们分别用这两个公式计算n!.
完整的代码如下：
#include "stdafx.h"
#include "math.h"
#define PI       3.1415926535897932384626433832795
#define E 2.7182818284590452353602874713527
struct bigNum
{
       double n; //科学计数法表示的尾数部分
       int    e; //科学计数法表示的指数部分,若一个bigNum为x,这x=n * 10^e
};
const double a1[]=
{     1.00, 1.0/12.0, 1.0/288, -139/51840,-571.0/2488320.0 };
const double a2[]=
{     1.0/12.0, -1.0/360, 1.0/1260.0 };
void printfResult(struct bigNum *p,char buff[])
{
       while (p->n >=10.00 )
       {p->n/=10.00; p->e++;}
       sprintf(buff,"%.14fe%d",p->n,p->e);
}
// n!=sqrt(2*PI*n)*(n/e)^n*(1+ 1/12/n+ 1/288/n^2 -139/51840/n^3-571/2488320/n^4+…)
void calcFac1(struct bigNum *p,int n)
{
       double logx,
              s,            //级数的和
              item;  //级数的每一项　　　
       int i;
       // x^n= pow(10,n * log10(x));
       logx= n* log10((double)n/E);
       p->e= (int)(logx);   p->n= pow(10.0, logx-p->e);
       p->n *= sqrt( 2* PI* (double)n);
      
      //求(1+ 1/12/n+ 1/288/n^2 -139/51840/n^3-571/2488320/n^4+…)
       for (item=1.0,s=0.0,i=0;i<sizeof(a1)/sizeof(double);i++)
       {
              s+= item * a1;
              item /= (double)n; //item= 1/(n^i)
       }
       p->n *=s;
}
//ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n^3 + 1/1260/n^5 -...)
void calcFac2(struct bigNum *p,int n)
{
       double logR;
       double s, //级数的和
       item;       //级数的每一项
       int i;
      
       logR=0.5*log(2.0*PI)+((double)n+0.5)*log(n)-(double)n;
      
       //　s= (1/12/n -1/360/n^3 + 1/1260/n^5)
       for (item=1/(double)n,s=0.0,i=0;i<sizeof(a2)/sizeof(double);i++)
       {
              s+= item * a2;
              item /= (double)(n)* (double)n; //item= 1/(n^(2i+1))
       }
       logR+=s;
      
       //根据换底公式，log10(n)=ln(n)/ln(10)
       p->e = (int)( logR / log(10));
       p->n = pow(10.00, logR/log(10) - p->e);
}
int main(int argc, char* argv[])
{
       struct bigNum r;
       char buff[32];
       int n;
       printf("n=?");
       scanf("%d",&n);
      
       calcFac1(&r,n);            //用第一种方法，计算n的阶乘
       printfResult(&r,buff);    //将结果转化一个字符串
       printf("%d!=%s by method 1/n",n,buff);
       calcFac2(&r,n);            //用第二种方法，计算n的阶乘
       printfResult(&r,buff);    //将结果转化一个字符串
       printf("%d!=%s by method 2/n",n,buff);
       return 0;
}
程序运行时间的测量
算法的好坏有好多评价指标，其中一个重要的指标是时间复杂度。如果两个程序完成一个同样的任务，即功能相同，处理的数据相同，那么运行时间较短者为优。操作系统和库函数一般都提供了对时间测量的函数，这么函数一般都会返回一个代表当前时间的数值，通过在运行某个程序或某段代码之前调用一次时间函数来得到一个数值，在程序或者代码段运行之后再调用一次时间函数来得到另一个数值，将后者减去前者即为程序的运行时间。
　在windwos平台(指windwow95及以后的版本，下同),常用的用于测量时间的函数或方法有三种：1.API函数GetTickCount或C函数clock,　2.API函数QueryPerformanceCounter，　3:汇编指令RDSTC
1．API函数GetTickCount：
函数原形：DWORD GetTickCount(VOID)；
该函数取回从电脑开机至现在的毫秒数，即每个时间单位为1毫秒。他的分辨率比较低，常用在测量用时较长程序，如果你的程序用时为100毫秒以上，可以使用这个函数.另一个和GetTickCount类似的函数是clock，该函数的回的时间的单位的是CLOCKS_PER_SEC，在windows95/2000操作系统，该值是1000,也就是说，在windows平台，这两个函数的功能几乎完全相同。
2.API函数QueryPerformanceCounter：
函数原形：BOOL QueryPerformanceCounter( LARGE_INTEGER *lpPerformanceCount);该函数取回当前的高分辨值performance计数器，用一个64bit数来表示。如果你的硬件不支持高分辨performance计数器，系统可能返加零。不像GetTickCount，每个计数器单位表示一个固定的时间1毫秒，为了知道程序确切的执行时间，你需要调用函数QueryPerformanceFrequency来得到每秒performance计数器的次数，即频率。
3.汇编指令RDTSC:
RDTSC 指令读取CPU内部的“时间戳(TimeStamp)"，它是一个64位无符号数计数器，这条指令执行完毕后，存储在EDX:EAX寄存中。该指令从intel奔腾CPU开始引入，一些老式的CPU不支持该指令，奔腾后期的CPU包括AMD的CPU均支持这条指令。和QueryPerformanceCounter类似，要想知道程序的确实的执行时间，必须知道CPU的频率，即平常所说的CPU的主频。不幸的是没有现成的函数可以得到CPU的频率。一种办法可行的办法延时一段指定时间，时间的测量可以用QueryPerformanceCounter来做，在这段时间的开始和结束调用RDTSC,得到其时钟数的差值，然后除以这段时间的的秒数就可以了。


使用特权

评论回复 赏
地板
gaoyang9992006|  楼主 | 2016-6-14 13:52 | 只看该作者

下面的代码给出使用3个函数封装和测试代码，用RDTSC指令来计时的代码参考了一个Ticktest的源代码，作者不详。

复制


getTime1，使用GetTickCount返回一个表示当前时间的值，单位秒。

getTime2，和getTime1类似，精度更高。

getTime3，返回一个64bit的一个计数器，欲转换为秒，需除以CPU频率。示例代码见函数test3.

#include "stdafx.h"

#include "windows.h"

#include "tchar.h"

double getTime1()

{

       DWORD t=GetTickCount();

       return (double)t/1000.00;

}

double getTime2() //使用高精度计时器

{     

       static LARGE_INTEGER s_freq;

       LARGE_INTEGER performanceCount;

       double t;

       if (s_freq.QuadPart==0)

       {

              if ( !QueryPerformanceFrequency( &s_freq))

                     return 0;

       }

      

       QueryPerformanceCounter( &performanceCount );

       t=(double)performanceCount.QuadPart / (double)s_freq.QuadPart;

       return t;

}

void test1()

{

       double t1,t2;

       t1=getTime1();

       Sleep(1000);

       t2=getTime1()-t1;

       printf("It take %.8f second/n",t2);

}

void test2()

{

       double t1,t2;

       t1=getTime2();

       Sleep(1000);

       t2=getTime2()-t1;

       printf("It take %.8f second/n",t2);

}

inline BOOL isNTOS() //检测是否运行在NT操作系统

{

       typedef BOOL (WINAPI *lpfnGetVersionEx) (LPOSVERSIONINFO);

       static int bIsNT=-1;

       if (bIsNT!=1)

              return (BOOL)bIsNT;

       // Get Kernel handle

       HMODULE hKernel32 = GetModuleHandle(_T("KERNEL32.DLL"));

       if (hKernel32 == NULL)

              return FALSE;

 #ifdef _UNICODE

    lpfnGetVersionEx lpGetVersionEx = (lpfnGetVersionEx) GetProcAddress(hKernel32, _T("GetVersionExW"));

 #else

    lpfnGetVersionEx lpGetVersionEx = (lpfnGetVersionEx) GetProcAddress(hKernel32, _T("GetVersionExA"));

 #endif

 

       if (lpGetVersionEx)

       {

              OSVERSIONINFO osvi;

              memset(&osvi, 0, sizeof(OSVERSIONINFO));

              osvi.dwOSVersionInfoSize = sizeof(OSVERSIONINFO);

              if (!lpGetVersionEx(&osvi))

                     bIsNT=FALSE;

              else

                     bIsNT=(osvi.dwPlatformId == VER_PLATFORM_WIN32_NT);

       }

       else

       {

              //Since GetVersionEx is not available we known that

              //we are running on NT 3.1 as all modern versions of NT and

              //any version of Windows 95/98 support GetVersionEx

              bIsNT=TRUE;

       }

       return bIsNT;

}

inline static BOOL checkRDSTC() //检测CPU是否支持RDSTC指令

{

       static int bHasRDSTC= -1;

       SYSTEM_INFO sys_info;

      

       if ( bHasRDSTC !=-1 )

              return (BOOL)bHasRDSTC;

       GetSystemInfo(&sys_info);

       if (sys_info.dwProcessorType==PROCESSOR_INTEL_PENTIUM)

       {

              try

              {

                     _asm

                     {

                     _emit 0x0f                           ; rdtsc    

                     _emit 0x31

                     }

              }

              catch (...)       // Check to see if the opcode is defined.

              {

                     bHasRDSTC=FALSE; return FALSE;

              }

              // Check to see if the instruction ticks accesses something that changes.

              volatile ULARGE_INTEGER ts1,ts2;

              _asm

              {

                     xor eax,eax

                     _emit 0x0f                                  ; cpuid

                     _emit 0xa2

                     _emit 0x0f                                  ; rdtsc

                     _emit 0x31

                     mov ts1.HighPart,edx

                     mov ts1.LowPart,eax

                     xor eax,eax

                     _emit 0x0f                                  ; cpuid

                     _emit 0xa2

                     _emit 0x0f                                  ; rdtsc

                     _emit 0x31

                     mov ts2.HighPart,edx

                     mov ts2.LowPart,eax

              }

              // If we return true then there's a very good chance it's a real RDTSC instruction!

              if (ts2.HighPart>ts1.HighPart)

                  bHasRDSTC=TRUE;

              else if (ts2.HighPart==ts1.HighPart && ts2.LowPart>ts1.LowPart)

                     bHasRDSTC=TRUE;

              else

              {

                     printf("RDTSC instruction NOT present./n");

                     bHasRDSTC=FALSE;

              }

       }

       else

              bHasRDSTC=FALSE;

       return bHasRDSTC;

}

//***********************************************

void getTime3( LARGE_INTEGER *pTime) //返加当前CPU的内部计数器

{

       if (checkRDSTC())

       {

              volatile ULARGE_INTEGER ts;

       //on NT don't bother disabling interrupts as doing

       //so will generate a priviledge instruction exception

              if (!isNTOS())

                     _asm cli

           //----------------     

        _asm

              {

                     xor eax,eax

            //-------------save rigister

                     push ebx

                     push ecx

                    

                     _emit 0x0f                                  ; cpuid - serialise the processor

                     _emit 0xa2

                    

                     //------------

                     _emit 0x0f                                  ; rdtsc

                     _emit 0x31

                    

                     mov ts.HighPart,edx

                     mov ts.LowPart,eax

                    

                     pop ecx

                     pop ebx

              }

        //-----------------

              if (!isNTOS())

                     _asm      sti

              //---------       

        pTime->QuadPart=ts.QuadPart;

       }

       else

           pTime->QuadPart=0;

}

// maxDetermainTime:最大测定时间,单位毫秒，在首次调用该函数时，

// 将花费maxDetermineTime的时间来测定CPU频率，以后的调用将直接返加静态变量的值

double GetCPUFrequency(DWORD maxDetermineTime )

{

       static double CPU_freq;

       LARGE_INTEGER period,t1,t2;

       register LARGE_INTEGER goal,current;

      

       if (CPU_freq>1000)      //this value have been initilization

              return CPU_freq;

       if (!QueryPerformanceFrequency(&period) || !checkRDSTC())

       {

              CPU_freq=-1.00;

              return CPU_freq;

       }

       QueryPerformanceCounter(&goal);

       goal.QuadPart += period.QuadPart * maxDetermineTime/1000;

       getTime3( &t1);  //开始计时

       do    //延时maxDetermineTime毫秒

       {

              QueryPerformanceCounter(¤t);

       } while(current.QuadPart<goal.QuadPart);

       getTime3(&t2);      //结束计时

      

       CPU_freq=double((t2.QuadPart-t1.QuadPart)*1000/maxDetermineTime);

      

       char buff[100];

       sprintf(buff,"Estimated the processor clock frequency =%gHz/n",CPU_freq);

    ::MessageBox(NULL,buff,"",MB_OK);

       return CPU_freq;

}

void test3()

{

       LARGE_INTEGER t,t1,t2;

       double f1,f2;

      

       GetCPUFrequency(100); //花费0.1秒时间计算CPU频率

       f1=getTime2();

       getTime3(&t1);

       Sleep(1000);

       getTime3(&t2);

       f2=getTime2();

       t.QuadPart=t2.QuadPart-t1.QuadPart;

       printf("It take %.8f second by getTime3/n",(double)t.QuadPart/GetCPUFrequency(100));

       printf("It take %.8f second by getTime2/n",f2-f1);

}

int main(int argc, char* argv[])

{

       test1();

       test2();

       test3();

       return 0;

}


使用特权

评论回复 赏
5楼
gaoyang9992006|  楼主 | 2016-6-14 13:52 | 只看该作者
入门篇之一
在《大数阶乘之计算从入门到精通－大数的表示》中，我们学习了如何表示和存储一个大数。在这篇**中，我们将讨论如何对大数做乘法运算，并给出一个可以求出一个大整数阶乘的所有有效数字的程序。
大整数的存储和表示已经在上一篇**做了详细的介绍。其中最简单的表示法是:大数用一个字符型数组来表示，数组的每一个元素表示一位十进制数字，高位在前，低位在后。那么，用这种表示法，如何做乘法运算呢？其实这个问题并不困难，可以使用模拟手算的方法来实现。回忆一下我们在小学时，是如何做一位数乘以多位数的乘法运算的。例如：2234*8。
  [img][/img]
  
我们将被乘数表示为一个数组A[], a[1],a[2],a[3],a[4]分别为2,2,3,4,a[0]置为0。
Step1: 从被乘数的个位a[4]起，取出一个数字4.
Step2: 与乘数8相乘，其积是两位数32，其个位数2作为结果的个位,存入a[4], 十位数3存入进位c。
Step3: 取被乘数的上一位数字a[3]与乘数相乘，并加上上一步计算过程的进位C,得到27，将这个数的个位7作为结果的倒数第二位，存入a[3]，十位数2存入进位c。
Step4：重复Step3,取a（i依次为4,3,2,1）与乘数相乘并加上c,其个位仍存入a,　十位数字存入c，直到i等于1为止。
Step5:将最后一步的进位c作为积的最高位a[0]。
这一过程其实和小学学过的多位数乘以一位数的珠算乘法一模一样，学过珠算乘法的朋友还有印象吗？
在计算大数阶乘的过程中，乘数一般不是一位数字，那么如何计算呢？我们可以稍作变通，将上次的进位加上本次的积得到数P,　将P除以10的余数做为结果的本位，将P除以10的商作为进位。当被乘数的所有数字都和乘数相乘完毕后，将进位C放在积的最前面即可。下面给出C语言代码。
一个m位数乘以n位数，其结果为m+n-1,或者m+n位，所以需首先定义一个至少m+n个元素的数组，并置前n位为0。
计算一个m位的被乘数乘以一个n位的整数k，积仍存储于数组a
void mul(unsigned char a[],unsigned long k,int m,int n)
{
     int i;
    unsigned long p;
    unsigned long c=0;
    for ( i=m+n-1; i>=n;i--)
    {
           p= a * k +c;
           a=(unsigned char)( p % 10);
           c= p / 10;
    }
   
    while (c>0)
    {
           a=(unsigned char)( c % 10);
           i--;
           c /=10;
    }
}
int main(int argc, char* argv[])
{
    int i;
    unsigned char a[]={0,0,0,2,3,4,5};
    mul(a,678,4,3);
    i=0;
    while ( a==0)
           i++;
    for (;i<4+3;i++)
        printf("%c",a+’0’); //由于数a（0<=a <=9）对应的可打印字任符为’0’到’9’,所以显示为i＋’0’
               return 0;
}
从上面的例子可知，在做乘法之前，必须为数组保留足够的空间。具体到计算n!的阶乘时，必须准备一个能容纳的n!的所有位数的数组或者内存块。即数组采有静态分配或者动态分配。前者代码简洁，但只适应于n小于一个固定的值，后者灵活性强，只要有足够的内存，可计算任意n的阶乘，我们这里讨论后一种情况，如何分配一块大小合适的内存。
n!有多少位数呢？我们给出一个近似的上限值：n! <(n+1)/2的n次方，下面是推导过程。
Caes 1: n是奇数，则中间的那个数mid= (n+1)/2, 除了这个数外，我们可以将1到n之间的数分成n/2组，每组的两个数为 mid-i和mid+i (i=1到mid-1),如1,2,3,4,5,6,7 可以分为数4，和3对数，它们是(3,5),(2,6)和(1,7),容易知道，每对数的积都于小mid*mid，故n!小于(n+1)/2 的n的次方。
Case 2: n 是个偶数，则中间的两个数(n－1)/2和(n+1)/2, 我们将(n+1)/2记做mid,则其它的几对数是(mid-2,mid+1),(mid-3)(mid+2)等等，容易看出,n!小于mid 的n次方。
由以上两种情况可知，对于任意大于1的正整数n, n!<(n+1)/2的n次方。
如果想求出n!更准确的上限，可以使用司特林公式，参见该系列**《阶乘之计算从入门到精通－近似计算之二》。

使用特权

评论回复 赏
6楼
gaoyang9992006|  楼主 | 2016-6-14 13:53 | 只看该作者
到此，我们已经解决大数阶乘之计算的主要难题，到了该写出一个完整程序的时候了，下面给出一个完整的代码。
#include "stdio.h"
#include "stdlib.h"
#include "memory.h"
#include "math.h"
#include "malloc.h"
void calcFac(unsigned long n)
{
   unsigned long i,j,head,tail;
int blkLen=(int)(n*log10((n+1)/2)); //计算n!有数数字的个数
   blkLen+=4;  //保险起见，多加4位
’
  if (n<=1)
   {       printf("%d!=0/n",n);     return;}
   char *arr=(char *)malloc(blkLen);        
   if (arr==NULL)
   {       printf("alloc memory fail/n"); return ;}
   memset(arr,0,sizeof(char)*blkLen);
   head=tail=blkLen-1;
   arr[tail]=1;
   
   for (i=2;i<=n;i++)
   {
           unsigned long c=0;
           for (j=tail;j>=head;j--)
           {
                    unsigned long prod=arr[j] * i +c;
                    arr[j]=(char)( prod % 10);
                    c= prod / 10;
           }
           while (c>0)
           {
                    head--;   
                    arr[head]=(char)(c % 10);
                    c/=10;
           }
   }
   printf("%d!=",n);
   for (i=head;i<=tail;i++)
           printf("%c",arr+'0');
   printf("/n");
   free(arr);
}
void testCalcFac()
{
    int n;
    while (1)
    {
            printf("n=?");
            scanf("%ld",&n);
            if (n==0)        
                   break;
            calcFac(n);
    }
}
int main(int argc, char* argv[])
{
    testCalcFac();
    return 0;
}

使用特权

评论回复 赏
7楼
gaoyang9992006|  楼主 | 2016-6-14 13:59 | 只看该作者
太多了，不贴了，有两万字限制。。。
大数阶乘.pdf (882.21 KB)




使用特权

评论回复 赏
返回列表
发新帖 我要提问
高级模式
B Color Image Link Quote Code 收费 Smilies
您需要登录后才可以回帖 登录 | 注册

本版积分规则

gaoyang9992006 Image Image
认证：西安公路研究院南京院
简介：主要工作从事监控网络与通信网络设计，以及从事基于嵌入式的通信与控制设备研发。擅长单片机嵌入式系统物联网设备开发，音频功放电路开发。

1965

主题

15945

帖子

208

粉丝