CSDN - 专家门诊 -
主  题:
擂台:超大整数高精度快速算法
作  者: gxqcn (GxQ)
信 誉 值: 100
所属论坛: 专题开发 数据结构与算法
问题点数: 100
回复次数: 156
发表时间: 2004-5-10 14:08:44

我用 C++ 开发了一套《超大整数高精度快速算法》,前后花了十余年。

实现了如下功能:

⊙ 高精度快速加法
⊙ 高精度快速减法
⊙ 高精度快速乘法
⊙ 高精度快速除法
⊙ 高精度快速同余
⊙ 高精度快速乘方
⊙ 高精度快速开方
⊙ 高精度快速乘方取模
⊙ 高精度快速计算 Fibonacci 数列
⊙ 高精度快速阶乘、双阶乘
⊙ 高精度快速求排列
⊙ 高精度快速求组合
⊙ 高精度快速求最大公约数(支持群组运算)
⊙ 高精度快速求最小公倍数(支持群组运算)
⊙ 高精度快速“等幂和”(支持群组运算)
⊙ 高精度快速任意进位制转换
⊙ 高精度计时器(有暂停、累计、复位等功能)
⊙ 强大而灵活的输出

该套算法效率非常高,举例如下:

在 P4 1.7GHz 256MB 的机器上:
生成所有不大于 1亿的素数仅需 0.828s;
精确计算 100,000! 需要 2.078s;
精确计算 200,000! 需要 6.935s。

(我的算法纯粹是数学上的优化,未内嵌一行汇编代码)

最新版下载地址:http://maths.diy.myrice.com/software.htm#02

因为 CSDN 是一个高手云集的地方,所以特借这块宝地广泛征集类似算法或软件,以期达到相互促进的目的。

注:
1、声明, 我的程序快只能说明我对此算法下了功夫,并不表明水平有多高,事实上,我的水平很有限(甚至还不会汇编);
2、并不要求所有功能均须强于我现有的软件,有一两项也请及时告诉我,谢谢!
3、算法或软件不限于自身原创,欢迎提供第三方链接信息,谢谢!

回复人: wlpwind(robin) ( ) 信誉:105 2004-5-10 16:20:39 得分:0
 

佩服。

Top
回复人: mysword(一怒拔剑) ( ) 信誉:100 2004-5-10 17:03:50 得分:0
 

没有源码


Top
回复人: aheadyes() ( ) 信誉:100 2004-5-10 17:17:39 得分:0
 

佩服++

Top
回复人: NowCan(((((( ★ )))))) ( ) 信誉:100 2004-5-10 18:24:55 得分:0
 

佩服,哪天和GMP比试一下。

Top
回复人: ZhangYv(has been killed) ( ) 信誉:129 2004-5-10 19:39:38 得分:0
 

你的 素数 和 N!是怎么计算出来的啊,没有使用一些数学方法吧?
求素数和求阶乘以前有讨论过一些数学方法,所以用这个标准来衡量是不太可靠的。在“以前的精华区”里有,可惜我刚才一点精华区发现帖子都没有了,心寒...不知道是不是BUG!

Top
回复人: ZhangYv(has been killed) ( ) 信誉:129 2004-5-10 19:47:29 得分:0
 

如果你的作品不是作为共享,不妨贴出代码让大家见识一下吧:)

精华区的BUG...我发消息去问了,希望那些帖子不要没掉.

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-5-10 19:58:09 得分:0
 

哈哈!老弟终于来了,我的 阶乘计算器 最后一个版本还没有完成,如果我的新版比你的算得更快,你还出新版吗?

Top
回复人: arya(行者) ( ) 信誉:100 2004-5-10 23:19:09 得分:0
 

有没有跟apfloat比较过?
http://www.apfloat.org/apfloat/


Top
回复人: programfanny() ( ) 信誉:100 2004-5-11 2:06:22 得分:0
 

佩服得了得! 你的高精度精到哪里?能无限位精度吗?我说的是具有自适应性能.

Top
回复人: shines(郭子) ( ) 信誉:100 2004-5-11 6:31:27 得分:0
 

你的质数的算法也太慢了:)

生成所有不大于 1亿的素数仅需 0.828s;

看我最近写的,运行结果:(还开着CSDN论坛的2个IE窗口,一个正在以100多KB的速度的BT下载,赛扬II 900, 384M内存)

PrimeLab Pro - Windows(R) CPU Identifier
-----------------------------------------

    Microsoft Windows 2000 Advanced Server Service Pack 4(Build 2195)
    -----------------------------------------------------------------
    TotalPhys: (383MB)392692KB, AvailPhys: 27336KB, FreePercent: 7%

    Pentium III (0.18 um) With 256 KB On-Die L2 Cache  - [CPU #0]
    -------------------------------------------------------------
        | Clock Frequency: 896MHz
        |
        + On-Chip Hardware Features
            | L1 Cache: 32KB
            | L2 Cache: 128KB

-----------------------------------------

Input A Number [2-4294967295]: 100000000

Number Range [2-99999999]:
Total have 5761455 Prime(s).
Used Time: 0.56384904 second(s).

Press any key to Exit ...


Top
回复人: shines(郭子) ( ) 信誉:100 2004-5-11 6:54:52 得分:0
 

上面的测试由于我的CPU太忙,安静的时候应该更快,不过我只是产生标记位而已,并没有生成质数表,问一下:生产质数是指把质数保存在数组吗?还是只保存标记位(按标记位查是否是质数)?如果要输出到数组,我还得扫描一下,可能需要一点时间

如果只是统计质数的个数可以更快一点

仔细一看是整数的高精度(其实写成整数的意义不大),
最好写成小数的,然后用来计算圆周率PI,虽然有不少的代码,但是自己写一个也非常有成就感

Top
回复人: languagec(C++)(各有所求) ( ) 信誉:100 2004-5-11 7:36:51 得分:0
 

是吗?能这么快!? 把代码贴出来比比看。

Top
回复人: programfanny() ( ) 信誉:100 2004-5-11 9:21:14 得分:0
 

同意楼上的观点  我们关心的是源代码! 算法思想也行啊.说说怎么实现的吧.

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-5-11 19:12:46 得分:0
 


非常高兴,小弟我第一次来宝地,就得到各位的捧场。

下面集中回复一些问题:

to NowCan:
    能与其它专业软件一比高下,这是我的目标,可惜能力不足,真希望有条件的朋友,
    将本人这套算法与一些专业数学软件进行对比,给出客观公正的评测数据,本人将不胜感激。
    如你提到的 GMP, arya(行者) 提到的 apfloat,以及著名的 Mathematica,。。。

to programfanny:
    该套算法具有“无限位精度”,即你所说的“自适应性能”。(我采用 STL 去动态管理内存分配)

to shines:
    看来我们还是本家呢。以前我的外号是“郭大侠”,并非指工夫有多高,而是说我太憨厚老实了:) 
    说直点,可能就是有点呆傻吧。
    
    关于“素数”的问题,我切入的时机非常晚,也就在“五一”前后,主要目的是为了开发更快的“阶乘”
    算法(整个“五一”基本上都在闭关修炼,谁让咱太笨了)。主要算法思想先后参考了 liangbch(宝宝)  
    和 mars22(三月瓜) 两位大侠的,然后自己又做了些优化工作,虽然速度已很快,但不敢奢望可以天下无
    敌,毕竟在此花的工夫并不多。从你给的数据来看,效率确实比我高。我的程序兼顾了用户的内存资源,
    否则速度可以再提高,但似乎意义已不大了。
    
to liangbch:
    宝宝大哥,你是我特别应该感谢的,你的阶乘计算器真的很快,并激发了我开发更快算法的激情,我们私
    底下也通过 e-mail 进行了多次交流,相互促进不少。衷心祝愿你的新版早点出来,到时可记得一定发给
    我一份。由于我将在秋季迎接我的 baby 诞生,所以如果现有版本未发现大的bug,将不会再更新它了。
    
    关于大数乘法,如果位数在一定范围,我采用的是普通的“硬乘法”;比较大的则采用分治法。由于曾开发过“等幂和计算器”PowCalc,积累了一些“群组运算”的算法技巧,所以使“阶乘”速度非常快。

另外,需说明的是:
1、我的程序用了多线程,虽然控制方便了,速度却略有下降(因为要经常判断是否需强制退出);
2、用 C++ 类,效率略低;
3、未用汇编优化;
4、我的 myrice 主页今天无法浏览(系统提示:“LYCOS主页服务免费空间系统维护中,暂停页面访问。系统会在5月13日后恢复访问。”),为方便大家下载软件,特将镜像站点列举如下:

⊙ http://maths.diy.myrice.com/
⊙ http://www.cccen.com/web/math/
⊙ http://www.freewebs.com/maths/
⊙ http://rd.bbs.xilu.com/(数学研发论坛)

(注:我的主页和论坛已被多家知名搜索引擎收录,甚至双双被 Google 收录进“World > Chinese Simplified > 科学 > 数学”网页目录中 )


Top
回复人: BlueSky2008(懒惰是程序员的美德) ( ) 信誉:160 2004-5-12 9:38:49 得分:0
 

最好做成一个高精度计算库,这样实用价值可能更高点。

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-5-12 15:34:32 得分:0
 

它的全称即为《超大整数完全精度快速算法库》,已封装成 dll,提供了 169 个公共接口,调用非常方便。

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-5-14 7:48:09 得分:0
 

怎么昨天一天没人响应?郁闷中。。。

Top
回复人: jp1984(毛咏洁) ( ) 信誉:100 2004-5-14 12:21:20 得分:0
 

这个板块更新确实慢的可以。。C++那块 很快。。
  用你的那个 统计质数的 把我机器算死机了,我算的是4000000000 以内质数
  
   依然万分佩服你。。你是数学专业的吧?

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-5-14 12:56:17 得分:0
 


    首先申明,我设计的那个 PrimeNumber 不是为了统计质数的,而是为了把里面的素数统统纠出来。由于是未采用边搜边存的方式,所以如果选择了输出到窗口或文件,将会生成一个很大的字符串,需要耗很大的资源,所以感觉好象当机了。

    我很喜欢数学,可惜我不是数学专业的,上天作弄人。

Top
回复人: junmayang(小鱼儿) ( ) 信誉:100 2004-5-14 12:57:13 得分:0
 

楼主你太强了,我写的算法算10000!需要五分钟。你的算法这么快,难以想象啊,便服!
你的阶乘算法是怎么设计的啊?

Top
回复人: Gulfing(一切皆有可能) ( ) 信誉:100 2004-5-14 14:41:53 得分:0
 

十年啊!

Top
回复人: fairyboy(善财童子) ( ) 信誉:100 2004-5-14 17:32:13 得分:0
 

我做过一个大整数类,
使用技术:
1,SSE2优化
2,C++内嵌汇编
3,预计算,整一堆算好的数,二分查找
4,整点指令和SSE2指令交错使用
其他技术没用,不好意思说出速度,说出来可能没人相信


Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-7 10:41:15 得分:0
 

首先,祝贺 CSDN 技术社区恢复正常!只是可惜丢失了一些帖子。

在 2004-06-01,HugeCalc 正式升级为 Ver 2.1.0.1,下面是其改进点:

	⊙ 解掉了“开平方”、“开立方”算法中的bug

	⊙ 解掉了“超大整数除以小整数”算法中的bug

	⊙ 改进素数搜索算法、并对乘法模块进行了优化

	⊙ 核心模块首次采用内嵌汇编优化

	⊙ PrimeNumber 实现了分段搜索存储功能

	⊙ PrimeNumber、Factorial、HugeCalc 中的 32bit 入参上限均扩大为 2^32-1

//=========================================================================

该版改进最大的是 PrimeNumber,在 P4 1.7GHz CPU、256MB RAM 的机器上,计算 π( 2^32 -1 ) 仅需 0.440s

下面的数据是计算阶乘的耗时记录(运行 Factorial.exe):


 10,000!   0.038s     35,660位
 20,000!   0.127s     77,338位
 40,000!   0.425s    166,714位
 80,000!   1.429s    357,507位
100,000!   2.154s    456,574位
200,000!   7.042s    973,351位
400,000!  22.726s  2,067,110位
--------2004.05.16

 10,000!   0.031s     35,660位
 20,000!   0.109s     77,338位
 40,000!   0.406s    166,714位
 80,000!   1.375s    357,507位
100,000!   2.063s    456,574位
200,000!   6.813s    973,351位
400,000!  21.875s  2,067,110位
--------2004.05.17

 10,000!   0.031s     35,660位
 20,000!   0.125s     77,338位
 40,000!   0.406s    166,714位
 80,000!   1.375s    357,507位
100,000!   2.078s    456,574位
200,000!   6.766s    973,351位
400,000!  21.734s  2,067,110位
800,000!  69.422s  4,375,040位
--------2004.05.18

 10,000!   0.031s     35,660位
 20,000!   0.109s     77,338位
 40,000!   0.406s    166,714位
 80,000!   1.344s    357,507位
100,000!   2.015s    456,574位
200,000!   6.594s    973,351位
400,000!  21.344s  2,067,110位
800,000!  68.641s  4,375,040位
--------2004.05.24

 10,000!   0.031s     2.8462… x 10^35,659
 20,000!   0.109s     1.8192… x 10^77,337
 40,000!   0.406s     2.0916… x 10^166,713
 80,000!   1.343s     3.0977… x 10^357,506
100,000!   2.000s     2.8242… x 10^456,573
200,000!   6.500s     1.4202… x 10^973,350
400,000!  21.000s     2.5344… x 10^2,067,109
800,000!  67.718s     5.6846… x 10^4,375,039
--------2004.05.25

 10,000!   0.031s     2.8462… x 10^35,659
 20,000!   0.109s     1.8192… x 10^77,337
 40,000!   0.390s     2.0916… x 10^166,713
 80,000!   1.328s     3.0977… x 10^357,506
100,000!   1.969s     2.8242… x 10^456,573
200,000!   6.438s     1.4202… x 10^973,350
400,000!  20.828s     2.5344… x 10^2,067,109
800,000!  66.921s     5.6846… x 10^4,375,039
--------2004.05.28

最新版下载地址:

⊙ http://maths.diy.myrice.com/software.htm#02
⊙ http://www.cccen.com/web/math/software.htm#02
⊙ http://www.freewebs.com/maths/software.htm#02

(现在,该套软件已被华军、天空、电脑之家、天天等下载站点收录)

Top
回复人: lyr311(老刘:别总在CSDN上逛!!!) ( ) 信誉:100 2004-6-7 20:55:45 得分:0
 

真是高手啊,佩服佩服!赶上Mathematica软件了!

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-7 21:51:30 得分:0
 

顶一下

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-9 9:57:46 得分:0
 

如果我的回复不丢,你应该看到了我的汇编代码,还是原来的意见,公开源代码

大家相互借鉴

另外希望能写出FFT/FNT的算法

否则,不能谈高速,我记的超过一定bit的乘法采用分治有效,再大的要FFT,再大就只有NTT了


而且,有一个浮点除法对不能采用 x /y = x * (1 / y) 算法的小整数(bit少的整数)很有效果

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-9 10:06:53 得分:0
 

嘿嘿,你的算法计时有问题啊

我用服务器算1000000! 显示1'33.374

任务管理器显示2'15
我是双P4 2.0 XEON HT,没有任何计算延迟的

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-9 11:20:52 得分:0
 

to yaos(夜夜秋雨孤灯下·无心人) :

    最新的一版仅仅将一个需要“同时得到商和余数”的函数改用汇编(可省去一次乘法指令),其它均未用汇编优化。

    关于算法原理,在适当的时候我会在适当的地方公开,也许就在本论坛,或是整理在我的主页中。

    FFT/FNT 需要浮点运算,这不是我的强项,况且该算法库的开发初衷是为了自己研究数论和组合数学,虽然需要超大整数,但一般也不是非常非常大;NTT 则需要先有模运算全面支持才行,现这块我未做太多的算法优化工作。

    关于计时问题,需要说明的是它是不包括输出时间的。你说的显示时间偏小,那时因为将 1000000! 输出到窗口或文件需要额外的时间消耗。请运行 Factorial.exe,将 “Output” 栏所有选项均不勾选,你会发现基本一致。

另:请下载 http://www.freewebs.com/maths/download/HugeCalc.rar (152 KB),里面的 HugeCalc.exe 采用了新的 RichEdit 控件,功能更强(但需用户已安装VC)。

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-10 23:12:13 得分:0
 

我采用apfloat库和自己的阶乘算法,已经有些接近擂主的速度了,超过liangbch(宝宝)的是毫无疑问的。

apfloat竟然没有以10进制的移位,我将改进一下,相信速度会有所提高。另外,使用apfloat采用硬乘的方法计算N!依然还是很慢的,但采用有效的阶乘算法以后,速度提升率是非常高的,而且N越大就越有优势,目前版本200000!已经和擂主的相差无几了,10000!大约是擂主的2倍时间,

我采用自己写的大整数,不过是基于2进制的,计算10000!也是擂主的2倍时间,乘法部分使用的是分治法,不过可惜的是基于2进制的,输出还要经过漫长的转换(比计算N!本身还长)

继续改进。。。

Top
回复人: SweetJerry() ( ) 信誉:100 2004-6-10 23:16:32 得分:0
 

蛮强的,不知道计算PI怎么样,100万位的.以前有个程序,100W位大概28秒,PIII800EB/256M,注意不是SuperPI,因为SuperPI要幔些.

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-10 23:19:27 得分:0
 

本来帖子我一直都在看(几乎每天),不想说话而已,今天无聊了说几句吧

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-10 23:46:49 得分:0
 

SweetJerry(): 你说的可能就是apfloat的Demo,aptest

比较速度是要在同样的计算方法的前提下才公平的,你说的可能是aptest的第一种计算方法Chudnovskys's binsplit,在我的赛扬II900上,同样使用高斯迭代法(SuperPI使用的是高斯迭代法),aptest需130秒,SuperPI是90秒

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-11 8:25:40 得分:0
 

to shines(郭靖的郭,英雄的雄,光辉的辉) :
    首先,恭喜恭喜你,愿你的程序能早日面市。
    apfloat库在大整数的高精度乘法方面有些独到之处,这可能是使你计算较大数的阶乘效率提升较快。
    我在设计阶乘、双阶乘、排列、组合等功能时,运用了大量的数学知识优化,未专注于对程序代码的技巧及高精度乘法算法的优化改进。

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-11 11:51:12 得分:0
 

除了使用apfloat库以外,我也没有对程序代码做什么优化,无非是采用了质因子分解再相乘的方法,其中采用了迭代相乘的方法(减少了很多次乘法计算),我不知道你所谓的数学知识的优化是指什么,不过你可以用你的方法套用一下apfloat,看看效率是否有所提高(当然,如果你所说的数学知识的优化不适合做简单的大数计算的话,那就另当别论)

如果你需要一个apfloat的通用模板我可以提供给你,VC++ console/或SDK/MFC的(请注明),你就不必去劳神去整理apfloat了

个人觉得apfloat还有提升的空间,毕竟它是一个较通用的库,不过提升的性能不会很多

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-11 12:25:38 得分:0
 

我用到了数论知识(素数分解)、不等式、组合数学等,也同时尽可能地减少乘法次数(也许即你所说的“迭代相乘”)。这些优化是独立的,可以与任何大整数的高精度计算配合使用。

在Ver1.2.0.1,我未采用素数分解法,但用了一种技巧性很高的算法(至少我是这么认为;该算法完全独创),其速度与liangbch(宝宝)高级算法基本持平,并略微胜出。该算法具有很高的理论价值,也许我会整理一下给大家。

可否将apfloat的通用模板VC++ console/MFC 发一份给我,先谢谢了!
我的E-Mail在个人主页或软件上都可找到。

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-11 12:40:06 得分:0
 

我待会就发份VC++ console的给你,VC++ SDK的还没有整理出来,其实是一样的,质数基于的界面不一样,一个是基于GUI,一个基于console,有个注意的问题,由于使用了Microsoft Macro Assembler (MASM)汇编,你必须有MASM编译器ml.exe,可以去微软直接下载MASM,或打vs sp5 + vs processer pack 5.0即可。

其实我也想到一个比较愚蠢的想法,我可以用你的大整数库,套用我的算法来试试看

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-11 12:53:40 得分:0
 

好哇。到时我将 HugeCalc.h 和 HugeCalc.lib 文件寄给你一份,就可完全调用 HugeCalc.dll 了。

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-11 12:58:29 得分:0
 

我发到你的163.com的信箱了,你收到信给我回信就可以了,我也是163.com的

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-11 13:51:00 得分:0
 

已发出,注意查收。

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-11 16:43:39 得分:0
 

我下载了http://www.apfloat.org/apfloat/ 的源代码,并编译通过。并利于它和另外两个软件 pifast 4.1 和 superpi 1.2 计算100万位圆周率,其速度如下:

  软件名称               速度(单位秒)
  superpi 1.2            70
  apfloat 方法1          24.55
  apfloat 方法2          111.39 
  apfloat 方法3          107.62
  pifast 4.1 方法1       18.66
  pifast 4.1 方法2       26.39

  注意:计算时间均不包括输出时间
  我以前只知道计算PI的中最出名的2个程序,superpi 和 pifast,但是这两个程序均不能得到源代码,今天居然看到了非常了不起的计算pi程序的源码,真是非常高兴.
  如果大家对计算圆周率感兴趣,请参考陈先生的个人主页:http://www.jason314.com

Top
回复人: NewStarSE(新星) ( ) 信誉:99 2004-6-11 17:06:32 得分:0
 

哎,擂主怎么还不公开你的算法?

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-11 18:49:29 得分:0
 

既然是擂台,没必要那么快公布源码嘛,而且光是有源码你也不一定看得懂啊,总要整理一下,把原理告诉你,所以需要一点时间的,再说公不公布源码也是人家的权力嘛

to gxqcn(GxQ):
  睡了一觉,终于收到你的email了,最近精神欠佳,睡了一觉好多了:)我稍后测试一下能否调起来,谢谢你的lib,你很客气:)

to liangbch(宝宝):
  记得我学编程第一个想写的程序就是PI的程序,当时是用的几何方法,用正N多变形逼近圆计算其边长的和,结果由于精度问题你可以知道,很快就算不下去了。apfloat我2,3年前就知道了,可由于水平问题以前没精力看(也没水平),现在看还算不太吃力,非常清晰,文档里说apfloat乘法用的是FFT,我还没有仔细看

  
  关于分治法我是再另外一个较为简单的计算PI的程序里面看到的,也叫做karatsuba(据说是这个人最早发现的),其实我还是一个菜鸟,只是积累的东西多了,笨鸟多飞嘛。前阵子刚开始研究N!的时候找到一个Integer大整数类,后来疯狂搜索才发现它是GNU里的libg++的整数类,在GNU里还有一个高精度整数库是gmp,不过由于其代码依赖gcc和linux太多,很难移植到Windows平台,但其中有不少各种CPU和平台的汇编优化,以及x86的MMX, SSE, SSE2等代码,但由于gcc下的汇编格式和MS的不太一样,移植比较麻烦。另外还有一个库叫libI++,用的是karatsuba,光看代码就感觉效率不太高,不过程序容易读懂

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-11 20:14:08 得分:0
 

to  shines:
   我和你的经历有几分类似:
   1.和你一样,最早得到的分治法的算法 是在一个计算圆周率的程序中看到的,现在我使用的算法主要基于此,但有几处改进:
   A.将1000进制的16位表示法换成 10^9 进制,用32位代码重写。
   B.在做减法时,先比较大小,尽量减少求补数的次数。
   C.将分治法扩展,任意长度乘法有选择的递归调用分治法核心,实现了任意长度大数的快速算法。
  2。我也写过计算圆周率的程序,关于其算法介绍最早看的就是:http://www.jason314.com,其后使用马庭公司 用嵌入式汇编写了一个计算圆周率的程序,包含在我的超级计算器中(该程序最早于2000.10月发布于本站,我也看过其他计算圆周率的程序,包括4行代码计算圆周率至800位。
  3。关于大数运算库,我看到的有好几个,基本都是在csdn上得到网址的,写的好的不多,大约在1年以前,我在业余数学天地上看到了 NTL, 并粗略的分析了一下,发现其大数乘法运算用的是分治法,就没再深入。
  4。GMP 和 apfloat 的介绍是 在本帖看到了,前些日子下了一个GMP,发现代码非常庞大,可跨越多种硬件平台,但需要在gcc下编译,由于时间有限,暂时没有深入。
  5。今天看了一下 apfloat ,觉得是个好东西,咱们可共同学习一下,交流经验,共同提高。



Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-11 20:16:48 得分:0
 

顺便告诉大家我发布在csdn上的几个程序:
  1。超级计算器: http://www.csdn.net/cnshare/soft/3/3229.shtm
  2. 快速阶乘计算器:http://www.csdn.net/cnshare/soft/19/19260.shtm


Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-11 21:03:00 得分:0
 

to liangbch(宝宝):

阶乘的地址我早就知道了,不过下载不了,你自己可以下载试试看。

提示如下:

服务器停机说明:

服务器调整、维护,软件频道暂停服务。

如果你想找地方提高给人下载,不如你发给我,我帮你传到我的空间去,我的email我用CSDN的短消息发给你 


Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-11 21:53:18 得分:0
 

to gxqcn(GxQ):
   采用你的CHugeInt初步比较的结果:

             apfloat        CHugeInt

10000!       0.0708 s       0.0664 s
20000!       0.1737 s       0.2201 s
40000!       0.5263 s       0.7461 s
80000!       1.9049 s       2.5347 s
100000!      2.3294 s       3.9508 s



Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-11 22:00:03 得分:0
 

以上均采用一致的方法,经过测试,可以看到在15000!以前,CHugeInt快过apfloat,但大于15000!以后apfloat则高了不少,不过可以看出,你的计算N!的方法优于我的方法,完全可以抵消CHugeInt的劣势

采用我的这个方法,200000!的时间和HugeCalc几乎一致,但400000!略有超越,期待你的测试结果

Top
回复人: wsmagic(飞翔的心) ( ) 信誉:100 2004-6-13 5:37:56 得分:0
 

我是shines,由于大家都不回贴,只好请马甲出来秀一下了

to liangbch(宝宝):

估计我们看到的分治法可能不是同一个程序,但估计非常类似

我最近才把他由Base 10000的10进制,改为Base 65536(2^16)和Base 2^32的2进制,
但我发现Base 2^32比Base 2^16还慢,我暂时没弄Base 10^9,应该是比Base 10000快
相似的是,我也做了和你B,C一样的工作:)

你的calcFac.zip我已经收到了,我已经上传到:
http://www.77studio.net/files/calcFac.zip

宝宝的3.0版也是非常的快

我把我的方法改进了一下,稍快了一点,累乘比较小的数,我原来说尽量减少乘法次数,其实并没有完全做到,稍后我把2和5拿出来,这样会减少一些计算量,最后以10移位一下就好了,应该会又进一步,不过分治法在bit较少的时候还是稍快一些

Top
回复人: wsmagic(飞翔的心) ( ) 信誉:100 2004-6-13 5:40:30 得分:0
 

另外说一句,宝宝的calcFac.zip里提高了1.0, 1.1, 1.2的源码以及2.0, 3.0的EXE,值得一看

地址见上

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-14 9:07:23 得分:0
 

我准备公开我的快速阶乘计算器 2.0版的源代码,有需要者报上名来。

Top
回复人: fireofhell(懒鬼) ( ) 信誉:99 2004-6-14 15:31:49 得分:0
 

huming81@citiz.net

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-14 17:46:33 得分:0
 

阶乘,是求一组数列的乘积,其效率的高低,一个是取决于高精度乘法算法,二是针对阶乘自身特点算法的优化。
前者,是一种广为习知的技术,虽然不同算法效率可天壤之别,但终究可以通过学习获得,甚至可用鲁迅先生的“拿来主义”;
后者,则有许多的发展空间,这里我将介绍一种由我完全独创的一种方法,欢迎大家讨论,如能起到“抛砖引玉”的效果,那就真正达到了目的。
    
我设计“阶乘”时,始终遵循如下原则:
一、参与高精度乘法算法的两数,大小应尽可能地相近;
二、尽可能将乘法转化为乘方;
三、尽可能地优先调用平方;
    
言归正转,下面以精确计算 1000! 为例,阐述该算法:
记 F1(n)=n*(n-1)*...*1;F2(n)=n*(n-2)*...*(2 or 1);Pow2(r)=2^r;
有 F1(2k+1) = F2(2k+1) * F2(2k)
          = Pow2(k) * F2(2k+1) * F1(k),
   F1(2k) = Pow2(k) * F2(2k-1) * F1(k),
及 Pow2(u) * Pow2(v) = Pow2(u+v),
∴1000! = F1(1000)
         = Pow2(500)*F2(999)*F1(500)
         = Pow2(750)*F2(999)*F2(499)*F1(250)
         = Pow2(875)*F2(999)*F2(499)*F2(249)*F1(125)
         = Pow2(937)*F2(999)*F2(499)*F2(249)*F2(125)*F1(62)
         = Pow2(968)*F2(999)*F2(499)*F2(249)*F2(125)*F2(61)*F1(31)
         = Pow2(983)*F2(999)*F2(499)*F2(249)*F2(125)*F2(61)*F2(31)*F1(15)
         = ...
如果你预存了某些小整数的阶乘(比如这里的F1(15)),则可提前终止分解,否则直至右边最后一项为F1(1)为至,
这样,我们将1000!转化为2的整数次幂与一些连续奇数的积(或再乘以一个小整数的阶乘);

再定义F2(e,f)=e*(e-2)*...*(f+2),这里仍用“F2”,就当是“函数重载”好了:),
则 F2(e) = F2(e,-1) = F2(e,f)*F2(f,-1) (e为奇数)
∴ F2(999) = F2(999,499)*F2(499,249)*F2(249,125)*F2(125,61)*F2(61,31)*F2(31),
   F2(499) = ____________F2(499,249)*F2(249,125)*F2(125,61)*F2(61,31)*F2(31),
   F2(249) = ________________________F2(249,125)*F2(125,61)*F2(61,31)*F2(31),
   F2(125) = ____________________________________F2(125,61)*F2(61,31)*F2(31),
   F2( 61) = ________________________________________________F2(61,31)*F2(31),
   F2( 31) = __________________________________________________________F2(31),
∴ F1(10000) = Pow2(983) * F2(999,499) \
                         * [F2(499,249)^2] * [F2(249,125)^3] \
                         * [F2(61,31)^4] * [F2(31)^5]
这样,我们又将阶乘转化为了乘方运算。上式实际上是个形如 a * b^2 * c^3 * d^4 * e^5 的式子;

我们再将指数转化为二进制,可得到如下公式:
a * b^2 * c^3 * d^4 * e^5 = (a*c*e)*[(b*c)^2]*[(d*e)^4],
即可转化成可利用高效的平方算法。

最后,我再提供大家一个“累乘比较小的数”的技巧,这是我自创的,相信会对大家有借鉴作用:
一、应采用“倒序”,比方说,应按 999 * 997 * ... 的顺序;
二、可预先设定一个数组,记录如果当前起始数组的最大值不大于某个值,则连乘多少个数不会溢出;
三、可利用“几何平均值≤算术平均值”定理,可推推导出“ k 个正整数连乘不大于某个值,其充分条件是其和不大于某个临界值 Sum(max)”,我们只需要事先计算出它们的对应关系并保存下来。

上述算法是 Ver1.2.0.1 的算法关键点,其效率已略高于liangbch(宝宝)的高级算法3.0版。

在最新的 HugeCalc ver2.1.0.1 中,采用的是更彻底的分解成质因数的方法,技巧性倒反不如上面的强(但却需要有高效的素数筛法),但里面的很多思想仍在延续。

好了,今天暂到此为止,我得回家吃饭了。

申明:以上为原创,未经许可,严禁剽窃!!!

//====================================

注:我家无法上网,只能在公司才能与大家交流;

to shines(无尽的等待WOW魔兽世界,有点不耐烦了) :
   我没有对apfloat进行测试,但从其文档中发现其用的是FFT+CRT,所以高精度的乘法效率比较高,等我以后完全有空时才可能去实现它了,并相信不会比它的效率低(又说大话了)。

to liangbch(宝宝):
   可以把你的源代码给我一份吗?谢谢!

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-14 19:40:45 得分:0
 

我透露以下我的算法中用到的一些技巧。

  1. 在我的阶乘计算器3.0中,使用了GXQ  所说的 三条原则 的第一条:
   “一、参与高精度乘法算法的两数,大小应尽可能地相近",这样做的原因是:总的除法次数最少。

  2. 在我的阶乘计算器4.2中版(只有设想,未完成 ),使用了GXQ  所说的 三条原则 的第一条:
    一、参与高精度乘法算法的两数,大小应尽可能地相近;
    二、尽可能将乘法转化为乘方;
    三、尽可能地优先调用平方;

  3。在我的阶乘计算器3.2中(只有设想,尚未来得及写),同样想到了下面的GXM 的2个原则。
   一、应采用“倒序”,比方说,应按 999 * 997 * ... 的顺序;
   二、可预先设定一个数组,记录如果当前起始数组的最大值不大于某个值,则连乘多少个数不会溢出;

  4.通过我不断地改进算法的过程中,提高速度的关键是,提高 大数乘法运算的速度,其他技巧对于提高总的运算速度收效甚微,我相信,如果使用了FNT算法,当n很大时,计算速度会有质的改变。

  5。我的各个版本速度(包括计划中的版本)提高的关键是:改进乘法运算的速度,其各个版本乘法基本特点是:
   1.0 版,算法最简单,最朴素,使用1位数乘以多位数。
   1.1 版,将 int64 / int 32 改为用汇编语言书写
   1.2 版, 在1.1 基础上稍作优化,将几个较小的数相乘,得到一个10^9以内的数,再和多位数相乘。
   2.0 版 ,先将几个一位数 相乘,得到一个16个长度的 多位数,在和 先前的计算结果相乘,使用多位数 和多位数 相乘算法,大大较少了 除法指令的次数。
   3.X 版,将几个数相乘,得到一个长度 接近 2^n的数,尽量使用 两个均为 2^n次方的数相乘,当n交大时,使用分治法。
   4.X版,使用分解质因数算法,这个算法最早不是我提出的,但我已经想出了一个很好的实现方法,但从某些方面将,这不是一个好算法,因为1:算法复杂,2:对计算结果要求为指定精度(如精确到100位十进制数),效率不高,而其他算法,可精确到任一精度,且精度较越低,速度越快。
   5.X版,和 3.X 版 算法基本相同,只是大数运算引入FNT 算法。

   6. 说明:1000! = F1(1000)= Pow2(500)*F2(999)*F1(500)=...,此算法我也曾经想到,但没有深入,原因是这个算法有点复杂。现在看来,这种算法应该是一种不彻底的分解质因数算法。

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-15 16:50:40 得分:0
 

to liangbch(宝宝):
    让“参与高精度乘法算法的两数,大小应尽可能地相近”,这样做的主要原因是可使总的乘法指令变少(当两正数积一定时,其差越小,则其和也越小,这里的“和”从某个侧面可反映分的段数和;另外,可以更有效地发挥某些“高精度乘法算法”的效率)。
    另外,你说分解质因数算法“对计算结果要求为指定精度”?不甚明白。我好象不存在这个问题呵,是不是我理解的偏差?

to all:
    对于昨天下午发的关于高精度阶乘快速算法,我已重新整理,并发布在我的个人主页,及 CSDN 的文档中心,内容和形式略有改动,欢迎大家浏览。URL:

http://dev.csdn.net/develop/article/29/29226.shtm

http://maths.diy.myrice.com/florilegium/factorial.htm




Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-15 19:06:47 得分:0
 

to GXQ:
   让“参与高精度乘法算法的两数,大小应尽可能地相近”,这样做的主要原因是(我的观点):虽然乘法指令可能减小了,但这还是其次,重要的除法指令减少了。下面的例子均使用10^9表示法,rad(基数)为1000000000.

   例1:有4个数,n1,n2,n3,n4,其长度均为5,其最高位约为30000, 故任意两数相乘,积长度为9。下面两种算法的乘法和除法的次数为:
   方法1:
    step 1: p1=n1 * n2 :共 5*5=25次乘法, 10次除法,p1长度为9
    step 2: p2=p1 * n3 :共 9*5=45次乘法, 14次除法,p2长度为13.5
    step 3: 积=p2 * n4 :共 14*5=70次乘法,19次除法,p3长度为18
    //---------------------------------------------------------
                         140次乘法,  43次除法

    方法2:
    step 1: p1=n1 * n2 :共 5*5=25次乘法, 10次除法,p1长度为9
    step 2: p2=p3 * n4 :共 5*5=25次乘法, 10次除法,p2长度为9
    step 3: 积=p1 * p2 :共 9*9=81次乘法, 18次除法,p3长度为18
    //---------------------------------------------------------
                              131次乘法, 38次除法
               		
     可以看出方法2使用的乘法和除法指令都比方法一少。

     例2:有4个数,n1,n2,n3,n4,其长度均为5,其最高位略小于10^9,故任意两个数长度的积为10,则下面两种算法的乘法和除法的次数为:
   
  方法1:
    step 1: p1=n1 * n2 :共 5*5=25次乘法, 10次除法,p1长度为10
    step 2: p2=p1 * n3 :共 10*5=50次乘法,15次除法,p2长度为15
    step 3: p3=p2 * n4 :共 15*5=75次乘法,20次除法,p3长度为20
    //---------------------------------------------------------
                        150次乘法, 45次除法

    方法2:
    step 1: p1=n1 * n2 :共 5*5=25次乘法,   10次除法, p1长度为10
    step 2: p2=p3 * n4 :共 5*5=25次乘法,   10次除法, p2长度为10
    step 3: p3=p1 * p2 :共 10*10=100次乘法,20次除法, p3长度为20
    //---------------------------------------------------------
                        150次乘法,  40次除法
    可以看到,在例2中,乘法的次数没有减少,而除法次数则减少了,因为除法指令远比乘法指令耗时,所以方法2可以明显提高速度。
  
  说明:如果使用10进制(包括10^4进制,10^9进制)表示一个多位数,计算其乘积的过程必须做除法.设n1,n2 分别为被乘数和乘数(为简化起见,设两数长度相等),从高位到低,n1可表示为 数组a[i] (0<=i<n), n2可表示为 b[i](0<=i<n),则积可表示为数组p[i],0<=i<2n, p[i]=a[i-1]*b[0] + a[i-2]*b[1] + a[i-3]*b[2] ...+ a[0] * b[i-1],p[i] 是一个大于 rad(基数)的数,因此,必须将p[i] % rad 得到本位,p[i] / rad 得到进位,由于求余 和 除法可1步完成,因此 两个 n 位 * n位的大数相乘,至少做 2n次除法运算。
   

 “对计算结果要求为指定精度”等解释:
   你的阶乘计算器及大数库为完全精度,但是在实际应用中,我们可能不需要全部有效数字,而是仅仅需要部分有效数字,同时需要更高的计算速度,固一些大数库和应用程序提供了任意精度的计算,如apfloat 和mathematica,如果我们需要一个可计算任意精度的阶乘计算器,(我的超级计算器就可以提供任意精度的计算,windows自带的计算器仅提供32位有效数字的计算), 使用分解质因数算法,计算超大数 (如1000万)而需要很低的精度 (如32位)则太慢了,因为在这种情况下,用筛法求素数可能比乘法计算更费时间,而使用我的普通算法,稍加改进,在计算n!,可使计算时间反比于要求的精度时。


Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-16 17:28:05 得分:0
 

to liangbch(宝宝):
    素数的生成应该不成为瓶颈,我生成出1000万内的所有素数也仅需要0.030s。如果是不要求高精度,阶乘计算转换成对数运算应该是效率较高的。

to all:
    在CSDN 的文档中心发布的《高精度快速阶乘算法》,昨天因审核缘故可能无法浏览,现在已通过审核,可以正常浏览了。URL:

http://dev.csdn.net/develop/article/29/29226.shtm

http://maths.diy.myrice.com/florilegium/factorial.htm(我的主页版)

    另外,为了与大家更好的交流,以及对 HugeCalc 的广泛测试,如需其对应的头文件及库文件者,请留下联系方式,我会尽量满足。

Top
回复人: BaiYunfeng(Kidan) ( ) 信誉:100 2004-6-17 21:01:16 得分:0
 

请问你用的什么方法生成1000万内的素数啊?

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-18 9:54:28 得分:0
 

嘿嘿

我曾经试图做过FNT

发现流程太简单了,比FFT简单多了
就是要实现的模算法有点复杂,没做

嘿嘿,记录在一张纸上了,不知道夹在什么书中,忘记了

就是有一个缺点,位数限制的太死

不如NTT

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-18 9:58:01 得分:0
 

另外,对于小于100000000的质数,现代的计算机,只需要几秒就能算出全部清单

我记得我发表在这里一个贴子,每个小时是1000亿整数以上的素数搜索速度,可能还高,怕写错了 :)



Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-18 10:01:59 得分:0
 

另外,对于乘幂,有减少乘法次数的优化算法

比如6次方

x * x = x ^ 2

x ^ 2 * x ^ 2 = x ^ 4

x ^ 4 * x ^ 2 = x ^ 6


采用树状分解后,能得到比二进制分解还快速的算法!!(上边的例子不足以说明)

Top
回复人: williamdog(啊狗) ( ) 信誉:100 2004-6-18 12:23:23 得分:0
 

我 菜狗 十分  佩服 你!

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-18 17:13:14 得分:0
 

to yaos(夜夜秋雨孤灯下·无心人):
    你说的那个方法需要比较大的临时空间,其效率将被打折扣。

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-18 17:21:00 得分:0
 

昨夜,我在家翻一本数学辞典,终于明白了某些术语,现将本贴出现的专业术语列举如下:

1、FFT 快速傅立叶变换
2、FNT 费马数变换
3、NTT 数论变换
4、CRT 中国剩余定理(孙子定理)

Top
回复人: BaiYunfeng(Kidan) ( ) 信誉:100 2004-6-18 21:21:56 得分:0
 

最少乘的次数用DP来算就是了。
讲讲FFT,谢谢。

Top
回复人: mxfeng(地势坤,君子以厚德载物-永远都吃亏的人) ( ) 信誉:100 2004-6-21 17:24:29 得分:0
 

各位,能告诉我你们研究这个干什么吗?谁给你们钱?


Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-21 22:55:38 得分:0
 

http://www.77studio.net/files/Factorial_0621.zip

里面包含apfloat的和CHugeInt的程序,CHugeInt的已经和HugeCalc的基本一样快了,可以代为测试一下

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-21 23:10:49 得分:0
 

研究这东西(研究本身,不单指研究本贴的内容)本来就不在乎经济利益这个东西的,研究的深入会推动很多相关的技术的进步,这就是价值

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-22 9:21:52 得分:0
 

to shines:

    非常感谢你的贡献!可以将HugeCalc与著名的apfloat进行对比测试。

    看来你已完成对apfloat的改造,将base改成了10^9;apfloat在大整数乘法方面用到了NTT+CRT(Fast Number Theoretic Transform, FNT,注意此时的FNT并非“费马数变换”),其效率非常高,所以对非常大的阶乘计算非常有利。
    我的HugeCalc尚未采用这些高效算法(一是代码量大,二是理论还没吃透),所以瓶颈出现在高精度乘法方面,这在开发它时已有预见。

    另外,发现在存盘时,都是9位一节,但HugeCalc没有直接提供这样的接口,可是你将输出自己改造了一下?:) 还有,你是直接调用 HugeCalc.dll 的接口函数 Factorial(),还是自己的阶乘算法?

    BTW,我非常赞同你上贴的观点,许多基础研究虽然不能马上带来经济效益,但它的发展往往可推动许多技术领域的变革,有时还是翻天覆地的,真希望我们的国家能真正重视起来,希望某些决策者不要那么短视。。。(哎!现在对国内的学术氛围感到悲观,太多的功利行为,压制、剽窃、。。。,不便多言)

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-22 11:13:15 得分:0
 

我现在正准备写一个 大数阶乘计算器(极速版),它可以计算任意大数(目前暂定义在36位10进制)的阶乘 ,速度极快,不论n的为多大,速度应该可以控制在1ms以内,最低精度为14位有效数字,最高精度和 n 的位数有关, 如果n位K位数,可精确到 2K+1 位有效数字,例如求9841245454564564133312314567!,可精确到 57位有效数字,这个算法的核心是对数的高精度计算,我目前正在 设计对数的算法(我自己独创的).
   如果大家对此有兴趣,可以一起写这个程序,到时看看谁的更好。




Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-22 14:54:38 得分:0
 

阶乘计算有一个非常好的近似公式:
    n!≈(√(2πn))*((n/e)^n)*(1+1/(12n))
其中π=3.1415926535897932384626433832795...;e=2.71828182845904520536...
对上式取对数,而后再取逆运算即可。

我猜想 liangbch(宝宝) 的大数阶乘可能用的是类似上面的算法,不知对否?

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-22 15:59:40 得分:0
 

是这类似的公式,我用到的公式是:
     n!≈(√(2πn))*((n/e)^n)* e^(1+1/(12n)),实际计算中,可泰勒公式将e^(1+1/(12n))展开前2项。


Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-22 16:50:16 得分:0
 

呵呵,这个我也研究过的,我给一个更加近似的公式,其实我不知道怎么推导的,就自己试了几个系数

这是更准确的公式:
n!≈(√(2πn))*((n/e)^n)* e^(1+1/(12n)+1/(288*n*n)-139/(51840*n*n*n)),

typedef unsigned __int64 Integer64;

Integer64 factorial_approximation(unsigned long n)
{
    // gamma(n) ==> sqrt(2*pi*n)*n^n*exp(-n)
    double dn = (double)n;
    double dbVal = 0.0;
    //ulFloorPower2 = (unsigned long) (sqrt(2 * PI * dn) * pow(dn, dn) * exp((-1.0) * dn));
    //ulFloorPower2 = (unsigned long) (pow(dn, dn) * exp((-1.0) * dn) * sqrt(2.0 * PI * dn));
    dbVal = (pow(dn, dn) * exp((-1.0) * dn) * sqrt(2.0 * PI * dn));
    dbVal = dbVal * (1.0 + 1.0 / (12.0 * dn) + 1.0 / (288.0 * dn * dn)
        - 139.0 / (51840.0 * dn * dn * dn));// - 401734.0 / (2121322203.0 * dn * dn * dn * dn));
        //391.0 / (2269322.4489795 * n^4)
        //2879629312.0 / (15205637551104.0 * n^4)
        //- 391.0 / (87178291591.0));
    Integer64 uint64_f = (Integer64) dbVal;
    return uint64_f;
}

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-22 17:04:14 得分:0
 

to gxqcn(GxQ):

    的确,我发现HugeCalc的输出方式不够灵活,好在还可以处理一下,自己写了个处理函数,不知道有没有bug,再我的测试基本上没有问题(输出格式上),匆匆写完就发上来了

    本来想把我的思路写上来的,我想跟你的应该很接近的,因为时间都差不多,我还研究了2,3种质数分解的方案,我的程序还可以优化,没有使用从后面开始相乘,没有使用不等式的知识,不过经验判断,也快不了多少,所以稍后再试了

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-22 17:07:51 得分:0
 

注意:上两位的近似公式都给错了!(Why?)

检验方法:令 n = 10,用Windows自带的计算器分别计算精确值,及用近似公式得到的值。

这里,我提供大家一个链接以供参考:http://mathworld.wolfram.com/Factorial.html


Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-22 17:16:58 得分:0
 

公式我是copy你们的,呵呵,要看正确的结果,调一下我的函数就知道了:)

应该是这样:
n!≈gamma(n)=n^n * e^(-n) * sqrt(2*PI*n) * (1+1/(12*n)+1/(288*n*n)-139/(51840*n*n*n))

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-22 23:01:03 得分:0
 

sorry,不小心将公式写错了,我给出的公式来自 算法导论影印版本,这个公式一般叫做斯特林公式:
   sqrt(2*PI*n) * (n/e)^n * exp(1/(12*n)) < n! < sqrt(2*PI*n) * (n/e)^n * exp(1/(12*n)+1)。
    我们令 n! 左边的为fac1,右边的为fac2,则对于任意n,有:
    误差<fac1-fac2
    = exp(1/(12*n)-exp(1/(12*n)+1)
    ≈1/12*n-1/(12n+1)
    ≈1/(144 *n^2)
    误差<1/(144 *n^2)
    所以该公式仅仅适用于n比较大的情况,当n较小时,有较大的误差。

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-23 2:50:50 得分:0
 

工程上一般用这个计算:

long  double  GammaLn(  long  double  x  )  
{  
   if  (  x  <  0  )  exit(  0  );  
   if  (  x  ==  0  )  return  0;  
   long  double  xbuf,  y,  tmp,  ser;  
   static  long  double  cof[6]  =  
   {  
       76.18009172947146,  -86.50532032941677,  24.01409824083091,  -1.231739572450155,  0.001208650973866179,  -0.000005395239384953  
   };  
   int  i;  
   y  =  xbuf  =  x;  
   tmp  =  xbuf  +  5.5;  
   tmp  -=  (  xbuf  +  0.5  )  *  log(  tmp  );  
   for  (  ser  =  1.000000000190015,  i  =  0;  i  <=  5;  i++  )  
   {  
       y++;  
       ser  +=  cof[i]  /  y;  
   }  
   return  (  -tmp  +  logl(  2.5066282746310005  *  ser  /  xbuf  )  );  
 
}  
 
//---------------------------------------------------------------------------  
long  double  Gamma(  long  double  x  )  
{  
   long  double  temp1,  temp2;  
   if  (  x  <  0  )  
   {  
       if  (  int(  x  )  ==  x  )  return  1.0e300;  
       temp1  =  fabs(  x  );  
       temp2  =  expl(  GammaLn(  temp1  )  );  
       return  -PI  /  (  x  *  temp2  *  sinl(  PI  *  x  )  );  
   }  
   else  
       return  (  expl(  GammaLn(  x  )  )  );  
}  
还有其他的办法,精度要高一些。

但是如果用VC编译的可执行文件,浮点数限制在1.0e300多一点,结果再大会溢出,其他编译器,结果可以求到1100!

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-23 2:54:37 得分:0
 

比如:
1100!= Gamma( 1101 ) = 5.34370848809264E2869

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-23 3:00:29 得分:0
 

Gammaln( 1e8 ) = 1742068066.10383

Top
回复人: geeksky(§illuSioN§) ( ) 信誉:100 2004-6-23 13:24:23 得分:0
 

同意楼上的,求阶乘用伽玛函数好

Top
回复人: freeia(后知后觉) ( ) 信誉:100 2004-6-23 13:55:49 得分:0
 

"我用 C++ 开发了一套《超大整数高精度快速算法》,前后花了十余年。"
就这么一点东西,花了楼主10余年的业余时间,能力严重不足。

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-23 14:47:30 得分:0
 

你说的不错。我说你没说错,并不表明我就认同了你对我的评价!

    首先我要告诉你,十余年里,前四年不得不学见到就头大的化学(什么无机化学、有机化学、物理化学、高分子化学。。。),而后又在工厂晃荡了近四年,期间没有上过网,更不知email、BBS为何物,周围没有可交流的同事,也找不到专业的书籍,一切全凭自创,而且全是业余时间。请见我在帮助文档 HugeCalc.chm 里写到其发展历史:
    1、1992-1996 大学期间,在 DOS 下用 BASIC 编程;
    2、1997-2000 工作期间,想方设法接触计算机,没有机会的日子就在自己的“人脑”里运行改进代码;
    3、2000-2002 进入电脑公司,在 Windows 下编程;利用业余时间将该套算法全面用 Microsoft Visual Basic 6.0 改写,并成功开发了PowCalc,被中国幻方研究者协会推荐使用;
    4、2003-2004 全面用 Microsoft Visual C++ 6.0(及最新的补丁版 VS6sp6) 改写该套算法,速度和性能得到大幅提升。

    开发这套算法的出发点是为了自己方便研究一些组合数学、数论问题。

    其实,该套算法,真正快速发展是在今年春节之后,在将它公开发布、与广大网友广泛交流之后,尤其在是发现了 CSDN 这块藏龙卧虎之地之后。

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-23 17:41:58 得分:0
 

佩服楼主的恒心!任何事情都是一样,没有恒心和毅力光凭聪明很难做出点什么的。
不知道楼主是否打算把它封装起来,建立一个类,然后for VC,for BCB,for Delphi等等,服务于对此有兴趣的人员。



Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-24 8:34:43 得分:0
 

非常感谢 kerbcurb,你的话让我联想到《世界需要热心肠》歌词所唱(一句知心的话语,也许赛过万钧雷霆):

    一个篱笆三个桩,一个好汉三人帮。谁也难免遇到险阻,谁也难免遇到忧伤。只要你我热心相助,懦夫也会变成金刚,只要你我热心相助,懦夫也会变成金刚。一句知心的话语,也许赛过万钧雷霆;一句亲切的呼唤,能有起死回生的力量。只要你我热心相助,懦夫也会变成金刚,只要你我热心相助,懦夫也会变成金刚。人生道路多曲折,人生道路又漫长。只要你我热心相助,懦夫也会变成金刚。为了一切都幸福,世界需要热心肠。

    不敢自诩好汉,但一样需要“你我热心相助”,这个世界需要热心肠。

    我现在已封装成一个类,具有很好的扩展性,可以方便地向有理分数,任意精度的实数推广。现在已for VC,由于对 BCB、Delphi等不了解,所以暂无法保证与之兼容。

Top
回复人: NowCan(((((( ★ )))))) ( ) 信誉:100 2004-6-24 12:37:18 得分:0
 

BCB我倒是很熟,Delphi我帮不了你了。

关于兼容问题,如果用静态库,VC和Borland的格式不同,根本无法使用。如果用DLL,好像又不能导出类,我觉得像这种纯计算的东西没必要用类。

Top
回复人: wenminghu(看到人生的希望) ( ) 信誉:100 2004-6-24 13:19:27 得分:0
 

还是建议楼主多多看看别人的代码,有时候自己埋头苦干是很幼稚的想法,你要走的路别人可能早已经走过了。可以借鉴别人的现有成果。

Top
回复人: wenminghu(看到人生的希望) ( ) 信誉:100 2004-6-24 13:22:14 得分:0
 

什么都自己写并不代表自己就水平很高,能把别人的长处学到手,为己所用,然后根据需要加以改进才是高手。
例如,所谓算法,可以找一个非线性数学方向加以研究,而不是低层次重复研究。
amd公司的主页有一个数学库算法的连接,别人是根据最新的SSE2, SSE3写的优化代码,倒不如和他比较一下,就知道自己工作的价值和意义了。

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-24 13:50:22 得分:0
 

我在我的超级计算器中实现了四则运算和指数、对数、三角函数、反三角函数的算法(四则运算可精确到几万位,函数运算设定为300位以内),但自认为对数和三角函数的算法速度不够快,因此一直希望得到更好的算法,昨天试了一下apfloat(其算法暂时没有弄懂),计算 对数 至1万位需要4-5秒(不包括常数的初始化),看了这个算法也不是顶级的算法,请问各位高手,有谁知道更好的对数算法或者网络资源。


Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-24 17:43:14 得分:0
 

HugeCalc bug 警报:

现象:6/2=2, 10/2=4, 14/2=6 9/3=2 15/3=4 ... 
几率:被除数、除数均为 CHugeInt 类型,且 除数<10^9 时。 
原因:因为我将主要精力放在超大整数的算法上,对一些常规整数的计算测试不足。
状态:Fixed

麻烦各位到我的主页下载最新版,已得到库文件的朋友,原 .h 和 .lib 文件仍可使用。

    上述 bug 是 HugeCalc 外部反馈回来的第一个bug! 发现者是一位北京航空航天大学博士研究生,通讯专业,在做编码解码的时候经常要用到大数的四则运算,试用了我的 HugeCalc 感觉非常好(原文:“库已经收到,确实很好用,速度很快,是我使用过的最好的库了”)。今天下午他给我反馈了上述 bug。

    为了尽可能的消除bug,特悬赏如下:如果发现我的软件出现bug,请及时反馈给我,可通过email,或直接在本论坛上发贴,并请详细描述bug现象,一经证实,本人将另开贴给首发现者送分,100分/bug,送完为止。


Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-24 20:56:19 得分:0
 

liangbch(宝宝)我对你的工作很感兴趣,是否可以公开部分源码?
kerbcurb@hotmail.com

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-25 0:40:04 得分:0
 

刚刚完成了一个可以计算ln(2)和ln(10)的程序,计算前者需要0.75s,后者需要0.62s,(CPU PIV 2.6G),
   算法:
   lg(2)=lg((1+1/3)/(1-1/3))
   lg(10) :自创的,算法如下:
     step 1. 使用反正切公式 计算 x1=ln(4/3)=ln((1+1/7)/(1-1/7)) 
     step 2. 使用反正切公式,计算 x2=ln(65610/65536) =ln((1+37/65573)/(1-37/65573))
     step 3. 计算 ln(10)=ln(65610/6561)=ln(65610/65536 * (4/3)^8 )
              =ln(65610/65536)+ 8*ln(4/3)
              =x2 + 8*x1 


Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-25 1:51:01 得分:0
 

我今天才看到liangbch(宝宝)6月15号寄给我的一个高精度N!的文档,我很少收mail,如果寄了记得提醒我,粗看了一下,写得还可以,跟你说的一样,用处不是很大,但是还是比较佩服作者的能力

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-25 8:19:25 得分:0
 

楼上的ln2计算,我晕了,多少位??

有逼近算法的,很快阿



Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-25 8:36:18 得分:0
 

对不起,忘了说精确到多少位了,是10000位。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-25 8:44:55 得分:0
 

慢死了 嘿嘿

你算法也许没问题,程序有问题吧

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-25 21:48:42 得分:0
 

请问 liangbch(宝宝) 和gxqcn(GxQ) 两位朋友1000位的整数乘以2000位的整数,大概需要多长时间,我的程序在P4 1.7 / 256下需要170~180个TickCount,我使用数组实现的,包括输出时间

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-26 8:36:44 得分:0
 

楼上,讨论你的问题首先要统一基于的进制,你是采用10^9进制,还是2^32进制?

再纠正一下gxqcn(GxQ)的一个错误的看法,apfloat本身就是基于10^9,base=10的时候,base=2的时候是2^30,我并没有修改对apfloat做任何的修改,值得注意的是HugeCalc计算到400000!的时间,我基于apfloat写的程序可以计算到950000!,在我的机器上

还有一点当n<20000左右(大概是这个值吧,程序改动过后具体值忘了)的时候,不使用分段相乘会更快一些,基本接近HugeCalc的速度,至少比现在的快一些,现在的是比较慢的,由于懒,没有写进去,所以当N小的时候比较慢一些

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-26 18:26:27 得分:0
 

shines(无尽的等待WOW魔兽世界,有点不耐烦了) 
采用10^9进制,我只是做了成法和加法,数组存储,10000位乘10000位数据包括输出3000个Tickcount,我不清楚TickCount与秒之间的转换,感觉在1秒多

Top
回复人: shines(郭子) ( ) 信誉:100 2004-6-26 18:39:42 得分:0
 

TickCount? 是GetTickCount()函数获得的吗?3000TickCount=3000ms=3s

你用的不会是Linux吧,Windows下,GetTickCount和timeGetTime都是以ms为单位的

位数大的时候最好用分治法,还有10000位是数字,还是10000个BYTE?

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-26 23:56:01 得分:0
 

谢谢 shines

10000位是普通十进制数字

我不会分治法,略微知道一些apfloat。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-6-27 9:14:12 得分:0
 

10000位整数用普通乘法效率低些,用FFT/FNT复杂些
分治好点

应该在0.001秒内完成乘法

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-28 11:35:18 得分:0
 

to kerbcurb():
    正如 yaos 所说,可在 0.001秒内完成计算并输出(测试计算:500!x900!,分别为 1135、2270 位)。

to liangbch(宝宝):
    如果需要同时得到 log(2)、log(10),后者可用如下方式得到,速度应更快:
        log(10)=log(8*( (1+1/9)/(1-1/9) )) = 3log(2) + log( (1+1/9)/(1-1/9) ),
你说呢?

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-28 13:11:21 得分:0
 

to gxqcn:
  单独计算ln(2)和ln(10)用我的方法,ln(10)更快一些,如果我同时计算ln(2)和ln(10),我会使用下面的方法:

  ln(10)=ln(65610/6561)=ln(65610/65536 * (4/3)^8 )
        =ln(65610/65536)+ 8*ln(4/3)
        =ln((1+37/65573)/(1-37/65573))+8*ln((1+1/7)/(1-1/7)) 

  ln(2) = ln(1024)/10= (ln(1000)+ln(1.024))/10
        =( 3ln(10)+ln(1+3/125))/10
        =( 3ln(10)+ln((1+3/253)(1-3/253)))/10

   这里迭代次数最慢的是  ln((1+1/7)/(1-1/7)),其次是ln((1+3/253)(1-3/253)),最快的是ln((1+37/65573)/(1-37/65573)),
   以精确到10000位,共需要 10000/log10(49)+10000/log10(7112)+10000/log10(3140846)=10051次迭代。
  
  而使用先计算ln(2),后计算ln(10)的算法,则需要
  计算 ln((1+1/3)(1-1/3)) 和 ln( (1+1/9)/(1-1/9) ),以精确到10000位计算,需要
  10000/log10(9)+10000/log10(81)=15718次迭代。



Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-28 13:14:16 得分:0
 

我又想了想,我的方法计算稍复杂,由于分子不为1,计算时需要做乘法,但具体到底哪个更快,需要程序来验证。

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-6-28 22:07:33 得分:0
 


我编了程序,对200以内的数进行搜索,得到比较好的计算对数的方法:
同时求2,3的对数,比较快的方法
	(4/3)^2* (9/8)^1 =2
	(4/3)^3* (9/8)^2 =3
同时求2,3,5的对数,比较快的方法
	(9/8)^3*(10/9)^2 * (16/15)^2 =2
	(9/8)^5*(10/9)^3 * (16/15)^3 =3
	(9/8)^7*(10/9)^5 * (16/15)^4 =5


Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-29 1:49:23 得分:0
 

18000位*18000位
351tickcount = 0.34秒
包括输入输出时间

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-6-29 17:31:41 得分:0
 

to kerbcurb():效率应至少还有十倍的提升空间。

to liangbch(宝宝):
    你作的工作很细致,由你给出的对数公式,其速度应有很大的提高,因为展开后分子均为1,分母也较大,且尽可能地共享数据。

另外,我今天提出了另一个问题,欢迎大家讨论,URL:
    http://community.csdn.net/Expert/topic/3131/3131084.xml?temp=.3031732

Top
回复人: kerbcurb() ( ) 信誉:100 2004-6-30 9:49:01 得分:0
 

我改进了算法,18000位*18000位约1秒,10^9进制,使用unsigned __int64类型存储数据

Top
回复人: booming(All) ( ) 信誉:73 2004-7-2 22:26:42 得分:0
 

不必用实际时间.
一般用时间复杂度比较好.

问一下,两个N位的整数相乘你的算法效率可以少于Tn=O(N^ln3)吗?
如果不行,这个软件没有创新.


Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-5 1:12:14 得分:0
 

楼上你说的是O(n^log2(3))吧,那FFT的O(n*log(n))算不算创新啊?可是fft很多程序都可以,虽然目前我还不会。。。

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-7-5 10:24:10 得分:0
 

to booming(All):
    请见我于 2004-06-07 10:41:14 提供的测试报告,当结果位数翻倍时,耗时基本要增加到3倍左右(经更严格的计算,其比值略还小于 log2(3) )。

    我现在正紧锣密鼓地开发下一版,其乘法算法复杂度将为 O(n*log(n)),预计在八月中旬正式发布,为此可能要搭上在此期间的所有双休日和零星业余时间。我将尽可能完成该计划,因为我的儿子将在八月十九号左右出生(注:因为做B超的医生是个熟人,所以提前知道了是个胖小子),否则新版将会delay,delay ...

Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-5 14:14:15 得分:0
 

恭喜楼主将要喜得贵子,偶GF都还没有出世,还是先把嫂子照顾好,这事可不能马虎啊,估计一辈子也就那么一两此,回家多干点家务活,程序是要写的,老婆也是要兼顾滴(叫照顾比较好)。

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-7-8 13:08:26 得分:0
 

这几天,shines一直在不遗余力的改进他的 大数 乘法算法,参照大家的建议,他已经实现了速度非常快的乘法子程序,参照帖子 315200,3121009。

   我好早就对大数乘法计算作了"研究”,1996年的一期 新浪潮 杂志 刊登了一个文章,内容就是讲述大数乘法的,作者号称比目前最好的算法快100多倍,我仔细看了代码,发现其算法并非很先进,经过 将核心代码用8086(那时我用的编译器是TURBO C 2.0,目标代码为16位)重写,速度提高了60倍,一时高兴,整理成文章,寄到杂志社,文章虽未发表,但却得到一年的杂志。

  关于FFT算法,我曾经尝试过,使用普通的FFT算法,即
   1。对被乘数进行FFT变换,每个复数表示4位10进制数
   2。对乘数进行FFT变换,每个复数表示4位10进制数
   3。计算内积并除以n(长度)
   4. 对积进行逆变换。

但有2个问题未能解决,
   1:使用的内存非常之多。n位十进制数用需要16n字节。 
   2。浮点数的累计误差使得当数的位数很长时,结果错误。
  最终发现,在当长度很小时,速度很慢,长度很大时又会因累计误差导致结果错误,使用FFT算法意义不大,只好作罢.
   FNT算法早有所耳闻,但不知其算法,借了一本 陈琦 60年代的书,《快速数论变换》
发现:这是一本纯数学的书,并未提出一个可直接编码的算法,而且难以读懂,只好作罢。因此最快的算法只能做到分治法,另一种手段就是用汇编优化代码,其核心就是 n * m位数需要mn次乘法,m+n次除法。

   我以前的阶乘计算器 基本算法 就是基于我给出的汇编代码,但是主要用c实现,没有完全用汇编优化,本来我想写一个汇编内核的乘法核心,但我这几天一直在琢磨 对数 、指数 、正弦、反正弦的快速算法,这个优化 乘法算法的计划就搁置了,因为基于同样的算法,所有我的乘法运算的速度应该和你的差不多,你在阶乘算法上稍稍想写办法,我相信你可以写出 一流速度 计算阶乘的程序。 我把 我,你、GXQ 的写的使用分治法计算阶乘的这类程序叫做 一流速度的程序, 如果我们在此基础上,再加上 FNT 算法,在计算乘法时根据乘数长度 来采用不同的乘法核心,包括汇编代码的硬乘法,分治法,FNT算法,那么,我们的算法就可称为超一流的算法了,也就是理论界速度最高的算法了。


Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-7-8 13:09:36 得分:0
 

应该是:"参照帖子 3150200,3121309"


Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-9 7:11:08 得分:0
 

讨论:一个基于10^9进制的大数乘法汇编代码,大家来批批
http://community.csdn.net/Expert/topic/3150/3150200.xml?temp=.4317591

汇编的大数运算核心算法
http://community.csdn.net/Expert/topic/3121/3121309.xml?temp=.86237612

今天改正一个小错误,终于放心了,关于FFT/FNT是我的下一个目标,FFT的浮点误差是可以修正的,但具体怎么修正我也不知道,但看介绍都是这么说的(修正的步骤是必须的),至于FNT,至少有一个apfloat可以参考,大不了,狠一点,全搬下来,然后自己消化(这是最后的没有办法才用的办法)

Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-9 7:12:07 得分:0
 

其实看论坛,很多人都懂FFT做乘法,懂的不妨开个讲座,偶来听听

Top
回复人: zairwolf(君子兰) ( ) 信誉:100 2004-7-9 11:21:41 得分:0
 

大家都这么强,我跟个帖子。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-11 20:00:53 得分:0
 

FFT关键是采用64比特的算法,中间结果80位,至于整数,要小于8比特才合适

这个能避免误差,适合几百万十进制数字的乘法了

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-11 20:12:39 得分:0
 

是避免误差!!!!

如果采用10进制,别超过100

Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-11 22:53:12 得分:0
 

在FFT中,误差a的计算公式为:
a <= 6*n^2*B^2*log(n)*e   (对于double双精度的浮点数,e约等于1.0E-16)
(B为基数,n为计算的长度(一般是2的幂次方))

为了是获得准确的值,误差a最好小于0.5,小于0.25更佳,但上面的公式是考虑最坏的情况下,实际的误差要小得多,所以一般一个复数可以代表3个或4个十进制数,即B=1000,或者B=10000

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-7-14 14:44:44 得分:0
 

我已完成了一个基于 NTT+CRT 的大整数高精度快速算法,
其核心模块参考了 apfloat,并作了适当改进:
1、改变数据存储顺序,使低位在前,高位在后,带来的好处是对齐容易,进位时不需移动数组;
2、改写了mod运算类,使之更高效、更方便。

测试的结果为,相对前一版,速度有很大提高,
比如,精确计算一百万的阶乘,在 P4 1.7GHz / 256MB RAM 的机器上需要 46.547s,
用 apfloat 则需 17.113s(用 shines(郭子) 提供的 FactorialApfloat.exe )

经分析,造成 apfloat 速度比我快的原因是,
其对耗时比较大的模块用 FPU (浮点汇编)进行了优化,
所以其大数运算速度基本上保持在我的2-3倍左右。

而我则未用任何汇编优化,也未调用浮点指令的函数。

看来,我下一步的努力的方向就在此了。。。


Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-15 1:44:16 得分:0
 

终于实现目前最快的分治法大数乘法及阶乘程序,我的大整数类(基数为10^9)目前只实现了乘法,乘方,左移。先贴个计算N!(阶乘)的计算时间(赛扬II 900),明天再详细比较,先睡了

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 5000! :

Calc N! used time: 0.00990573 second(s).
Output used time : 0.00843068 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 10000! :

Calc N! used time: 0.03084414 second(s).
Output used time : 0.01426829 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 20000! :

Calc N! used time: 0.10197665 second(s).
Output used time : 0.02734593 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 40000! :

Calc N! used time: 0.32939969 second(s).
Output used time : 0.05489971 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 80000! :

Calc N! used time: 1.10303013 second(s).
Output used time : 0.12778188 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 100000! :

Calc N! used time: 1.64062807 second(s).
Output used time : 0.14348751 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 160000! :

Calc N! used time: 3.60637567 second(s).
Output used time : 0.23549390 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 200000! :

Calc N! used time: 5.43043068 second(s).
Output used time : 0.30177243 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Calculate 400000! :

Calc N! used time: 17.82116330 second(s).
Output used time : 0.63812887 second(s).

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-15 1:49:39 得分:0
 

用时大概是HugeCalc用时的75%左右

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-15 8:00:20 得分:0
 

嘿嘿,GMP
AMD750 
10000! 计算 9ms 输出24ms

Top
回复人: pomelowu(羽战士) ( ) 信誉:100 2004-7-15 14:02:39 得分:0
 

gz

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-7-15 14:46:13 得分:0
 

我刚刚对阶乘计算器3.0 输出算法作了改进(计算过程没有改),输出速度从0.70s 提高到0.25s(输出格式和GXQ的完全相同,输出到txt文件),以下是100000!的输出结果。如有感兴趣者,可尝试重写一个算法,看看能否超过我?

 please input a num, 0:exit
n=?100000
计算过程耗时1.619099s
输出过程耗时0.024702
please input a num, 0:exit
n=?100000
计算过程耗时1.621863s
输出过程耗时0.024941
please input a num, 0:exit
n=?100000
计算过程耗时1.621124s
输出过程耗时0.025613
please input a num, 0:exit
n=?

Top
回复人: wsmagic(飞翔的心) ( ) 信誉:100 2004-7-15 15:50:21 得分:0
 

to yaos(等待WoW前玩私服吧·无心人):
  GMP是基于2^32的吧,你输出之前要转换成10进制才能显示(这个很耗时间的,N小的时候还能忍受,试试把N算大一点看看),我们的输出时间是写到文件的用时,数据都是基于10^9存储的,根本没有转换的过程

  如果GMP的是基于2^32或2^16的,你试试计算200000!看看(还不够大就试试1000000!),一定吓死你,即使它的进制转换算法很快!

Top
回复人: wsmagic(飞翔的心) ( ) 信誉:100 2004-7-15 15:53:27 得分:0
 

我是shines, 不知道你的输出时间是不是包括转换进制的时间,还是包括了其他诸如显示输出或输出到文本的时间,

to liangbch(宝宝):
  你的输出的确很快啊,你的CPU是多少的啊。。。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-15 16:55:15 得分:0
 

说句玩笑:

如果我用GMP计算1000000!的阶乘,计算时间确实要吓死我,不过不是太多,而是太少
嘿嘿

我为了能上网,安装的是红旗,不知道能不能安装GMP,晚上我试一下
估计输出时间为10秒,计算时间小于10秒

等我的好消息

另:最近在培训学生打字,领导不让上网,连机器都让给了学生,嘿嘿,惨,已经10天没玩游戏了
哈哈

打算开帖子讨论FNT,谁有兴趣??

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-15 16:59:28 得分:0
 

还要声明:

我一直拿GMP做例子,不是否定各位大哥的功劳,是证明2进制好,嘿嘿
我可能不能做到GMP的水平,连各位大哥都赶不上呢,嘿嘿

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-7-15 18:05:48 得分:0
 

to wsmagic(飞翔的心) :
   我的时间是输出总时间,包括数转化成字符串,每10个数字加1个空格,每50个字符加1个回车换行等格式化操作这些计算时间,和打开文件,写、关闭文件 这些I/O时间。
   使用的技巧:
   1.使用了相当快的数转化为字符串算法,过去的方法是“fprintf(fp,"%d",arr[i]);”
   2.文件写入使用fwrite函数,一次写64k左右.
硬件配置:PIV 2.6G,80M硬盘。
  稍候请收我的邮件。


Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-7-15 18:29:48 得分:0
 

Gxq,shines:
  我给你们发送了邮件,请查收。

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-7-15 18:48:44 得分:0
 

我又测试了输出时间,发现其并非和文件大小成正比,而是有一个基本的开销,可能是打开文件比较耗时,所有当n较大时,才能反映其真实的性能。
  值        文件大小   输出时间     平均每M字节耗时
                    (包括转化时间)  
  10000!   40K        0.018         0.45s  
  100000!  500K       0.025         0.050s
  200000!  1065K      0.029         0.027
  400000!  2261K      0.048         0.021

下面是抓取得屏幕内容:
please input a num, 0:exit
n=?10000
计算过程耗时0.029097s
输出过程耗时0.018407
please input a num, 0:exit
n=?100000
计算过程耗时1.629555s
输出过程耗时0.024947
please input a num, 0:exit
n=?200000
计算过程耗时5.807814s
输出过程耗时0.028452
please input a num, 0:exit
n=?400000
计算过程耗时21.038139s
输出过程耗时0.051946
please input a num, 0:exit
n=?100000
计算过程耗时1.632746s
输出过程耗时0.024607
please input a num, 0:exit
n=?200000
计算过程耗时5.796285s
输出过程耗时0.028658
please input a num, 0:exit
n=?400000
计算过程耗时20.979330s
输出过程耗时0.044438
please input a num, 0:exit
n=?

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-7-15 20:17:43 得分:0
 

我终于找到了不用 FPU 速度会慢一半的依据:
Using the FPU the modular multiplication takes about 33 clock cycles on the Pentium whereas using the integer unit it would take about 52 clock cycles.

我的 FNT 在 CarryCRT 过程中用的是更有效的手段,可惜我对 FPU 不熟,否则一定可以比 apfloat 速度还快!不过,现在速度已完全可与之保持一个小常数的比值,也就是说,在算法复杂度同一水平。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-15 21:36:15 得分:0
 

#include <stdio.h>
#include <gmp.h>

int main(void)
{
unsigned long n;
mpz_t f;

mpz_init(f);

printf("Please Input a Number: ");
scanf("%d", &n);
if (n <= 1)
  printf("The %d! is 1", n);
else
  {
  mpz_fac_ui(f, n);
  printf( "The %d! is ");
  gmp_printf( "%Zd\n", f);
  }

return 0;
}

1000000! 计算过程23秒,总过程95秒,输出到文件不知道为什么不行,是输出到屏幕的时间

P4 2.0 256DDR,嘿嘿


Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-15 21:51:04 得分:0
 

嘿嘿
FPU有隐含的处理时间,用了你就知道了
节约的时间有限

上边的GMP的输出估计是因为数值太大的原因,实际的转化过程应小

Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-15 22:23:53 得分:0
 

to liangbch(宝宝), gxqcn(GxQ):
  你们的信我都收到了,今天有点困,先睡了

to yaos:

你可以用手表或者手机估计一下转换的时间,看到开始输出就停止计时

我的分治法(赛扬II 900):
Please input a number [0 = exit]:
N = ?1000000

Calc N! used time: 83.26649840 second(s).

Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-16 8:37:57 得分:0
 

to yaos:
  昨天晚上我还担心你自己写的阶乘算法会有些慢,原来GMP自己写了阶乘的函数mpz_fac_ui(f, n);,总的来说GMP真的有点慢,计算时间都有些不太快(不知其算法怎么样,不好说,相信不一定是最好的),转换的过程就自然不必说了,一定是长过计算时间的

to liangbch(宝宝):
  我上面问“不知道你的输出时间是不是包括转换进制的时间,还是包括了其他诸如显示输出或输出到文本的时间,”是问yaos的,不过现在知道他是包括转换进制和显示屏幕的时间,
  我的输出时间和你是一样的,都是包括转换成字符串,格式化,打开文件,写、关闭文件 这些I/O时间

你的输出真是快到不行了,看来充分利用缓冲区是一个很好的技巧,“2.文件写入使用fwrite函数,一次写64k左右.”,看我的测试结果,看看10000的时候,够快吧,怎么比你上面测试的还快啊:(

please input a num, 0:exit
n=?10000
计算过程耗时0.058362s
输出过程耗时0.002710
please input a num, 0:exit
n=?100000
计算过程耗时3.389233s
输出过程耗时0.020861
please input a num, 0:exit
n=?200000
计算过程耗时12.100887s
输出过程耗时0.041514

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 10:55:06 得分:0
 

有点不公平啊,GMP输出到屏幕,你们输出到文件

嘿嘿,windows系统下,输出到屏幕一个近5M的文本,恐怕100秒不够吧
不很明白linux的屏幕输出原理,肯定比windows快很多,但是肯定比文件操作慢

用手表也不准确,输出前,可能有隐含操作

我找GMP转换函数直接转换吧

嘿嘿

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 11:21:25 得分:0
 

shines啊,是全部95秒啊,计算时间才23秒啊

刚做输出到字符中,105秒,郁闷啊,嘿嘿,转换时间105-23=82,倒!!
linux下top查看的,99%以上CPU时间,还是比较准确的,不过内存占用一直在变,有动态内存管理,哎

下边是他的原代码,计算fact的函数

/* mpz_fac_ui(result, n) -- Set RESULT to N!.

Copyright 1991, 1993, 1994, 1995, 2000, 2001, 2002 Free Software Foundation,
Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"


/* Enhancements:

   Data tables could be used for results up to 3 or 4 limbs to avoid
   fiddling around with small quantities.

   The product accumulation might be worth splitting out into something that
   could be used elsewhere too.

   The tree of partial products should be done with TMP_ALLOC, not mpz_init.
   It should be possible to know a maximum size at each level.

   Factors of two could be stripped from k to save some multiplying (but not
   very much).  The same could be done with factors of 3, perhaps.

   The prime factorization of n! is easy to determine, it might be worth
   using that rather than a simple 1 to n.  The powering of primes could do
   some squaring instead of multiplying.  There's probably other ways to use
   some squaring too.  */


/* for single non-zero limb */
#define MPZ_SET_1_NZ(z,n)       \
  do {                          \
    mpz_ptr  __z = (z);         \
    ASSERT ((n) != 0);          \
    PTR(__z)[0] = (n);          \
    SIZ(__z) = 1;               \
  } while (0)

/* for single non-zero limb */
#define MPZ_INIT_SET_1_NZ(z,n)                  \
  do {                                          \
    mpz_ptr  __iz = (z);                        \
    ALLOC(__iz) = 1;                            \
    PTR(__iz) = __GMP_ALLOCATE_FUNC_LIMBS (1);  \
    MPZ_SET_1_NZ (__iz, n);                     \
  } while (0)

/* for src>0 and n>0 */
#define MPZ_MUL_1_POS(dst,src,n)                        \
  do {                                                  \
    mpz_ptr    __dst = (dst);                           \
    mpz_srcptr __src = (src);                           \
    mp_size_t  __size = SIZ(__src);                     \
    mp_ptr     __dst_p;                                 \
    mp_limb_t  __c;                                     \
                                                        \
    ASSERT (__size > 0);                                \
    ASSERT ((n) != 0);                                  \
                                                        \
    MPZ_REALLOC (__dst, __size+1);                      \
    __dst_p = PTR(__dst);                               \
                                                        \
    __c = mpn_mul_1 (__dst_p, PTR(__src), __size, n);   \
    __dst_p[__size] = __c;                              \
    SIZ(__dst) = __size + (__c != 0);                   \
                                                        \
  } while (0)


void
mpz_fac_ui (mpz_ptr result, unsigned long int n)
{
#if SIMPLE_FAC
  /* Be silly.  Just multiply the numbers in ascending order.  O(n**2).  */
  unsigned long int k;
  mpz_set_ui (result, 1L);
  for (k = 2; k <= n; k++)
    mpz_mul_ui (result, result, k);
#else

  /* Be smarter.  Multiply groups of numbers in ascending order until the
     product doesn't fit in a limb.  Multiply these partial product in a
     balanced binary tree fashion, to make the operand have as equal sizes
     as possible.  When the operands have about the same size, mpn_mul
     becomes faster.  */

  unsigned long  k;
  mp_limb_t      p, p1, p0;

  /* Stack of partial products, used to make the computation balanced
     (i.e. make the sizes of the multiplication operands equal).  The
     topmost position of MP_STACK will contain a one-limb partial product,
     the second topmost will contain a two-limb partial product, and so
     on.  MP_STACK[0] will contain a partial product with 2**t limbs.
     To compute n! MP_STACK needs to be less than
     log(n)**2/log(BITS_PER_MP_LIMB), so 30 is surely enough.  */
#define MP_STACK_SIZE 30
  mpz_t mp_stack[MP_STACK_SIZE];

  /* TOP is an index into MP_STACK, giving the topmost element.
     TOP_LIMIT_SO_FAR is the largets value it has taken so far.  */
  int top, top_limit_so_far;

  /* Count of the total number of limbs put on MP_STACK so far.  This
     variable plays an essential role in making the compututation balanced.
     See below.  */
  unsigned int tree_cnt;

  /* for testing with small limbs */
  if (MP_LIMB_T_MAX < ULONG_MAX)
    ASSERT_ALWAYS (n <= MP_LIMB_T_MAX);

  top = top_limit_so_far = -1;
  tree_cnt = 0;
  p = 1;
  for (k = 2; k <= n; k++)
    {
      /* Multiply the partial product in P with K.  */
      umul_ppmm (p1, p0, p, (mp_limb_t) k);

#if GMP_NAIL_BITS == 0
#define OVERFLOW (p1 != 0)
#else
#define OVERFLOW ((p1 | (p0 >> GMP_NUMB_BITS)) != 0)
#endif
      /* Did we get overflow into the high limb, i.e. is the partial
	 product now more than one limb?  */
      if OVERFLOW
	{
	  tree_cnt++;

	  if (tree_cnt % 2 == 0)
	    {
	      mp_size_t i;

	      /* TREE_CNT is even (i.e. we have generated an even number of
		 one-limb partial products), which means that we have a
		 single-limb product on the top of MP_STACK.  */

	      MPZ_MUL_1_POS (mp_stack[top], mp_stack[top], p);

	      /* If TREE_CNT is divisable by 4, 8,..., we have two
		 similar-sized partial products with 2, 4,... limbs at
		 the topmost two positions of MP_STACK.  Multiply them
		 to form a new partial product with 4, 8,... limbs.  */
	      for (i = 4; (tree_cnt & (i - 1)) == 0; i <<= 1)
		{
		  mpz_mul (mp_stack[top - 1],
			   mp_stack[top], mp_stack[top - 1]);
		  top--;
		}
	    }
	  else
	    {
	      /* Put the single-limb partial product in P on the stack.
		 (The next time we get a single-limb product, we will
		 multiply the two together.)  */
	      top++;
	      if (top > top_limit_so_far)
		{
		  if (top > MP_STACK_SIZE)
		    abort();
		  /* The stack is now bigger than ever, initialize the top
		     element.  */
		  MPZ_INIT_SET_1_NZ (mp_stack[top], p);
		  top_limit_so_far++;
		}
	      else
		MPZ_SET_1_NZ (mp_stack[top], p);
	    }

	  /* We ignored the last result from umul_ppmm.  Put K in P as the
	     first component of the next single-limb partial product.  */
	  p = k;
	}
      else
	/* We didn't get overflow in umul_ppmm.  Put p0 in P and try
	   with one more value of K.  */
	p = p0;
    }

  /* We have partial products in mp_stack[0..top], in descending order.
     We also have a small partial product in p.
     Their product is the final result.  */
  if (top < 0)
    MPZ_SET_1_NZ (result, p);
  else
    MPZ_MUL_1_POS (result, mp_stack[top--], p);
  while (top >= 0)
    mpz_mul (result, result, mp_stack[top--]);

  /* Free the storage allocated for MP_STACK.  */
  for (top = top_limit_so_far; top >= 0; top--)
    mpz_clear (mp_stack[top]);
#endif
}


Top
回复人: shines(郭子) ( ) 信誉:100 2004-7-16 12:44:09 得分:0
 

汇报一下,目前我着眼于FFT的研究,FFT的代码优化是个严重的问题

Top
回复人: chenzhengzhanglu(陈正) ( ) 信誉:100 2004-7-16 13:47:42 得分:0
 

佩服大家!

有点同意   wenminghu(看到人生的希望)  的话!

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 14:34:30 得分:0
 

哪个FFT?

浮点的不好处理,俺原来做过一次,纳闷问什么变换一次,再逆过来,有比较大的误差


Top
回复人: wsmagic(飞翔的心) ( ) 信誉:100 2004-7-16 14:40:51 得分:0
 

to yaos:
  就是快速傅立叶变换啊。

  我用printf("%d", data[i]);这样的方式,在Windows下输出10000000!的结果大致是50多秒,赛扬900,应该有理由相信Linux下屏幕输出会快一点,如果使用STL下的 cout << 输出那个才叫慢呢,虽然屏幕好像没那么花眼,但是也柔和得太慢了:)

to chenzhengzhanglu(陈正):
  我想大家都在努力啊

我在这个贴子里也只能以这个ID出现了,因为shines回复已经超过30次了(汗,这个蝈蝈俊,我想骂死他的)

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 15:23:49 得分:0
 

哦,我曾经做过FPU代码优化,没成功
不知SSE2能行不,嘿嘿,用汇编速度至少提高3-10倍

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-7-16 15:30:11 得分:0
 

这正是我的新开发的FNT偏慢的瓶颈所在。
尽管我采用了更合理的数据结构,和更优化的流程。
看来,又有一块硬骨头需要啃下去了。。。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 15:46:51 得分:0
 

SSE2专门有乘加指令计算
a*b+c*d

Top
回复人: liangbch(宝宝) ( ) 信誉:105 2004-7-16 15:51:35 得分:0
 

to gxq:
   影响 FNT 算法最重要是的模乘法(a * b % p)的速度如何, 你的模乘法是如何实现的,是使用double型求商再使用乘法-减法余数, 还是使用unsigned __int64 型直接求余数,还是使用整数汇编指令求余,还是使用整形汇编指令以乘代除法求余数。
  我觉得你可以在这个方面作一些尝试,看看那个更快。



Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-7-16 16:00:10 得分:0
 

我用的是“unsigned __int64 型直接求余数”,
可否告诉我后两者如何实现,因为你的汇编实在太强了。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 16:08:40 得分:0
 

__asm
  {
  mov eax, a
  mov edx, b
  mul edx
  div p
  }

输出 eax商 ed余数

Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-7-16 16:21:23 得分:0
 

哦,明白了,在这个特定的前提下(a、b 均小于 p),
a * b % p 最多只需调一次 div 指令,不必判断高位是否要做 div 的问题。
多谢提醒。

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 17:42:59 得分:0
 

介于karatsuba和FFT中间的乘法:
T(n) = O(n*2*sqrt(2*logn)*logn)

Top
回复人: yaos(等待WoW前玩私服吧·无心人) ( ) 信誉:100 2004-7-16 18:37:39 得分:0
 

介于karatsuba和FFT中间的乘法:
T(n) = O(n*2^sqrt(2*logn)*logn)


Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-07-21 12:26:00 得分:0
 
喜讯:HugeCalc Ver 3.0.0.1 (beta) 发布了!
=========================================

它采用了“快速数论算法”,速度得以再次提升!
里面的参数参考了apfloat,同时改写和优化了许多核心代码。

计算一百万的阶乘,在 P4 1.7GHz / 256MB RAM 的机器上需要 13.500s,
速度已超过 apfloat!

有兴趣者可下载:http://www.freewebs.com/maths/download/HugeCalc.rar 
正式版将于 2004-08-08 正式发布(现仅需完善 help 文档了)。

//=================================================================

据说,论坛有回复 30 条限制(对楼主居然也不网开一面?),
所以,如果我一旦发言受限,将会立即结帖,并新开贴讨论。

原本打算来者有分,可惜我的级别太低,只能提供 100 分的给分权限,
到时,我将尽可能地公平给分,如果有不满意的,还请多多包涵!

(在本贴申明,主要是节省一次宝贵的发言机会;同时也消除一些不必要/无意义的误会)
Top
回复人: gxqcn(GxQ) ( ) 信誉:100 2004-07-21 16:02:00 得分:0
 
test report:

P4 CPU 1.70GHz / 256MRAM / WinXP
================================

   10,000!   0.031s     0.28462... x 10^35,660
   20,000!   0.109s     0.18192... x 10^77,338
   40,000!   0.390s     0.20916... x 10^166,714
   80,000!   1.328s     0.30977... x 10^357,507
  100,000!   1.969s     0.28242... x 10^456,574
  200,000!   6.438s     0.14202... x 10^973,351
  400,000!  20.828s     0.25344... x 10^2,067,110
  800,000!  66.921s     0.56846... x 10^4,375,040
--------2004.05.28


HugeCalc 3.0.0.1 (beta)
   10,000!   0.031s     0.28462... x 10^35,660
   20,000!   0.078s     0.18192... x 10^77,338
   40,000!   0.219s     0.20916... x 10^166,714
   80,000!   0.562s     0.30977... x 10^357,507
  100,000!   0.671s     0.28242... x 10^456,574
  200,000!   1.640s     0.14202... x 10^973,351
  400,000!   4.015s     0.25344... x 10^2,067,110
  800,000!   8.859s     0.56846... x 10^4,375,040
1,000,000!  13.469s     0.82639... x 10^5,565,709
--------2004.07.21


apfloat 2.40 (shines(郭子) 提供的 FactorialApfloat.exe,2004-06-21)
   10,000!   0.039s
   20,000!   0.089s
   40,000!   0.225s
   80,000!   0.721s
  100,000!   0.830s
  200,000!   2.118s
  400,000!   5.134s
  800,000!  11.411s
1,000,000!  17.120s


Mathematica 4.0.0.0
   10,000!   0.031s
   20,000!   0.141s
   40,000!   0.406s
   80,000!   1.094s
  100,000!   1.422s
  200,000!   4.046s
  400,000!  10.016s
  800,000!  24.719s
1,000,000!  44.219s
Top