回 帖 发 新 帖 刷新版面

主题:[原创]关于求任意数阶乘的算法

//factorial.h
/********************************************************************
*  文件名:   factorial.h
*  文件描述: 求阶乘
*  创建人:   陈泽丹 ,2006年3月23日       (QQ:82314038) 
*  版本号:   1.0
*  修改记录:
********************************************************************/
/*============================================================
    自定义类型(ElemType) :                        double
    类名:            Factorial
    -    x:            ElemType
    -    next:        Factorial* 
    //执行进位的操作    
    -    carry(Factorial* &c, ElemType xx):        ElemType 
    //显示结果的重载函数
    -    display(Factorial* t):                    void
    //格式化显示结果
    -    format(ElemType xx,int count):            void
    //构造函数防止next越界
    +    Factorial():x(0),next(NULL){}
    //设置类的初始值
    +    set(ElemType xx):                        void
    //执行乘法运算
    +    M(ElemType xx):                            void
    //显示结果
    +    display():                                void
    //显示结果的位数
    +    digit():                                double 
作  者:陈泽丹 2006/3/23
============================================================*/
/* --------------------------------------------------------------------
补充:
1,    以下CARRY与DIGIT的更新需同步,
2,    以下CARRY与DIGIT值越大,效率越高风险越大
3,    以下CARRY与DIGIT值越低,效率越差风险越小
4,    为合理利用空间,本程序采用动态分配的方案,故若阶乘数太大,请注意你的
    内存问题。
(解释:由于进位的堆积,如某数据段一万乘一万,进位后又刚好与下一个本身
就很大数据段的数据堆积一起,于是“突围了”。。。,所以越后面越有数据
段范围崩溃危险,改良方法是每一个数据段“递归”或改“类乘类(而非本例的'类乘
实数')”,不过这就不是此文三言两语说得清的了。另外,即使数据段定位再小也有风
险,另忘了数值小点的数乘大数还是可成大数的。)
----------------------------------------------------------------------*/

#ifndef HEADER_FACTORIAL
#define HEADER_FACTORIAL

#include <iostream.h>
#define        ElemType    double
#define        CARRY        1000        //进位的条件
#define        DIGIT        3            //数据段的位数

class Factorial
{
private:
    ElemType x;
    Factorial* next;        
    ElemType carry(Factorial* &c, ElemType xx);
    void display(Factorial* t);    
    void format(ElemType xx,int count);
public:
    Factorial():x(0),next(NULL){}
    void set(ElemType xx);
    void M(ElemType xx);
    void display();
    double digit();
};

#endif
/*======================================================================*/
//factorial.cpp
//factor类实现部份

#include "factorial.h"

void Factorial::set(ElemType xx)
{
    x=xx;
    Factorial *t=this;
    while (t != NULL)
    {
        if ( t->x > CARRY) 
        { 
            t->x=carry(t->next,t->x); 
        }
        t=t->next;
    }
}

void Factorial::M(ElemType xx)
{
    Factorial *t=this;
    while (t!=NULL)
    {
        t->x*=xx;
        t=t->next;
    }
    t=this;
    while (t != NULL)
    {
        if ( t->x > CARRY) 
        { 
            t->x=carry(t->next,t->x); 
        }
        t=t->next;
    }
}

ElemType Factorial::carry(Factorial* &c, ElemType xx)
{
    if ( c==NULL) 
    { 
        c= new Factorial;
        c->next=NULL;
    }
    int temp=xx/CARRY;
    c->x=c->x+temp;
    return xx-temp*CARRY;
}


void Factorial::display()
{
    display(this);
}
void Factorial::display(Factorial *t)
{
    if (t!=NULL)
    {
        display(t->next);
        if (t->next!=NULL)    { format(t->x,DIGIT); }
        else cout<<t->x;
        if (t != this) cout<<",";
    }
}
void Factorial::format(ElemType xx, int count)
{
    if (count>0)
    {
        int temp=xx;
        format(xx/10,count-1);
        cout<<temp%10;
    }
}
double Factorial::digit()
{
    Factorial *t=this;
    double i=0;
    while(t->next!=NULL)
    {
        i+=1;
        t=t->next;
    }
    int d=0;
    int temp=t->x;
    while(temp)
    {
        d++;
        temp=temp/10;
    }
    return i*DIGIT+d;
}

/*=================================================================*/
//Main.cpp
//主函数实现部份

#include "factorial.h"
#include "iostream.h"

void main()
{

    Factorial c;
    ElemType start, end;
    cout<<"本程序实现连乘效果(阶乘是连乘的特例)"<<endl;
    cout<<"请输入起始数:";
    cin>>start;
    c.set(start);
    cout<<"您输入的的起始数为";
    c.display();
    cout<<endl;
    cout<<"请输入终止数: ";
    cin>>end;
    for (int i=start+1;i<=end; i++) { c.M(i); }
    c.display();
    cout<<endl;
    cout<<"共有"<<c.digit()<<"位"<<endl;

}

回复列表 (共17个回复)

11 楼

哈哈,在这个论坛上竟然有人认识我,不胜荣幸。

计算阶乘的各个版本我的电脑上是有的。但我准备对3.0以上的各个程序重写,等写好了再发布,所以在上一个回复中没有提到。新版不但对硬乘法作彻底的优化(估计速度可提高20%),对大数乘法引入FFT算法,以取代现在4.x系列的FNT(快速数论变换)算法。
   除此之外,我准备进一步充实这个程序级,增加迷你版(用汇编语言编写,分为dos版和win32版,执行文件非常小),极速版(使用斯特林公式或者硬乘法,精度稍低(比如8-128位十进制数),但速度极快,估计硬乘法的版本也要比windows自带的计算器快上10000倍,采用斯特林公式就更是快了许多).
  
   再进一步,我可能会发布一篇文章,名字暂定为:阶乘之计算从入门到精通。内容将设计到
    1)各种乘法的(硬乘法,分治法,FFT)算法。
    2)目前所看到的各种阶乘程序(包括我的各个版本的程序)的算法和性能分析。
    3)程序设计中的算法(代码)优化方法和技巧,这一部分会涉及到汇编语言.
    文中将附带大量的代码和程序,但最终经过优化过的程序可能不提供,这项工作将是一项浩繁的工作,文章的最终发表可能需要等上较长的时间。

12 楼

哈哈,在这个论坛上竟然有人认识我,不胜荣幸。

计算阶乘的各个版本我的电脑上是有的。但我准备对3.0以上的各个程序重写,等写好了再发布,所以在上一个回复中没有提到。新版不但对硬乘法作彻底的优化(估计速度可提高20%),对大数乘法引入FFT算法,以取代现在4.x系列的FNT(快速数论变换)算法。
   除此之外,我准备进一步充实这个程序级,增加迷你版(用汇编语言编写,分为dos版和win32版,执行文件非常小),极速版(使用斯特林公式或者硬乘法,精度稍低(比如8-128位十进制数),但速度极快,估计硬乘法的版本也要比windows自带的计算器快上10000倍,采用斯特林公式就更是快了许多).
  
   再进一步,我可能会发布一篇文章,名字暂定为:阶乘之计算从入门到精通。内容将设计到
    1)各种乘法的(硬乘法,分治法,FFT)算法。
    2)目前所看到的各种阶乘程序(包括我的各个版本的程序)的算法和性能分析。
    3)程序设计中的算法(代码)优化方法和技巧,这一部分会涉及到汇编语言.
    文中将附带大量的代码和程序,但最终经过优化过的程序可能不提供,这项工作将是一项浩繁的工作,文章的最终发表可能需要等上较长的时间。

13 楼

目前我看到的程序中,速度最好的程序有3个,1。我的程序。2。郭先强的程序。3。北大刘楚雄写的一个程序。 程序3速度最快,其大数乘法运算采用FFT算法,这是最核心的,技术含量最高的部分,并非作者独立编写,而是直接使用了 ooura 的FFT 源代码,可从在http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html下载.

14 楼

waiting...

15 楼

/*
  写一个计算阶乘的程序吧,这个程序可以计算到40万的阶乘。这个程序可能是最简单的计算阶乘的程序,而且速度也不慢。如果计算小于1万的阶乘,下面的代码可以不用循环
   while (carry>0)
  {
     pBeg--;
     *pBeg=(WORD)(carry % RAD);
     carry /= RAD;
  }
*/

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

#define PI  3.1415926535897932384626433832795
#define RAD 10000
typedef unsigned long DWORD;
typedef unsigned short WORD;

//用stirling 公式估算结果长度,稍微算得大一点
DWORD calcResultLen(DWORD n,DWORD rad)
{
    double r=0.5*log(2*PI) + (n+0.5)*log(n)-n;
    return (DWORD)(r/log(rad)+2);
}

int main(int argc, char* argv[])
{
    DWORD i,n,carry,prod,len;
    WORD *buff,*pBeg,*pEnd,*p;
    //----------------------------
    printf("please input n:");
    scanf("%d",&n);
    
    if (n==0)
    { printf("%d!=1",n); return 0;}

    //---计算并分配所需的存储空间
    len=calcResultLen(n,RAD);
    buff=(WORD*)malloc( sizeof(WORD)*len);
    if (buff==NULL)
    return 0;
 
    //以下代码计算n! 
     pBeg=pEnd=buff+len-1;
     for (*pEnd=1,i=2;i<=n;i++)
     {
         for (carry=0,p=pEnd;p>=pBeg;p--)
    {
        prod=(DWORD)(*p) * i +carry;
        *p=(WORD)(prod % RAD);
        carry=prod / RAD;
    }
    while (carry>0)
    {
       pBeg--;
       *pBeg=(WORD)(carry % RAD);
       carry /= RAD;
    }
    }
    //显示计算结果
    printf("%d!=%d",n,*pBeg);
    for (p=pBeg+1;p<=pEnd;p++)
       printf("%04d",*p);
    free(buff);//释放分配的内存
    return 0;
}

16 楼

liangbch原来这么早就研究过了啊~我的程序只能达到900ms左右,估计和最快的差一两个数量级吧,向你学习~

你给的那个最快的fft的网页打不开了~~现在要继续提高速度,要从哪方面着手,我的方法是把阶乘拆成素因子的方法~你提到的分治的乘法,会比FFT快吗,复杂度是O(n^1.x)的吗?

17 楼

我的博客(blog.csdn.net/liangbch)开通了,内有大数阶乘之计算系列文章

我来回复

您尚未登录,请登录后再回复。点此登录或注册