前面的估算因子是近似值，η=（2B2+1）/（2B1+1）或者直接写为B2/B1，这里B2=2B2+1，若p为预估因子，则B2就是真因子的比值，若知道了B2的精确值则真因子就可以丝毫不差的求出，但我们预先是不可能知道，无法知道小数点后几十几百位，只能试验B的值的整数部分，我的程序设定B=3～2001，所以前面的预估因子和F的值只能是接近实际，而且还差很多，由于是从最高位向低位试验所以已经很节约程序了，但B的值要从3到2001没有快速乘法除法程序很费力，前面的只用了几秒钟分解的整数也是预先知道比值约为1.9，所以从3开始很快就得到结果。若不知道的话就从2001开始向下算到3那就费力了，我还分解了一个61位的整数那个也是预先知道两个因子的比值的整数部分为1042，从1043开始算的约几分钟才出来结果。
所以没有大整数的快速乘法除法程序是不行的。
前面估计的因子与实际的差距好大，如何更进一步接近实际呢？以至达到仅末尾5位以内与实际不同？其实我们可以用如下的简单方法，是个经验公式：
设p<q，p已经接近实际，pq+r=M，则q^2-M=M1，设m1=√M1，n=q-m1，则2n+1接近实际或略大于实际，最多迭代2次就可以非常接近理想结果即只有末尾几位与实际不同，略小于或略多于实际。
其实原理容易证明的，我们都知道（n+1）^2-n2=2n+1，据此就可以推导证明出来的 ，我没有证明只是个经验。
需要注意的是，此法是向上接近实际的，若p是大于实际的，那么此法的结果就更大，越来越远离实际了，所以必须注意不能迭代多次无穷迭代下去就算不完了，电脑就死机了。
我的电脑打不开本网站，我是用手机登录的，所以我的分解大整数的程序无法上传，怎么弄？
欢迎感兴趣的朋友指教，欢迎探讨！
举个例子：14585113=997*14629，而14585113/983=14837.3479，14837^2=220136569，220136569-14585113=205551456，√205551456=14337，14837-14337=500，2*500+1=1001，1001-997=4，比983更接近实际，1001>997所以不能再迭代下去了，再算下去就更远了。
这个法也是复杂，更简单的是这个（也向上接近实际的有可能迭代2次就远远大于实际了）：
设q1=（p+q）/2，q2=（q1）^2-M，则p1=q1-√q2，即为接近实际的预测因子。
举例：（14629+997）/2=7813，7813^2-14585113=46457856，√46457856=6816，7813-6816=997，997*14629=14585113.
此法是常用的法，费马分解法是以此为基础发展演变出来的。
实际费马法就是选择的p和q的值是远远大于实际的，而且是p+1，2，……，这样实验的，所以因子是不会超过实际因子的，但pq不是随机选择的，是用改进算法得到一个系列的数据库，这些方法有多种，所以，某些有效的方法可以分解很大的整数，对于大因子的积用到这种方法是有效的，如果是小因子用常规法反而更快。所以，某些RSA密码也不是绝对安全的。
本帖最后由 ysr 于 2021-10-24 08:31 编辑
此图就是大衍求一术的过程，是秦九韶弄出来的，与矩阵变换求乘法的逆元的原理是本质上一样的，比外国早了500年。
这个原理我没有推导，应该正确，我验证了一下，也是需要分类输出的，比如7/4，结果是上排为2，1，下排为5，1，实际答案是5-2=3，3*7/4=5……1，还要注意的是当某数除以1的情况，也就是下排先变成1的时候，商为某数减1而不是某数除以1仍然得某数，比如3/1则商为2余数为1.
还有一种情况就是比如23/61，最后结果是，上排为3，1，下排为5，7，实际答案是3+5=8，8*23=184，184/61=3……1。
如何分类，还不清楚，有待研究。
看上去简单，过程不明显，不容易理解，原理可以参考矩阵变换法求乘法的逆元。
比如：61/23得到，上排38，1，下排5，1，实际答案是38+5=43，验证：43*61-1=2623-1=2622，而2622/23=114，成立。
用C会比较好，C有个大数运算库，能省略求证的时候对大数进行计算的复杂处理
以下代码是我在试图分解一个RSA的公钥
#include <nbtheory.h>
#include <stdio.h>
#include <windows.h>
#include <sstream>
#include <iostream>
#pragma comment (lib,"./cryptlib.lib")
using namespace CryptoPP;
using namespace std;
Integer mod_exp(Integer a,Integer b,Integer c)
{
Integer res, t;
res = 1 % c;
t = a % c;
while (b.NotZero())
{
if (b % 2==1)
{
res = res * t%c;
}
t = t * t%c;
b = b / 2;
}
return res;
}
Integer exgcd(Integer a, Integer b, Integer *x, Integer *y)
{
if (b == 0) {
(*x) = 1, (*y) = 0;
return a;
}
Integer r = exgcd(b, a%b, x, y);
Integer t = *y;
*y =(*x)- a / b * (*y);
*x = t;
return r;
}
Integer searchUtil(Integer target, Integer l, Integer r)
{
if (l + 1 == r)
return l;
else if (l + 1 < r)
{
Integer m = l + (r - l) / 2;
if (m*m == target)
return m;
else if (m*m < target)
return searchUtil(target, m, r);
else if (m*m > target)
return searchUtil(target, l, m);
}
}
Integer searchSqrtFloor(Integer target)
{
return searchUtil(target, 1, target);
}
unsigned int main()
{
Integer divv("141200934203348100398237203092056603810497665953866503938436185453188653157534841542476799951117838232625775007972565407157639296613520419503480411605352582268797113156766586228582255481135809096970885597928003710586870460226020496368765994413857405308613414568674463104615847526186185331402709548604450804313");
Integer divs("11876433034222735702746321410377879376272094689220724132957900598176502768763874640725383992846504809710949618722994122810563152826799311453851843016297619");
while (1)
{
if (divs % 10 == 7)
{
divs = divs -4;
}
else
{
divs = divs-2;
}
if (IsPrime(divs))
{
cout << dec << divs << endl;
if (divv%divs != 0)
{
continue;
}
else
{
break;
}
}
}
cout <<"find it:"<< dec << divs << endl;
}
本帖最后由 xyaoy 于 2021-10-25 16:59 编辑
技术员 发表于 2012-6-5 16:35
支持蜘蛛侠！请继续努力。开根大数算法有难度，如果你给我算法，我可以帮你试试。
求最近的整数平方根
Integer searchUtil(Integer target, Integer l, Integer r)
{
if (l + 1 == r)
return l;
else if (l + 1 < r)
{
Integer m = l + (r - l) / 2;
if (m*m == target)
return m;
else if (m*m < target)
return searchUtil(target, m, r);
else if (m*m > target)
return searchUtil(target, l, m);
}
}
Integer searchSqrtFloor(Integer target)
{
return searchUtil(target, 1, target);
}
xyaoy 发表于 2021-10-25 08:58
求最近的整数平方根
Integer searchUtil(Integer target, Integer l, Integer r)
{
谢谢沟通指导！我不会VC，用VB程序的确速度太慢。
你这个平方根是用的迭代法吗？据说是迭代法计算速度快。
欢迎指导！谢谢！