真正的高精度科学计算器 - 第11页 - 开发应用 - 数学研发论坛

云梦发表于 2013-1-16 18:56:27

感觉是特殊问题特殊处理，不可能处理一般情况
mathe 发表于 2013-1-16 17:35 http://bbs.emath.ac.cn/images/common/back.gif
一般情况当然很好处理，也就不需要什么特殊的方法。针对不同的问题采取不同的方法才是灵活可取的。如果不灵活的话，那么解这个方程要求的精度太高了。

云梦发表于 2013-1-16 19:00:32

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

wayne 发表于 2013-1-16 23:46:50

针对不同的问题采取不同的方法才是灵活可取的。
云梦发表于 2013-1-16 18:56 http://bbs.emath.ac.cn/images/common/back.gif
大型软件如Mathematica，MATLAB都不敢说一条命令通吃所有问题，
你这个小小的计算器也不大可能点一下按钮，或者输入方程，敲一个命令，就给出图片中方程解的超高的精度。
所以我宁愿相信，你特地为这个方程写了一段代码，放在计算器里面，
然后说我的计算器可以完美的解决这种问题，而其他的数学软件不能。这有什么意义呢。
========
先不谈意义，就题论题，该问题总得有具体的算法实现吧，
不采用mathe的处理方式，我实在是想不出其他的方法，
能否共享一下此题的算法实现

https://bbs.emath.ac.cn/data/attachment/forum/month_1301/13011416222e3d36e2373bf1e9.jpg

云梦发表于 2013-1-17 07:59:30

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

云梦发表于 2013-1-17 08:14:01

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

云梦发表于 2013-1-17 08:17:46

如果你能看懂我的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

云梦发表于 2013-1-17 08:20:21

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 ' ok

Ax = Qbx
Bx = Qby
hsfx = "-"
Call Addition 'Ln
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

云梦发表于 2013-1-17 08:20:57

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 发表于 2013-1-17 08:26:30

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

这也可看做是一种特殊的优化。

打个比方，要算1000万个连续“7”构成的大整数的平方，如果能敏感的抓住其特征，转化为
(7//9xx(10^1000-1))^2=(10^2000-2xx10^1000+1)xx49//81
这样，将比普通的通用算法效率高很多。
我在 HugeCalc 里就大量应用了类似优化，
北大的刘楚雄老师曾以该特殊值（因为好输入）测试 HugeCalc，
结果速度之快令其大吃一惊！

云梦发表于 2013-1-17 09:01:25

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

页: 1 2 3 4 5 6 7 8 9 10 [11] 12 13

数学研发论坛's Archiver