RosettaCodeData/Task/Pathological-floating-point.../Nim/pathological-floating-point...

116 lines
3.6 KiB
Nim

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)