#!/usr/bin/perl -w

# dwhite@nist.gov 3/23/04  - based on BASIC code
# in "Astronomical Algorithms" 2nd ed. by Jean Meeus, pg 17-19

use strict;
use vars qw( $opt_v $x $j $pi );

$opt_v = shift;
if ((defined $opt_v) && ($opt_v ne "-v")) {
    die "Usage: $0 [-v]\n\tthis calculates the internal accuracy\n\tof the programming language.  -v is verbose\n";
}

# log info about this machine
# print "$0 : you appear to be using \n";
# print "\tVENDOR   \t\"$ENV{VENDOR}\"\n";
# print "\tMACHTYPE \t\"$ENV{MACHTYPE}\"\n";
# print "\tOSTYPE   \t\"$ENV{OSTYPE}\"($^O)\"\n";
# print "$0 : you appear to be using \n\tVENDOR   \t\"$ENV{VENDOR}\"\n\tMACHTYPE \t\"$ENV{MACHTYPE}\"\n\tOSTYPE   \t\"$ENV{OSTYPE}\"($^O)\"\n\tHOSTTYPE \t$ENV{HOSTTYPE}\" \n";

# start of code block 1
$x = 1 ;
$j = 0 ;
$x *= 2 ;
if (defined $opt_v) { print "Testing significant bits, significant digits...\n"; }
while (($x + 1) != $x) {
if (defined $opt_v) { print "\t$j\t$x\n"; }
$j++ ;
$x *= 2 ;
}
if (defined $opt_v) { print "\t$j\t$x\n"; }
print "\n", $j , " significant bits in mantissa of floating number\n" , int($j * 0.30103) , " significant digits in a decimal number (", $j * 0.30103 , ")\n above is only for SIMPLE ARITHMETICS, not trig functions!\n";
# end of code block 1

# simple pi check
$pi = atan2(1,1) * 4 ;
print "\narctan(1)*4 [aka pi] is $pi \n check this against     3.141592653589793238462643383279502\n";

print "\nthe 2nd column here should NOT list diverging numbers...\n";
print "\$x=1.0/3.0;  for(\$j=1;\$j<31;\$j++) { \$x = (9*\$x+1)*\$x-1 ; }\n";
# start of code block 2
$x = 1.0/3.0 ;
for($j=1;$j<31;$j++) {
    $x = (9*$x+1)*$x-1 ;
    if (defined $opt_v) { print "$j\t$x\n"; }
    if ((!defined $opt_v) && ($j % 6 == 1)) { print "$j\t$x\n"; }
}
print "However, they probably will diverge on your machine.\n";
# end of code block 2

# another simple check pg. 18
print "\nsquare 1.0000001 27 times...\n";
$x = 1.0000001 ;
for($j=0;$j<27;$j++) {
    $x *= $x ;
}
print "674530.4707 is the expected result to 10 sigdigs,\n$x is your calculated result\n";

# more quick checks pg. 19
print "\n\$x = 4.34 ; \$j = int(100*(\$x-int(\$x))); ";
$x = 4.34 ;
$j = int(100*($x-int($x)));
print "\n\$j should be 34 and you calculated \$j to be $j \n";

# the order of the addition in these 2 tests can be a factor
print "\n2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 0 and you get ";
print 2 +0.2 +0.2 +0.2 +0.2 +0.2 -3 , "\n";

print "0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0 and you get ";
print 0.2 +0.2 +0.2 +0.2 +0.2 +2 -3 , "\n";

print "2 + (5 * 0.2) - 3 = 0 and you get ";
print 2 + (5*0.2) -3 , "\n";

# start of code block 3
for($j=0;$j<=100;$j+=0.1) {
    $x = $j;
}
print "\nfor(\$j=0;\$j<=100;\$j+=0.1) { \$x = \$j; } ";
print "\n$x should equal 100\n";
# end of code block 3

# page 20
# another simple check - try big numbers or 2**x-1
# 255, 65535, 16777215, 4294967295, 1099511627775 = 2**(8*x)
print "\nSimple test that may fail with large numbers:\n";
my @a = ( 255, 65535, 16777215, 4294967295, 1099511627775, 281474976710655 );
for ($j=0;$j<=$#a;$j++) {
    my $b = $a[$j]/10 ; my $c = 10 * $b ;
    print "\$a = $a[$j]; \$b = \$a/10 ; \$c = 10 * \$b ; \$a -\$c = ", $a[$j]-$c , " = 0\n";
}

# page 20 
print "\nsqrt(25)-5 = ", sqrt(25)-5 , " : should be 0\n";
print "sqrt(25)-int(sqrt(25)) = ", sqrt(25)-int(sqrt(25)) , " : should be 0\n";


print "\n";
exit;

__END__

original code is

# start of code block 1
10 x=1
20 j=0
30 x=x*2
40 if x+1 <> x then 60
50 goto 80
60 j=j+1
70 goto 30
80 print j, j*0.30103
90 end
# end of code block 1

# start of code block 2
10 x=1/3
20 for j=1 to 30
30 x=(9*x+1)*x-1
40 print j,x
50 next j
60 end
# end of code block 2

# start of code block 3
10 for i=0 to 100 step 0.1
20 u=i
30 next i
40 print u
50 end
# end of code block 3