真正的高精度科学计算器

云梦

感觉是特殊问题特殊处理，不可能处理一般情况
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位，为特殊人群服务。