#! /usr/bin/python3
# File paranoia.py drives the Paranoia test. Execute this fie.
from __future__ import division # use // for integer division
import math
# ~~~~~~~ 095-open_comment
#
# William Kahan's Floating Point Test "Paranoia" recast in Python
#
# Jerome Coonen, v. 0.8 2020
#
# For an extended discussion of this program, please visit
# the Paranoia Study at www.arithmazium.org.
#
# Versions of Paranoia can be found on www.netlib.org:
# BASIC -- original source by William Kahan, 1982
# Pascal -- first translation, by Brian A. Wichmann, with
# David M. Gay and David G. Hough, 1986
# C -- adapted from Pascal by Thos Sumner, with
# David M. Gay, 1986
# Modula-2 -- adapted from Pascal by K. Y. Tan, 1986
# FORTRAN - single and double precision versions, 1986
# Forth -- adapted from C by Krishna Myneni, 2009
#
# BASIC version (C) 1983 by:
# Professor William M. Kahan,
# EECS Dept., Soda Hall, Room 441
# University of California
# Berkeley, California 94720 USA
#
# You may copy this program freely if you acknowledge its source.
# Please send comments on the Python version to
# Jerome Coonen at jcoonen@gmail.com
#
# ~~~~~~~ 100-initializations
# ----START CODE----
# Number of trials in half a dozen random and sequential tests.
NUM_TRIALS = 32
# Set all three flags to False for classic Paranoia flow.
verbose = False
ultra_verbose = True
hasty = True # When True, skips all user interaction
# ~~~~~~~ 110-globals
# Declarations of global variables other than CONSTANTS
err_failure = 0 # Indexes into the count and type lists below
err_serious = 1
err_defect = 2
err_flaw = 3
error_count = [0, 0, 0, 0]
error_type = ["Failure", "Serious", "Defect", "Flaw"]
fpecount = 0 # int count of exceptioins detected
milestone = 0 # int progress through the code
page_num = 1 # int vestige of glass tty pages
pow_err_cnt = 0 # int for testing x**y
save_pow_err_cnt = 0 # int for testing x**y
min_sqrt_error = math.nan
max_sqrt_error = math.nan
rand_seed = math.nan
rand_next = math.nan
# ~~~~~~~ 115-flags
flags = {
"mult_guard_digit": False,
"div_guard_digit": False,
"add_sub_guard_digit": False,
"mult_rounding": "other", # rounded and chopped are the alternatives
"div_rounding": "other",
"add_sub_rounding": "other",
"sqrt_rounding": "other",
"uses_sticky_bit": False,
"IEEE": False
}
# ~~~~~~~ 150-print_msgs
# Opening messages -- printed as execution begins.
def print_msgs(s):
"""Print a list of strings."""
for m in s:
print(m)
return
# ~~~~~~~ 155-instructions
def instructions():
"""Self-explanatory user instructions"""
instr = [
"Lest this program stop prematurely, i.e. before displaying\n",
" `END OF TEST',\n",
""
"try to persuade the computer NOT to terminate execution when an",
"error like Over/Underflow or Division by Zero occurs, but rather",
"to persevere with a surrogate value after, perhaps, displaying some",
"warning. If persuasion avails naught, don't despair but run this",
"program anyway to see how many milestones it passes, and then",
"amend it to make further progress.\n"
]
print_msgs(instr)
return
# ~~~~~~~ 160-heading
def heading():
"""Gives user a place to report issues."""
head = [
"Users are invited to help debug and augment this program so it will",
"cope with unanticipated and newly uncovered arithmetic pathologies.\n",
"Please send suggestions and interesting results to",
"\tJerome Coonen",
"\tjcoonen_at_gmail_dot_com\n\n",
"\t*(C) 1983 Prof. William M. Kahan",
"\t* You may copy this program freely if you acknowledge its source."
]
print_msgs(head)
return
# ~~~~~~~ 165-charcteristics
def characteristics():
"""Explain whtat the program reveals"""
chars = [
"Running this program should reveal these characteristics:",
" Radix B = 1, 2, 4, 8, 10, 16, 100, 256 ...",
" PRECISION = number of significant digits carried.",
" ULP_OF_ONE_PLUS = B/B^PRECISION = One Ulp",
" (Unit in the Last Place) of 1.000xxx .",
" ULP_OF_ONE_MINUS = 1/B^PRECISION = One Ulp",
" of numbers a little less than 1.0 .",
" Presence of guard digits for Mult., Div. and Subt.",
" Whether arithmetic is chopped, correctly rounded, or",
" something else for Mult., Div., Add/Subt. and Sqrt.",
" Whether a sticky bit is used correctly in rounding.",
" underflow_threshold = an underflow threshold.",
" min_positive and too_tiny_x tell whether underflow is",
" abrupt, gradual, or fuzzy.",
" overflow_threshold = an overflow threshold, roughly.",
" pos_big_B tells, roughly, whether Infinity is represented.",
" Comparisons are checked for consistency with subtraction",
" and for contamination with pseudo-zeros.",
" Sqrt is tested. Y^X is tested for integer X.",
" Extra-precise subexpressions are revealed but NOT YET tested.",
" Decimal-Binary conversion is NOT YET tested for accuracy."
]
print_msgs(chars)
return
# ~~~~~~~ 170-history
def history():
"""Present a bit of historical context and set user expectations."""
hist = [
"The program attempts to discriminate among",
" FLAWs, like lack of a sticky bit,",
" Serious DEFECTs, like lack of a guard digit, and",
" FAILUREs, like 2+2 == 5 .",
"Failures may confound subsequent diagnoses.",
"\nThe diagnostic capabilities of this program go beyond an earlier",
"program called `MACHAR', which can be found at the end of the",
"book `Software Manual for the Elementary Functions' (1980) by",
"W. J. Cody and W. Waite. Although both programs try to discover",
"the radix, precision, and range (over/underflow thresholds)",
"of the arithmetic, this program tries to cope with a wider variety",
"of pathologies, and to say how well the arithmetic is implemented.",
"\nThe program is based upon a conventional radix representation for",
"floating-point numbers, but also allows logarithmic encoding",
"(radix B=1) as used by certain early WANG machines.\n"
]
print_msgs(hist)
return
# ~~~~~~~ 200-utilities
# ----START UTILITIES----
def notify(s):
"""When the impossible happens, send up a flare."""
print("{:s} test appears to be inconsistent...".format(s))
print(" PLEASE NOTIFY THE AUTHOR!\n")
return
# ~~~~~~~ 210-maths
def fabs(x):
return math.fabs(x)
def floor(x):
return math.floor(x)
def log(x):
"""Return the natural log of the argument."""
return math.log(x)
def pow(x, y):
"""Return x**y."""
return math.pow(x, y)
def sqrt(x):
return math.sqrt(x)
# ~~~~~~~ 220-alt_input
def alt_input(s):
"""Local version of input allows hasty execution."""
if hasty:
print(s)
reply = "y"
else:
reply = input(s)
return reply
# ~~~~~~~ 222-pause
def pause():
"""Wait until the user hits RETURN."""
global milestone, page_num
s = alt_input("\nTo continue, press RETURN ")
print("\nDiagnosis resumes after milestone Number ", milestone,
" Page: ", page_num, "\n")
milestone += 1 # the only incremental change to milestone
page_num += 1
return
# ~~~~~~~ 225-test_cond
def test_cond(k, v, t):
"""Respond to test v of severity k, with error message t."""
if not v:
bad_cond(k, t)
print(".\n")
return
# ~~~~~~~ 230-bad_cond
def bad_cond(k, t):
"""Record and report error t of severity k.
Basic 920, for example.
"""
error_count[k] += 1
print(error_type[k] + ": " + t, end="")
return
# ~~~~~~~ 235-rand_frac
def rand_frac():
"""Generate random fractions.
Compute x = (rand_next + rand_seed)^5
rand_next = x - floor(x) + 0.000005 * x
Globals rand_next and rand_seed are initialized when needed.
Returns the new value of rand_next
"""
# ~~~~~~~ 237-rand_frac2
global rand_next, rand_seed
x = rand_next + rand_seed
y = x * x
y = y * y
x = x * y
y = x - math.floor(x)
rand_next = y + x * 0.000005
return rand_next
# ~~~~~~~ 240-print_if_err_cnt_positive
def print_if_err_cnt_positive():
'''Print an error message from a global count.'''
global pow_err_cnt
if pow_err_cnt > 0:
print("Similar discrepancies have occurred {:d} times."
.format(pow_err_cnt))
return
# ~~~~~~~ 250-find_big_two_to_the_nth
# =============================================
# Computations of fundamental system constants like radix and precisioon.
# =============================================
def find_big_2_to_nth():
"""Compute 2**n satisfying |((2**n + 1)-2**n)-1| >= 1."""
x = ONE
while True:
x = x + x
y = x + ONE
z = y - x
y = z - ONE
if MINUS_ONE + fabs(y) >= ZERO: break
# ... now x is just big enough that |((x+1)-x)-1| >= 1 ...
return x
# ~~~~~~~ 255-find_radix_from_2_nth
def find_radix_from_big_2_to_nth(big):
"""The radix is marked by a change in the 'tens' place.
Args:
big: smallest value satisfying |((big + 1) - big) - 1| >= 1
Returns:
radix, which may not be integral in log arithmetic
"""
# ~~~~~~~ 257-find_radix_from_2_nth2
y = ONE
while True:
b = big + y
y = y + y
b = b - big
if b != ZERO: break
return b
# ~~~~~~~ 260-find_precision_big_B_to_nth
def find_precision_big_B_to_nth(b):
"""Compute the number of B-digits in the arithmetic and the power
of B sufficient to have the ones place fall off the right.
Args:
b: the global radix B, accepted as an argument
Returns:
precision: number of B digits in arithmetic
power of B such that the low-order digit is the B's place
"""
# ~~~~~~~ 262-find_precision_big_B_to_nth2
big_b = ONE
precision = ZERO
while True:
precision = precision + ONE
big_b = big_b * b
y = big_b + ONE
if y - big_b != ONE: break
return precision, big_b
# ~~~~~~~ 265-find_ulp_of_one_plus
def find_ulp_of_one_plus(guess):
"""Compute u such at 1+u is the next larger value than 1."""
x = FOUR / THREE # rounds in ulps of 1
third = x - ONE
sixth = ONE_HALF - third
x = sixth + sixth
x = fabs(x - third) # sums two rounding errors
if x < guess:
x = guess
# Now x = (unknown no.) ulps of 1+.
# ~~~~~~~ 267-find_ulp_of_one_plus2
while True:
u = x
x = ONE_HALF * u + THIRTY_TWO * u * u
x = x + ONE
x = x - ONE
if u <= x or x <= ZERO: break
return u
# ~~~~~~~ 270-find_ulp_of_one_minus
def find_ulp_of_one_minus(guess):
"""Compute u such at 1-u is the next smaller value than 1."""
x = TWO / THREE
sixth = x - ONE_HALF
third = sixth + sixth
x = third - ONE_HALF
x = fabs(x + sixth)
if (x < guess):
x = guess
# Now x == (unknown no.) ulps of 1-.
# ~~~~~~~ 272-find_ulp_of_one_minus2
while True:
u = x
y = ONE_HALF * u + THIRTY_TWO * u * u
y = ONE_HALF - y
x = ONE_HALF + y
y = ONE_HALF - x
x = ONE_HALF + y
if u <= x or x <= ZERO: break
return u
# ~~~~~~~ 275-compute_A_and_inverse
def compute_A_and_inverse(radix):
"""Return value for constants A and ONE_OVER_A, based on radix.
a is 2 if radix is 1 (logarithmic) or a power of 2.
a is 10 if radix is a power of 10. Otherwise a is the radix.
There are two tests that "a" behaves well in division.
Args:
radix: passed from B
Returns:
a, 1 / a
"""
# ~~~~~~~ 277-compute_A_and_inverse2
a = max(TWO, radix) # Default, picking up log radix=1
for z in [TWO, NINE + ONE]:
y = radix
while True:
x = y
y = y / z
if floor(y) != y: break
if x == ONE:
a = z # Have found a power of 2 or 10
break
# ~~~~~~~ 279-compute_A_and_inverse3
one_over_a = ONE / a
x = a
for x in [a, B]:
y = ONE / x
z = x * y - ONE_HALF
test_cond(err_failure, z == ONE_HALF, "x * (1/x) differs from 1")
return a, one_over_a
# ~~~~~~~ 280-compute_C_and_inverse
def compute_C_and_inverse():
"""Return untility constant C = 1 / B**k, safely above the underflow
threshold, and its reciprocal.
C is a useful starting point for investigation of underflow.
It has the form 1 / B**k with k chosen so that C is farther than
a factor of 1 / B**PRECISION of underflowing to zero.
C is a normal number. It's reciprocal is a large power of B.
Returns:
C, 1 / C
Basic 4370-4400
"""
# ~~~~~~~ 282-compute_C_and_inverse2
d = ULP_OF_ONE_MINUS # default value for integer precision
if PRECISION != floor(PRECISION):
# Compute 1 / B**ceiling(PRECISION)
d = ONE_OVER_B
for _ in range(floor(PRECISION)):
d = d * ONE_OVER_B
# ~~~~~~~ 283-compute_C_and_inverse3
# d = 1 / B**k for k = PRECISION or floor(PRECISION)+1
# Compute 1/d * 1/d**2 * 1/d**4 * ... until the multiplier
# tapers off (due to underflow). This is the quick route to
# a very tiny 1 / B**k in ever-inceasing steps.
y = ONE
z = d
while True:
c = y
y = z
z = y * y
if y <= z or z + z <= z: break
# ~~~~~~~ 284-compute_C_and_inverse4
# Finally, step toward the minimum normal number by multiples of d.
# c stays one step behind y, so that at termination
# 1 >> c > B**(min+PRECISION) > y > B**mnin > z (underflows)
y = c
z = y * d
while True:
c = y
y = z
z = y * d
if y <= z or z + z <= z: break
return c, ONE / c
# ~~~~~~~ 285-compute_safe_ulps_of_one
def compute_safe_ulps_of_one():
"""Return a value useful in underflow and square root testing.
Given the constant C, a power of B not far above the underflow
threshold, find x such that (1+x)*C > C. Try x = 1 ulp,
but back off to B ulps of one. C is large enough that
underflow is not an issue, but if multiplication
lacks a guard digit, a single ulp might be lost.
Returns:
1 or B ulps of one
Basic: 4450-4460
"""
# ~~~~~~~ 287-compute_safe_ulps_of_one2
x = ULP_OF_ONE_PLUS
s = ONE + x
d = C * s
if d <= C:
# Fall back to B digits. A later test determines
# whether multiplication is just too awful.
x = B * ULP_OF_ONE_PLUS
return x
# ~~~~~~~ 290-compute_H_and_inverse
def compute_H_and_inverse():
"""Return min(1/B, 1/2) and reciprocal, noting H is a fraction."""
if B < TWO:
h_inv = TWO
else:
h_inv = B
return ONE / h_inv, h_inv
# ~~~~~~~ 295-standardize_fractional_precision
def standardize_fractional_precision(ulp, b):
"""Return the precision standardized so that any fraction part
is in units of 1/240. The factors of 240 = 3*5*16 ensures
that fractional octal, decimal, or hex digits will be handled
properly.
Args:
ulp: unit in the last place of a value just less than 1
b: radix
Returns:
sss
"""
# ~~~~~~~ 297-standardize_fractional_precision2
# Intuitively, ulp = b**(-precisioon), so the negative of the
# ratio of the logs retrieves the precision. But there may be
# extra bits in octal or hex or even decimal, making up a partial
# lowest-order digit. The factor 240 is a multiple of 8, 10, and
# 16 with a bit extra. The value y rounds off noise from the logs.
# TODO: resolve the apparent redundancy between the big factor 240
# and the check for being with a quarter.
x = - TWOFORTY * log(ulp) / log(b)
y = floor(ONE_HALF + x)
if (fabs(x - y) * FOUR < ONE):
x = y
p = x / TWOFORTY
y = floor(ONE_HALF + p)
if (fabs(p - y) * TWOFORTY < ONE_HALF):
p = y
return p
# ~~~~~~~ 300-test_for_extra_precise_subepressions
# =============================================
# Test for extra precison.
# =============================================
def test_for_extra_precise_subepressions():
"""Print messages if wider intermediate values detected.
"""
# ~~~~~~~ 302-test_for_extra_precise_subepressionsc2
# (4/3 -1) - 1/4 = 1/12 with rounding error in values bigger than 1.
# 3 * rounded value - 1/4 should be some rounding error in 1+.
# z2 should be an ulp of 1+, unless the computation of x below enjoys
# extra precision, allowing it to be even smaller.
x = fabs(((FOUR / THREE - ONE) - ONE / FOUR) * THREE - ONE / FOUR)
doit = True
while doit:
z2 = x
x = (ONE + (ONE_HALF * z2 + THIRTY_TWO * z2 * z2)) - ONE
doit = (z2 > x) and (x > ZERO)
# ~~~~~~~ 304-test_for_extra_precise_subepressions4
# Now look at ulps of 1- by starting with rounding noise in 3/4.
# x, y, and z are used in the three following loops.
x = y = z = fabs((THREE / FOUR - TWO / THREE) * THREE - ONE / FOUR)
while True:
z1 = z
z = (ONE / TWO
- ((ONE / TWO - (ONE_HALF * z1 + THIRTY_TWO * z1 * z1))
+ ONE / TWO)) + ONE / TWO
if z1 <= z or z <= ZERO: break
# ~~~~~~~ 306-test_for_extra_precise_subepressions6
# BUG in BASIC and C version -- the two loops below are not nested,
# but simply in sequence with the z1 loop just above.
# See lines 1730 and 1740 in the BASIC, the latter line should branch
# back to itself.
while True:
y1 = y
y = (ONE_HALF - ((ONE_HALF - (ONE_HALF * y1 + THIRTY_TWO * y1 * y1))
+ ONE_HALF)) + ONE_HALF
print("1730: y1 = {:.17e} y = {:.17e}".format(y1, y))
if y1 <= y or y <= ZERO: break
while True:
x1 = x
x = ((ONE_HALF * x1 + THIRTY_TWO * x1 * x1)
- ONE_MINUS_ULP) + ONE_MINUS_ULP
print("1740: x1 = {:.17e} x = {:.17e}".format(x1, x))
if x1 <= x or x <= ZERO: break
# ~~~~~~~ 308-test_for_extra_precise_subepressions8
if x1 != y1 or x1 != z1:
bad_cond(err_serious, "Disagreements among the values x1, y1, z1,")
print("respectively {:0.7e}, {:0.7e}, {:0.7e},".format(x1, y1, z1))
print("are symptoms of inconsistencies introduced")
print("by extra-precise evaluation of arithmetic subexpressions.")
notify("Possibly some part of this")
if (x1 == ULP_OF_ONE_MINUS or y1 == ULP_OF_ONE_MINUS
or (z1 == ULP_OF_ONE_MINUS)):
print("That feature is not tested further by this program.\n")
# ~~~~~~~ 310-test_for_extra_precise_subepressions10
elif z1 != ULP_OF_ONE_MINUS or z2 != ULP_OF_ONE_PLUS:
if z1 >= ULP_OF_ONE_MINUS or z2 >= ULP_OF_ONE_PLUS:
bad_cond(err_failure, "")
notify("PRECISION")
print("\tULP_OF_ONE_MINUS = {:0.7e}, z1 - ULP_OF_ONE_MINUS = {:0.7e}".format(ULP_OF_ONE_MINUS,
z1 - ULP_OF_ONE_MINUS))
print("\tULP_OF_ONE_PLUS = {:0.7e}, z2 - ULP_OF_ONE_PLUS = {:0.7e}".format(ULP_OF_ONE_PLUS,
z2 - ULP_OF_ONE_PLUS))
# ~~~~~~~ 312-test_for_extra_precise_subepressions12
else:
if (z1 <= ZERO or z2 <= ZERO):
print("Because of unusual B = {:0.f}".format(B))
print(", or exact rational arithmetic a result")
print("z1 = {:0.7e}, or z2 = {:0.7e} ".format(z1, z2))
notify("of an extra-precision")
# ~~~~~~~ 314-test_for_extra_precise_subepressions14
if z1 != z2 or z1 > ZERO:
x = z1 / ULP_OF_ONE_MINUS
y = z2 / ULP_OF_ONE_PLUS
if y > x:
x = y
q = - log(x)
print("Some subexpressions appear to be calculated extra")
print("precisely with about {:.g} extra B-digits, i.e.",
(q / log(B)))
print("roughly {:.g} extra significant decimals."
.format(q / log(10.)))
print("That feature is not tested further by this program.")
# ~~~~~~~ 320-test_normalized_subtraction
# =============================================
# Tests for normalization, guard digit, rounding, use of a sticky bit.
# =============================================
def test_normalized_subtraction():
"""Test for renormalization on magnitude subtraction.
Computes 1.0 via cancellation of two large values. Then adds
a tiny increment that would be rounded off were the value 1.0
not normalized after the subtraction.
"""
# ~~~~~~~ 322-test_normalized_subtraction2
if (B >= TWO):
x = BIG_B_NTH / (B * B)
y = x + ONE
z = y - x
t = z + ULP_OF_ONE_PLUS
x = t - z
if x == ULP_OF_ONE_PLUS:
print("Subtraction appears to be normalized, as it should be.")
else:
test_cond(err_failure, False,
"Subtraction is not normalized x=y,x+z != y+z!")
return
# ~~~~~~~ 325-does_mult_have_guard_digit
def does_mult_have_guard_digit():
"""Test for a guard digit in mulitiplicaiton.
Run three tests, 1*X and X*1 for X with nonzero low-order
digit, and a product known to require one extra digit, before
normalization.
Returns:
boolean answer to the question
"""
# ~~~~~~~ 327-does_mult_have_guard_digit2
print("\nChecking for guard digit in *, /, and -.")
y = ONE_MINUS_ULP * ONE
z = ONE * ONE_MINUS_ULP
x = ONE_MINUS_ULP - ONE_HALF
y = (y - ONE_HALF) - x
z = (z - ONE_HALF) - x
x = ONE + ULP_OF_ONE_PLUS
t = x * B
r = B * x
x = t - B
x = x - B * ULP_OF_ONE_PLUS
t = r - B
t = t - B * ULP_OF_ONE_PLUS
x = x * (B - ONE)
t = t * (B - ONE)
gmult = (x == ZERO) and (y == ZERO) and (z == ZERO) and (t == ZERO)
test_cond(err_serious, gmult, "* lacks a Guard Digit, so 1*x != x")
return gmult
# ~~~~~~~ 330-test_mult_low_digits
def test_mult_low_digits():
"""Try two tests of the form (1 + tiny)**2."""
z = B * ULP_OF_ONE_PLUS
x = ONE + z
y = fabs((x + z) - x * x) - ULP_OF_ONE_PLUS
# ~~~~~~~ 332-mullowdig2
x = ONE - ULP_OF_ONE_PLUS
z = fabs((x - ULP_OF_ONE_PLUS) - x * x) - ULP_OF_ONE_MINUS
test_cond(err_failure, y <= ZERO and z <= ZERO,
"* gets too many final digits wrong.\n")
# ~~~~~~~ 335-does_div_have_guard_digit
def does_div_have_guard_digit():
"""Check for guard digit in division.
Returns:
boolean reply to question
"""
y = ONE - ULP_OF_ONE_PLUS
x = ONE + ULP_OF_ONE_PLUS
z = ONE / y
y = z - x
# ~~~~~~~ 337-does_div_have_guard_digit2
x = ONE / THREE
z = THREE / NINE
x = x - z
# ~~~~~~~ 339-does_div_have_guard_digit3
t = NINE / TWENTY_SEVEN
z = z - t
gdiv = (x == ZERO) and (y == ZERO) and (z == ZERO)
test_cond(err_defect, gdiv,
"Division lacks a Guard Digit, so error can exceed 1 ulp\n"
+ "or 1/3 and 3/9 and 9/27 may disagree")
# ~~~~~~~ 341-does_div_have_guard_digit6
y = ONE_MINUS_ULP / ONE
x = ONE_MINUS_ULP - ONE_HALF
y = (y - ONE_HALF) - x
x = ONE + ULP_OF_ONE_PLUS
t = x / ONE
x = t - x
# ~~~~~~~ 343-does_div_have_guard_digit7
# JTC: gdiv now includes all 5 tests not just last 4 like original.
gdiv &= (x == ZERO) and (y == ZERO)
test_cond(err_serious, gdiv, "Division lacks a Guard Digit, so x/1 != x")
return gdiv
# ~~~~~~~ 345-does_add_sub_have_guard_digit
def does_add_sub_have_guard_digit():
"""Check for guard digit in add and subtract.
Returns:
boolean reply to question
"""
x = ONE / (ONE + ULP_OF_ONE_PLUS)
y = x - ONE_HALF - ONE_HALF
test_cond(err_serious, y < ZERO, "Computed value of 1/1.000..1 >= 1")
# ~~~~~~~ 347-does_add_sub_have_guard_digit2
x = ONE - ULP_OF_ONE_PLUS
y = ONE + B * ULP_OF_ONE_PLUS
z = x * B
t = y * B
r = z / B
s = t / B
x = r - x
y = s - y
test_cond(err_failure, x == ZERO and y == ZERO,
"* and/or / gets too many last digits wrong")
# ~~~~~~~ 349-does_add_sub_have_guard_digit3
y = ONE - ULP_OF_ONE_MINUS
x = ONE - ONE_MINUS_ULP
y = ONE - y
t = B - ULP_OF_ONE_PLUS
z = B - B_MINUS_ULP
t = B - t
gaddsub = ((x == ULP_OF_ONE_MINUS) and (y == ULP_OF_ONE_MINUS)
and (z == ULP_OF_ONE_PLUS) and (t == ULP_OF_ONE_PLUS))
test_cond(err_serious, gaddsub,
"- lacks Guard Digit, so cancellation is obscured")
return gaddsub
# ~~~~~~~ 350-how_does_mult_round
def how_does_mult_round():
"""Check type of rounding in multiplication.
Returns:
"rounded", "chopped", "other"
"""
# ~~~~~~~ 352-how_does_mult_round2
rmult = "other"
y2 = ONE + ULP_OF_ONE_PLUS
y1 = ONE - ULP_OF_ONE_PLUS
x = ONE_AND_HALF - ULP_OF_ONE_PLUS
y = ONE_AND_HALF + ULP_OF_ONE_PLUS
z = (x - ULP_OF_ONE_PLUS) * y2
t = y * y1
z = z - x
t = t - x
# ~~~~~~~ 354-how_does_mult_round4
x = x * y2
y = (y + ULP_OF_ONE_PLUS) * y1
x = x - ONE_AND_HALF
y = y - ONE_AND_HALF
# ~~~~~~~ 356-how_does_mult_round6
if ((x == ZERO) and (y == ZERO) and (z == ZERO) and (t <= ZERO)):
x = (ONE_AND_HALF + ULP_OF_ONE_PLUS) * y2
y = ONE_AND_HALF - ULP_OF_ONE_PLUS - ULP_OF_ONE_PLUS
z = ONE_AND_HALF + ULP_OF_ONE_PLUS + ULP_OF_ONE_PLUS
t = (ONE_AND_HALF - ULP_OF_ONE_PLUS) * y1
# ~~~~~~~ 358-how_does_mult_round8
x = x - (z + ULP_OF_ONE_PLUS)
s1 = y * y1
s = z * y2
t = t - y
y = (ULP_OF_ONE_PLUS - y) + s1
z = s - (z + ULP_OF_ONE_PLUS + ULP_OF_ONE_PLUS)
# ~~~~~~~ 360-how_does_mult_round10
s1 = (y2 + ULP_OF_ONE_PLUS) * y1
y1 = y2 * y1
s1 = s1 - y2
y1 = y1 - ONE_HALF
# ~~~~~~~ 362-how_does_mult_round12
if ((x == ZERO) and (y == ZERO) and (z == ZERO) and (t == ZERO)
and (s1 == ZERO) and (y1 == ONE_HALF)):
rmult = "rounded"
print("Multiplication appears to round correctly. ")
# ~~~~~~~ 364-how_does_mult_round14
elif ((x + ULP_OF_ONE_PLUS == ZERO) and (y < ZERO) and (z + ULP_OF_ONE_PLUS == ZERO)
and (t < ZERO) and (s1 + ULP_OF_ONE_PLUS == ZERO)
and (y1 < ONE_HALF)):
rmult = "chopped"
print("Multiplication appears to chop. ")
# ~~~~~~~ 366-how_does_mult_round16
else:
print("* is neither chopped nor correctly rounded. ")
if (rmult == "rounded"
and not flags["mult_guard_digit"]):
notify("Multiplication")
# ~~~~~~~ 368-how_does_mult_round18
else:
print("* is neither chopped nor correctly rounded. ")
return rmult
# ~~~~~~~ 370-how_does_div_round
def how_does_div_round():
"""Check for rounding in division.
Returns:
"rounded", "chopped", "other"
"""
# ~~~~~~~ 372-how_does_div_round2
rdiv = "other"
y2 = ONE + ULP_OF_ONE_PLUS
y1 = ONE - ULP_OF_ONE_PLUS
z = ONE_AND_HALF + ULP_OF_ONE_PLUS + ULP_OF_ONE_PLUS
x = z / y2
t = ONE_AND_HALF - ULP_OF_ONE_PLUS - ULP_OF_ONE_PLUS
y = (t - ULP_OF_ONE_PLUS) / y1
z = (z + ULP_OF_ONE_PLUS) / y2
# ~~~~~~~ 374-how_does_div_round4
x = x - ONE_AND_HALF
y = y - t
t = t / y1
z = z - (ONE_AND_HALF + ULP_OF_ONE_PLUS)
t = (ULP_OF_ONE_PLUS - ONE_AND_HALF) + t
# ~~~~~~~ 376-how_does_div_round6
if not ((x > ZERO) or (y > ZERO) or z > ZERO or t > ZERO):
x = ONE_AND_HALF / y2
y = ONE_AND_HALF - ULP_OF_ONE_PLUS
z = ONE_AND_HALF + ULP_OF_ONE_PLUS
x = x - y
t = ONE_AND_HALF / y1
y = y / y1
t = t - (z + ULP_OF_ONE_PLUS)
# ~~~~~~~ 378-how_does_div_round8
y = y - z
z = z / y2
y1 = (y2 + ULP_OF_ONE_PLUS) / y2
z = z - ONE_AND_HALF
y2 = y1 - y2
y1 = (ONE_MINUS_ULP - ULP_OF_ONE_MINUS) / ONE_MINUS_ULP
# ~~~~~~~ 380-how_does_div_round10
if ((x == ZERO) and (y == ZERO) and (z == ZERO) and (t == ZERO)
and (y2 == ZERO) # BUG in C version: duplicate test
and (y1 - ONE_HALF == ONE_MINUS_ULP - ONE_HALF)):
rdiv = "rounded"
print("Division appears to round correctly. ")
if not flags["div_guard_digit"]:
notify("Division")
# ~~~~~~~ 382-how_does_div_round12
elif ((x < ZERO) and (y < ZERO) and (z < ZERO) and (t < ZERO)
and (y2 < ZERO) and (y1 - ONE_HALF < ONE_MINUS_ULP - ONE_HALF)):
rdiv = "chopped"
print("Division appears to chop. ")
# ~~~~~~~ 384-how_does_div_round15
if rdiv == "other":
print("/ is neither chopped nor correctly rounded. ")
return rdiv
# ~~~~~~~ 385-how_does_add_sub_round
def how_does_add_sub_round():
"""Check rounding on add/subtract.
Returns:
"rounded", "chopped", "other"
"""
# ~~~~~~~ 387-how_does_add_sub_round2
round_add_sub = "other"
x = ONE - ULP_OF_ONE_MINUS * ULP_OF_ONE_MINUS
y = ONE + ULP_OF_ONE_PLUS * (ONE - ULP_OF_ONE_PLUS)
z = ONE_MINUS_ULP - ONE_HALF
x = (x - ONE_HALF) - z
y = y - ONE
# ~~~~~~~ 389-how_does_add_sub_round4
if (x == ZERO) and (y == ZERO):
round_add_sub = "chopped"
print("Add/Subtract appears to be chopped. ")
# ~~~~~~~ 391-how_does_add_sub_round6
while flags["add_sub_guard_digit"]: # JTC TODO jibe with Basic
x = (ONE_HALF + ULP_OF_ONE_PLUS) * ULP_OF_ONE_PLUS
y = (ONE_HALF - ULP_OF_ONE_PLUS) * ULP_OF_ONE_PLUS
x = ONE + x
y = ONE + y
x = (ONE + ULP_OF_ONE_PLUS) - x
y = ONE - y
if x != ZERO or y != ZERO: break
# ~~~~~~~ 393-how_does_add_sub_round8
x = (ONE_HALF + ULP_OF_ONE_PLUS) * ULP_OF_ONE_MINUS
y = (ONE_HALF - ULP_OF_ONE_PLUS) * ULP_OF_ONE_MINUS
x = ONE - x
y = ONE - y
x = ONE_MINUS_ULP - x
y = ONE - y
if x != ZERO or y != ZERO: break
# ~~~~~~~ 395-how_does_add_sub_round10
if round_add_sub == "chopped":
# ***JTC: BASIC 2690 sets R3 = F1 - O2*R3 triggering 3
# messages when boch chopped and rounded are detected.
round_add_sub = "other"
else:
round_add_sub = "rounded"
print("Addition/Subtraction appears to round correctly. ")
# ~~~~~~~ 397-how_does_add_sub_round12
if not flags["add_sub_guard_digit"]:
# ***JTC: BUG, annot happen, guard digit is checked above.
notify("Add/Subtract")
break # Execute this block just once, with early exit possible.
# ~~~~~~~ 399-how_does_add_sub_round14
if round_add_sub == "other": # pick up all the odd cases
print("Addition/Subtraction neither rounds nor chops.")
return round_add_sub
# ~~~~~~~ 400-rounding_symmetry_test
def rounding_symmetry_test():
"""Test sign symmetry in simple rounding case.
Returns:
True if passes, False otherwise
"""
# ~~~~~~~ 402-rounding_symmetry_test2
verdict = True
x = ONE + ONE_HALF * (ONE + ONE_HALF)
y = (ONE + ULP_OF_ONE_PLUS) * ONE_HALF
z = x - y
t = y - x
s = z + t
if (s != ZERO):
verdict = False
bad_cond(err_flaw, "(x - y) + (y - x) is non zero!\n")
return verdict
# ~~~~~~~ 405-has_sticky_bit
def has_sticky_bit():
"""Create a rounding situation designed to require a sticky bit.
Returns:
True if sticky detected, False otherwise
"""
# ~~~~~~~ 407-has_sticky_bit2
verdict = False
x = (ONE_HALF + ULP_OF_ONE_MINUS) * ULP_OF_ONE_PLUS
y = ONE_HALF * ULP_OF_ONE_PLUS
z = ONE + y
t = ONE + x
# ~~~~~~~ 409-has_sticky_bit4
if (z - ONE <= ZERO) and (t - ONE >= ULP_OF_ONE_PLUS):
z = t + y
y = z - x
# ~~~~~~~ 411-has_sticky_bit6
if (z - t >= ULP_OF_ONE_PLUS) and (y - t == ZERO):
x = (ONE_HALF + ULP_OF_ONE_MINUS) * ULP_OF_ONE_MINUS
y = ONE_HALF * ULP_OF_ONE_MINUS
z = ONE - y
t = ONE - x
# ~~~~~~~ 413-has_sticky_bit8
if (z - ONE == ZERO) and (t - ONE_MINUS_ULP == ZERO):
z = (ONE_HALF - ULP_OF_ONE_MINUS) * ULP_OF_ONE_MINUS
t = ONE_MINUS_ULP - z
q = ONE_MINUS_ULP - y
# ~~~~~~~ 415-has_sticky_bit10
if ((t - ONE_MINUS_ULP == ZERO)
and (ONE_MINUS_ULP - ULP_OF_ONE_MINUS - q == ZERO)):
z = (ONE + ULP_OF_ONE_PLUS) * ONE_AND_HALF
t = (ONE_AND_HALF + ULP_OF_ONE_PLUS) - z + ULP_OF_ONE_PLUS
x = ONE + ONE_HALF / B
y = ONE + B * ULP_OF_ONE_PLUS
z = x * y
# ~~~~~~~ 417-has_sticky_bit12
if (t == ZERO) and (x + B * ULP_OF_ONE_PLUS - z == ZERO):
if (B != TWO):
x = TWO + ULP_OF_ONE_PLUS
y = x / TWO
if (y - ONE) == ZERO:
verdict = True
else:
verdict = True
# ~~~~~~~ 419-has_sticky_bit14
return verdict
# ~~~~~~~ 425-test_mult_commutivity
def test_mult_commutivity(count):
"""Run some random tests to check that x*y = y*x. """
global rand_next, rand_seed
rand_seed = sqrt(3.0)
rand_next = ONE / THREE
# ~~~~~~~ 425a-test_mult_commutivity2
i = 0
doit = True
while doit:
x = rand_frac()
y = rand_frac()
yx = y * x
xy = x * y
xyyx = xy - yx
i = i + 1
doit = (i < count) and (xyyx == ZERO)
# ~~~~~~~ 425b-test_mult_commutivity2
if (xyyx == ZERO):
# One rather less random test lays out the operands, products,
# and the computed difference.
x = ONE + ONE_HALF / THREE
y = (ULP_OF_ONE_PLUS + ULP_OF_ONE_MINUS) + ONE
xy = x * y
yx = y * x
xyyx = ((ONE + ONE_HALF / THREE)
* ((ULP_OF_ONE_PLUS + ULP_OF_ONE_MINUS) + ONE)
- ((ULP_OF_ONE_PLUS + ULP_OF_ONE_MINUS) + ONE)
* (ONE + ONE_HALF / THREE))
# ~~~~~~~ 427-test_mult_commutivity2
# In the usual case, count random tests and 1 extra test are run.
# xyyx is the final arbiter of success.
if xyyx != ZERO:
print("x * y == y * x fails for x = {:0.17e} y = {:0.17e}"
.format(x, y))
print(" x * y = {:0.17e} y * x = {:0.17e}".format(xy, yx))
print(" and x * y - y * x = {:0.17e}".format(xyyx))
bad_cond(err_defect, "x * y == y * x trial fails.")
else:
print(" No failures found in {:d} integer pairs.".format(count + 1))
# ~~~~~~~ 430-sqrt_single_bound_test
# =============================================
# Tests of square root accuracy and monotonicity.
# =============================================
def sqrt_single_bound_test(tlist):
"""Test x = sqrt(x*x), raise error if not, and track max and min
signed error in units of ulp.
Uses an artful decomposition of x into xa + xb to facilitate test for
equality.
***TODO: elaborate that this decomposition is exact -- because the values
of x all have trailing 0s?
Args:
tlist - a list of tuples of:
x - value to be tested
ulp - 1/2 ulp of argument x
error_kind - level of alarm to raise
Returns:
True if all tests pass, False otherwise
"""
# ~~~~~~~ 432-sqrt_single_bound_test2
global min_sqrt_error, max_sqrt_error
success = True
for (x, ulp, error_kind) in tlist:
xb = x * ONE_OVER_B
xa = x - xb
err = ((sqrt(x * x) - xb) - xa) / ulp
# ~~~~~~~ 434-sqrt_single_bound_test4
if err != ZERO:
success = False
if err < min_sqrt_error:
min_sqrt_error = err
if err > max_sqrt_error:
max_sqrt_error = err
bad_cond(error_kind, "\n")
print("sqrt( {:0.17e}) - {:0.17e} = {:0.17e}".format(x * x, x, ulp * err))
print("\tinstead of correct value 0 .")
return success
# ~~~~~~~ 433-test_sqrt_int
def test_sqrt_int():
"""Test sqrt(x**2) for consecutive integer in a B-ade."""
x = TWO # Default start point for log arithmetic
y = B
if B != ONE: # ***TODO: log arithmetic issue
while True:
x = y
y = B * y
if y - x >= NUM_TRIALS: break
ulp = x * ULP_OF_ONE_PLUS # ***TODO -- why was Py _MINUS, with others _PLUS?
for _ in range(NUM_TRIALS):
x = x + ONE # when NUM_TRIALS is size of a binade, x hits higher B**k
if not sqrt_single_bound_test([(x, ulp, err_defect)]): break
return
# ~~~~~~~ 435-test_sqrt_monotonicity
def test_sqrt_monotonicity():
"""Test for monotonicity around the radix B."""
i = - 1
x = B_MINUS_ULP
y = B
z = B + B * ULP_OF_ONE_PLUS # B plus 1 ulp
not_monot = False
monot = False
# ~~~~~~~ 435a-test_sqrt_monotonicity2
while not (not_monot or monot):
i = i + 1
# CASE: y = B
x = sqrt(x)
q = sqrt(y)
z = sqrt(z)
if x > q or q > z:
not_monot = True
# ~~~~~~~ 435b-test_sqrt_monotonicity4
else:
q = floor(q + ONE_HALF)
# B = q*q only if radix B is 16 or 4
if (i <= 0) and (B != q * q):
monot = True # exit early unless B is 16 or 4
# ~~~~~~~ 435c-test_sqrt_monotonicity6
elif i <= 0:
# CASE: y = sqrt(B), if B = 4, 16 or other square
y = q
x = y - ULP_OF_ONE_PLUS
z = y + ULP_OF_ONE_PLUS
# ~~~~~~~ 435d-test_sqrt_monotonicity8
elif i == 1:
# CASE: 7 = 1/sqrt(B), if B = 4, 16, or other square
y = y * ONE_OVER_B
x = y - ULP_OF_ONE_MINUS
z = y + ULP_OF_ONE_MINUS
else:
# Have survived three cases
monot = True
# ~~~~~~~ 435e-test_sqrt_monotonicity10
if monot:
print("sqrt has passed a test for Monotonicity.")
else:
bad_cond(err_defect, "")
print("sqrt(x) is non-monotonic for x near {:0.7e} .".format(y))
# ~~~~~~~ 440-update_sqrt_errors
def update_sqrt_errors(down, up):
"""Update the sqrt global min and max errors."""
global min_sqrt_error, max_sqrt_error
if up > max_sqrt_error:
max_sqrt_error = up
if down < min_sqrt_error:
min_sqrt_error = down
return
# ~~~~~~~ 446-num_theory_sqrt_test
def num_theory_sqrt_test(x, D, y, z2):
"""Test utility."""
# ~~~~~~~ 446a-num_theory_sqrt_test1
test = 0
if x - B >= z2 - B and x - z2 <= BIG_B_NTH - z2:
test = 1
if ultra_verbose:
print("\t\tnum_theory_sqrt_test x = {:.6f} D = {:.6f}".format(x, D))
X2 = sqrt(x * D)
if ultra_verbose:
print("\t\t X2 = {:.6f} = sqrt( {:.6f} )"
.format(X2, x * D))
y2 = (X2 - z2) - (y - z2)
X2 = X8 / (y - ONE_HALF)
X2 = X2 - ONE_HALF * X2 * X2
# ~~~~~~~ 446b-num_theory_sqrt_test2
down_err = (y2 + ONE_HALF) + (ONE_HALF - X2)
up_err = y2 - X2
if ultra_verbose:
print("\t\t y2 = {:.5e} X2 = {:.5e}".format(y2, X2))
print("\t\t dn = {:.5e} up = {:.5e}"
.format(down_err, up_err))
update_sqrt_errors(down_err, up_err)
if ultra_verbose:
print("\tnew X2 = {:.17e} y2 = {:.17e}".format(X2, y2))
elif verbose:
print("\thasty exit from num_theory_sqrt_test()")
return test
# ~~~~~~~ 448-num_theory_sqrt_value
def num_theory_sqrt_value(x, D, Q, z, z1):
"""Test utility."""
x = z1 * Q
x = floor(ONE_HALF - x / B) * B + x
Q = (Q - x * z) / B + x * x * (D / B)
z = z - TWO * x * D
# ~~~~~~~ 448a-num_theory_sqrt_value1
if z <= ZERO:
if ultra_verbose:
print("\t\t\tnum_theory_sqrt_value: z = {:.17} <= ZERO".format(z))
z = - z
z1 = - z1
# ~~~~~~~ 448b-num_theory_sqrt_value2
D = B * D
if ultra_verbose:
print("\tnew D = {:.17e} Q = {:.17e}".format(D, Q))
return x, D, Q, z, z1
# ~~~~~~~ 450-test_power_result
# =============================================
# Tests of X**K for integer X and K.
# =============================================
def test_power_result(x, z, q):
"""Check that x equals z**q, usually for exact integer cases."""
global pow_err_cnt
y = pow(z, q)
if ultra_verbose:
print("\t{:.17e} = {:.17e} ^ {:.17e}".format(y, z, q))
if y != x:
# ~~~~~~~ 452-test_power_result2
if pow_err_cnt <= 0:
if z == ZERO and q <= ZERO:
print("WARNING: computing")
else:
bad_cond(err_defect, "computing")
print("\t({:0.17e}) ^ ({:0.17e})".format(z, q))
print("\tyielded {:0.17e};\n".format(y))
print("\twhich compared unequal to correct {:0.17e} ;"
.format(x))
print("\t\tthey differ by {:0.17e} .".format(y - x))
pow_err_cnt += 1 # ... count discrepancies.
return
# ~~~~~~~ 455-test_integer_powers
def test_integer_powers(x, z, i, max_i):
"""Test z**k for k = i .. max_i."""
while True:
q = float(i)
if ultra_verbose:
print("test_int_pow with x, z, q = {:.17e} {:.17e} {:.4f}"
.format(x, z, q))
test_power_result(x, z, q)
# ~~~~~~~ 457-test_integer_powers2
i += 1
x = z * x
# Stay within the range of exactly representable integers.
if x >= BIG_B_NTH or i > max_i: break
return
# ~~~~~~~ 520-tiny_powers_of_B
# =============================================
# Tests for underflow behavior.
# =============================================
def tiny_powers_of_B(c):
"""Return tiny values, stepping toward zero by factors of radix or 2.
Starting with a value like C = 1/B**k close to the underflow
threshold, walk down by factors of H = min(1/B, 1/2). Compute
less_tiny, less_tiny, and z differing by factors of H.
Value less_tiny avoids anomalies. Value z typically falls to
zero. The factor H protects against B=1.
Args:
c: 1/B**k, expected to be the constant C or comparable
Returns:
less_tiny: next tiniest positive value, by a factor of H
tiny: tiniest positive value
z: too tiny value, that underflows to zero or a pseudo-zero
Basic 4430-4440
"""
# ~~~~~~~ 521-tiny_powers_of_B2
tiny = c
z = tiny * H
while True:
less_tiny = tiny
tiny = z
z = z * H
# Terminate if z pins at a nonzero minimal value or
# if z falls to an epsilon or 0.
if tiny <= z or z + z <= z: break
print("***Approaching underflow by powers of the radix:")
print("*** less_tiny = {:.7e}, tiny = {:.7e},"
+ "z = {:.7e}".format(less_tiny, tiny, z))
return less_tiny, tiny, z
# ~~~~~~~ 522-tiny_values_and_difference
def tiny_values_and_difference(one_plus):
"""Compute the tiniest values from a starting value with nonzero
low-order bits. This function is parallel to tiny_powers_of_B(),
but looks for nonzero differences of tiny values, which are
evidence of gradual underflow.
If multiplication is too inaccurate, raise an error and return default
values, shown in parens under Returns.
Args:
one_plus: 1 plus enough ulps that they shouldn't be lost when
multiplied by a power of B
# ~~~~~~~ 523-tiny_values_and_difference2
Returns:
tiny: smallest well-behaved value (0) *** why not tiny_B?
tinier: first value to exhibit underflow (tiny_B)
too_tiny: underflows to 0 or an epsilon value (too_tiny_B)
u_threshold: smallest well-behaved value, if underflow degrades
into nonzero values, or 0 if underflow flushes to 0 (tiny_B)
tiny_delta: difference between two values caused by underflow (0)
save_tiny1: value of tiny where tiny_delta computed (0)
save_tiny2: recomputed tiny where tiny_delta computed (0)
Basic 4450-4520
"""
# ~~~~~~~ 524-tiny_values_and_difference3
d = C * one_plus
# 1 + SAFE_ULPS_OF_ONE should be enough ulps of 1 that d > C.
if d <= C:
# But if not, we punt with default values.
bad_cond(err_failure,
"multiplication gets too many last digits wrong.")
pause()
return 0, tiny_B, too_tiny_B, tiny_B, 0, 0, 0
# ~~~~~~~ 525-tiny_values_and_difference4
else:
tinier = d
too_tiny = tinier * H
u_threshold = ZERO
tiny_delta = ZERO # indicates no value computed yet
save_tiny1 = ZERO
save_tiny2 = ZERO
while True:
tiny = tinier
tinier = too_tiny
# This loop mirrors the computation in tiny_powers_of_B(),
# except that the launch point is C*s with nonzero low-order
# bits, versus C, which is just a power of the radix.
if (tiny_delta + tiny_delta) <= tiny_delta:
if ultra_verbose:
print("\tUnderflow loop: within if, tinier = {:.7e}"
.format(tinier))
save_tiny2 = tinier * ONE_OVER_H
tiny_delta = fabs(tiny - save_tiny2) # may underflow to 0
save_tiny1 = tiny
if (u_threshold == ZERO) and (tiny != save_tiny2):
# Test for 0 captures first time through, in case
# tiny_delta underflows to 0.
u_threshold = tiny
# Placing this statement after the (tiny*H)/H test above
# emphasizes that only nonzero values of tinier qualify for
# values degraded by underflow.
too_tiny = too_tiny * H
if ultra_verbose:
print("\tUnderflow loop: too_tiny = {:.7e}".format(too_tiny))
if (tinier <= too_tiny
or too_tiny + too_tiny <= too_tiny): break
return (tiny, tinier, too_tiny, u_threshold,
tiny_delta, save_tiny1, save_tiny2)
# ~~~~~~~ 526-does_tiny_value_misbehave
def does_tiny_value_misbehave(z):
"""Check for abnormal behavior of a value presumed
to be tiny or zero.
Given a value z, check that (z + z)/z is 2.0 or very close, and that
1.0*z = z*1.0 = z/1.0 = z. This is especially interesting near the
underflow threshold, but the function should work unless z is so
huge that doubling it overflows.
Args:
z - test value
Returns:
boolean reply to the question
Basic 4640-4680
"""
# ~~~~~~~ 527-does_tiny_value_misbehave2
reply = False
if z == ZERO: # nothing complex to test in the easy case
return reply
print("Since comparison denies z = 0, evaluating ", end="")
print("(z + z) / z should be safe.")
z_quo = ONE
try:
z_quo = (z + z) / z
except ZeroDivisionError:
reply = True
error_count[err_serious] = error_count[err_serious] + 1
print("But the division triggered an exception.")
print("This is a VERY SERIOUS DEFECT!")
print("What the machine gets for (z + z) / z is {:0.17e} ."
.format(z_quo))
if fabs(z_quo - TWO) < B * ULP_OF_ONE_PLUS:
print("This is O.K., provided Over/Underflow", end="")
print(" has NOT just been signaled.")
else:
reply = True
if z_quo < ONE or z_quo > TWO:
error_count[err_serious] = error_count[err_serious] + 1
print("This is a VERY SERIOUS DEFECT!")
else:
error_count[err_defect] = error_count[err_defect] + 1
print("This is a DEFECT!")
# ~~~~~~~ 528-does_tiny_value_misbehave3
r1 = z * ONE
r2 = ONE * z
r3 = z / ONE
if (z == r1) and (z == r2) and (z == r3):
if reply:
pause()
else:
reply = True
bad_cond(err_defect, "What prints as z = ", end="")
print("{:0.17e}\n\tcompares different from ".format(z))
if z != r1:
print("z * 1 = {:0.17e} ".format(r1))
if (z != r2) and (r2 != r1):
print("1 * z == {:0.17e}".format(r2))
if z != r3:
print("z / 1 = {:0.17e}".format(r3))
if r2 != r1:
error_count[err_defect] = error_count[err_defect] + 1
bad_cond(err_defect, "Multiplication does not commute!\n")
print("\tComparison alleges that 1 * z = {:0.17e}"
.format(r2))
print("\tdiffers from z * 1 = {:0.17e}".format(r1))
pause()
return reply
# ~~~~~~~ 529-is_gradual_underflow_IEEE
def is_gradual_underflow_IEEE(t):
"""Test gradual underflow for IEEE compatibility.
Args:
t - tiniest nonzero number discovered
Returns:
boolean reply to question
Basic 5150
"""
# ~~~~~~~ 530-is_gradual_underflow_IEEE2
print("Underflow is gradual; it incurs Absolute Error =")
print("(roundoff in underflow_threshold) < min_positive.")
# Test for IEEE rounding at the very bottom of the range.
# The test looks for a double rounding, first to a
# higher precision, then the target. The computation is just
# under half an ulp, but an extra rounding will push it to a
# half-way case that will round up.
# ~~~~~~~ 531-is_gradual_underflow_IEEE2
# y = t/C * (3/2 + u)
# x = 1/C * (1 + u)
# y/x = t * (3/2 + u) * (1 - u + u**2 - ...)
# = t * (1 + 1/2 - u/2 + u**2/2 - ...)
# = t, barely rounding down
y = t * ONE_OVER_C
y = y * (ONE_AND_HALF + ULP_OF_ONE_PLUS)
x = ONE_OVER_C * (ONE + ULP_OF_ONE_PLUS)
y = y / x
return (y == t)
# ~~~~~~~ 532-test_tiny_differences
def test_tiny_differences():
""""With no gradual underflow, look for x != z but x-z = 0.
This function is called after cases 1, 2, and 3 of underflow
triage, where case 4 captures gradual underflow. The idea is
to show how the nonzero differences between tiny numbers are
not representable in the arithmetic.
Usually, base_tiny == underflow_threshold but if it's smaller,
then r * underflow_threshold is the geometric mean of the two.
In the typical non-gradual-underflow cases when tiny_x and
underflow_threshold differ by a factor of 1 or H, the else
clause below is taken.
Basic 5000-5080
"""
# ~~~~~~~ 533-test_tiny_differences2
print("")
try:
r = sqrt(base_tiny / underflow_threshold)
except ZeroDivisionError:
print("base_tiny / underflow_threshold failed!")
r = H + H # fall through to else case below
if r <= H:
z = r * underflow_threshold
x = z * (ONE + r * H * (ONE + H))
else:
z = underflow_threshold
x = z * (ONE + H * H * (ONE + H))
# ~~~~~~~ 533a-test_tiny_differences3
if (x != z) and (x - z == ZERO):
bad_cond(err_flaw, "")
print("x = {:0.17e}\n\tis not equal to z = {:0.17e} ."
.format(x, z))
w = x - z
print("yet x - z yields {:0.17e} .".format(w))
print(" Should this NOT signal Underflow, ")
print("this is a SERIOUS DEFECT\nthat causes ")
print("confusion when innocent statements like");
print(" if (x == z) ... else")
print(" ... (f(x) - f(z)) / (x - z) ...")
print("encounter Division by ZERO although actually")
try:
print("x / z = 1 + {:g} .".format((x / z - ONE_HALF) - ONE_HALF))
except ZeroDivisionError:
print("x / z fails!")
return
# ~~~~~~~ 534-test_for_narrow_range
def test_for_narrow_range():
"""Test that the 5th power of half an ulp of 1 doesn't underflow.
Basic 5170-5210
"""
# ~~~~~~~ 535-test_for_narrow_range2
y2 = ULP_OF_ONE_MINUS * ULP_OF_ONE_MINUS
y = y2 * y2
y2 = y * ULP_OF_ONE_MINUS
if y2 <= underflow_threshold:
# ***ISSUE: Should this be > underflow_threshold, given message?
if y > min_positive:
bad_cond(err_defect, "")
i = 5
else:
bad_cond(err_serious, "")
i = 4
print("Range is too narrow; ULP_OF_ONE_MINUS^{:d} Underflows."
.format(i))
return
# ~~~~~~~ 536-test_extreme_underflow
def test_extreme_underflow(ut):
"""Test the behavior of underflow_threshold**2.
The code computes
ut+ = ut * (1 + eps)
ut++ = ut * B * (1 + eps)
and then determines where the ultra-small product lies on the line
---------- 0 ---------- ut+ ---------- ut++ ----------
serious OK defect serious
Args:
ut: expect underflow_threshold, but works with any tiny value
Basic 5230-5300
"""
# ~~~~~~~ 537-test_extreme_underflow2
# Compute log(ut) base 1/H, where H=min(1/B, 1/2), so it's
# log base 2, 8, 10 or 16. Use the 240 factor to pick up
# fractional digits base 8, 10, or 16.
y = (-floor(ONE_HALF - TWOFORTY * log(ut) / log(ONE_OVER_H))
/ TWOFORTY)
y2 = y + y
# ~~~~~~~ 538-test_extreme_underflow3
print("Since underflow occurs below the threshold")
print("threshold = ({:0.17e}) ^ ({:0.17e})\nonly underflow "
.format(ONE_OVER_H, y), end="")
print("should afflict the expression\n\t({:0.17e}) ^ ({:0.17e});"
.format(ONE_OVER_H, y2))
print("actually calculating yields:", end="")
# Python has no underflow exception, so no try/except here.
ultra = pow(ONE_OVER_H, y2)
# Can't happen -- bad_cond(err_serious, " trap on underflow.\n")
print(" {:0.17e} .".format(ultra))
# ~~~~~~~ 539-test_extreme_underflow4
# Note use of SAFE_ULPS_OF_ONE versus ULP_OF_ONE_PLUS to ensure
# multiplication doesn't lose the low bits, for lack of a
# guard digit or worse.
if ultra < ZERO or ultra > (B + B * SAFE_ULPS_OF_ONE) * ut:
bad_cond(err_serious, "this is not between 0 and underflow\n")
print(" threshold = {:0.17e} .".format(ut))
elif not (ultra > ut * (ONE + SAFE_ULPS_OF_ONE)):
print("This computed value is O.K.")
else:
bad_cond(err_defect, "this is not between 0 and underflow\n")
print(" threshold = {:0.17e} .".format(ut))
return
# ~~~~~~~ 540-test_for_pseudo_zero
def test_for_pseudo_zero(pseudo):
"""Test for a nonzero value that violates basic numerical laws.
Args:
pseudo: tiny value that might misbehave
Basic 4580-4630
"""
# ~~~~~~~ 541-test_for_pseudo_zero2
if (pseudo == ZERO):
return
print("")
# Test pseudo for "phoney-zero" behavior, violating either
# pseudo < tiny_x or pseudo < pseudo + pseudo
# These are the loop termination tests in tiny_values_and_difference().
if (pseudo <= ZERO):
bad_cond(err_failure, "Positive expressions can underflow to an " +
"allegedly negative value\n")
print("pseudo that prints out as: {:g} .".format(pseudo))
x = -pseudo
if x <= ZERO:
print("But -pseudo, which should be", end="")
print("positive, isn't; it prints out as {:g} ."
.format(x))
# ~~~~~~~ 542-test_for_pseudo_zero3
else:
bad_cond(err_flaw, "Underflow can stick at an allegedly positive\n")
print("value pseudo that prints out as {:g} ."
.format(pseudo))
# One final test for strangeness around the underflow threshold.
# The test is made just for the side-effects, not the reulst value.
discard = does_tiny_value_misbehave(pseudo)
return
# ~~~~~~~ 548-e_squared
def e_squared():
"""Compute exp(2) from the Taylor series for exp(x).
We compute exp(1/2) and square it twice.
exp(1/2) = 1 + 1/2 + 1/8 + 1/48 + 1/(48*2*4) + 1/(48*4*4*5) + ...
= 13/8 + 1/48 * (1 + 1/8 + 1/80 + 1/960 + ...)
Compute the sum as a two-piece hi/lo accumulation and then
evaluate 13/8 + sum/48.
Basic 5320-5350
"""
# ~~~~~~~ 550-e_squared2
# Initialize the loop values so that the first term is 1, addied
# into the initial sum of zero.
x = ZERO
k = 2
term = TWO * THREE
x_lo = ZERO
if verbose:
print("Loop to compute 1/48th the sum of terms 3 to k")
print("of the Taylor series of exp(1/2).")
while True:
old_x = x
k = k + 1
term = term / (k + k)
inc = x_lo + term
x = old_x + inc
# After adding the next term into the sum, capture x_lo, the
# amount of inc that is not captured by x.
x_lo = (old_x - x) + inc
if verbose:
print("\texp2 x={:.6e} x_lo={:.6e} term={:.6e}"
.format(x, x_lo, term))
# Exit the loop when the value x is no longer changing.
if x <= old_x: break
# ~~~~~~~ 550a-e_squared3
e12 = (ONE_AND_HALF + ONE / EIGHT) + x / (ONE_AND_HALF * THIRTY_TWO)
e = e12 * e12
e2 = e * e
if verbose:
print("exp(1/2) = {:.12f} exp(1) = {:.12f} exp(2) = {:.12f}"
.format(e12, e, e2))
return e2
# ~~~~~~~ 554-test_x_to_xp1_over_xm1
def test_x_to_xp1_over_xm1(e2):
""" Compare x^(x+1)/(x-1) to exp(2) for x above and below 1.
Args:
e2: computed value of exp(2)
Basic 5360-5440
"""
print("Testing x^((x + 1) / (x - 1)) vs. " +
"exp(2) = {:0.12f} as x -> 1.".format(exp2))
# BUG: x runs away from 1, perhaps "x near 1"?
# Run tests below 1, then above 1, with a start value
# and a next value one ulp away.
if verbose:
print("Looping on values just below and just above 1.")
for x, next_x in [(ONE_MINUS_ULP, ONE_MINUS_ULP - ULP_OF_ONE_MINUS),
(ONE + ULP_OF_ONE_PLUS,
ONE + ULP_OF_ONE_PLUS + ULP_OF_ONE_PLUS)]:
# ~~~~~~~ 554a-test_x_to_xp1_over_xm12
test_cnt = 1
while True:
z = x - ONE_OVER_B
z = (x + ONE) / (z - (ONE - ONE_OVER_B))
dev = pow(x, z) - exp2
if verbose:
print("\tx = {:.17f} z = {:.6e} deviation of".format(x, z)
+ " x^z = {:.3e}".format(dev))
# ~~~~~~~ 554b-test_x_to_xp1_over_xm13
if (fabs(dev) > TWOFORTY * ULP_OF_ONE_PLUS):
delta_x = (x - ONE_OVER_B) - (ONE - ONE_OVER_B)
bad_cond(err_defect, "Calculated")
print(" {:0.17e} for".format(pow(x, z)))
print("\t(1 + ({:0.17e}) ^ ({:0.17e});".format(delta_x, z))
print("\tdiffers from correct value by {:0.17e} ."
.format(dev))
print("\tThis much error may spoil financial")
print("\tcalculations involving tiny interest rates.")
return
# ~~~~~~~ 554c-test_x_to_xp1_over_xm14
else:
# Jump by 1 ulp, then 2, then 4, etc.
z = (next_x - x) * TWO + next_x
x = next_x
next_x = z
z = ONE + (x - ONE_MINUS_ULP) * (x - ONE_MINUS_ULP)
# *** BUG: P has (z > ONE), which means test done once only.
# Keep testing until too close to 1 or count too big.
if (z <= ONE) and (test_cnt < NUM_TRIALS):
test_cnt += 1
else:
break
# BUG: removed unneeded error flag if N == 0:
print("Accuracy seems adequate.")
return
# ~~~~~~~ 556-extreme_z_to_q
def extreme_z_to_q():
"""Test A and 1/A to powers reaching 1/C and C.
A is 2 or 10. C is 1/B**k, tiny but well away from underflow.
C is an integer power of 1/A, so try the four possibilities of
A and 1/A producing C and 1/C.
Basic 5460-5490
"""
global pow_err_cnt # var used by power tests
print("Testing powers z^q at four nearly extreme values.")
pow_err_cnt = 0 # this test resets the global value to 0 without restore.
# ~~~~~~~ 556a-test_x_to_xp1_over_xm14
q = floor(ONE_HALF - log(C) / log(A)) # C is a power of B and thus A
for z in [A, ONE_OVER_A]:
x = ONE_OVER_C
test_power_result(x, z, q)
q = -q
x = C
test_power_result(x, z, q)
print_if_err_cnt_positive()
if pow_err_cnt == 0:
print(" ... no discrepancies found.")
return
# ~~~~~~~ 560-big_powers_of_B
def big_powers_of_B():
"""Searching for the overflow threshold,
return the NEGATIVES of the largest powers ofB.
With a strategy similar to underflow, use three values in
step to get up to and just beyond the overflow limit.
Returns:
bigger: value that saturates like infinity or triggers an exception
big: down a factor of 1/H
less_big: down a factor of 1/H from big
no_except: boolean indicating no exception triggered
Basic 5530-5530
TODO: Resolve that we use negatives because of possible 2's complement
representation, that could allow the largest magnitude to be
negative.
"""
# ~~~~~~~ 560a-big_powers_of_B2
big = -ONE_OVER_C
bigger = ONE_OVER_H * big
no_except = True
try:
while True:
less_big = big
big = bigger
bigger = ONE_OVER_H * big
if verbose:
print("\tbigger = {:.12e} big = {:.12e}"
.format(bigger, big))
if bigger >= big: break # remember, these values < 0
except OverflowError: # Python does not catch floating point overflow
no_except = False
bigger = big
if verbose:
print("\tEncountered overflow exception.")
return bigger, big, less_big, no_except
# ~~~~~~~ 564-no_trap_overflow_threshold
def no_trap_overflow_threshold(big, less_big):
""" Determine the threshold when there is no trap on overflow.
Typically, the search for the threshold will result in bigger_B and
big_B pinned at Infinity or some huge finite value. Value less_big_B
is guaranteed to be a finite, unexceptional value. This routine either
returns big_B as a finite saturating value or fills out less_big_B
with significant digits to reach right up to big_B.
Args:
big: power of 1/H that may be infinite or saturating finite value
less_big: tame value a nominal factor of H smaller than big
Returns:
oflo: the overflow threshold
"""
# ~~~~~~~ 564a-no_trap_overflow_threshold2
oflo = -less_big # default start
pos_big = -big
y = less_big * (ONE_OVER_H * ULP_OF_ONE_PLUS - ONE_OVER_H)
z = y + ((ONE - ONE_OVER_H) * ULP_OF_ONE_PLUS) * less_big
if z < pos_big:
oflo = z
if verbose:
print("Try full-precision overflow threshold.")
elif y < pos_big:
oflo = y
if verbose:
print("Try reduced precision overflow threshold.")
if pos_big - oflo < pos_big:
oflo = pos_big
if verbose:
print("Finite saturating power of B overflow threshold.")
return oflo
# ~~~~~~~ 568-trap_overflow_threshold
def trap_overflow_threshold(big):
"""Determine the threshold when overflow traps.
The difference from the non-trap case is that the middle value, big,
returned from the saerch is the greatest power of the radix, so the
task is to fill it out to max value without increasing the exponent.
Argunment:
big: biggest finite power of the radix
Returns:
oflo: overflow threshold
"""
oflo = big * (ONE_OVER_H * ULP_OF_ONE_PLUS - ONE_OVER_H)
oflo = oflo + ((ONE - ONE_OVER_H) * ULP_OF_ONE_PLUS) * big
return oflo
# ~~~~~~~ 572-test_sqrt_tiny
def test_sqrt_tiny(tlist):
"""Check sqrt(x)**2 for values near the underflow threshold."""
for z in tlist:
if z != ZERO:
sqrt_z = sqrt(z)
y = sqrt_z * sqrt_z
if ((y / (ONE - B * SAFE_ULPS_OF_ONE) < z)
or y > (ONE + B * SAFE_ULPS_OF_ONE) * z):
if sqrt_z > ULP_OF_ONE_MINUS:
bad_cond(err_serious, "")
else:
bad_cond(err_defect, "")
print("Comparison alleges that what prints as z = {:0.17e}"
.format(z))
print(" is too far from sqrt(z) ^ 2 = {:0.17e} .".format(y))
return
# ~~~~~~~ 576-test_sqrt_big
def test_sqrt_big(tlist):
"""Check sqrt(x)**2 for values near the overflow threshold."""
for z in tlist:
sqrt_z = sqrt(z)
x = (ONE - B * SAFE_ULPS_OF_ONE) * sqrt_z
y = sqrt_z * x
if ((y < (ONE - TWO * B * SAFE_ULPS_OF_ONE) * z) or (y > z)):
if x < BIG_B_NTH:
bad_cond(err_serious, "")
else:
bad_cond(err_defect, "")
print("Comparison alleges that z = {:0.17e}".format(z))
print(" is too far from sqrt(z) ^ 2 ({:0.17e}) .".format(y))
# ~~~~~~~ 580-check_range_balance
def check_range_balance():
"""Check balnce of overflow and underflow thresholds."""
x = underflow_threshold * overflow_threshold
y = B * B
if verbose:
print("\tThe product of the overflow and"
+ "underflow thresholds is {:.4e}".format(x))
if x * y < ONE or x > y:
if x * y < ULP_OF_ONE_MINUS or x > y / ULP_OF_ONE_MINUS:
bad_cond(err_defect, "Badly")
else:
bad_cond(err_flaw, "")
print((" unbalanced range; underflow_threshold * overflow_threshold"
+ "= {:0.17e}\n\tis too far from 1.").format(x))
return
# ~~~~~~~ 584-test_unit_ratios
def test_unit_ratios(tlist):
"""Try x/x for interesting values."""
for x in tlist:
y = x
# Beware values near over/underflow thresholds.
try:
diff = (y / x - ONE_HALF) - ONE_HALF
except ZeroDivisionError:
print(" x / x traps when x = {:g}", x)
if diff != ZERO:
if diff == (-ULP_OF_ONE_MINUS and
(fabs(x) < 1 or fabs(x) > B)):
bad_cond(err_flaw, "")
else:
bad_cond(err_serious, "")
print(" x / x differs from 1 when x = {:0.17e}", x)
print(" instead, x / x - 1/2 - 1/2 = {:0.17e} .", diff)
return
# ~~~~~~~ 588-test_div_by_zero
def test_div_by_zero():
"""Try 1/0 and 0/0."""
my_zero = ZERO
print("")
print("What message and/or values does Division by ZERO produce?")
print("This can interupt your program. You can ", end="")
print("skip this part if you wish.")
s = alt_input("Do you wish to compute 1 / 0?\n")
if s.startswith("N") or s.startswith("n"):
print("O.K. skipping")
else:
print("\tTrying to compute 1 / 0 produces ...")
try:
print(" \t{:0.7e} .".format(ONE / my_zero))
except ZeroDivisionError:
print("\tComputing 1 / 0 triggered an exception.")
s = alt_input("\nDo you wish to compute 0 / 0?\n")
if s.startswith("N") or s.startswith("n"):
print("O.K. skipping")
else:
print("\tTrying to compute 0 / 0 produces ...")
try:
print("\t{:0.7e} .", ZERO / my_zero)
except ZeroDivisionError:
print("\tComputing 0 / 0 triggered an exception.")
return
# ~~~~~~~ 700-start_execution
# PARANOIA EXECUTION BEGINS
verbose |= ultra_verbose # if it's ultra, it's verbose, too
# Small floating point constants built from coerced integers.
ZERO = float(0)
ONE = float(1)
TWO = ONE + ONE
THREE = TWO + ONE
FOUR = THREE + ONE
FIVE = FOUR + ONE
EIGHT = FOUR + FOUR
NINE = THREE * THREE
TWENTY_SEVEN = NINE * THREE
THIRTY_TWO = FOUR * EIGHT
TWOFORTY = FOUR * FIVE * THREE * FOUR
MINUS_ONE = -ONE
ONE_HALF = ONE / TWO
ONE_AND_HALF = ONE + ONE_HALF
# ~~~~~~~ 710-start_messages
milestone = 0 # ==============================
instructions()
pause()
heading()
pause()
characteristics()
pause()
history()
pause()
# ~~~~~~~ 715-zero_one
milestone = 7 # ==============================
# Basic 880-950
print("Program is now RUNNING tests on small integers:")
test_cond(err_failure, (ZERO + ZERO == ZERO)
and (ONE - ONE == ZERO)
and (ONE > ZERO)
and (ONE + ONE == TWO),
"0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2")
# ~~~~~~~ 717-zero_one2
z = -ZERO
if (z != 0.0):
error_count[err_failure] += 1
print("Comparison alleges that Minus Zero, obtained by setting")
print(" X = 0 and then Z = -X , is nonzero!")
# The next two lines artificially set two global constants for the test.
ULP_OF_ONE_PLUS = 0.001
B = 1
discard = does_tiny_value_misbehave(z)
# ~~~~~~~ 720-small_ints
test_cond(err_failure, (THREE == TWO + ONE)
and (FOUR == THREE + ONE)
and (FOUR + TWO * (-TWO) == ZERO)
and (FOUR - THREE - ONE == ZERO),
"3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0")
test_cond(err_failure, (MINUS_ONE == (0 - ONE))
and (MINUS_ONE + ONE == ZERO)
and (ONE + MINUS_ONE == ZERO)
and (MINUS_ONE + fabs(MINUS_ONE) == ZERO)
and (MINUS_ONE + MINUS_ONE * MINUS_ONE == ZERO),
"-1 != 0-1, -1+1 != 0, 1+(-1) != 0, "
+ "(-1)+abs(-1) != 0, or -1+(-1)*(-1) != 0")
test_cond(err_failure, ONE_HALF + MINUS_ONE + ONE_HALF == ZERO,
"1/2 + (-1) + 1/2 != 0")
# ~~~~~~~ 725-med_ints
milestone = 10 # ==============================
test_cond(err_failure, (NINE == THREE * THREE)
and (TWENTY_SEVEN == NINE * THREE) and (EIGHT == FOUR + FOUR)
and (THIRTY_TWO == EIGHT * FOUR)
and (THIRTY_TWO - TWENTY_SEVEN - FOUR - ONE == ZERO),
"9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0")
test_cond(err_failure, (FIVE == FOUR + ONE)
and (TWOFORTY == FOUR * FIVE * THREE * FOUR)
and (TWOFORTY / THREE - FOUR * FOUR * FIVE == ZERO)
and (TWOFORTY / FOUR - FIVE * THREE * FOUR == ZERO)
and (TWOFORTY / FIVE - FOUR * THREE * FOUR == ZERO),
"5 != 4+1, 240 != 4*5*3*4, 240/3 != 80, 240/4 != 60, or 240/5 != 48")
if (error_count[err_failure] == 0):
print("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.")
# ~~~~~~~ 730-radix
# Basic 1160-1230
print("Searching for radix B and PRECISION.")
big_2_to_nth = find_big_2_to_nth()
B = find_radix_from_big_2_to_nth(big_2_to_nth)
if B < TWO: # TODO: log case
B = ONE
print("B = {:0.4f} .\n".format(B))
# ~~~~~~~ 735-precision
# Compute PRECISION and BIG_B_NTH = B**PRECISION barely too big
# to satisfy (BIG_B_NTH + 1) - BIG_B_NTH == 1 .
# Basic 1240-1270
if B == ONE: # defer to the secondary calculation below
PRECISION = ZERO
BIG_B_NTH = big_2_nth
else:
PRECISION, BIG_B_NTH = find_precision_big_B_to_nth(B)
# ~~~~~~~ 737-precision2
guess_ulp_of_one_minus = ONE / BIG_B_NTH
guess_ulp_of_one_plus = B * guess_ulp_of_one_minus
print("Closest relative separation found is ULP_OF_ONE_MINUS = {:0.7e} ."
.format(guess_ulp_of_one_minus))
# ~~~~~~~ 740-recalc_radix
print("Recalculating radix and precision")
ULP_OF_ONE_PLUS = find_ulp_of_one_plus(guess_ulp_of_one_plus)
ULP_OF_ONE_MINUS = find_ulp_of_one_minus(guess_ulp_of_one_minus)
if (ULP_OF_ONE_MINUS == guess_ulp_of_one_minus):
print("confirms closest relative separation ULP_OF_ONE_MINUS ")
else:
print("gets better closest relative separation ULP_OF_ONE_MINUS = "
+ "{:0.7e} .".format(ULP_OF_ONE_MINUS))
# ~~~~~~~ 742-recalc_radix2
BIG_B_NTH = ONE / ULP_OF_ONE_MINUS
ONE_MINUS_ULP = (ONE_HALF - ULP_OF_ONE_MINUS) + ONE_HALF
save_b = B
B = floor(0.01 + ULP_OF_ONE_PLUS / ULP_OF_ONE_MINUS)
B_OVER_TWO = B / TWO
ONE_OVER_B = ONE / B
if (B == save_b):
print("B confirmed.")
else:
print("MYSTERY: recalculated B = {:0.7e} .".format(B))
test_cond(err_defect, B <= EIGHT + EIGHT,
"B is too big: roundoff problems")
test_cond(err_flaw, (B == TWO) or (B == 10)
or (B == ONE), "B is not as good as 2 or 10")
# ~~~~~~~ 745-cmp_ulps_of_one
milestone = 20 # ==============================
test_cond(err_failure, ONE_MINUS_ULP - ONE_HALF < ONE_HALF,
"(1-ULP_OF_ONE_MINUS)-1/2 < 1/2 is FALSE, prog. fails?")
# BUG in C version. The point is to try two cases of x.
for x in [ONE_MINUS_ULP, ONE + ULP_OF_ONE_PLUS]:
y = x - ONE_HALF
z = y - ONE_HALF
test_cond(err_failure, (x != ONE)
or (z == ZERO), "Comparison is fuzzy,x=1 but x-1/2-1/2 != 0")
# ~~~~~~~ 750-the_constants
milestone = 25 # ==============================
# ... B_MINUS_ULP = nextafter(B, 0)
B_MINUS_ULP = B - ONE
B_MINUS_ULP = (B_MINUS_ULP - ULP_OF_ONE_PLUS) + ONE
# A = 2 if B = 2**k or arithmetic is logarithmic; it's 10 if B=10
# C = 1/B**k not too close to underflow, useful for over/underflow analysis.
# H = min(1/B, 1/2).
A, ONE_OVER_A = compute_A_and_inverse(B)
C, ONE_OVER_C = compute_C_and_inverse()
H, ONE_OVER_H = compute_H_and_inverse()
SAFE_ULPS_OF_ONE = compute_safe_ulps_of_one() # 1 or B ulps of one
# ~~~~~~~ 755-frac_prec
# Purify Integers was the original comment. The following test assumes
# B=1 for log arithmetic, before any fractional part of the precision
# is "purified".
if (B != ONE):
PRECISION = standardize_fractional_precision(ULP_OF_ONE_MINUS, B)
if (PRECISION != floor(PRECISION)) or (B == ONE):
print("PRECISION cannot be characterized by an Integer number")
print("of significant digits but, by itself, this is a minor flaw.")
if B == ONE:
print("logarithmic encoding has precision characterized solely "
+ "by ULP_OF_ONE_MINUS.")
else:
print("The number of significant digits of the B is {:0.4f} ."
.format(PRECISION))
# ~~~~~~~ 760-check_prec
test_cond(err_serious, ULP_OF_ONE_PLUS * NINE * NINE * TWOFORTY < ONE,
"PRECISION worse than 5 decimal figures ")
# ~~~~~~~ 765-first_prec_tests
milestone = 30 # ==============================
test_for_extra_precise_subepressions()
pause()
milestone = 35 # ==============================
test_normalized_subtraction()
# ~~~~~~~ 770-guard_digit
flags["mult_guard_digit"] = does_mult_have_guard_digit()
test_mult_low_digits()
flags["div_guard_digit"] = does_div_have_guard_digit()
flags["add_sub_guard_digit"] = does_add_sub_have_guard_digit()
# ~~~~~~~ 772-guard2
if (ONE_MINUS_ULP != ONE and ONE_MINUS_ULP - ONE >= ZERO):
bad_cond(err_serious, "comparison alleges (1-ULP_OF_ONE_MINUS) < 1 although")
print(" subtraction yields (1-ULP_OF_ONE_MINUS) - 1 = 0 , thereby vitiating")
print(" such precautions against division by zero as")
print(" ... if (x == 1.0) {.....} else {.../(x-1.0)...}")
# ~~~~~~~ 774-guard3
if flags["mult_guard_digit"]:
print(" * appears to have a guard digit, as it should.")
if flags["div_guard_digit"]:
print(" / appears to have a guard digit, as it should.")
if flags["add_sub_guard_digit"]:
print(" - appears to have a guard digit, as it should.")
milestone = 40 # ==============================
pause()
# ~~~~~~~ 775-rounding
print("Checking rounding on multiply, divide and add/subtract. ")
flags["mult_rounding"] = how_does_mult_round()
milestone = 45 # ==============================
flags["div_rounding"] = how_does_div_round()
# ~~~~~~~ 777-rounding2
test_cond(err_failure, ONE_OVER_B * B - ONE_HALF == ONE_HALF,
"B * ( 1 / B ) differs from 1")
# ~~~~~~~ 779-rounding3
milestone = 50 # ==============================
test_cond(err_failure, (ONE_MINUS_ULP + ULP_OF_ONE_MINUS) - ONE_HALF == ONE_HALF
and (B_MINUS_ULP + ULP_OF_ONE_PLUS) - ONE == B - ONE,
"Incomplete carry-propagation in Addition")
# ~~~~~~~ 782-rounding4
flags["add_sub_rounding"] = how_does_add_sub_round()
# ~~~~~~~ 785-sticky
# The boolean set_sticky is used to indicate the presence of a sticky bit
# below. It is cleared if the first simple test fails.
set_sticky = rounding_symmetry_test()
# ~~~~~~~ 790-sticky2
if (flags["mult_guard_digit"] and flags["div_guard_digit"]
and flags["add_sub_guard_digit"]
and flags["mult_rounding"] == "rounded"
and flags["div_rounding"] == "rounded"
and flags["add_sub_rounding"] == "rounded"
and floor(B_OVER_TWO) == B_OVER_TWO):
print("Checking for sticky bit.")
set_sticky &= has_sticky_bit()
flags["uses_sticky_bit"] = set_sticky
# ~~~~~~~ 795-round_summary
if flags["uses_sticky_bit"]:
print("Sticky bit apparently used correctly.")
else:
print("Sticky bit used incorrectly or not at all.")
test_cond(err_flaw, (flags["mult_guard_digit"]
and flags["div_guard_digit"]
and flags["add_sub_guard_digit"]
and flags["mult_rounding"] != "other"
and flags["div_rounding"] != "other"
and flags["add_sub_rounding"] != "other"),
"lack(s) of guard digits or failure(s) to correctly round or chop\n"
+ "(noted above) count as one flaw in the final tally below")
# ~~~~~~~ 800-mul_cmmute
milestone = 60 # ==============================
print("")
print("Does Multiplication commute? ", end="")
print("Testing on {:d} random pairs.".format(NUM_TRIALS))
test_mult_commutivity(NUM_TRIALS)
# ~~~~~~~ 805-sqrt
milestone = 70 # ==============================
print("\nRunning test of square root(x).")
test_cond(err_failure, (ZERO == sqrt(ZERO))
and (- ZERO == sqrt(- ZERO))
and (ONE == sqrt(ONE)), "Square root of 0.0, -0.0 or 1.0 wrong")
min_sqrt_error = ZERO # the global sqrt err bounds
max_sqrt_error = ZERO
ok = sqrt_single_bound_test([
(B, ULP_OF_ONE_PLUS, err_serious),
(ONE_OVER_B, ONE_OVER_B * ULP_OF_ONE_MINUS, err_serious),
(BIG_B_NTH, ONE, err_serious), # BUG: Missing from C
(ULP_OF_ONE_MINUS, ULP_OF_ONE_MINUS * ULP_OF_ONE_MINUS, err_serious),
])
if not ok:
pause()
# ~~~~~~~ 810-sqrt_int
print("Testing if sqrt(x * x) == x for {:d} Integers x.".format(NUM_TRIALS))
test_sqrt_int()
print("Test for sqrt monotonicity.")
test_sqrt_monotonicity()
milestone = 80 # ==============================
if verbose:
print("\tAfter monotonicity test:")
print("\t\tmin_sqrt_error = {:0.6e}".format(min_sqrt_error))
print("\t\tmax_sqrt_error = {:0.6e}".format(max_sqrt_error))
# ~~~~~~~ 810a-sqrt_int1
print("Test sqrt for accuracy.")
# TODO: explain this bump - scaling?
min_sqrt_error = min_sqrt_error + ONE_HALF
max_sqrt_error = max_sqrt_error - ONE_HALF
y = (sqrt(ONE + ULP_OF_ONE_PLUS) - ONE) / ULP_OF_ONE_PLUS
up_err = (y - ONE) + ULP_OF_ONE_PLUS / EIGHT
down_err = y + ULP_OF_ONE_PLUS / EIGHT
update_sqrt_errors(down_err, up_err)
if verbose:
print("\tAfter 1 + ulp test:")
print("\t\tmin_sqrt_error = {:0.17e}".format(min_sqrt_error))
print("\t\tmax_sqrt_error = {:0.17e}".format(max_sqrt_error))
# ~~~~~~~ 810b-sqrt_int2
y = (((sqrt(ONE_MINUS_ULP) - ULP_OF_ONE_PLUS) - (ONE - ULP_OF_ONE_PLUS))
/ ULP_OF_ONE_MINUS)
if verbose:
print("sqrt(1-ulp) = {:0.17e}".format(sqrt(ONE_MINUS_ULP)))
print("y = {:0.17e}".format(y))
up_err = y + ULP_OF_ONE_MINUS / EIGHT
down_err = (y + ONE) + ULP_OF_ONE_MINUS / EIGHT
update_sqrt_errors(down_err, up_err)
if verbose:
print("\tAfter 1 - ulp test:")
print("\t\tmin_sqrt_error = {:0.17e}".format(min_sqrt_error))
print("\t\tmax_sqrt_error = {:0.17e}".format(max_sqrt_error))
# ~~~~~~~ 810c-sqrt_int3
ulp = ULP_OF_ONE_PLUS
x = ulp
for (x, ulp) in [(ULP_OF_ONE_PLUS, ULP_OF_ONE_PLUS),
(ULP_OF_ONE_PLUS
* floor(EIGHT / (NINE * sqrt(ULP_OF_ONE_PLUS))),
ULP_OF_ONE_PLUS),
(-ULP_OF_ONE_MINUS, ULP_OF_ONE_MINUS)]:
if verbose:
print("\t3-term sqrt loop: x = {:.17e} ulp = {:.17e}".format(x, ulp))
y = sqrt((x + ULP_OF_ONE_MINUS + x) + ONE_MINUS_ULP)
y = ((y - ULP_OF_ONE_PLUS) - ((ONE - ULP_OF_ONE_PLUS) + x)) / ulp
z = ((ULP_OF_ONE_MINUS - x) + ONE_MINUS_ULP) * ONE_HALF * x * x / ulp
down_err = (y + ONE_HALF) + z
up_err = (y - ONE_HALF) + z
update_sqrt_errors(down_err, up_err)
if verbose:
print("\tAfter ulp of 1+ cases")
print("\t\tmin_sqrt_error = {:0.17e}".format(min_sqrt_error))
print("\t\tmax_sqrt_error = {:0.17e}".format(max_sqrt_error))
# ~~~~~~~ 815-sqrt_num_theory
milestone = 85 # ==============================
sqrt_rnd_odd = False
sqrt_anomalous_int = False
tests = 0
flags["sqrt_rounding"] = "other"
while B != ONE:
print("Testing whether sqrt is rounded or chopped.")
D = floor(ONE_HALF + pow(B, ONE + PRECISION - floor(PRECISION)))
# ... == B^(1 + fract) if (PRECISION == Integer + fract.
x = D / B
y = D / A
if x != floor(x) or y != floor(y):
sqrt_anomalous_int = True
break
# ~~~~~~~ 815c-sqrt_num_theory3
x = ZERO
z2 = x
y = ONE
y2 = y
# **** B = 8 #**********************************************
# *** D = 16 #**********************************************
z1 = B - ONE
FourD = FOUR * D
# ~~~~~~~ 815f-sqrt_num_theory6
while True:
if (y2 > z2):
Q = B
y1 = y
while True:
x1 = fabs(Q + floor(ONE_HALF - Q / y1) * y1)
Q = y1
y1 = x1
if ultra_verbose:
print("\t\t\tloop inner x1 = {:.6e} Q = {:.6e} y1 = {:.6e}"
.format(x1, Q, y1))
if x1 <= ZERO: break
if (Q <= ONE):
z2 = y2
z = y
# ~~~~~~~ 815i-sqrt_num_theory9
y = y + TWO
x = x + EIGHT
y2 = y2 + x
if y2 >= FourD:
y2 = y2 - FourD
if ultra_verbose:
print("\t\t\tloop x = {:.6e} x1 = {:.6e}".format(x, x1))
print("\t\t\t y = {:.6e} y1 = {:.6e} y2 = {:.6e}".format(y, y1, y2))
print("\t\t\t z = {:.6e} z1 = {:.6e} z2 = {:.6e}".format(z, z1, z2))
if y >= D: break
# ~~~~~~~ 815l-sqrt_num_theory12
X8 = FourD - z2
# ***JTC: fix global Q used by num_theory_sqrt_value
Q = (X8 + z * z) / FourD
X8 = X8 / EIGHT
if ultra_verbose:
print("Loop finds Q = {:.17e} X8 = {:.17e}".format(Q, X8))
if Q != floor(Q):
sqrt_anomalous_int = True
break
# ~~~~~~~ 815n-sqrt_num_theory14
while True:
x = z1 * z
x = x - floor(x / B) * B
if x != ONE:
z1 = z1 - ONE
if x == ONE or z1 <= ZERO: break
if z1 <= ZERO and x != ONE:
sqrt_anomalous_int = True
break
# ~~~~~~~ 815p-sqrt_num_theory16
if z1 > B_OVER_TWO:
z1 = z1 - B
while True:
x, D, Q, z, z1 = num_theory_sqrt_value(x, D, Q, z, z1)
if ULP_OF_ONE_PLUS * D >= ONE_MINUS_ULP: break
if D * B - D != BIG_B_NTH - D:
sqrt_anomalous_int = True
break
# ~~~~~~~ 815r-sqrt_num_theory18
z2 = D
y = D + (ONE + z) * ONE_HALF
x = D + z + Q
tests += num_theory_sqrt_test(x, D, y, z2)
y = D + (ONE - z) * ONE_HALF + D
x = D - z + D
x = x + Q + x
tests += num_theory_sqrt_test(x, D, y, z2)
x, D, Q, z, z1 = num_theory_sqrt_value(x, D, Q, z, z1)
if D - z2 != BIG_B_NTH - z2:
sqrt_anomalous_int = True
break
# ~~~~~~~ 815t-sqrt_num_theory20
y = (D - z2) + (z2 + (ONE - z) * ONE_HALF)
x = (D - z2) + (z2 - z + Q)
tests += num_theory_sqrt_test(x, D, y, z2)
y = (ONE + z) * ONE_HALF
x = Q
tests += num_theory_sqrt_test(x, D, y, z2)
break # exit ficttitious loop in one pass
# ~~~~~~~ 815u-sqrt_num_theory22
if verbose:
print("\tAfter all sqrt test computations:")
print("\t\tmin_sqrt_error = {:0.17e}".format(min_sqrt_error))
print("\t\tmax_sqrt_error = {:0.17e}".format(max_sqrt_error))
if tests == 0 or sqrt_anomalous_int:
bad_cond(err_failure, "Anomalous arithmetic with Integer < ")
print("B^PRECISION = {:0.7e}".format(BIG_B_NTH))
print(" fails test whether sqrt rounds or chops.")
sqrt_rnd_odd = True
else:
if min_sqrt_error >= ZERO and max_sqrt_error <= ZERO:
flags["sqrt_rounding"] = "rounded"
print("Square root appears to be correctly rounded.")
else:
if (max_sqrt_error + ULP_OF_ONE_PLUS > ULP_OF_ONE_PLUS - ONE_HALF
or min_sqrt_error > ONE_HALF
or min_sqrt_error + B < ONE_HALF):
sqrt_rnd_odd = True
else:
flags["sqrt_rounding"] = "chopped"
print("Square root appears to be chopped.")
# ~~~~~~~ 815w-sqrt_num_theory24
if sqrt_rnd_odd:
print("Square root is neither chopped nor correctly rounded.")
print("Observed errors run from {:0.7e} "
.format(min_sqrt_error - ONE_HALF), end="")
print("to {:0.7e} ulps.".format(ONE_HALF + max_sqrt_error))
test_cond(err_serious, max_sqrt_error - min_sqrt_error < B * B,
"sqrt gets too many last digits wrong")
# ~~~~~~~ 820-pow_small_int
milestone = 90 # ==============================
pause()
print("Testing powers z^i for small Integers z and i.")
pow_err_cnt = 0
# ... test powers of zero.
# If(-1)^k is invalid, replace MINUS_ONE by ONE.
for (z, i) in [(-ZERO, 0), (MINUS_ONE, -4)]:
x = ONE
test_integer_powers(x, z, 0, 3)
# For higher odd powers z is the starting value of x.
test_integer_powers(z, z, 1023, 1026)
print_if_err_cnt_positive()
save_pow_err_cnt = pow_err_cnt
pow_err_cnt = 0
m = int(floor(TWO * log(BIG_B_NTH) / log(A)))
for z in [A, ONE_OVER_A]: # Executes 2x if A = 1/A = 1
test_integer_powers(z, z, 1, m)
# ~~~~~~~ 825-pow_primes
milestone = 100 # ==============================
# Powers of B have been tested, sonext try a few primes.
# The simpler loop skips the trick of stepping by 2 through the odd
# numbers from 3, filtering the further multiples of 3.
m = NUM_TRIALS
for z in [3, 5, 7, 11, 13, 17, 19, 23]:
test_integer_powers(z, z, 1, m)
if pow_err_cnt > 0:
print("Errors like this may invalidate financial calculations")
print("\tinvolving interest rates.")
print_if_err_cnt_positive()
pow_err_cnt += save_pow_err_cnt
if pow_err_cnt == 0:
print("... no discrepancies found.")
if pow_err_cnt > 0:
pause()
else:
print("")
# ~~~~~~~ 900-uflow1
milestone = 110 # ==============================
# UNDERFLOW AND THE BEHAVIOR OF TINY VALUES
# Starts at Basic 4330
print("Seeking three values related to Underflow:")
print(" underflow_threshold = smallest well-behaved value")
print(" min_positive = smallest nonzero value discovered")
print(" too_tiny_x = a value underflowed to zero or a pseudo-zero")
# First step: compute tiny powers of the radix. H and C are two useful
# constants.
# H = min(1/B, 1/2); ONE_OVER_H = max(B, 2)
# C = 1/B**k, not too close to the underflow thereshold; ONE_OVER_C = B**k
less_tiny_B, tiny_B, too_tiny_B = tiny_powers_of_B(C)
# ~~~~~~~ 900a-uflow1a
# Constant SAFE_ULPS_OF_ONE = B OR 1 * ULP_OF_ONE_PLUS should have
# enough units in the last place of 1 to leave something nonzero
# in multiplication, even if there's no guard digit.
one_plus_safe = ONE + SAFE_ULPS_OF_ONE
# Compute a rich set of values analogous to the tiniest powers of B, but with
# nonzero low-order bits and attention to the differences of small values.
(less_tiny_x, tiny_x, too_tiny_x, underflow_threshold, add_sub_tiny,
add_sub_hi, add_sub_lo) = tiny_values_and_difference(one_plus_safe)
print("***less_tiny_x, tiny_x, too_tiny_x = {:.7e} {:.7e} {:.7e}"
.format(less_tiny_x, tiny_x, too_tiny_x))
print("***underflow_threshold, add_sub_tiny, add_sub_hi = {:.7e} {:.7e} {:.7e}"
.format(underflow_threshold, add_sub_tiny, add_sub_hi))
# ~~~~~~~ 902-uflow2
# The two sets of tiny values characterize the underflow behavior.
# The following comments mirror BASIC comments 4530-4560. The relational ">~"
# may be loosely read "greater than or kind of equal to". It refers to the odd
# behavior of some arithmetics near 0, where some "pseudo zero" values behave
# like zero in some contexts and nonzero in others.
#
# 1 >> C = 1/B^k >= less_tiny_B > tiny_B >~ too_tiny_B >~ 0
# where the tiny values differ by a factor of H
#
# 1 >> d = (1+SAFE_ULPS_OF_ONE)*C >= underflow_threshold
# >= add_sub_hi >= less_tiny_x > tiny_x >~ too_tiny_x >~ 0
# where underflow_threshold = d * H^k is the first NONZERO value to
# violate (x*H)/H = x else underflow_threshold = 0
# and where add_sub_hi = underflow_threshold * H^k is the first to satisfy
# add_sub_tiny <- |(x*H)/H - x| > 0 else add_sub_hi = less_tiny_x
test_for_pseudo_zero(too_tiny_x)
# ~~~~~~~ 903-uflow3
milestone = 120 # ==============================
# less_tiny_B and less_tiny_x result from similar loops, one starting
# with C, the other with (1 + SAFE_ULPS_OF_ONE) * C.
# Check to see who is the tiniest of them all, but scale up
# to avoid faulty comparisons as on CDC 6000.
# If the number encoding is such that what "looks like" the smallest power
# of the radix is interpreted as zero, then the smallest nonzero values
# will not be power of B. This test catches that case. The value h_fact
# records that the smallest power of B is too big by nearly a factor of H.
if ONE_OVER_C * less_tiny_B > ONE_OVER_C * less_tiny_x:
h_fact = H
min_positive = tiny_x
else:
h_fact = ONE
min_positive = tiny_B
# ~~~~~~~ 904-uflow4
# Compared here are min_positive, the smallest nonzero value
# gotten by successively multiplying by H, and add_sub_tiny,
# which is |x - (x*H)/H| for a tiny x.
if (add_sub_tiny != ZERO) and (add_sub_tiny != min_positive):
bad_cond(err_defect, "")
if (add_sub_tiny < min_positive):
print("Products underflow at a higher", end="")
print(" threshold than differences.")
if too_tiny_x == ZERO: # Only if there are NO pseudo-zeros.
# Don't admit add_sub_tiny if it could be a pseudo-zero.
min_positive = add_sub_tiny
else:
print("Difference underflows at a higher", end="")
print(" threshold than products.")
print("Smallest strictly positive number found is min_positive = {:g} ."
.format(min_positive))
# ~~~~~~~ 905-uflow5
# Now consolidate the three triples of tiny results into the values
# from the loop computation. Break into several cases...
# This defines base_tiny, which may be needed as a final result in some
# anomalous arithmetics.
if does_tiny_value_misbehave(min_positive):
base_tiny = less_tiny_B # Back up by a factor of 1/H, for CDC 7600
else:
base_tiny = min_positive
# ~~~~~~~ 906-uflow6
# The Basic code switches on four cases depending on whether two values are
# zero or nonzero.
# add_sub_tiny: the difference between two tiny unequal values may
# underflow to zero
# underflow_threshold: intended to be the smallest normal number
# unaffected by underflow, but will not have been set yet for
# arithmetic flushes underflow abruptly to zero
#
# -----------CASES---- 4 3 2 1
# add_sub_tiny non0 non0 0 0
# underflow_threshold non0 0 non0 0
#
# In high-quality arithmetic, both values are nonzero, so it's just a matter
# of checking for gradual underflow and its quality.
# If on the other hand, either value is zero, a specific action is taken and
# then all three cases fall into an explicit check for differences of zero
# arising from unequal values.
# ~~~~~~~ 907-uflow7
if add_sub_tiny and underflow_threshold:
# BASIC case 4
if ((add_sub_hi == underflow_threshold) and (add_sub_tiny == min_positive)
and (fabs(underflow_threshold - add_sub_tiny / SAFE_ULPS_OF_ONE)
<= add_sub_tiny)):
flags["IEEE"] = is_gradual_underflow_IEEE(min_positive)
# ~~~~~~~ 908-uflow8
else:
# BASIC cases 1-3. Case 3 falls through to a final test.
if underflow_threshold == ZERO: # Cases 1-2
if add_sub_tiny:
# Case 2
underflow_threshold = add_sub_hi # last good behavior
bad_cond(err_failure,
"Underflow confuses Comparison, which alleges that")
print("add_sub_hi == add_sub_lo while denying that |add_sub_hi"
+ " - add_sub_lo| == 0; these values")
print("print out as add_sub_hi = {:0.17e}, "
+ "add_sub_lo ={:0.17e} .".format(add_sub_hi, add_sub_lo))
print("|add_sub_hi - add_sub_lo| = {:0.17e} ."
.format(fabs(add_sub_hi - add_sub_lo)))
# BUG: C lacks the following message.
print("and add_sub_hi / add_sub_lo| = 1.0 + {:0.17e} ."
.format((add_sub_hi / add_sub_lo - ONE_HALF) - ONE_HALF))
# ~~~~~~~ 909-uflow9
else:
# Case 1
underflow_threshold = base_tiny # Safe value computed above
# A sanity check determines that the result from
# tiny_x_and_difference() matches expectation, based on the
# loop on powers of B. The factor one_plus_safe accounts for
# the difference between the tiny_B and tiny_x loops.
# h_fact accounts for the case that what ought to be the
# tiniest power of B in the representation is given the value
# zero. If these don't match, the threshold is backed up a
# factor of B.
if ((ONE_OVER_C * add_sub_hi)
!= ((ONE_OVER_C * less_tiny_B) * one_plus_safe * h_fact)):
underflow_threshold = less_tiny_B
bad_cond(err_failure,
"Either accuracy deteriorates as numbers\n")
print("approach a threshold = {:0.17e}"
.format(underflow_threshold))
print(" coming down from {:0.17e}".format(C))
print(" or else multiplication gets too many ", end="")
print("last digits wrong.")
pause()
# Cases 1-3 fall through to a final test for zero differences.
test_tiny_differences()
# ~~~~~~~ 910-uflow10
print("The Underflow threshold is {:0.17e}, below which"
.format(underflow_threshold))
print("calculation may suffer larger Relative error than ", end="")
print("merely roundoff.")
test_for_narrow_range()
milestone = 130 # ==============================
test_extreme_underflow(underflow_threshold)
# ~~~~~~~ 930-exp2
milestone = 140 # ==============================
print("")
# ...calculate exp2 == exp(2) == 7.389056099...
exp2 = e_squared()
test_x_to_xp1_over_xm1(exp2)
# ~~~~~~~ 932-z_to_q
milestone = 150 # ==============================
# Test z^q at 4 nearly extreme values.
extreme_z_to_q()
print("")
# ~~~~~~~ 934-oflow1
milestone = 160 # ==============================
pause()
print("Searching for Overflow threshold:")
print("This may generate an error.")
bigger_B, big_B, less_big_B, no_oflow_except = big_powers_of_B()
# In a system that saturates at INFINITY, only less_big_B will be finite.
overflow_threshold = less_big_B # first approxomation, to be refined
print("Can z = -y overflow?")
print("Trying it on big_B = {:0.17e} .".format(big_B))
v = -big_B
pos_big_B = v
if overflow_threshold - big_B == overflow_threshold + pos_big_B:
print("Seems O.K.")
else:
print("finds a ")
bad_cond(err_flaw, "-(-big_B) differs from big_B.")
if bigger_B != big_B:
bad_cond(err_serious, "")
print("overflow past {:0.17e}\n\tshrinks to {:0.17e} ."
.format(big_B, bigger_B))
# ~~~~~~~ 938-oflow2
if no_oflow_except:
overflow_threshold = no_trap_overflow_threshold(big_B, less_big_B)
else:
overflow_threshold = trap_overflow_threshold(big_B)
print("Computed overflow_threshold = {:0.17e} .".format(overflow_threshold))
if no_oflow_except:
print("Overflow saturates at pos_big_B = {:0.17e} .".format(pos_big_B))
else:
print("There is no saturation value because " +
"the system traps on overflow.")
# ~~~~~~~ 934a-oflow1a
print("----Starting overflow watch.----")
print("No Overflow should be signaled for overflow_threshold * 1 = ", end="")
v = overflow_threshold * ONE
print("{:0.17e}".format(v))
print(" nor for overflow_threshold / 1 = ", end="")
v = overflow_threshold / ONE
print("{:0.17e} .".format(v))
print("Any overflow signal during this watch is a DEFECT.")
print("----End of overflow watch.---")
# ~~~~~~~ 938-oflow_comparisions
milestone = 170 # ==============================
if not (-overflow_threshold < overflow_threshold and
-pos_big_B < pos_big_B and
-underflow_threshold < overflow_threshold and
underflow_threshold < overflow_threshold):
bad_cond(err_failure, "Comparisons involving ")
print("+-{:.6e}, +-{:.6e}\nand +-{:.6e} are confused by Overflow."
.format(overflow_threshold, pos_big_B, underflow_threshold))
# ~~~~~~~ 942-tiny_sqrt_tests
milestone = 175 # ==============================
print("")
test_sqrt_tiny([underflow_threshold, min_positive, too_tiny_x])
milestone = 180 # ==============================
# Peel off historical case from the more typical else clause.
if B == 2 and PRECISION == 56 and too_tiny_B != ZERO and -ZERO != ZERO:
print("FAILURE: Attempts to evaluate SQR(overflow_threshold)")
print("in Double Precision n BASIC on the IBM PC display the word")
print("`Overflow' and then disable the Keyboard! This is disastrous.")
else:
test_sqrt_big([overflow_threshold, pos_big_B])
# ~~~~~~~ 950-range_unit_div_zero
milestone = 190 # ==============================
pause()
print("Checking balance of over/underflow thresholds.")
check_range_balance()
milestone = 200 # ==============================
print("Checking x/x for various x, some extreme.")
test_unit_ratios([ONE_MINUS_ULP,
ONE + ULP_OF_ONE_PLUS,
overflow_threshold,
underflow_threshold,
B])
milestone = 210 # ==============================
test_div_by_zero()
# ~~~~~~~ 960-final_report
milestone = 220 # ==============================
pause()
print("")
msg = ["FAILUREs encountered =",
"SERIOUS DEFECTs discovered =",
"DEFECTs discovered =",
"FLAWs discovered ="
]
for i in range(0, 4):
if error_count[i]:
print("The number of {:>29s} {:d}.".format(msg[i], error_count[i]))
# ~~~~~~~ 960a-final_report1
print("")
if (error_count[err_failure] or error_count[err_serious] or
error_count[err_defect] or error_count[err_flaw]):
if (error_count[err_failure] + error_count[err_serious]
+ error_count[err_defect] == 0) and (error_count[err_flaw] > 0):
print("The arithmetic diagnosed seems ", end="")
print("Satisfactory though flawed.")
# ~~~~~~~ 960b-final_report2
elif ((error_count[err_failure] + error_count[err_serious] == 0)
and error_count[err_defect] > 0):
print("The arithmetic diagnosed may be Acceptable", end="")
print(" despite inconvenient Defects.")
# ~~~~~~~ 960c-final_report3
elif error_count[err_failure] or error_count[err_serious]:
print("The arithmetic diagnosed has", end="")
print(" unacceptable Serious Defects.")
if error_count[err_failure]:
print("Potentially fatal FAILURE may have spoiled this", end="")
print(" program's subsequent diagnoses.")
# ~~~~~~~ 960d-final_report4
else:
print("No failures, defects nor flaws have been discovered.")
if not (flags["mult_rounding"] == "rounded"
and flags["div_rounding"] == "rounded"
and flags["add_sub_rounding"] == "rounded"
and flags["sqrt_rounding"] == "rounded"):
print("The arithmetic diagnosed seems Satisfactory.")
else:
if (flags["uses_sticky_bit"] and
(B - TWO) * (B - NINE - ONE) == ZERO):
print("Rounding appears to conform to ", end="")
print("IEEE standard ", end="")
if ((B == TWO) and
((PRECISION - FOUR * THREE * TWO)
* (PRECISION - TWENTY_SEVEN
- TWENTY_SEVEN + ONE) == ZERO)):
print("754", end="")
else:
print("854", end="")
if flags["IEEE"]:
print(".")
else:
print(",\nexcept for possibly Double Rounding", end="")
print(" during Gradual Underflow.")
print("The arithmetic diagnosed appears to be Excellent!")
# ~~~~~~~ 960e-final_report5
if fpecount:
print("\nA total of {:d} floating point exceptions were registered."
.format(fpecount))
print("END OF TEST.")
exit(0)