真正的高精度科学计算器

云梦

感觉是特殊问题特殊处理，不可能处理一般情况 mathe 发表于 2013-1-16 17:35

一般情况当然很好处理，也就不需要什么特殊的方法。针对不同的问题采取不同的方法才是灵活可取的。如果不灵活的话，那么解这个方程要求的精度太高了。

云梦

我的计算器精度虽然不高，但在绝大多数情况下均能保持这个精度，复杂和简单的运算正是数学函数以及表达式的组成元素。 例如：如果不用数学软件，计算不完全伽马函数及正则不完全伽马函数，在0-10^14内的实数你可否能保持高的精度及正确结果？

wayne

针对不同的问题采取不同的方法才是灵活可取的。 云梦 发表于 2013-1-16 18:56

大型软件如Mathematica，MATLAB都不敢说一条命令通吃所有问题， 你这个小小的计算器也不大可能点一下按钮，或者输入方程，敲一个命令，就给出图片中方程解的超高的精度。 所以我宁愿相信，你特地为这个方程写了一段代码，放在计算器里面， 然后说我的计算器可以完美的解决这种问题，而其他的数学软件不能。这有什么意义呢。 ======== 先不谈意义，就题论题，该问题总得有具体的算法实现吧， 不采用mathe的处理方式，我实在是想不出其他的方法， 能否共享一下此题的算法实现

云梦

不需要什么特殊代码，只是利用了该方程的特点,x=a+s*10^-r。这样就可以把整数和小数分离，只计算小数部分。83楼也就是利用了近似的特点得到了近似值。不过他的x值只是我程序中x的初始值罢了。 下面就需要用牛顿法进行迭代，在迭代时只计算小数部分，这样不就得到较精确的结果了吗?其实这只是一个小技巧。技巧编程可以解决很多不容易解决的问题。

云梦

把方程整理一下，变成以下格式，即可以求出方程的初始值。 把左式与右式的差进行求导，得到下面最后一个公式。

云梦

如果你能看懂我的VB程序，你就分析一下，吧。 NameF.Text = "F(x) =" Dim Ii, xChen0, Xp0 As Single Dim cv As String * 2000 xChen0 = ychen Xp0 = xP Call SJbz Select Case ychen Case 64 xP = 24 xChen = 80 Case 128 xP = 45 xChen = 160 Case 256 xP = 85 xChen = 300 End Select Xspd = 0 Call SJbz Wx = Ax An = SToo(0) Bn = SToo(1) Ax = An Bx = Bn hsfx = "-" Call Addition If Left(Ax.St, 1) <> 0 And An.bz = "" And Bn.bz = "" Then 'A -B <> 0 '初始值 x0 If Left(Xy.St, 1) = "1" Or Left(Xy.St, 1) = "3" Then 'C=0 If Ax.bz = "" Then '如果A>B A-B交换 An = SToo(1) Bn = SToo(0) End If Call Dj01 '计算初始值 x0 Ax = SToo(3) If Xspd = 0 Then For Ii = 0 To 100 Call Dj2012 'x 带入原方程 'SToo(4) = Ax Call Dj2012F '导数 SToo(5) = Ax Cx = SToo(4) Ax = SToo(5) Call multx_ 'F(x)/F'(x) Bx = Cx Ax = SToo(3) hsfx = "-" Call Addition If Left(Bx.St, 1) = "0" Or Bx.zs < -xChen - 20 Then SToo(3) = Ax Call Dj2012 SToo(4) = Ax Exit For End If SToo(3) = Ax Next Else '仅计算小数部分 '(2by)^(1/2)*d^b Ax = Bn Bx = Bn hsfx = "+" Call Addition '2b Cx = SToo(3) Call mult_ '2by Ax = Cx Call Lnx Cx.bz = "" Cx.St = "5" Cx.zs = -1 Call mult_ Ax = Cx Call Exp '(2by)^(1/2) Px = Ax Ax = Dn 'd Call Lnx Cx = Bn 'b Call mult_ Ax = Cx Call Exp 'd^b Cx = Px Call mult_ Ax = Cx SToo(4) = Ax '2b/a*(b^2-a^2)^(1/2)d^b*y Ax = Ajn(1) Cx = Bn Call mult_ Ax = An Call multx_ Ax = Cx Bx = Cx hsfx = "+" Call Addition Px = Ax Ax = Dn Call Lnx Cx = Bn Call mult_ Ax = Cx Call Exp Cx = Px Call mult_ '*d^b Ax = SToo(3) Call mult_ Ax = SToo(3) Bx = Cx hsfx = "+" Call Addition End If Else '数值验证 ,=0 Call SJbz Call Dj2012 End If Else '如果 A=B 根x=A Ax = An SToo(3) = An End If '结果显示处理 Form1.Height = 5900 If Left(Xy.St, 1) = "1" Then '显示根 If Xspd = 0 Then Ax = SToo(3) xChen = xChen - fn0 Call sc(Ax.bz, Ax.St, Ax.zs) Else '只计算小数部分时 For Ii = xChen + 20 To Bn.zs + 1 Step -1 If Mid(Bn.St, Ii, 1) <> "0" Then Bn.St = Left(Bn.St, Ii) End If Next Call sc(Ax.bz, Ax.St, Ax.zs) Ax = SToo(3) cv = Left(Bn.St, Bn.zs + 1) + "." + IIf(Ii - Bn.zs = 0, "", Mid(Bn.St, Bn.zs + 1, Ii - Bn.zs)) + " + " + Left(Ax.St, 1) + "." + Mid(Ax.St, 2, ychen - 1) + " E" + CStr(Ax.zs) '在下窗口按精度显示 Text11.Text = cv End If Else '显示验证结果 Ax = SToo(4) Call sc(Ax.bz, Ax.St, Ax.zs) End If xChen = xChen0 xP = Xp0

云梦

Function Dj2012() '原方程 Yn1 = Ax Cx = Ax Call mult_ Xx = Cx 'x^2 Ax = An Cx = An Call mult_ Yx = Cx 'A^2 Ax = Bn Cx = Bn Call mult_ Rx = Cx 'B^2 Ax = Xx Bx = Yx hsfx = "-" 'x^2-a^2 Call Addition Call Lnx Ajn(10) = Ax 'ln(x^2-a^2) Cx.bz = "" Cx.St = "5" Cx.zs = -1 Call mult_ Ax = Cx Ajn(7) = Ax 'Ln(x^2 - a^2)^(1/2)) Call Exp Ajn(2) = Ax '(x^2 - a^2)^(1/2) Ax = Yn1 Bx = Ajn(2) hsfx = "-" Call Addition Ajn(4) = Ax 'x-(x^2 - a^2)^(1/2) Ax = Xx Bx = Rx 'x^2-b^2 hsfx = "-" Call Addition Call Lnx Cx.bz = "" Cx.St = "5" Cx.zs = -1 Call mult_ Ax = Cx Ajn(8) = Ax 'Ln(x^2 - b^2)^(1/2)) Call Exp Ajn(3) = Ax '(x^2 - b^2)^(1/2) Ax = Yn1 Bx = Ajn(3) hsfx = "-" Call Addition Ajn(5) = Ax 'x-(x^2 - b^2)^(1/2) Call Lnx Ax = Ajn(4) Call Lnx Qbx = Ax 'Ln(x-(x^2-a^2)^(1/2)) Ax = Ajn(5) Call Lnx Qby = Ax 'Ln(x-(x^2-b^2)^(1/2)) Ax = Bn Cx = An Call multx_ 'a/b Ajn(6) = Cx Ax = Cx Call Lnx 'ln(a/b) Ajn(9) = Ax Bx = Ajn(8) hsfx = "+" Call Addition Bx = Ajn(7) hsfx = "-" Call Addition Call Exp Wn = Ax '[a(x^2-b^2)^(1/2)/(b(x^2-a^2)^(1/2)] ok Ax = Qbx Bx = Qby hsfx = "-" Call Addition 'Ln[a(x-(x^2-b^2)^(1/2))/(b(x-(x^2-a^2)^(1/2))] Bx = Ajn(9) hsfx = "-" Call Addition 'b/a Cx = Yn1 Call mult_ Ax = Cx Call Exp Bx = Ax Ax = Wn hsfx = "-" Call Addition SToo(4) = Ax End Function

云梦

Function Dj2012F() '求导数 '1) Ax = Ajn(7) Bx = Ajn(8) hsfx = "+" Call Addition Bx = Ax Ax = Ajn(9) hsfx = "-" Call Addition Call Exp Cx = Yn1 Call mult_ Ajn(14) = Cx 'ax/(b(x^2-a^2)^(1/2)(x^2-b^2)^(1/2)) ok '2) Ax = Ajn(8) Bx = Ajn(7) hsfx = "-" Call Addition Bx = Ajn(10) hsfx = "-" Call Addition Bx = Ajn(9) hsfx = "+" Call Addition Call Exp Cx = Yn1 Call mult_ Ajn(15) = Cx 'ax(x^2-b^2)/(b(x^2-a^2)^(3/2)) ok '3) Ax = Ajn(3) Cx = Yn1 Call multx_ Bx = Cx Ax.bz = "" Ax.St = "1" Ax.zs = 0 hsfx = "-" Call Addition '1-x/(x^2-b^2)^(1/2) Tx = Ax '4) Ax = Qby Bx = Ax hsfx = "+" Call Addition Bx = Ax Ax = Qbx hsfx = "-" Call Addition Bx = Ajn(9) hsfx = "-" Call Addition Call Exp Cx = Tx Call mult_ 'b(x-(x^2-a^2)^(1/2))(1-x/(x^2-b^2)^(1/2))/a/(x-(x^2-b^2)^(1/2)) ok Tx = Cx '5) Ax = Ajn(2) Cx = Yn1 Call multx_ Bx = Cx Ax.bz = "" Ax.St = "1" Ax.zs = 0 hsfx = "-" Call Addition '1-x/(x^2-a^2)^(1/2) Cx = Ax Ax = Ajn(5) Call multx_ Ax = Ajn(6) Call multx_ Ax = Cx Bx = Tx hsfx = "-" Call Addition Cx = Ajn(5) Call mult_ Ax = Ajn(4) Call multx_ Ax = Ajn(6) Call mult_ Ax = Yn1 Call mult_ Tx = Cx '6) Ax = Qbx Bx = Qby hsfx = "-" Call Addition Ax = Ax Bx = Ajn(9) hsfx = "-" Call Addition Bx = Tx hsfx = "+" Call Addition Tx = Ax '7) Ax = Qbx Bx = Qby hsfx = "-" Call Addition Bx = Ajn(9) hsfx = "-" Call Addition Cx = Yn1 Call mult_ Ax = Cx Call Exp Cx = Tx Call mult_ Bx = Cx Ax = Ajn(14) hsfx = "-" Call Addition Bx = Ajn(15) hsfx = "-" Call Addition End Function

gxqcn

不需要什么特殊代码，只是利用了该方程的特点,x=a+s*10^-r。这样就可以把整数和小数分离，只计算小数部分。83楼也就是利用了近似的特点得到了近似值。不过他的x值只是我程序中x的初始值罢了。 下面就需要用牛顿法进 ... 云梦 发表于 2013-1-17 07:59

这也可看做是一种特殊的优化。 打个比方，要算1000万个连续“7”构成的大整数的平方，如果能敏感的抓住其特征，转化为

$(7//9xx(10^1000-1))^2=(10^2000-2xx10^1000+1)xx49//81$

这样，将比普通的通用算法效率高很多。 我在 HugeCalc 里就大量应用了类似优化， 北大的刘楚雄老师曾以该特殊值（因为好输入）测试 HugeCalc， 结果速度之快令其大吃一惊！

云梦

当然，计算器毕竟就是计算器，离不开按钮，但是每个按钮下的文章就无法想象了。如果这个方程经常被使用，当然可以放在按钮下，解决了很多不易记忆的命令行和不必要的重复劳动，这正是本计算器的目的之所在。 唯一最大的弱点就是计算速度不令人满意,HugeCalc 我也想利用，可小数转换成整数很麻烦，且计算中不需要那么多的整数，能保证一定精度就可以了。大整数运算需要消耗大量的存储资源，用VB调用HugeCalc中的函数虽然基本运算速度大大提高，可毕竟涉及的是大量的小数，因此转换消费的时间降低了效率。 作为计算器具有128位精度其实已经大大超出了常规计算器，且丰富的函数也是现计算器中没有的，更高的精度会降低运算速度，所以常用的是64位和128位，256位以上似乎没有什么意义，但仍保留到1024位，为特殊人群服务。