Category : Files from Magazines
Archive   : PCTJ8804.ZIP
Filename : ACCURACY.PAS
(c) 1987, PC Tech Journal and Ziff Communications Co.
written by Jim Roberts.
Uncomment the correct constant declarations for your compiler.
Insert correct version number in COMPIL variable,
name of hardware system & math coprocessor in MACHIN variable.
}
program ACCURACY ;
{$N+}
type
{for Turbo 4.0}
real = extended;
accarray = array[1..5] of real ;
{for Turbo}
const
COMPIL : string[18] = 'Turbo Pascal 4.0' ;
MACHIN : string[18] = 'Your hardware' ;
MINERR : real = 1.0E-17 ;
LOGMIN : real = 17.0 ;
N = 10 ;
Y : real = 1.0 ;
STEP : real = 0.2 ;
ITER : integer = 20 ;
ITERTRIG : integer = 5 ;
LOG10E : real = 0.43429448190325182765 ;
PI : real = 3.14159265358979323846 ;
PIO2 : real = 1.57079632679489661923 ;
ROOT2 : real = 1.4142135623730950488 ;
ROOT3 : real = 1.7320508075688772935 ;
SQRTO2 : real = 0.7071067811865475244 ;
{for Pascal-2
const
COMPIL = 'Oregon Pascal-2 2.2A' ;
MACHIN = 'Your hardware' ;
MINERR = 1.0E-17 ;
LOGMIN = 17.0 ;
N = 10;
Y = 1.0 ;
STEP = 0.2 ;
ITER = 20 ;
ITERTRIG = 10 ;
LOG10E = 0.43429448190325182765 ;
PI = 3.14159265358979323846 ;
PIO2 = 1.57079632679489661923 ;
ROOT2 = 1.4142135623730950488 ;
ROOT3 = 1.7320508075688772935 ;
SQRTO2 = 0.7071067811865475244 ;
}
var
i, j, k, l, m, ntest : integer ;
a,b,c : array[1..N,1..N] of real;
sum, X : real ;
xx, zz, quot : real ;
a0, a1, d0, d1, frac : real ;
p, p2 : real ;
th, err, logerr, diverr, val, funct: accarray ;
testerr : array[1..10] of real ;
toterr : real ;
procedure filla;
var
i,j : integer;
begin
for i := 1 to N do
for j := 1 to N do
if i <> j then a[i,j] := Y
else a[i,j] := X + Y;
end;
procedure fillb;
var
i,j : integer;
f,d : real ;
begin
f := X + N*Y ;
d := 1.0/(X*f) ;
for i := 1 to N do
for j := 1 to N do
if i <> j then b[i,j] := -Y*d
else b[i,j] := (-Y+f)*d ;
end;
procedure fillc;
var
i,j : integer;
begin
for i := 1 to N do
for j := 1 to N do
c[i,j] := 0.0;
end;
procedure mult;
var
i, j, k : integer;
begin
for i := 1 to N do
for j := 1 to N do
begin
sum := 0.0;
for k := 1 to N do
sum := sum + a[i,k] * b[k,j];
c[i,j] := sum;
end;
end;
procedure sumit ;
var
i, j : integer;
begin
sum := 0.0;
for i := 1 to N do c[i,i] := c[i,i] - 1.0 ;
for i := 1 to N do
for j := 1 to N do
sum := sum + abs(c[i,j]);
end;
function osgn(n:integer) : real ; {negative if n odd}
begin
if n = 2*(n div 2) then osgn := 1.0
else osgn := -1.0 ;
end;
procedure header ;
begin
writeln('ACCURACY: double precision reals tester: ', COMPIL, '; ',MACHIN) ;
writeln(' V 1.7 (c) 1987, PC Tech Journal and Ziff Communications Co.');
writeln(' written by Jim Roberts.');
writeln('Test 1 checks multiplication and addition, then division and subtraction.');
writeln('Test 2 measures the accuracy of the Turbo trig functions (there is no tan()).');
writeln('Test 3 finds the truncation error in some exponential and sqrt identities.');
writeln('ACCURACY is the rounded negative log of error. Program may exit abnormally.');
writeln('NOTE: an increase of 1 in the rating means a factor of TEN less accurate.');
writeln('Interpretation <0.0 - 0.5 => Excellent 1.0 - 1.5 => Fair');
writeln(' of RATING: 0.5 - 1.0 => Good 1.5 < => Poor');
writeln ;
writeln(' TESTS ACCURACY RATING ');
end;
procedure arith ;
{TEST 1: well-conditioned combinatorial matrix times its inverse.}
var
i, k, l, m : integer ;
begin
zz := 0.30 ; {factor used to control decrease of condition of matrix}
for l := 1 to 5 do
begin
xx := zz*(3-l) ;
X := exp(xx/LOG10E); {slowly decreases conditioning}
filla ; fillb ; fillc ;
mult ; sumit ;
err[l] := sum/sqr(N) ; {error is average absolute error per element}
if err[l] > MINERR then logerr[l] := -ln(err[l]) * LOG10E
else logerr[l] := LOGMIN;
testerr[1] := testerr[1] + (LOGMIN - logerr[l]) ;
end ;
testerr[1] := testerr[1]/5.0 ;
write('#1a: 10x10 matrix ');
for i := 1 to 5 do write(logerr[i]:5:1);
writeln(' ',testerr[1]:6:2) ;
{ infinite product for 1-x: run in reverse to test division }
sum := 0.0 ;
for l := 1 to 5 do
begin
xx := (1 - l) / 4.0 ;
zz := exp((xx-2.0)/LOG10E); {increases number of factors for convergence}
xx := 1.0 - zz ; { zz ÷ .01 => loss of 2 sig figs }
{The following formula for the number of factors is designed to give
sufficient accuracy, while avoiding underflow in the powers of xx.
It gives a more uniform computation from compiler to compiler.}
m := 12+l ;
quot := 1.0 ;
for k := 1 to m do
begin
quot := quot / (1.0 + xx) ;
xx := xx * xx ;
end ;
err[l] := abs(1.0 - quot/zz)*0.01 ;
{ factor of 0.01 to compensate for cancellation errors }
if err[l] > MINERR then diverr[l] := -ln(err[l]) * LOG10E
else diverr[l] := LOGMIN;
sum := sum + (LOGMIN - diverr[l]) ;
logerr[l] := diverr[l] ; {needed for later average}
end ;
sum := sum / 5.0 ;
write('#1 : infinite product ');
for i := 1 to 5 do write(diverr[i]:5:1);
writeln(' ',sum:6:2) ;
{ test continued fraction for tangent against exact values for five angles:
this is a test of division and subtraction, not of the tangent.}
th[1] := PI/12.0 ;
th[2] := PI/6.0 ;
th[3] := PI/4.0 ;
th[4] := PI/3.0 ;
th[5] := 5.0*PI/12.0 ;
val[1] := 2.0 - ROOT3 ;
val[2] := 1.0 / ROOT3 ;
val[3] := 1.00 ;
val[4] := ROOT3 ;
val[5] := 2.0 + ROOT3 ;
sum := 0.0 ;
m := 8 ; { this number of iterations gives sufficient accuracy }
for l := 1 to 5 do
begin
a0 := 2.0 * m + 1.0 ;
p2 := th[l] ;
p := sqr(p2) ;
d0 := a0 - p / (a0 + 2.0) ;
for k := 1 to m do
begin
a1 := a0 - 2.0 ;
d1 := a1 - p / d0 ;
a0 := a1 ;
d0 := d1 ;
end ;
frac := p2 / d0 ;
funct[l] := frac ;
end ;
for l := 1 to 5 do
begin
err[l] := abs(1.0 - val[l]/funct[l]) ;
if err[l] > MINERR then diverr[l] := -ln(err[l]) * LOG10E
else diverr[l] := LOGMIN;
sum := sum + (LOGMIN - diverr[l]);
end ;
sum := sum / 5.0 ;
write('#1 : continued fraction ');
for i := 1 to 5 do write(diverr[i]:5:1);
writeln(' ',sum:6:2) ;
for i := 1 to 5 do
begin
logerr[i] := 0.5*(logerr[i] + diverr[i]) ;
testerr[2] := testerr[2] + (LOGMIN - logerr[i]);
end ;
testerr[2] := testerr[2]/5.0 ;
write('#1b: division average ');
for i := 1 to 5 do write(logerr[i]:5:1);
writeln(' ',testerr[2]:6:2) ;
end;
procedure trig ; {TEST 2: first, errors in some sine identities }
var
i, j, k, l : integer ;
begin
for l := 1 to 5 do logerr[l] := 0.0 ;
for j := 1 to ITERTRIG do
begin
p := j - 1 ;
th[1] := PI/12.0 + p*PI ;
th[2] := PI/6.0 + p*PI ;
th[3] := PI/4.0 + p*PI ;
th[4] := PI/3.0 + p*PI ;
th[5] := 5.0*PI/12.0 + p*PI ;
val[1] := osgn(j-1)*ROOT2*(ROOT3-1.0)*0.25 ;
val[2] := osgn(j-1)*0.5 ;
val[3] := osgn(j-1)*SQRTO2 ;
val[4] := osgn(j-1)*0.5*ROOT3 ;
val[5] := osgn(j-1)*ROOT2*(ROOT3+1.0)*0.25 ;
for l := 1 to 5 do funct[l] := sin(th[l]) ;
for l := 1 to 5 do
begin
err[l] := abs(1.0 - val[l]/funct[l]) ;
if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
else logerr[l] := logerr[l] + LOGMIN ;
end;
end;
for l := 1 to 5 do logerr[l] := logerr[l] / ITERTRIG ;
for l := 1 to 5 do testerr[3] := testerr[3] + (LOGMIN - logerr[l]);
testerr[3] := testerr[3]/5.0 ;
write('#2a: sine function ');
for i := 1 to 5 do write(logerr[i]:5:1);
writeln(' ',testerr[3]:6:2) ;
{ compare sin()/cos() with exact values for tan() }
for l := 1 to 5 do logerr[l] := 0.0 ;
for j := 1 to ITERTRIG do
begin
p := j - 1 ;
th[1] := PI / 12.0 + (j-1)*PI ;
th[2] := PI / 6.0 + (j-1)*PI ;
th[3] := PI / 4.0 + (j-1)*PI ;
th[4] := PI / 3.0 + (j-1)*PI ;
th[5] := 5.0 * PI / 12.0 + (j-1)*PI ;
val[1] := 2.0 - ROOT3 ;
val[2] := 1.0 / ROOT3 ;
val[3] := 1.00 ;
val[4] := ROOT3 ;
val[5] := 2.0 + ROOT3 ;
for l := 1 to 5 do funct[l] := sin(th[l])/cos(th[l]) ;
for l := 1 to 5 do
begin
err[l] := abs(1.0 - val[l]/funct[l]) ;
if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
else logerr[l] := logerr[l] + LOGMIN ;
end ;
end;
for l := 1 to 5 do logerr[l] := logerr[l] / ITERTRIG ;
for l := 1 to 5 do
testerr[4] := testerr[4] + (LOGMIN - logerr[l]) ;
testerr[4] := testerr[4]/5.0 ;
write('#2b: tangent|sin()/cos()');
for i := 1 to 5 do write(logerr[i]:5:1) ;
writeln(' ',testerr[4]:6:2) ;
{ compare arctan() with tan() for consistency }
for l := 1 to 5 do logerr[l] := 0.0 ;
for j := 1 to ITER do
begin
for l := 1 to 5 do th[l] := (5*j+l-5)*PIO2/(5*ITER+1) ;
for l := 1 to 5 do val[l] := sin(th[l])/cos(th[l]) ;
for l := 1 to 5 do funct[l] := arctan(val[l]) ;
for l := 1 to 5 do
begin
err[l] := abs(1.0 - th[l]/funct[l]) ;
if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
else logerr[l] := logerr[l] + LOGMIN ;
end;
end;
for l := 1 to 5 do logerr[l] := logerr[l] / ITER ;
for l := 1 to 5 do
testerr[5] := testerr[5] + (LOGMIN - logerr[l]);
testerr[5] := testerr[5]/5.0 ;
write('#2c: arctan function ');
for i := 1 to 5 do write(logerr[i]:5:1) ;
writeln(' ',testerr[5]:6:2) ;
end;
procedure transc ; {TEST 3: ln() and exp() for consistency}
var
i, j, l : integer ;
begin
for l := 1 to 5 do logerr[l] := 0.0 ;
for j := 1 to ITER do
begin
for l := 1 to 5 do th[l] := (5*j+l-5)*STEP ;
for l := 1 to 5 do val[l] := exp(th[l]) ;
for l := 1 to 5 do funct[l] := ln(val[l]) ;
for l := 1 to 5 do
begin
err[l] := abs(1.0 - th[l]/funct[l]) ;
if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
else logerr[l] := logerr[l] + LOGMIN;
end;
end;
for l := 1 to 5 do logerr[l] := logerr[l] / ITER ;
for l := 1 to 5 do
testerr[6] := testerr[6] + (LOGMIN - logerr[l]);
testerr[6] := testerr[6]/5.0 ;
write('#3a: ln() & exp() ') ;
for i := 1 to 5 do write(logerr[i]:5:1) ;
writeln(' ',testerr[6]:6:2) ;
end;
procedure roots ; { sqrt() and sqr() identities }
var
i, j, l : integer ;
begin
for l := 1 to 5 do logerr[l] := 0.0 ;
for j := 1 to ITER do
begin
for l := 1 to 5 do th[l] := (5*j+l-5)*STEP ;
for l := 1 to 5 do val[l] := sqrt(th[l]) ;
for l := 1 to 5 do funct[l] := sqr(val[l]) ;
for l := 1 to 5 do
begin
err[l] := abs(1.0 - th[l]/funct[l]) ;
if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
else logerr[l] := logerr[l] + LOGMIN;
end;
end;
for l := 1 to 5 do logerr[l] := logerr[l] / ITER ;
for l := 1 to 5 do
testerr[7] := testerr[7] + (LOGMIN - logerr[l]) ;
testerr[7] := testerr[7]/5.0;
write('#3b: sqrt() & sqr() ') ;
for i := 1 to 5 do write(logerr[i]:5:1) ;
writeln(' ',testerr[7]:6:2) ;
end;
begin { program ACCURACY begins here }
header ;
for i := 1 to 10 do testerr[i] := 0.0 ;
arith ;
trig ;
transc ;
roots ;
ntest := 7 ;
toterr := 0.0;
for i := 1 to ntest do toterr := toterr + testerr[i] ;
toterr := toterr/ntest ;
writeln('Overall rating: ',toterr:6:2);
end.
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/