Computer Arithmetic Emulation For Rounding-Error Analysis

William Schwartz

Cerro Coso Community College

Introduction

1. Applications.
2. Package features.
3. How to use the package.
4. Simple examples.

Applications

1. You are suspicious of an algorithm's performance.
2. Two different output devices give different results using the 'same' algorithm.
3. Development of new algorithms.
4. Checking an algorithm over its input-parameter space by sampling the input-parameter space.
5. Teaching.

Package Features

IEEE roundings
1. round-to-zero (choping) 
2. round-to-nearest (regular rounding:  e.g., less than .5, round down; 
		greater than or equal to .5, round up)
3. round-to-positive-infinity
4. round-to-negative-infinity

Real and complex arithmetic for the following: addition, subtraction, multiplication, division, 
	exponentiation (for exponentiation, only square roots; integer powers of negative numbers; 
	integer powers of complex numbers), absolute value, and sign; for real numbers, all the trig 
	functions, their inverses, (e^x),  ln(x), and log(base,x) ; logical operators And,Or, Negation, 
	Equal, Unequal, SameQ, UnsameQ, less than, greater than, less than or equal, and greater 
	than or equal.

For binary operations the software converts the computer numbers to rationals, applies the exact 
	binary operation, and then converts the result to a computer number using the rounding selected.

For unary operations the software converts the computer number to rational, applies the exact 
	unary operation, and then converts the result to a computer number using the rounding selected.

Conversion of real constants (Integer, Rational, Real) to computer numbers is done automatically;
	only numbers represented by variables must be converted to computer numbers by the user.

One can use any integer base, 2 or greater.

The package checks for underflow and overflow. 

CN[{1, 1, 0, 0, 1, 1, 0, 0, 1, 1}, -3, -1]

C$[CN[{1, 1, 0, 0, 0, 0, 0, 0, 0, 0}, 2, 1], CN[{1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 3, 1]]

How To Use The Package

FUNCTION comp1(a, b, x, y, z )

	INTEGER 		x,	 		z
	REAL 					a, 			b
	COMPLEX 		y

		comp1 = (21/10)*((x/y)/2) ** 6 + 3 ** 4 + (y ** x)/x ** z) + y/(34/10) + y/(3/4) + Sin (z)/y + 1000/(a - b)

		RETURN
		END

Needs["BS$REA`"]

ShowParameters
SetDigits[60]

a=1000000000001/10000000000000;
b=1000000000000/10000000000000;
x=3;
y=3+4*I;
z=5;
TrueValue=(21/10)*((((x/y)/2))^6)+(3^4)+(y^x)/(x^z)+y/(34/10)+y/(3/4)+Sin[z]/y+1000/(a-b);
{TrueValue,"          ",N[TrueValue,20]}
{a,b,x,y,z}=Map[ConvertToCN,{a,b,x,y,z}];
CAE$Value=(21/10)*((((x/y)/2))^6)+(3^4)+(y^x)/(x^z)+y/(34/10)+y/(3/4)+Sin[z]/y+1000/(a-b);
CAE$Value=CN$ToRational[CAE$Value];
{CAE$Value,"          ", N[CAE$Value,20]}
N[TrueValue-CAE$Value]

Example 1: Overflow

SetDigitsND[16]

a=(2^124);
b=(2^(-5));
{a,b}=Map[ConvertToCN,{a,b}]

CAE$Value=a/b
CN$ToRational[CAE$Value]

a=-1;
b=0;
{a,b}=Map[ConvertToCN,{a,b}]

CAE$Value=a/b
CN$ToRational[CAE$Value]

Example 2: Do Loop Increments

n=100;
inc=8/9;
index=0;
Do[
		index=index+inc
(*EndDo*),{i,n}];
TrueIndex=index;

SetDigitsND[8]
inc=ConvertToCN[8/9];
index=ConvertToCN[0];
Do[
		index=index+inc
(*EndDo*),{i,n}];
CAE$Index=index;
Grid[{{"true value","CAE value"},{N[TrueIndex,6],N[CN$ToRational[CAE$Index],6]}},
	TableDirections->{Column},Spacings->{4,.3}]

Example 3: Quadratic Formula

root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a))                 root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a))

a=1/8;
b=211/12;
c=11/20;
TrueRoot1=N[((-b-(Sqrt[(b^2)-4 a c]))/(2 a)),10];
TrueRoot2=N[((-b+(Sqrt[(b^2)-4 a c]))/(2 a)),10];

SetBaseND[2]
SetRoundingND["RoundToZero"]

SetDigitsND[2];
a1=a;b1=b;c1=c;
{a,b,c}=Map[ConvertToCN,{a1,b1,c1}];
CAE$Root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a));
CAE$Root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a));
Print[2,"           ",TableForm[{{TrueRoot1,TrueRoot2},N[{CN$ToRational[CAE$Root1],CN$ToRational[CAE$Root2]},10]},TableDirections->{Column,Row},
	TableHeadings->{{"true root","CAE root"},{"root1","root2"}},TableAlignments->{Center}]]

SetDigitsND[4]
{a,b,c}=Map[ConvertToCN,{a1,b1,c1}];
CAE$Root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a));
CAE$Root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a));
Print[4,"           ",TableForm[{{TrueRoot1,TrueRoot2},N[{CN$ToRational[CAE$Root1],CN$ToRational[CAE$Root2]},10]},TableDirections->{Column,Row},
	TableHeadings->{{"true root","CAE root"},{"root1","root2"}},TableAlignments->{Center}]]

SetDigitsND[8]
{a,b,c}=Map[ConvertToCN,{a1,b1,c1}];
CAE$Root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a));
CAE$Root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a));
Print[8,"           ",TableForm[{{TrueRoot1,TrueRoot2},N[{CN$ToRational[CAE$Root1],CN$ToRational[CAE$Root2]},10]},TableDirections->{Column,Row},
	TableHeadings->{{"true root","CAE root"},{"root1","root2"}},TableAlignments->{Center}]]

SetDigitsND[16]
{a,b,c}=Map[ConvertToCN,{a1,b1,c1}];
CAE$Root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a));
CAE$Root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a));
Print[16,"           ",TableForm[{{TrueRoot1,TrueRoot2},N[{CN$ToRational[CAE$Root1],CN$ToRational[CAE$Root2]},10]},TableDirections->{Column,Row},
	TableHeadings->{{"true root","CAE root"},{"root1","root2"}},TableAlignments->{Center}]]

SetDigitsND[32]
{a,b,c}=Map[ConvertToCN,{a1,b1,c1}];
CAE$Root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a));
CAE$Root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a));
Print[32,"           ",TableForm[{{TrueRoot1,TrueRoot2},N[{CN$ToRational[CAE$Root1],CN$ToRational[CAE$Root2]},10]},TableDirections->{Column,Row},
	TableHeadings->{{"true root","CAE root"},{"root1","root2"}},TableAlignments->{Center}]]

SetDigitsND[45]
{a,b,c}=Map[ConvertToCN,{a1,b1,c1}];
CAE$Root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a));
CAE$Root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a));
Print[45,"           ",TableForm[{{TrueRoot1,TrueRoot2},N[{CN$ToRational[CAE$Root1],CN$ToRational[CAE$Root2]},10]},TableDirections->{Column,Row},
	TableHeadings->{{"true root","CAE root"},{"root1","root2"}},TableAlignments->{Center}]]

SetDigitsND[46]
{a,b,c}=Map[ConvertToCN,{a1,b1,c1}];
CAE$Root1=((-b-(Sqrt[(b^2)-4 a c]))/(2 a));
CAE$Root2=((-b+(Sqrt[(b^2)-4 a c]))/(2 a));
Print[46,"           ",TableForm[{{TrueRoot1,TrueRoot2},N[{CN$ToRational[CAE$Root1],CN$ToRational[CAE$Root2]},10]},TableDirections->{Column,Row},
	TableHeadings->{{"true root","CAE root"},{"root1","root2"}},TableAlignments->{Center}]]

Example 4: Runga Kutta Order 4

y'(t)=-y-(t^2)+3, y(0)=1/2 over the interval [0,12]

y(t)=-(1/2)(e^(-t))(1-2(e^t)-4t (e^t)+2(t^2) (e^t))

SUBROUTINE RK4( t0, y0, h, steps, t, y )

	INTEGER 		i				steps
	REAL 					h,			k1,					k2,					k3,					k4,					t,					t0, 				y,				y0
	DIMENSION					t(1),				y(1)

		t(1)=t0
		y(1)=y0
		DO 10 i=1,steps
				k1=h*f(	t(i),  y(i)  )
				k2=h*f(	t(i)+h/2,  y(i)+k1/2  )
				k3=h*f(	t(i)+h/2,  y(i)+k2/2  )
				k4=h*f(	t(i)+h,  y(i)+k3  )
				t(i+1)=t0+i*h
				y(i+1)=y(i)+(k1 + 2*k2 + 2*k3 + k4)/6
10	ENDDO
		RETURN
		END

FUNCTION f(t,y)
	REAL 					f				t, 			y
		f=-y-(t**2)+3
		RETURN
		END

Clear[f]
f[t_,y_]:=-y-t^2+3

Clear[RK4]
RK4[f_,t0_,y0_,h_,steps_]:=

		Module[

				{
				i,t={t0},y={y0}
				},


				Do[
						k1=h*f[  t[[i]], y[[i]]  ];
						k2=h*f[  t[[i]]+h/2,  y[[i]]+k1/2];
						k3=h*f[  t[[i]]+h/2,  y[[i]]+k2/2];
						k4=h*f[  t[[i]]+h,  y[[i]]+k3];
						t=Append[t,t0+i*h];
						y=Append[y,y[[i]]+(k1+2*k2+2*k3+k4)/6],
				(*EndDo*) {i,steps}];

				Return[Transpose[{t,y}]];

		(*EndModule*)];
(*EndRK4*)

code for computing and displaying results

Clear[f]
f[t_,y_]:=-y-t^2+3
y0=1/2;
Clear[y]
y[t_]=y[t]/.Flatten[DSolve[{y'[t]==-y[t]-t^2+3,y[0]==y0},y,t]];

Clear[CompareRK4]
CompareRK4[y_,f_,t0_,y0_,h_,steps_]:=

		Module[

				{
				hCN,len1,len2,res,res1,res2,t0CN,y0CN
				},


				res1=RK4[f,t0,y0,h,steps];
				len1=Length[res1];

				{t0CN,y0CN,hCN}=Map[ConvertToCN,{t0,y0,h} ];

				res2=RK4[f,t0CN,y0CN,hCN,steps];
				res2=Map[CN$ToRational,res2,{2}];
				len2=Length[res2];

				res=MapThread[{#1[[1]],0,Abs[#1[[1]]-#2[[1]]],y[ #1[[1]] ],#1[[2]],#2[[2]],
					Abs[y[ #1[[1]] ]-#1[[2]]], Abs[#2[[2]]-#1[[2]]] }&,{res1,res2}];

				res=Map[{IntFracS[#[[1]],5,4],IntFracS[#[[2]],5,4],IntFracS[#[[3]],5,6],
					IntFracS[#[[7]],5,10],IntFracS[#[[8]],5,10]}&,SetPrecision[res,30]];

				Return[
						TableForm[res,TableAlignments->{Center,Center},
							TableHeadings->{None,{"t","|t-(Subscript[t, a])|","|t-(Subscript[t, c])|","|y-(Subscript[y, a])|","|(Subscript[y, a])-(Subscript[y, c])|"}},TableSpacing->{1,7}]
			(*EndReturn*)];

		(*EndModule*)];
(*EndCompareRK4*)

Clear[CompareRK4Dec]
CompareRK4Dec[y_,f_,t0_,y0_,h_,steps_]:=

		(* 
		for step size 1/4. 
		*)

		Module[

				{
				hCN,len1,len2,res,res1,res2,t0CN,y0CN
				},


				res1=RK4[f,t0,y0,h,steps];
				len1=Length[res1];

				{t0CN,y0CN,hCN}=Map[ConvertToCN,{t0,y0,h} ];

				res2=RK4[f,t0CN,y0CN,hCN,steps];
				res2=Map[CN$ToRational,res2,{2}];
				len2=Length[res2];

				res=MapThread[{#1[[1]],0,Abs[#1[[1]]-#2[[1]]],y[ #1[[1]] ],#1[[2]],#2[[2]],
					Abs[y[ #1[[1]] ]-#1[[2]]], Abs[#2[[2]]-#1[[2]]] }&,{res1,res2}];

				res=Map[{IntFracS[#[[1]],5,4],IntFracS[#[[2]],5,4],IntFracS[#[[3]],5,6],
					IntFracS[#[[7]],5,10],IntFracS[#[[8]],5,10]}&,SetPrecision[res,30]];

				(* for step size 1/4. *)

				res1={res[[1]]};
				res=Join[res1,Table[res[[i]],{i,5,49,4}]];

				Return[
						TableForm[res,TableAlignments->{Center,Center},
							TableHeadings->{None,{"t","|t-(Subscript[t, a])|","|t-(Subscript[t, c])|","|y-(Subscript[y, a])|","|(Subscript[y, a])-(Subscript[y, c])|"}},TableSpacing->{1,7}]
			(*EndReturn*)];

		(*EndModule*)];
(*EndCompareRK4Dec*)

Clear[IntFracS]
IntFracS[num_,IntDig_,FracDig_]:=

		(*
		given the number in num, the number is converted to a string containing IntDig 
			integer digits and FracDig decimal digits: IntFracS[33.4545, 3, 8] would 
			produce "033.45450000" .
		*)

		Module[

				{exp1,IntBuf,NumList,str1},


				{NumList,exp1}=RealDigits[FractionalPart[Abs[num]],10,FracDig+3];
				IntBuf=StringJoin[Map[ToString,Table[0,{IntDig}]]];
				str1=StringTake[IntBuf<>ToString[IntegerPart[Abs[num]]],{-IntDig,-1}]<>"."<>
							Apply[StringJoin,Map[ToString,Take[Join[Table[0,{-exp1}],NumList],FracDig]]];
				If[Sign[num]===-1,
						If[StringTake[str1,1]==="0",
								str1=StringReplacePart[str1,"-",{1,1}]
						(*EndIf*)]
				(*EndIf*)];
				Return[str1];

		(*EndModule*)];
(*EndIntFracS*)

Results

SetDigitsND[24]
SetRoundingND["RoundToZero"]

Below, t and y represent, respectively, time and the solution of the differential equation. 
	The subscripts blank, a, and c represent, respectively, a true value for the situation/problem, 
	a true value for the algorithm, a computer arithmetic emulation value.

24 digits, step size 4, number of steps 3

CompareRK4[y,f,0,y0,4,3]

24 digits step size 1/4, number of steps 48 (not all times are shown)

CompareRK4Dec[y,f,0,y0,1/4,48]

SetDigitsND[42]

42 digits, step size 4, number of steps 3

CompareRK4[y,f,0,y0,4,3]

42 digits, step size 1/4, number of steps 48 (not all times are shown)

CompareRK4Dec[y,f,0,y0,1/4,48]

SetDigitsND[41]

41 digits, step size 1/4, number of steps 48 (not all times are shown)

CompareRK4Dec[y,f,0,y0,1/4,48]

The End