Rosetta Code

Pathological floating point problems: Difference between revisions

Page Discussion

Pathological floating point problems (view source)

Revision as of 11:33, 12 August 2019

21,383 bytes added ,  4 years ago
→‎{{header|Visual Basic .NET}}: Added VB.NET
Line 2,790:
u(100) : 100
</pre>
 
=={{header|Visual Basic .NET}}==
'''Compiler:''' Roslyn Visual Basic (language version >=15.3, e.g. with Visual Studio 2015)
{{works with|.NET Framework|4.6.2}}
{{works with|.NET Core|2.1}}
{{libheader|BigRationalLibrary|1.0.0}}
 
.NET has three built-in non-integral types, Single, Double, and Decimal.
 
Single and Double are 32-bit and 64-bit IEEE floating-point numbers, respectively, and so are susceptible to the rounding errors that this task is designed to demonstrate.
 
Decimal is a fixed-point number intended for use in financial calculations where a high number of significant digits are required and magnitudes are reasonably small. It has a 128-bit size and is stored as a 1-bit sign, 96-bit integer, and scaling factor (from 1 to 10^28) that determines the position of the decimal point. Decimal math in .NET is notoriously slow.
 
.NET Framework and Core do have the arbitrary-precision BigInteger structure, but do not have an arbitrary-precision rational type, though Microsoft did create a prototype BigRational that never made it into the library. Its source has been released on GitHub and a NuGet package was published for it in 2011 (at [https://www.nuget.org/packages/BigRationalLibrary/]). The type is fully functional and is sufficient for these tasks, so this program has been designed to optionally use it (and must be compiled with a reference to the library to do so).
 
The program defines three compiler constants that toggle:
* USE_BIGRATIONAL: whether to use BigRationalLibrary; when False, a mock type and functions are defined that allows the code to be compiled without a reference to the library.
* BANDED_ROWS: whether to change the console colors to format output tables with alternating white-on-black and black-on-white rows.
* INCREASED_LIMITS: whether to additionally display n = 1000 for Wrong Convergence Sequence and up to year 40 for The Chaotic Bank Society.
 
Because of operator overloading, the implementations are visually very similar. Integral literals are used where possible as they implicitly convert to the other types.
 
The non-floating-point implementations calculate e using the reciprocals of factorials formula because Math.E is a double and is not sufficiently precise.
 
Because BigRational.ToString() returns a string containing its exact numerator and denominator, the BigRational results are converted to Decimal for display.
 
The following sections source code must be located in a single file.
 
Main() procedure and output formatting:
 
<lang vbnet>Imports System.Globalization
Imports System.Numerics
Imports Numerics
 
#Const USE_BIGRATIONAL = True
#Const BANDED_ROWS = True
#Const INCREASED_LIMITS = True
 
#If Not USE_BIGRATIONAL Then
' Mock structure to make test code work.
Structure BigRational
Overrides Function ToString() As String
Return "NOT USING BIGRATIONAL"
End Function
Shared Narrowing Operator CType(value As BigRational) As Decimal
Return -1
End Operator
End Structure
#End If
 
Module Common
Public Const FMT_STR = "{0,4} {1,-15:G9} {2,-24:G17} {3,-32} {4,-32}"
Public ReadOnly Property Headings As String =
String.Format(CultureInfo.InvariantCulture,
FMT_STR,
{"N", "Single", "Double", "Decimal", "BigRational (rounded as decimal)"})
 
<Conditional("BANDED_ROWS")>
Sub SetConsoleFormat(n As Integer)
If n Mod 2 = 0 Then
Console.BackgroundColor = ConsoleColor.Black
Console.ForegroundColor = ConsoleColor.White
Else
Console.BackgroundColor = ConsoleColor.White
Console.ForegroundColor = ConsoleColor.Black
End If
End Sub
 
Function FormatOutput(n As Integer, x As (sn As Single, db As Double, dm As Decimal, br As BigRational)) As String
SetConsoleFormat(n)
Return String.Format(CultureInfo.CurrentCulture, FMT_STR, n, x.sn, x.db, x.dm, CDec(x.br))
End Function
 
Sub Main()
WrongConvergence()
Console.WriteLine()
ChaoticBankSociety()
 
Console.WriteLine()
SiegfriedRump()
 
SetConsoleFormat(0)
End Sub
End Module</lang>
 
===Task 1: Converging sequence===
Somewhat predictably, Single fairs the worst and Double slightly better. Decimal lasts past iteration 20, and BigRational remains exactly precise.
 
<lang vbnet>Module Task1
Iterator Function SequenceSingle() As IEnumerable(Of Single)
' n, n-1, and n-2
Dim vn, vn_1, vn_2 As Single
vn_2 = 2
vn_1 = -4
 
Do
Yield vn_2
vn = 111 - (1130 / vn_1) + (3000 / (vn_1 * vn_2))
vn_2 = vn_1
vn_1 = vn
Loop
End Function
 
Iterator Function SequenceDouble() As IEnumerable(Of Double)
' n, n-1, and n-2
Dim vn, vn_1, vn_2 As Double
vn_2 = 2
vn_1 = -4
 
Do
Yield vn_2
vn = 111 - (1130 / vn_1) + (3000 / (vn_1 * vn_2))
vn_2 = vn_1
vn_1 = vn
Loop
End Function
 
Iterator Function SequenceDecimal() As IEnumerable(Of Decimal)
' n, n-1, and n-2
Dim vn, vn_1, vn_2 As Decimal
vn_2 = 2
vn_1 = -4
 
' Use constants to avoid calling the Decimal constructor in the loop.
Const i11 As Decimal = 111
Const i130 As Decimal = 1130
Const E000 As Decimal = 3000
 
Do
Yield vn_2
vn = i11 - (i130 / vn_1) + (E000 / (vn_1 * vn_2))
vn_2 = vn_1
vn_1 = vn
Loop
End Function
 
#If USE_BIGRATIONAL Then
Iterator Function SequenceRational() As IEnumerable(Of BigRational)
' n, n-1, and n-2
Dim vn, vn_1, vn_2 As BigRational
vn_2 = 2
vn_1 = -4
 
' Same reasoning as for Decimal.
Dim i11 As BigRational = 111
Dim i130 As BigRational = 1130
Dim E000 As BigRational = 3000
 
Do
Yield vn_2
vn = i11 - (i130 / vn_1) + (E000 / (vn_1 * vn_2))
vn_2 = vn_1
vn_1 = vn
Loop
End Function
#Else
Iterator Function SequenceRational() As IEnumerable(Of BigRational)
Do
Yield Nothing
Loop
End Function
#End If
 
<Conditional("INCREASED_LIMITS")>
Sub IncreaseMaxN(ByRef arr As Integer())
ReDim Preserve arr(arr.Length)
arr(arr.Length - 1) = 1000
End Sub
 
Sub WrongConvergence()
Console.WriteLine("Wrong Convergence Sequence:")
 
Dim displayedIndices As Integer() = {3, 4, 5, 6, 7, 8, 20, 30, 50, 100}
IncreaseMaxN(displayedIndices)
 
Dim indicesSet As New HashSet(Of Integer)(displayedIndices)
 
Console.WriteLine(Headings)
 
Dim n As Integer = 1
' Enumerate the implementations in parallel as tuples.
For Each x In SequenceSingle().
Zip(SequenceDouble(), Function(sn, db) (sn, db)).
Zip(SequenceDecimal(), Function(a, dm) (a.sn, a.db, dm)).
Zip(SequenceRational(), Function(a, br) (a.sn, a.db, a.dm, br))
If n > displayedIndices.Max() Then Exit For
 
If indicesSet.Contains(n) Then
Console.WriteLine(FormatOutput(n, x))
End If
 
n += 1
Next
End Sub
End Module</lang>
 
{{out}}
Note that "Wrong" is not a heading for the first column--a weird optical effect, that would be.
<pre>Wrong Convergence Sequence:
N Single Double Decimal BigRational (rounded as decimal)
3 18.5 18.5 18.5 18.5
4 9.37837791 9.378378378378379 9.378378378378378378378378378 9.378378378378378378378378378
5 7.80114746 7.8011527377521688 7.80115273775216138328530259 7.8011527377521613832853025937
6 7.15434647 7.1544144809753334 7.154414480975249353527890606 7.1544144809752493535278906539
7 6.80583048 6.8067847369248113 6.806784736923632983941755925 6.8067847369236329839417565963
8 6.57857943 6.592632768721792 6.592632768704438392741992887 6.5926327687044383927420027764
20 100 98.349503122165359 6.04355210719488789087813234 6.0435521101892688677774773641
30 100 99.999999999998934 101.88552052291609961584734802 6.006786093031205758530554048
50 100 100 100.00000000000000000000000068 6.0001758466271871889456140207
100 100 100 100.0 6.0000000193194779291040868034
1000 100 100 100.0 6.0000000000000000000000000000</pre>
 
===Task 2: The Chaotic Bank Society===
Decimal appears to be doing well by year 25, but begins to degenerate ''the very next year'' and soon overflows.
 
<lang vbnet>Module Task2
Iterator Function ChaoticBankSocietySingle() As IEnumerable(Of Single)
Dim balance As Single = Math.E - 1
Dim year As Integer = 1
 
Do
balance = (balance * year) - 1
Yield balance
year += 1
Loop
End Function
Iterator Function ChaoticBankSocietyDouble() As IEnumerable(Of Double)
Dim balance As Double = Math.E - 1
Dim year As Integer = 1
 
Do
balance = (balance * year) - 1
Yield balance
year += 1
Loop
End Function
 
Iterator Function ChaoticBankSocietyDecimal() As IEnumerable(Of Decimal)
' 27! is the largest factorial decimal can represent.
Dim balance As Decimal = CalculateEDecimal(27) - Decimal.One
Dim year As Integer = 1
 
Do
balance = (balance * year) - Decimal.One
Yield balance
year += 1
Loop
End Function
 
#If USE_BIGRATIONAL Then
Iterator Function ChaoticBankSocietyRational() As IEnumerable(Of BigRational)
' 100 iterations is precise enough for 25 years.
Dim balance As BigRational = CalculateEBigRational(100) - BigRational.One
Dim year As Integer = 1
 
Do
balance = (balance * year) - BigRational.One
Yield balance
year += 1
Loop
End Function
#Else
Iterator Function ChaoticBankSocietyRational() As IEnumerable(Of BigRational)
Do
Yield Nothing
Loop
End Function
#End If
 
Function CalculateEDecimal(terms As Integer) As Decimal
Dim e As Decimal = 1
Dim fact As Decimal = 1
 
For i As Integer = 1 To terms
fact *= i
e += Decimal.One / fact
Next
 
Return e
End Function
 
#If USE_BIGRATIONAL Then
Function CalculateEBigRational(terms As Integer) As BigRational
Dim e As BigRational = 1
Dim fact As BigInteger = 1
 
For i As Integer = 1 To terms
fact *= i
e += BigRational.Invert(fact)
Next
 
Return e
End Function
#End If
 
<Conditional("INCREASED_LIMITS")>
Sub IncreaseMaxYear(ByRef year As Integer)
year = 40
End Sub
 
Sub ChaoticBankSociety()
Console.WriteLine("Chaotic Bank Society:")
Console.WriteLine(Headings)
 
Dim maxYear As Integer = 25
IncreaseMaxYear(maxYear)
 
Dim i As Integer = 0
For Each x In ChaoticBankSocietySingle().
Zip(ChaoticBankSocietyDouble(), Function(sn, db) (sn, db)).
Zip(ChaoticBankSocietyDecimal(), Function(a, dm) (a.sn, a.db, dm)).
Zip(ChaoticBankSocietyRational(), Function(a, br) (a.sn, a.db, a.dm, br))
If i >= maxYear Then Exit For
Console.WriteLine(FormatOutput(i + 1, x))
i += 1
Next
End Sub
End Module</lang>
 
{{out}}
<pre>Chaotic Bank Society:
N Single Double Decimal BigRational (rounded as decimal)
1 0.718281865 0.71828182845904509 0.7182818284590452353602874714 0.7182818284590452353602874713
2 0.43656373 0.43656365691809018 0.4365636569180904707205749428 0.4365636569180904707205749427
3 0.309691191 0.30969097075427054 0.3096909707542714121617248284 0.3096909707542714121617248281
4 0.238764763 0.23876388301708218 0.2387638830170856486468993136 0.2387638830170856486468993124
5 0.193823814 0.1938194150854109 0.1938194150854282432344965680 0.1938194150854282432344965623
6 0.162942886 0.16291649051246537 0.1629164905125694594069794080 0.1629164905125694594069793739
7 0.140600204 0.14041543358725761 0.1404154335879862158488558560 0.1404154335879862158488556174
8 0.124801636 0.12332346869806088 0.1233234687038897267908468480 0.1233234687038897267908449393
9 0.123214722 0.10991121828254791 0.1099112183350075411176216320 0.1099112183350075411176044541
10 0.232147217 0.099112182825479067 0.0991121833500754111762163200 0.0991121833500754111760445416
11 1.55361938 0.090234011080269738 0.0902340168508295229383795200 0.0902340168508295229364899583
12 17.6434326 0.082808132963236858 0.0828082022099542752605542400 0.0828082022099542752378795006
13 228.364624 0.076505728522079153 0.0765066287294055783872051200 0.0765066287294055780924335089
14 3196.10474 0.071080199309108139 0.0710928022116780974208716800 0.0710928022116780932940691248
15 47940.5703 0.066202989636622078 0.0663920331751714613130752000 0.0663920331751713994110368720
16 767048.125 0.059247834185953252 0.0622725308027433810092032000 0.0622725308027423905765899521
17 13039817 0.0072131811612052843 0.0586330236466374771564544000 0.0586330236466206398020291865
18 234716704 -0.87016273909830488 0.0553944256394745888161792000 0.0553944256391715164365253585
19 4.45961728E+09 -17.533092042867793 0.0524940871500171875074048000 0.0524940871442588122939818127
20 8.91923415E+10 -351.66184085735586 0.0498817430003437501480960000 0.0498817428851762458796362544
21 1.8730392E+12 -7385.898658004473 0.0475166030072187531100160000 0.0475166005887011634723613427
22 4.12068642E+13 -162490.77047609841 0.0453652661588125684203520000 0.0453652129514255963919495414
23 9.47757884E+14 -3737288.7209502636 0.0434011216526890736680960000 0.0433998978827887170148394524
24 2.27461897E+16 -89694930.302806318 0.0416269196645377680343040000 0.0415975491869292083561468582
25 5.68654735E+17 -2242373258.570158 0.0406729916134442008576000000 0.0399387296732302089036714552
26 1.47850229E+19 -58301704723.824112 0.0574977819495492222976000000 0.0384069715039854314954578354
27 3.99195623E+20 -1574146027544.251 0.5524401126378290020352000000 0.0369882306076066503773615576
28 1.11774772E+22 -44076088771240.031 14.468323153859212056985600000 0.0356704570129862105661236148
29 3.24146835E+23 -1278206574365962 418.58137146191714965258240000 0.0344432533766001064175848308
30 9.72440521E+24 -38346197230978864 12556.441143857514489577472000 0.0332976012980031925275449256
31 3.01456563E+26 -1.1887321141603448E+18 389248.67545958294917690163200 0.0322256402380989683538926936
32 9.64661E+27 -3.8039427653131035E+19 12455956.614706654373660852224 0.0312204876191669873245661979
33 3.18338125E+29 -1.2553011125533242E+21 411046567.28531959433080812339 0.0302760914325105817106845313
34 1.08234959E+31 -4.268023782681302E+22 13975583286.700866207247476195 0.0293871087053597781632740646
35 3.78822341E+32 -1.4938083239384556E+24 489145415033.53031725366166682 0.0285488046875922357145922644
36 1.36376043E+34 -5.3777099661784406E+25 17609234941206.091421131820006 0.0277569687533204857253215188
37 5.04591372E+35 -1.989752687486023E+27 651541692824624.38258187734022 0.0270078438728579718368961961
38 1.91744716E+37 -7.5610602124468873E+28 24758584327335725.538111338928 0.0262980671686029298020554545
39 ∞ -2.9488134828542859E+30 965584788766093294.9863422182 0.0256246195755142622801627278
40 ∞ -1.1795253931417144E+32 38623391550643731798.453688728 0.0249847830205704912065091156</pre>
 
===Task 3: Rump's example===
Because Decimal is fixed point, the polynomial as displayed in the task cannot be evaluated for the specified input because of intermediate values (b<sup>8</sup>, for instance) being too large.
 
Possibly due to JIT optimizations (''and nasal demons resulting from attempting to ensure that floating-point numbers are (un)equal''), the outputs for Single and Double are identical (both equal to that for Double in debug builds) for optimized (release) builds targeting x86 or AnyCPU-prefer-x86. Since the the machine the author used is x64, the Single was likely stored (with extra precision as the CLR specification allows) in a 64-bit register.
 
Each of these three tasks are susceptible to this phenomenon; the outputs of this program for floating-point arithmetic should not be expected to be the same on different systems.
 
<lang vbnet>Module Task3
Function SiegfriedRumpSingle(a As Single, b As Single) As Single
Dim a2 = a * a,
b2 = b * b,
b4 = b2 * b2,
b6 = b4 * b2
 
' Non-integral literals must be coerced to Single using the type suffix.
Return 333.75F * b6 +
(a2 * (
11 * a2 * b2 -
b6 -
121 * b4 -
2)) +
5.5F * b4 * b4 +
a / (2 * b)
End Function
 
Function SiegfriedRumpDouble(a As Double, b As Double) As Double
Dim a2 = a * a,
b2 = b * b,
b4 = b2 * b2,
b6 = b4 * b2
 
' Non-integral literals are Doubles by default.
Return 333.75 * b6 +
(a2 * (
11 * a2 * b2 -
b6 -
121 * b4 -
2)) +
5.5 * b4 * b4 +
a / (2 * b)
End Function
 
Function SiegfriedRumpDecimal(a As Decimal, b As Decimal) As Decimal
Dim a2 = a * a,
b2 = b * b,
b4 = b2 * b2,
b6 = b4 * b2
 
' The same applies for Decimal.
Return 333.75D * b6 +
(a2 * (
11 * a2 * b2 -
b6 -
121 * b4 -
2)) +
5.5D * b4 * b4 +
a / (2 * b)
End Function
 
#If USE_BIGRATIONAL Then
Function SiegfriedRumpRational(a As BigRational, b As BigRational) As BigRational
' Use mixed number constructor to maintain exact precision (333+3/4, 5+1/2).
Dim c1 As New BigRational(333, 3, 4)
Dim c2 As New BigRational(5, 1, 2)
 
Dim a2 = a * a,
b2 = b * b,
b4 = b2 * b2,
b6 = b4 * b2
 
Return c1 * b6 +
(a2 * (
11 * a2 * b2 -
b6 -
121 * b4 -
2)) +
c2 * b4 * b4 +
a / (2 * b)
End Function
#Else
Function SiegfriedRumpRational(a As Integer, b As Integer) As BigRational
Return Nothing
End Function
#End If
 
Sub SiegfriedRump()
Console.WriteLine("Siegfried Rump Formula:")
Dim a As Integer = 77617
Dim b As Integer = 33096
 
Console.Write("Single: ")
Dim sn As Single = SiegfriedRumpSingle(a, b)
Console.WriteLine("{0:G9}", sn)
Console.WriteLine()
 
Console.Write("Double: ")
Dim db As Double = SiegfriedRumpDouble(a, b)
Console.WriteLine("{0:G17}", db)
Console.WriteLine()
 
Console.WriteLine("Decimal:")
Dim dm As Decimal
Try
dm = SiegfriedRumpDecimal(a, b)
Catch ex As OverflowException
Console.WriteLine("Exception: " + ex.Message)
End Try
Console.WriteLine($" {dm}")
Console.WriteLine()
 
Console.WriteLine("BigRational:")
Dim br As BigRational = SiegfriedRumpRational(a, b)
Console.WriteLine($" Rounded: {CDec(br)}")
Console.WriteLine($" Exact: {br}")
End Sub
End Module</lang>
 
{{out|note=debug build}}
<pre>Siegfried Rump Formula
Single: -6.27969943E+29
 
Double: -2.3611832414348226E+21
 
Decimal:
Exception: Value was either too large or too small for a Decimal.
0
 
BigRational:
Rounded: -0.8273960599468213681411650955
Exact: -54767/66192</pre>
 
{{out|note=release build}}
<pre>Siegfried Rump Formula:
Single: -2.36118324E+21
 
Double: -2.3611832414348226E+21
 
Decimal:
Exception: Value was either too large or too small for a Decimal.
0
 
BigRational:
Rounded: -0.8273960599468213681411650955
Exact: -54767/66192</pre>
 
=={{header|zkl}}==
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/264651"