How to resolve the algorithm Pathological floating point problems step by step in the Nim programming language

Published on 12 May 2024 09:40 PM

How to resolve the algorithm Pathological floating point problems step by step in the Nim programming language

Table of Contents

Problem Statement

Most programmers are familiar with the inexactness of floating point calculations in a binary processor. The classic example being: In many situations the amount of error in such calculations is very small and can be overlooked or eliminated with rounding. There are pathological problems however, where seemingly simple, straight-forward calculations are extremely sensitive to even tiny amounts of imprecision. This task's purpose is to show how your language deals with such classes of problems.

A sequence that seems to converge to a wrong limit. Consider the sequence:

As   n   grows larger, the series should converge to   6   but small amounts of error will cause it to approach   100.

Display the values of the sequence where   n =   3, 4, 5, 6, 7, 8, 20, 30, 50 & 100   to at least 16 decimal places.

The Chaotic Bank Society   is offering a new investment account to their customers. You first deposit   $e - 1   where   e   is   2.7182818...   the base of natural logarithms. After each year, your account balance will be multiplied by the number of years that have passed, and $1 in service charges will be removed. So ...

What will your balance be after   25   years?

Siegfried Rump's example.   Consider the following function, designed by Siegfried Rump in 1988.

Demonstrate how to solve at least one of the first two problems, or both, and the third if you're feeling particularly jaunty.

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Pathological floating point problems step by step in the Nim programming language

Source code in the nim programming language

import math, strutils, strformat
import decimal
import bignum

####################################################################################################
# Utilities.

proc format[T: DecimalType|Rat](val: T; intLen, fractLen: Positive): string =
  ## Format a decimal or a rational with "intLen" integer digits and "fractLen"
  ## fractional digits.
  let s = when T is DecimalType: ($val).split('.')
          else: ($val.toFloat).split('.')

  result = s[0].align(intLen) & '.'
  if s[1].len < fractLen:
    result.add s[1].alignLeft(fractLen, '0')
  else:
    result.add s[1][0..<fractLen]


proc `^`(a: Rat; b: Natural): Rat =
  ## Missing exponentiation operator for rationals.
  ## Adaptation of operator for floats.
  case b
  of 0: result = newRat(1)
  of 1: result = a
  of 2: result = a * a
  of 3: result = a * a * a
  else:
    var (a, b) = (a.clone, b)
    result = newRat(1)
    while true:
      if (b and 1) != 0:
        result *= a
      b = b shr 1
      if b == 0: break
      a *= a


####################################################################################################
# Task 1.

proc v[T: float|DecimalType|Rat](n: Positive): seq[T] =
  ## Return the "n" first values for sequence "Vn".
  var (v1, v2) = when T is float: (2.0, -4.0)
                 elif T is Rat: (newRat(2), newRat(-4))
                 else: (newDecimal(2), newDecimal(-4))

  result.add default(T)   # Dummy value to start at index one.
  result.add v1
  result.add v2
  for _ in 3..n:
    # Need to change evaluation order to avoid a bug with rationals.
    result.add 3000 / (result[^1] * result[^2]) - 1130 / result[^1] + 111


setPrec(130)  # Default precision is not sufficient.

let vfloat = v[float](100)
let vdecimal = v[DecimalType](100)
let vrational = v[Rat](100)

echo "Task 1"
echo "  n          v(n) float            v(n) decimal           v(n) rational"
for n in [3, 4, 5, 6, 7, 8, 20, 30, 50, 100]:
  echo &"{n:>3}    {vfloat[n]:>20.16f}   {vdecimal[n].format(3, 16)}   {vrational[n].format(3, 16)}"


####################################################################################################
# Task 2.

proc balance[T: float|DecimalType|Rat](): T =
  ## Return the balance after 25 years.
  result = when T is float: E - 1
           elif T is DecimalType: exp(newDecimal(1)) - 1
           else: newInt("17182818284590452353602874713526624977572470") /
                 newInt("10000000000000000000000000000000000000000000")

  var n = when T is float: 1.0 else: 1
  while n <= 25:
    result = result * n - 1
    n += 1

echo "\nTask 2."
echo "Balance after 25 years (float):     ", (&"{balance[float]():.16f}")[0..17]
echo "Balance after 25 years (decimal):   ", balance[DecimalType]().format(1, 16)
echo "Balance after 25 years: (rational): ", balance[Rat]().format(1, 16)


####################################################################################################
# Task 3.

const
  A = 77617
  B = 33096

proc rump[T: float|DecimalType|Rat](a, b: T): T =
  ## Return the value of the Rump's function.
  let C1 = when T is float: 333.75
           elif T is Rat: newRat(333.75)
           else: newDecimal("333.75")
  let C2 = when T is float: 5.5
           elif T is Rat: newRat(5.5)
           else: newDecimal("5.5")
  result = C1 * b^6 + a^2 * (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) + C2 * b^8 + a / (2 * b)

echo "\nTask 3"

let rumpFloat = rump(A.toFloat, B.toFloat)
let rumpDecimal = rump(newDecimal(A), newDecimal(B))
let rumpRational = rump(newRat(A), newRat(B))

echo &"f({A}, {B}) float =    ", rumpFloat
echo &"f({A}, {B}) decimal =  ", rumpDecimal.format(1, 16)
echo &"f({A}, {B}) rational = ", rumpRational.format(1, 16)


  

You may also check:How to resolve the algorithm Mertens function step by step in the SETL programming language
You may also check:How to resolve the algorithm Numerical and alphabetical suffixes step by step in the Julia programming language
You may also check:How to resolve the algorithm Bitmap/Histogram step by step in the Vedit macro language programming language
You may also check:How to resolve the algorithm Hello world/Text step by step in the Mathcad programming language
You may also check:How to resolve the algorithm Run-length encoding step by step in the XPL0 programming language