Rosetta Code

Pathological floating point problems: Difference between revisions

Page Discussion

Pathological floating point problems (view source)

Revision as of 14:18, 12 August 2019

17,700 bytes added ,  4 years ago
→‎{{header|C sharp}}: Added C#
Line 632:
v_100 : 6.000000019319477929060000000000000000000000000000000000000000000000000000000000
</pre>
 
=={{header|C sharp}}==
{{trans|Visual Basic .NET}}
'''Compiler:''' Roslyn C# (language version >=7.3)
{{works with|.NET Framework|4.6.2}}
{{works with|.NET Core|2.1}}
{{libheader|BigRationalLibrary|1.0.0}}
See VB.NET entry for details (Single is float in C#; Double is double, and Decimal is decimal).
 
The following sections source code must be located in a single file.
 
<lang csharp>#define USE_BIGRATIONAL
#define BANDED_ROWS
#define INCREASED_LIMITS
 
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Globalization;
using System.Linq;
using System.Numerics;
using Numerics;
 
using static Common;
using static Task1;
using static Task2;
using static Task3;
 
#if !USE_BIGRATIONAL
// Mock structure to make test code work.
struct BigRational
{
public override string ToString() => "NOT USING BIGRATIONAL";
public static explicit operator decimal(BigRational value) => -1;
}
#endif
 
static class Common
{
public const string FMT_STR = "{0,4} {1,-15:G9} {2,-24:G17} {3,-32} {4,-32}";
public static string Headings { get; } =
string.Format(
CultureInfo.InvariantCulture,
FMT_STR,
new[] { "N", "Single", "Double", "Decimal", "BigRational (rounded as Decimal)" });
 
[Conditional("BANDED_ROWS")]
static void SetConsoleFormat(int n)
{
if (n % 2 == 0)
{
Console.BackgroundColor = ConsoleColor.Black;
Console.ForegroundColor = ConsoleColor.White;
}
else
{
Console.BackgroundColor = ConsoleColor.White;
Console.ForegroundColor = ConsoleColor.Black;
}
}
 
public static string FormatOutput(int n, (float sn, double db, decimal dm, BigRational br) x)
{
SetConsoleFormat(n);
return string.Format(CultureInfo.CurrentCulture, FMT_STR, n, x.sn, x.db, x.dm, (decimal)x.br);
}
 
static void Main()
{
WrongConvergence();
 
Console.WriteLine();
ChaoticBankSociety();
 
Console.WriteLine();
SiegfriedRump();
 
SetConsoleFormat(0);
}
}</lang>
 
===Task 1: Converging sequence===
'''See''' [[#VB.NET Task 1]]
 
<lang csharp>static class Task1
{
public static IEnumerable<float> SequenceSingle()
{
// n, n-1, and n-2
float vn, vn_1, vn_2;
vn_2 = 2;
vn_1 = -4;
 
while (true)
{
yield return vn_2;
vn = 111f - (1130f / vn_1) + (3000f / (vn_1 * vn_2));
vn_2 = vn_1;
vn_1 = vn;
}
}
 
public static IEnumerable<double> SequenceDouble()
{
// n, n-1, and n-2
double vn, vn_1, vn_2;
vn_2 = 2;
vn_1 = -4;
 
while (true)
{
yield return vn_2;
vn = 111 - (1130 / vn_1) + (3000 / (vn_1 * vn_2));
vn_2 = vn_1;
vn_1 = vn;
}
}
 
public static IEnumerable<decimal> SequenceDecimal()
{
// n, n-1, and n-2
decimal vn, vn_1, vn_2;
vn_2 = 2;
vn_1 = -4;
 
// Use constants to avoid calling the Decimal constructor in the loop.
const decimal i11 = 111;
const decimal i130 = 1130;
const decimal E000 = 3000;
 
while (true)
{
yield return vn_2;
vn = i11 - (i130 / vn_1) + (E000 / (vn_1 * vn_2));
vn_2 = vn_1;
vn_1 = vn;
}
}
 
#if USE_BIGRATIONAL
public static IEnumerable<BigRational> SequenceRational()
{
// n, n-1, and n-2
BigRational vn, vn_1, vn_2;
vn_2 = 2;
vn_1 = -4;
 
// Same reasoning as for decimal.
BigRational i11 = 111;
BigRational i130 = 1130;
BigRational E000 = 3000;
 
while (true)
{
yield return vn_2;
vn = i11 - (i130 / vn_1) + (E000 / (vn_1 * vn_2));
vn_2 = vn_1;
vn_1 = vn;
}
}
#else
public static IEnumerable<BigRational> SequenceRational()
{
while (true) yield return default;
}
#endif
 
static void IncreaseMaxN(ref int[] arr)
{
int[] tmp = new int[arr.Length + 1];
arr.CopyTo(tmp, 0);
tmp[arr.Length] = 1000;
arr = tmp;
}
 
public static void WrongConvergence()
{
Console.WriteLine("Wrong Convergence Sequence:");
 
int[] displayedIndices = { 3, 4, 5, 6, 7, 8, 20, 30, 50, 100 };
IncreaseMaxN(ref displayedIndices);
 
var indicesSet = new HashSet<int>(displayedIndices);
 
Console.WriteLine(Headings);
 
int n = 1;
// Enumerate the implementations in parallel as tuples.
foreach (var x in SequenceSingle()
.Zip(SequenceDouble(), (sn, db) => (sn, db))
.Zip(SequenceDecimal(), (a, dm) => (a.sn, a.db, dm))
.Zip(SequenceRational(), (a, br) => (a.sn, a.db, a.dm, br)))
{
if (n > displayedIndices.Max()) break;
 
if (indicesSet.Contains(n))
Console.WriteLine(FormatOutput(n, x));
 
n++;
}
}
}</lang>
 
{{out}}
<pre>Wrong Convergence Sequence:
N Single Double Decimal BigRational (rounded as Decimal)
3 18.5 18.5 18.5 18.5
4 9.37837982 9.378378378378379 9.378378378378378378378378378 9.378378378378378378378378378
5 7.80116463 7.8011527377521688 7.80115273775216138328530259 7.8011527377521613832853025937
6 7.15456009 7.1544144809753334 7.154414480975249353527890606 7.1544144809752493535278906539
7 6.80883026 6.8067847369248113 6.806784736923632983941755925 6.8067847369236329839417565963
8 6.62275314 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===
'''See''' [[#VB.NET Task 2]]
 
<lang csharp>static class Task2
{
public static IEnumerable<float> ChaoticBankSocietySingle()
{
float balance = (float)(Math.E - 1);
 
for (int year = 1; ; year++)
yield return balance = (balance * year) - 1;
}
 
public static IEnumerable<double> ChaoticBankSocietyDouble()
{
double balance = Math.E - 1;
 
for (int year = 1; ; year++)
yield return balance = (balance * year) - 1;
}
 
public static IEnumerable<decimal> ChaoticBankSocietyDecimal()
{
// 27! is the largest factorial decimal can represent.
decimal balance = CalculateEDecimal(27) - 1;
 
for (int year = 1; ; year++)
yield return balance = (balance * year) - 1;
}
 
#if USE_BIGRATIONAL
public static IEnumerable<BigRational> ChaoticBankSocietyRational()
{
// 100 iterations is precise enough for 25 years.
BigRational brBalance = CalculateERational(100) - 1;
 
for (int year = 1; ; year++)
yield return brBalance = (brBalance * year) - 1;
}
#else
public static IEnumerable<BigRational> ChaoticBankSocietyRational()
{
while (true) yield return default;
}
#endif
 
public static decimal CalculateEDecimal(int terms)
{
decimal e = 1;
decimal fact = 1;
for (int i = 1; i <= terms; i++)
{
fact *= i;
e += decimal.One / fact;
}
 
return e;
}
 
#if USE_BIGRATIONAL
public static BigRational CalculateERational(int terms)
{
BigRational e = 1;
BigRational fact = 1;
for (int i = 1; i < terms; i++)
{
fact *= i;
e += BigRational.Invert(fact);
}
 
return e;
}
#endif
 
[Conditional("INCREASED_LIMITS")]
static void InceaseMaxYear(ref int year) => year = 40;
 
public static void ChaoticBankSociety()
{
Console.WriteLine("Chaotic Bank Society:");
Console.WriteLine(Headings);
 
int maxYear = 25;
InceaseMaxYear(ref maxYear);
 
int i = 0;
foreach (var x in ChaoticBankSocietySingle()
.Zip(ChaoticBankSocietyDouble(), (sn, db) => (sn, db))
.Zip(ChaoticBankSocietyDecimal(), (a, dm) => (a.sn, a.db, dm))
.Zip(ChaoticBankSocietyRational(), (a, br) => (a.sn, a.db, a.dm, br)))
{
if (i >= maxYear) break;
Console.WriteLine(FormatOutput(i + 1, x));
i++;
}
}
}</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.91923497E+10 -351.66184085735586 0.0498817430003437501480960000 0.0498817428851762458796362544
21 1.87303933E+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===
'''See''' [[#VB.NET Task 3]]
 
<lang csharp>static class Task3
{
public static float SiegfriedRumpSingle(float a, float b)
{
float
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);
}
 
public static double SiegfriedRumpDouble(double a, double b)
{
double
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 * b * b
- b6
- 121 * b4
- 2)
+ 5.5 * b4 * b4
+ a / (2 * b);
}
 
public static decimal SiegfriedRumpDecimal(decimal a, decimal b)
{
decimal
a2 = a * a,
b2 = b * b,
b4 = b2 * b2,
b6 = b4 * b2
;
 
// The same applies for decimal.
return
333.75m * b6
+ a2 * (
11 * a2 * b * b
- b6
- 121 * b4
- 2)
+ 5.5m * b4 * b4
+ a / (2 * b);
}
 
#if USE_BIGRATIONAL
public static BigRational SiegfriedRumpRational(BigRational a, BigRational b)
{
// Use mixed number constructor to maintain exact precision (333+3/4, 5+1/2).
var c1 = new BigRational(33375, 100);
var c2 = new BigRational(55, 10);
 
return c1 * BigRational.Pow(b, 6)
+ (a * a * (
11 * a * a * b * b
- BigRational.Pow(b, 6)
- 121 * BigRational.Pow(b, 4)
- 2))
+ c2 * BigRational.Pow(b, 8)
+ a / (2 * b);
}
#else
public static IEnumerable<BigRational> SiegfriedRumpRational()
{
while (true) yield return default;
}
#endif
 
public static void SiegfriedRump()
{
Console.WriteLine("Siegfried Rump Formula");
int a = 77617;
int b = 33096;
 
Console.Write("Single: ");
float sn = SiegfriedRumpSingle(a, b);
Console.WriteLine("{0:G9}", sn);
Console.WriteLine();
 
Console.Write("Double: ");
double db = SiegfriedRumpDouble(a, b);
Console.WriteLine("{0:G17}", db);
Console.WriteLine();
 
Console.WriteLine("Decimal:");
decimal dm = 0;
try
{
dm = SiegfriedRumpDecimal(a, b);
}
catch (OverflowException ex)
{
Console.WriteLine("Exception: " + ex.Message);
}
Console.WriteLine($" {dm}");
Console.WriteLine();
 
Console.WriteLine("BigRational:");
BigRational br = SiegfriedRumpRational(a, b);
Console.WriteLine($" Rounded: {(decimal)br}");
Console.WriteLine($" Exact: {br}");
}
}</lang>
 
{{out}}
Note that the output for Single is slightly different from VB.
<pre>Siegfried Rump Formula
Single: -6.338253E+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>
 
=={{header|Clojure}}==
Line 2,875 ⟶ 3,375:
End Module</lang>
 
<!--Anchors for the C# section to link to-->
===Task 1: Converging sequence===
==={{anchor|VB.NET Task 1}}Task 1: Converging sequence===
Somewhat predictably, Single fairs the worst and Double slightly better. Decimal lasts past iteration 20, and BigRational remains exactly precise.
 
Line 3,002 ⟶ 3,503:
1000 100 100 100.0 6.0000000000000000000000000000</pre>
 
==={{anchor|VB.NET Task 2}}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.
 
Line 3,153 ⟶ 3,654:
40 ∞ -1.1795253931417144E+32 38623391550643731798.453688728 0.0249847830205704912065091156</pre>
 
==={{anchor|VB.NET Task 3}}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.
 
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/264652"