数学研发论坛

 找回密码
 欢迎注册
查看: 409|回复: 5

[讨论] 多项式除法 算法可应用于多精度大数计算吗?

[复制链接]
发表于 2018-4-24 16:16:09 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
闲逛网络,偶尔发现篇“多项式除法及求模”,用到了“多项式求逆元”及系数反转方法,非常巧妙!

现在的问题是:多精度大数计算中的除法或模运算,可以借鉴上述算法吗?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-4-24 16:45:27 | 显示全部楼层
这个感觉和Montgomery算法的思想很类似呀
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-4-28 14:13:22 | 显示全部楼层
它与 Montgomery 算法思想确实有点类似。

可惜,我揣摩了几天,似乎不适合应用于大数计算中。
因为有理代数运算,系数允许负数、分数;而大数运算得到的结果必须为小于进制的非负整数。
两者的鸿沟还是蛮大的。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-4-28 14:58:30 | 显示全部楼层
假设我们要计算 \(123456 \div 789 =?\)

如果将它们看作10进制的代数式,则:
  1. PolynomialQuotient[x^5 + 2 x^4 + 3 x^3 + 4 x^2 + 5 x + 6,  7 x^2 + 8 x + 9, x]
复制代码

结果为:\[\frac{706}{2401} + \frac{36 x}{343} + \frac{6 x^2}{49} + \frac{x^3}{7}\]

采用先求多项式逆元的算法,流程如下:
  1. AR = 1 + 2 x + 3 x^2 + 4 x^3 + 5 x^4 + 6 x^5;
  2. BR = 7 + 8 x + 9 x^2;
  3. 1/7
  4. PolynomialMod[% (2 - % BR), x^2]
  5. PolynomialMod[% (2 - % BR), x^4]
  6. PolynomialMod[% AR, x^4]
复制代码

结果为:
\begin{align*}
&\frac{1}{7} \\
&\frac{1}{7} - \frac{8 x}{49} \\
&\frac{1}{7} - \frac{8 x}{49}+ \frac{x^2}{343} + \frac{496 x^3}{2401} \\
&\frac{1}{7} + \frac{6 x}{49}+ \frac{36 x^2}{343} + \frac{706 x^3}{2401} \\
\end{align*}
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-4-28 21:17:15 | 显示全部楼层
注意到分母7、49、343、2401都是\(7^n\);除式\(7 x^2 + 8 x + 9\)的最高幂项系数

验证确实只与最高幂项系数有关
  1. PolynomialQuotient[4 x^5 + 5 x^4 + 6 x^3 + 7 x^2 + 8 x + 9,  a x^2 + 2 x + 3, x]
复制代码


那么如何进一步改良,,,,,,

-----------

设想1,调整除数进制使最高幂项系数为1

-----------

设想2,除数被除数均乘上一个整数k,使得除数最高幂项系数为1
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-4-28 23:38:38 | 显示全部楼层
单次试验
  1. {200,100}//RandomInteger[{10^#,10^(#+1)-1}]&/@#&//(Print[#];#)&//Ceiling[10^(IntegerLength[#[[2]]]+1)/#[[2]]]#&//(Print[#];#)&//FromDigits[IntegerDigits[#],x]&/@#&//{#[[1]]/#[[2]],PolynomialQuotient[#[[1]],#[[2]],x]}&//#/.x->10&//{#[[1]]//Round,#[[2]],#[[2]]-#[[1]]//Round//Abs//{#,IntegerLength[#]}&}&
复制代码


多次试验
  1. Table[{200,100}//RandomInteger[{10^#,10^(#+1)-1}]&/@#&//Ceiling[10^(IntegerLength[#[[2]]]+1)/#[[2]]]#&//FromDigits[IntegerDigits[#],x]&/@#&//{#[[1]]/#[[2]],PolynomialQuotient[#[[1]],#[[2]],x]}&//#/.x->10&//#[[2]]-#[[1]]&//Round//Abs//IntegerLength,{100}]
复制代码

误差小于准确值长度一半,需要一次牛顿法修正,,,,
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2018-10-18 09:06 , Processed in 0.064949 second(s), 16 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表