double error

[Jack]

this is only an example of accumulated errors in double floating-point

Code: (Select All)

_Title "MatMul"

$Console:Only

_Dest _Console



Dim As Double a(3, 3), b(3, 3), c(3, 3), d(3, 3)

Dim As Long i



a(0, 0) = 1#: a(0, 1) = 1#: a(0, 2) = 1#: a(0, 3) = 1#

a(1, 0) = 2#: a(1, 1) = 4#: a(1, 2) = 8#: a(1, 3) = 16#

a(2, 0) = 3#: a(2, 1) = 9#: a(2, 2) = 27#: a(2, 3) = 81#

a(3, 0) = 4#: a(3, 1) = 16#: a(3, 2) = 64#: a(3, 3) = 256#



'b() = Inverse of a()

b(0, 0) = 4#: b(0, 1) = -3#: b(0, 2) = 4 / 3#: b(0, 3) = -1 / 4#

b(1, 0) = -13 / 3#: b(1, 1) = 19 / 4#: b(1, 2) = -7 / 3#: b(1, 3) = 11 / 24#

b(2, 0) = 3 / 2#: b(2, 1) = -2#: b(2, 2) = 7 / 6#: b(2, 3) = -1 / 4#

b(3, 0) = -1 / 6#: b(3, 1) = 1 / 4#: b(3, 2) = -1 / 6#: b(3, 3) = 1 / 24#



MatMul d(), a(), a()

MatMul c(), d(), b()

printm c()

Print

For i = 1 To 4

    MatMul d(), c(), c()

    MatMul c(), d(), b()

    printm c()

    Print

Next



Sub MatMul (ans() As Double, a() As Double, b() As Double)

    Dim As Long i, j, k, n, m, p

    If (UBound(a, 2) = UBound(b)) Then 'if valid dims

        n = UBound(a, 2)

        m = UBound(a)

        p = UBound(b, 2)



        For i = 0 To m

            For j = 0 To p

                For k = 0 To n

                    ans(i, j) = 0

        Next k, j, i



        For i = 0 To m

            For j = 0 To p

                For k = 0 To n

                    ans(i, j) = ans(i, j) + (a(i, k) * b(k, j))

        Next k, j, i

    Else

        Print "invalid dimensions"

    End If

End Sub



Sub printm (a() As Double)

    Dim As Long i, j, m, p

    m = UBound(a)

    p = UBound(a, 2)

    For i = 0 To m

        For j = 0 To p

            Print Using "####.##############"; a(i, j);

        Next j

        Print

    Next i

End Sub

output

Code: (Select All)

1.00000000000000   1.00000000000000   1.00000000000000   1.00000000000000
   2.00000000000011   4.00000000000000   8.00000000000011  15.99999999999997
   3.00000000000091   9.00000000000000  26.99999999999955  81.00000000000000
   4.00000000000000  16.00000000000000  64.00000000000182 256.00000000000000

   1.00000000000631   0.99999999999397   1.00000000000310   0.99999999999937
   2.00000000007003   3.99999999992497   8.00000000004081  15.99999999999153
   3.00000000030832   8.99999999965075  27.00000000019281  80.99999999995953
   4.00000000090040  15.99999999894135  64.00000000059481 255.99999999987494

   1.00000001287511   0.99999998735953   1.00000000618997   0.99999999878051
   2.00000017282946   3.99999983016824   8.00000008319284  15.99999998360744
   3.00000082591305   8.99999918818139  27.00000039771248  80.99999992162930
   4.00000253669350  15.99999750625284  64.00000122176243 255.99999975924175

   1.00003194036101   0.99996927234497   1.00001492581766   0.99999707177334
   2.00043128527955   3.99958508810619   8.00020154182073  15.99996046042222
   3.00206469592649   8.99801368628778  27.00096484445567  80.99981071149432
   4.00634687076126  15.99389407272611  64.00296593190251 255.99941812706993

   1.07888506286064   0.92429206751581   1.03674375970633   0.99279446622859
   3.06536455898004   2.97754344879331   8.49623463210287  15.90268727430023
   8.10050050155542   4.10492428973066  29.37575483774754  80.53410914096605
  19.67932001084409   0.95217113012768  71.30324807562465 254.56781665058406

[Pete]

If you want a computer you can count on, buy a calculator... or an abacus. Mine has self aware beads. I know that because they head of to New Orleans every year for Mardi Gras.

Hoping one of these days I can figure out how you handled roots with Dec-Float. String math is too slow for those calculations when carried out to several decimal places.

Pete