迭代法求n次方根：欲与POWER()比高低（另附：小数转分数、求最大公因数自定义函数）

aoe1981 · 发表于 2020-2-3 20:14

本帖最后由 aoe1981 于 2020-2-4 09:04 编辑

本帖题目虽长，但主旨只有一个：

1.利用迭代法求n次方根

至于：
2.小数转分数
3.求最大公因数
顶多算是主旨1探究过程中的“副产品”。

迭代求根主要基于以下迭代函数：
f(x)=x*(n-1)/n+C/(n*x^(n-1))
其中：
x^n=C
或：
C^(1/n)=x
即：
从乘方运算的角度：C是幂，n是指数，x是底数；
从开方运算的角度：C是被开方数，n是次数，x是方根；
所以混为一谈，是因为：乘方运算和开方运算其实是一回事。
乘方运算：y=x^n
开方运算：y=x^(1/n)
二者都可以统一到“幂函数”，乘方运算的定义域根据n判断，开方运算的定义域根据1/n判断，这一点至关重要，后面还有涉及。
迭代的过程很简单：
①猜测一个初值，比如：x0=1；
②x1=f(x0)；
③x0=x1；
④反复循环②③步，直至abs(x1-x0)小于事先指定的精度退出。

有了迭代函数，求n次方根表面上看很简单，初级版（问题多多）的自定义函数如下：

Option Explicit
'初级欠完善版备份（x0=1）
Public Function FangGen1(C As Double, Optional n As Double = 2, Optional e As Double = 0.000001, Optional pd As Boolean = False)
If C < 0 And (n And 1) = 0 Then FangGen = "负数无偶数次方根": Exit Function
If C = 1 And n = 0 Then FangGen = "1的0次方根有无穷多个": Exit Function
If C <> 1 And n = 0 Then FangGen = "非1的0次方根不存在": Exit Function
Dim x0 As Double, x1 As Double
x1 = (n - 1) / n + C / n
Do
x0 = x1
x1 = x0 * (n - 1) / n + C / (n * x0 ^ (n - 1))
Loop Until Abs(x1 - x0) < e
FangGen = x1
End Function

复制代码

只要稍微对比使用下工作表幂函数：=POWER()，就会发现：上面的自定义函数简直弱爆了，就是闹着玩的，多此一举。

深入研究发现：
1.幂函数的指数n或1/n小于0时，底数C≠0；
2.幂函数的指数n或1/n表成“最简分数”后：分母为奇数时（包含整数），底数C可以取遍实数集R；分母为偶数时，底数C只能取正实数集R+和0。

POWER()函数：
（一）可以计算形如：
=(-2)^(0.2)或=POWER(-2,0.2)
=(-2)^(-0.2)或=POWER(-2,-0.2)

（二）不能直接计算：
=(-2)^(0.8)或=POWER(-2,0.8)
=(-2)^(-0.8)或=POWER(-2,-0.8)

但是可以套用POWER()函数间接计算：
=((-2)^4)^0.2或=POWER(POWER(-2,4),0.2)
=((-2)^-4)^0.2或=POWER(POWER(-2,-4),0.2)
=((-2)^4)^-0.2或=POWER(POWER(-2,4),-0.2)

（三）不能直接计算：
=(-2)^(0.6)或=POWER(-2,0.6)
=(-2)^(-0.6)或=POWER(-2,-0.6)
同上，可以类似套用POWER()函数。

（备注：我在Excel 2013中测试）
大抵是：0.2=1/5，0.8=4/5，0.6=3/5，所谓“最简分数”，POWER()函数只认分子是1的。

上面的初级版的自定义开方函数对于底数C为负时，完全不能计算，因此我又设计了下面的加强版（很费劲、费劲、费劲！！！）：

Option Explicit
'迭代法求n次方根：欲与POWER()比高低
Public Function FangGen(C As Double, Optional n As Double = 2, Optional e As Double = 0.000001)
If C = 1 And n = 0 Then FangGen = "1的0次方根有无穷多个": Exit Function
If C <> 1 And n = 0 Then FangGen = "非1的0次方根不存在": Exit Function
If C = 0 And n < 0 Then FangGen = "0不能开负数次方": Exit Function
Dim pd As Boolean, n1 As Double, i As Long, gys As Double, fz As Double, fm As Double
If C < 0 Then
n1 = Abs(n)
Do Until Int(n1) = n1
n1 = n1 * 10
i = i + 1
Loop
gys = GCD_vba(n1, 10 ^ i)
fz = 10 ^ i / gys
fm = n1 / gys
If (fm And 1) = 0 Then '偶分母
FangGen = "错误": Exit Function
Else '奇分母
C = C * -1
If fm = 1 Then '分母为1（整数）
If (fz And 1) = 1 Then pd = True '奇整数
Else '非1奇分母
pd = True
End If
End If
End If
'******以下大致推断一个初值x0******
Dim C1 As Double, mws As Long, dws As Long, dws0 As Double, ds As Long
n1 = Abs(n)
C1 = Int(Abs(C))
ds = InStr(1, C1, "E+")
If ds > 0 Then
mws = Mid(C1, ds + 2) + 1
Else
mws = Len(C1 & "")
End If
For i = 0 To Int(n1) - 1
dws0 = (mws + i) / Int(n1)
If Int(dws0) = dws0 Then dws = dws0: Exit For
Next i
'******以上大致推断一个初值x0******
Dim x0 As Double, x1 As Double
x0 = 0.55 * 10 ^ dws
x1 = x0 * (n1 - 1) / n1 + C / (n1 * x0 ^ (n1 - 1))
Do
x0 = x1
x1 = x0 * (n1 - 1) / n1 + C / (n1 * x0 ^ (n1 - 1))
Loop Until Abs(x1 - x0) < e
If n < 0 Then FangGen = 1 / x1 Else FangGen = x1
If pd Then FangGen = -1 * FangGen
End Function
Public Function GCD_vba(x As Double, y As Double) As Double '最大公因数
Dim bcs#, cs#, q#, r#
If x < y Then bcs = y: cs = x Else bcs = x: cs = y
If cs = 0 Then GCD_vba = bcs: Exit Function
Do
q = Int(bcs / cs)
r = bcs - cs * q
If r = 0 Then
GCD_vba = cs
Exit Do
Else
bcs = cs
cs = r
End If
Loop
End Function

复制代码

这段代码繁琐了许多，但有许多收获：
我的自定义求n次方根函数：FangGen(被开方数C, 次数n, 精度e)，可以直接计算上面的类型（一）、（二）:

对比1.jpg

但是依旧不能计算类型（三），不过，这不是我的代码和算法的原因，因为vba本身不支持：

不能计算.jpg

奇怪的是，系统自带的计算器却可以计算类型（三）：

能计算1.jpg

至此，是不是能说：我的自定义函数FangGen()比POWER()强大呢？

对于上面所列，的确如此。

对比2.jpg

但是，当面对特别大（小）的底数、指数时，POWER()只要数据不超出最大范围且不违反数学法则，总能计算出结果，我的FangGen()却不一定。主要有四种组合：
大底数、大指数，各取正负，共4种；
大底数、小指数，各取正负，共4种；
小底数、大指数，各取正负，共4种；
小底数、小指数，各取正负，共4种；
以上总计16种情况，我测试了几种，已经发现了这一问题，没有完全测试。

正常的：

正常.jpg

出错的：

对此，我是做过改进的。不是算法本身的问题，关键是：迭代函数的初值x0并不好猜！！！

一个好的初值x0，不仅能减少迭代次数，而且还能避免迭代计算过程中产生“数据溢出”的错误。上图FangGen()不能计算的情况，只要能估计出一个比较接近真实值的初值，是可以计算出结果的。但这又好像是“正确的废话”！不过，我的确做了一些有价值的尝试。

一个整数部分为mws（幂位数）位的C，开n或1/n次方，其方根的位数dws（底位数）是可以计算求得的，这算是我上面代码的核心独创之处。相关部分如下：

'******以下大致推断一个初值x0******
Dim C1 As Double, mws As Long, dws As Long, dws0 As Double, ds As Long
n1 = Abs(n)
C1 = Int(Abs(C))
ds = InStr(1, C1, "E+")
If ds > 0 Then
mws = Mid(C1, ds + 2) + 1
Else
mws = Len(C1 & "")
End If
For i = 0 To Int(n1) - 1
dws0 = (mws + i) / Int(n1)
If Int(dws0) = dws0 Then dws = dws0: Exit For
Next i
'******以上大致推断一个初值x0******

复制代码

原理后述。
但，光有这一步，还远远不够，我只能取平均数，如图：

平均数.jpg

为了简便，我在代码中这样处理：

x0 = 0.55 * 10 ^ dws

复制代码

如果vba中能根据“数据溢出错误”调整代码走向，至少可以遍历：0.05、0.15、0.25、0.35、0.45、0.55、0.65、0.75、0.85、0.95，说不定会极大可能的试出合适的初值，但似乎这在vba中是无法实现的。

POWER()函数的强大之处在于：或许使用了更为先进的算法。

当然，最强大的还是系统自带的计算器。

哈哈。

附件如下：

迭代方根.zip (37.7 KB, 下载次数: 7)

（附件已更新至6楼补充，本楼文字介绍部分中的偏差不再纠正，可阅读6楼）

aoe1981 · 发表于 2020-2-3 20:21

下面是小数转分数自定义函数及效果图：

Public Function GCD_vba(x As Double, y As Double) As Double '最大公因数
Dim bcs#, cs#, q#, r#
If x < y Then bcs = y: cs = x Else bcs = x: cs = y
If cs = 0 Then GCD_vba = bcs: Exit Function
Do
q = Int(bcs / cs)
r = bcs - cs * q
If r = 0 Then
GCD_vba = cs
Exit Do
Else
bcs = cs
cs = r
End If
Loop
End Function
Public Function XsToFs(a As Double, Optional n As Integer = 1) As String '小数转分数
Dim zs As Long, fz As Double, fm As Double, gys As Double, pd As Boolean, i As Long
If a < 0 Then a = a * -1: pd = True
Do Until InStr(1, a, ".") = 0
a = a * 10
i = i + 1
Loop
a = Int(a)
gys = GCD_vba(a, 10 ^ i)
fm = 10 ^ i / gys
fz = a / gys
zs = fz \ fm
If n <> 1 Then '假分数
XsToFs = "0 " & fz & "/" & Abs(fm) '保留整数部分0是为了在表格中能够正常引用计算
Else '真分数
fz = fz - zs * fm
If fz = 0 Then XsToFs = zs & " " Else XsToFs = zs & " " & fz & "/" & Abs(fm)
End If
If pd Then XsToFs = "-" & XsToFs
End Function

复制代码

效果图中：对于真分数、假分数，保留整数部分的0，输出为：0 3/5，是为了保证输出结果被引用时能够正确计算，如果输出为： 3/5，则会被当成日期值进行错误的计算。

aoe1981 · 发表于 2020-2-3 20:31

本帖最后由 aoe1981 于 2020-2-3 21:08 编辑

以下是求最大公因数的自定义函数：

Public Function GCD_vba(x As Double, y As Double) As Double '最大公因数
Dim bcs#, cs#, q#, r#
If x < y Then bcs = y: cs = x Else bcs = x: cs = y
If cs = 0 Then GCD_vba = bcs: Exit Function
Do
q = Int(bcs / cs)
r = bcs - cs * q
If r = 0 Then
GCD_vba = cs
Exit Do
Else
bcs = cs
cs = r
End If
Loop
End Function

复制代码

基本算理是：辗转相除法。

请注意：当有负数时，工作表函数GCD()会出错，我的自定义函数GCD_vba会返回结果；当有小数时，两者返回结果存在差别；工作表函数GCD支持不定项参数，即可以一次性求两个数、三个数、四个数、……的最大公因数，我的自定义函数只能一次处理两个数的最大公因数，更多的需要套用函数使用，如图：

最大公因数.jpg

上图已更新到1楼附件中。

aoe1981 · 发表于 2020-2-3 21:23

顶楼举例补充说明：
工作表中支持POWER()函数套用的情况：

支持.jpg

在vba中不支持：

二者有差异啊。

tongkj · 发表于 2020-2-3 22:23

好功能

aoe1981 · 发表于 2020-2-4 08:57

本帖最后由 aoe1981 于 2020-2-4 09:03 编辑

发现我的自定义函数可以计算顶楼谈到的类型（三）：
=(-2)^(0.6)
=(-2)^(-0.6)

当出现#VALUE!的错误时，往往是可以计算的，情况大约有以下几点（值得讨论）：

1.如何判断整数？

我常用：Int(n1) = n1
但是似乎：InStr(1, n1, ".") = 0更靠谱，int在处理大数或循环小数时，会多判断显示位数。

这种情况似乎解决了。

2.如何判断偶数？

我学了：
If (fm And 1) = 0 Then '偶分母
If (fm Mod 2) = 0 Then '偶分母

但上二者面对：fm=1.06666666666667E+15这样的大数时，会出现“溢出错误”，我改用：
i = InStr(1, fm, "E+")
If i = 0 Then ws = Right(fm, 1) Else ws = Mid(fm, i - 1, 1)
If ws Mod 2 = 0 Then '偶分母或：
If (ws And 1) = 0 Then '偶分母

似乎解决了……但是如果是：fm=1.06666666666667E+35，取到的7其实并不是个位，这样处理是有问题的，木有办法呀……

3.还是初值x0估计偏差较大造成的。

经改造，类型（三）可以计算了，如图：

可以计算了.jpg

顶楼附件同步更新，介绍内容中的偏差不再纠正。

jurdy001 · 发表于 2020-2-9 10:36

aoe1981 发表于 2020-2-3 20:21
下面是小数转分数自定义函数及效果图：

受用厉害老师太强了

		自动登录	找回密码
密码			免费注册

[讨论] 迭代法求n次方根：欲与POWER()比高低（另附：小数转分数、求最大公因数自定义函数）

评分

最新热点 /1