擂台：超大整数高精度快速算法

gxqcn 2004-05-10 02:08:44
我用 C++ 开发了一套《超大整数高精度快速算法》，前后花了十余年。

实现了如下功能：

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

该套算法效率非常高，举例如下：

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

（我的算法纯粹是数学上的优化，未内嵌一行汇编代码）

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

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

注：
1、声明, 我的程序快只能说明我对此算法下了功夫，并不表明水平有多高，事实上，我的水平很有限（甚至还不会汇编）；
2、并不要求所有功能均须强于我现有的软件，有一两项也请及时告诉我，谢谢！
3、算法或软件不限于自身原创，欢迎提供第三方链接信息，谢谢！
...全文
4128 158 打赏 收藏 转发到动态 举报
写回复
用AI写文章
158 条回复
切换为时间正序
请发表友善的回复…
发表回复
不羡仙z 2021-07-20
  • 打赏
  • 举报
 回复

100000!我只用了50ms啊

zhuzebin 2004-08-19
  • 打赏
  • 举报
 回复
Mark!
gxqcn 2004-07-21
  • 打赏
  • 举报
 回复
喜讯：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 分的给分权限，
到时，我将尽可能地公平给分，如果有不满意的，还请多多包涵！

（在本贴申明，主要是节省一次宝贵的发言机会；同时也消除一些不必要/无意义的误会）
wsmagic 2004-07-21
  • 打赏
  • 举报
 回复
本贴有续集，欢迎继续

请大家转至：
http://community.csdn.net/Expert/topic/3197/3197332.xml?temp=.4934198
gxqcn 2004-07-21
  • 打赏
  • 举报
 回复
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
shines77 2004-07-16
  • 打赏
  • 举报
 回复
汇报一下，目前我着眼于FFT的研究，FFT的代码优化是个严重的问题
yaos 2004-07-16
  • 打赏
  • 举报
 回复
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
}
yaos 2004-07-16
  • 打赏
  • 举报
 回复
有点不公平啊，GMP输出到屏幕，你们输出到文件

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

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

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

嘿嘿
shines77 2004-07-16
  • 打赏
  • 举报
 回复
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
yaos 2004-07-16
  • 打赏
  • 举报
 回复
介于karatsuba和FFT中间的乘法：
T(n) = O(n*2^sqrt(2*logn)*logn)
yaos 2004-07-16
  • 打赏
  • 举报
 回复
介于karatsuba和FFT中间的乘法：
T(n) = O(n*2*sqrt(2*logn)*logn)
gxqcn 2004-07-16
  • 打赏
  • 举报
 回复
哦，明白了，在这个特定的前提下（a、b 均小于 p），
a * b % p 最多只需调一次 div 指令，不必判断高位是否要做 div 的问题。
多谢提醒。
yaos 2004-07-16
  • 打赏
  • 举报
 回复
__asm
{
mov eax, a
mov edx, b
mul edx
div p
}

输出 eax商 ed余数
gxqcn 2004-07-16
  • 打赏
  • 举报
 回复
我用的是“unsigned __int64 型直接求余数”，
可否告诉我后两者如何实现，因为你的汇编实在太强了。
liangbch 2004-07-16
  • 打赏
  • 举报
 回复
to gxq:
影响 FNT 算法最重要是的模乘法(a * b % p)的速度如何， 你的模乘法是如何实现的，是使用double型求商再使用乘法-减法余数， 还是使用unsigned __int64 型直接求余数，还是使用整数汇编指令求余，还是使用整形汇编指令以乘代除法求余数。
我觉得你可以在这个方面作一些尝试，看看那个更快。

yaos 2004-07-16
  • 打赏
  • 举报
 回复
SSE2专门有乘加指令计算
a*b+c*d
gxqcn 2004-07-16
  • 打赏
  • 举报
 回复
这正是我的新开发的FNT偏慢的瓶颈所在。
尽管我采用了更合理的数据结构，和更优化的流程。
看来，又有一块硬骨头需要啃下去了。。。
yaos 2004-07-16
  • 打赏
  • 举报
 回复
哦，我曾经做过FPU代码优化，没成功
不知SSE2能行不，嘿嘿，用汇编速度至少提高３－１０倍
wsmagic 2004-07-16
  • 打赏
  • 举报
 回复
to yaos:
就是快速傅立叶变换啊。

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

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

我在这个贴子里也只能以这个ID出现了，因为shines回复已经超过30次了（汗，这个蝈蝈俊，我想骂死他的）
yaos 2004-07-16
  • 打赏
  • 举报
 回复
哪个FFT？

浮点的不好处理，俺原来做过一次，纳闷问什么变换一次，再逆过来，有比较大的误差
加载更多回复（138）
汇编优化
csdn上果然藏龙卧虎呀，今天看到一系列汇编优化的讨论贴，认真看看对指令优化、SSE2的理解会很有帮助，把链接附上： 千分求汇编优化：UInt96x96To192(...)http://community.csdn.net/Expert/TopicView3.asp?id=5505130擂台：超大整数高精度快速算法http://topic.csdn.net/T/20040510/14/30497...
洛谷专题训练 ——【算法1-1】模拟与高精度
洛谷题单【算法1-1】模拟与高精度 ACM-ICPC在线模板 题单链接： 【算法1-1】模拟与高精度 下面的这一坨都是洛谷题单上的东东 题单简介 恭喜大家完成了第一部分语言入门，相信大家已经可以使用 C++ 写出一些简单程序了。 各位读者有听说过“建模”一词吗？所谓“建模”，就是把事物进行抽象，根据实际问题来建立对应的数学模型。“抽象”并不意味着晦涩难懂；相反，它提供了大量的便利。计算机很难直接去...
【算法1-1】洛谷模拟与高精度题目与解析
【算法1-1】洛谷 模拟与高精度 题目与解析 前言（引用自洛谷） 各位读者有听说过“建模”一词吗？所谓“建模”，就是把事物进行抽象，根据实际问题来建立对应的数学模型。“抽象”并不意味着晦涩难懂；相反，它提供了大量的便利。计算机很难直接去解决实际问题，但是如果把实际问题建模成数学问题，就会大大地方便计算机来“理解”和“解决”。 举个生活中常见的例子：我们拿到了某次数学考试的成绩单，现在需要知道谁考得最好。当然不能把成绩单对着电脑晃一晃，然后问“谁考得最好？”。需要通过一种途径让计算机来理解这个问题。这个问题可
阶乘优化算法
阶乘算法优化 2016年11月28日 10:17:37 阅读数：3574 原文地址：http://blog.csdn.net/yxnk/article/details/1665052 我的感言：首先，有一个概念上的认知，即根据阶乘定义而来的常规算法，如果是long int型只能正确计算到12左右的阶乘，如果用double型只能正确计算170左右的阶乘，当然这些只是大概，需要结合实际平台进行...
局部排序之擂台算法
n个数里找出前m个数（m<<n） 引子 每年十一月各大IT公司都不约而同、争后恐后地到各大高校进行全国巡回招聘。与此同时，网上也开始出现大量笔试面试题；网上流传的题目往往都很精巧，既能让考查基础知识，又在平淡中隐含了广阔的天地供优秀学生驰骋。 从10亿个浮点数中找出最大的1万个。 这是一道似易实难的题目，一般同学最容易中的陷阱就是没有重视这个“亿”字。因...
数据结构与算法

33,008

社区成员

35,326

社区内容

发帖
与我相关
我的任务
社区描述
数据结构与算法相关内容讨论专区
社区管理员
  • 数据结构与算法社区
加入社区
  • 近7日
  • 近30日
  • 至今

加载中

查看更多榜单
社区公告
暂无公告

试试用AI创作助手写篇文章吧

+ 用AI写文章