Perhaps you have a book containing trig functions, and you compare the answer to the value in the table. How do you know the table is correct? It was probably generated by computer program. How do you know its computer program is more accurate than TRIG.Cos? You could check the result with your pocket calculator. Are you going to take the word of a $29 piece of plastic containing a 4-bit processor over a real 16- or 32-bit computer?
When it comes right down to it, it becomes a question of trust. You will accept the math tables in the handbook because it is published by John Wiley & Sons, and therefore can not contain any errors. Perhaps you trust the pocket calculator because it was made by a well known manufacturer with a reputation you can trust. Your confidence may be boosted by the fact that the math tables and the calculator both agree. For one reason or another, you put your faith in something.
Listing 74 shows the program I used to test the TRIG.Cos function with input values from -5.0 degrees to +370.0 degrees in 1 degree steps. That's almost 400 values, but it only takes the computer a few seconds to generate them all. (You can easily test more values by changing the loop increment to 0.1 degree or 0.001 degree if you like.)
The program creates a data file listing the 376 input angles and the function result of each angle on the same line. At 50 lines per page, that's almost 8 pages of data. You can quickly scan the data, looking at angles that are multiples of 30 or 45 degrees and verify that those angles give the expected results. Then you can look for obviously wild values. If you don't find any, that gives you some confidence, but it doesn't prove all the values are correct. What are we going to do?
The trick is to try to obtain the result two different ways and compare the results. If they agree, the answer is almost certainly correct. If they disagree, you need to find a third method to independently check the answer.
One of the reasons for using the identity Cos(THETA) = Sin(THETA+PI_OVER_TWO) for the Meridian and Alsys versions of the TRIG package was to use a different way computing the cosine than the DEC version of the TRIG package used.
The DEC version is the one I trust. Digital Equipment Corporation has supplied scientific run-time libraries to FORTRAN programs for years, and I'm sure they have discovered an outstanding, accurate, method for computing all the trig functions. The Cos function in the TRIG package I use on the DEC machine simply calls Cosd function in the VAX/VMS math library.
The THETA+PI_OVER_TWO method I used in the Meridian and Alsys versions is mathematically inferior because it introduces round-off error. Suppose I want to find the cosine of an angle THETA near -90 degrees. Perhaps the angle THETA is exactly -1.57 radians. The cosine of that angle is 7.9632686e-4 (according to my pocket calculator). Suppose the value I use to approximate PI_OVER_TWO is 1.5708 radians. When I compute Sin(-1.57 + 1.5708) I get 7.9999986e-4, an error of 3.6729994e-6! That bothers mathematicians. It only bothers engineers if they know the angle THETA was measured to an accuracy better than 3.673 microradians (about 0.00021 degrees).
I ran the Cos_Test program on Meridian, Alsys, and DEC machines, expecting to get slightly different results. I wanted to know the magnitude of the difference. I also wanted to know what input angle gives the poorest results.
Cos_Test always creates an output file called COS.DAT. I could have modified the source code every time I moved the Cos_Test from one computer to another, or I could have made the program ask the user what file to put the output data in, but I decided it simpler to let it put the data in COS.DAT and use an operating system command to rename it to DECCOS.DAT, MERCOS.DAT, or ALSCOS.DAT after it was created.
Figure 40 shows some of the results I got with the Meridian version of the cosine routine. Eight pages of cosine data isn't very interesting reading, so I've shown a few segments of the output just to give you a sample of what it looks like. The DEC and Alsys results are almost identical, often differing by 1 in the least significant digit.
If a person compares two eight page printouts, and doesn't notice any major errors, that gives you some measure of confidence. People get careless sometimes, so you can't really be sure the printouts match. You could have more confidence if a computer meticulously compared both listings, line by line.
MS DOS has a file comparison utility that I hoped to used to compare them. When I tried it I found it wasn't totally successful. So many lines differed by one digit that the utility program got discouraged and gave up. That meant I had to write a program to read two files and find the differences. I could have made the program compute the average difference and standard deviation, but I was more interested in the maximum error.
The COS_DIF program is shown in Listing 75. It asks the user for two file names and compares them, line by line. The first entry on each line (the angle tested) should be identical. If it isn't, then one of the files has been corrupted. The program ends with an error message telling where the error was detected. Otherwise, it goes through the whole file, keeping track of the maximum errors in the positive and negative directions. When it is finished it prints the maximum differences. When I used it to compare the Meridian results to the DEC "truth", I got the results shown in Figure 41. Results comparing Alsys to DEC outputs were similar.
Now I have confidence in the Meridian and Alsys versions of the TRIG.COS function, because I have compared them to a standard I can trust (a routine in the DEC math library).
In a perfect world you would always get exactly what you started with. Unfortunately round-off errors give you a result that is nearly equal to the original value, which makes comparison more difficult. You can take two approaches to when comparing the results. (1) You can establish an acceptable error, and declare any difference between the original value and the reconstructed value smaller than this threshold value to be OK. (2) You can measure the difference between the original and reconstructed values, and keep track of the maximum positive and negative errors.
The first method is most useful for situations where you know how accurate you need to be. For example, suppose the COORDINATES package is going to be used in a program that guides a missile to a target. If the missile warhead has a lethal radius of X feet, then the COORDINATES package must not introduce any inaccuracy resulting in a miss distance greater than X feet. It would probably be a good idea to be conservative, and set the threshold at X/2 feet, or X/10 feet, or whatever is appropriate to the application.
The second method is better for situations where you want to know what is the best you can do. For example, you could use the second method to find the maximum error introduced by the COORDINATES package. If the calculations can be off by as much as Y feet, then you will have to design the warhead to have a lethal radius of at least Y feet. It would probably be a good idea to be conservative, and make it 2Y feet or 10Y feet.
Both methods have weakness, especially in situations where there is one really wild point and many points that just barely fail. Suppose there is one input condition that gives absolutely crazy results (perhaps a location near the origin where a division by a value near 0 occurs), and several other points that are moderately over threshold. The first method will tell you many points failed, but won't tell you the magnitude of the one spectacular failure. You might look at four or five points, see they are all just barely out of tolerance, and decide it's close enough. The second method will tell you the magnitude of the failure at that one awful point, but not tell you that there were a dozen other points that were 20% over the threshold. You might look at the one wild point, decide it is a pathological case that can't occur in normal operation, and think everything is OK.
The even functions, such as square root, are partially symmetrical. You can square 2 to get 4 and take the square root of 4 to get 2 again, but the method doesn't work when the original value is -2. This doesn't mean you can't use the method on partially symmetrical functions, you just have to be more careful. When testing the square root, for example, you could test it with positive values expecting to get the original value back. Then test a second time with all negative values, expecting to get the absolute value back. When testing trigonometric functions, you may want to test one quadrant at a time because inverse functions return mirror images of the input angle in certain quadrants. (The Arcsine of the Sine of 100 degrees is 80 degrees).
Inverse testing may leave some gaps in the test suite. The square root test described above won't tell you what will happen if you try to take the square root of a negative number because all the input values were created by squaring a positive or negative number, the result of which can never be negative. Furthermore, if you are testing a 32-bit integer square root, the test cases will be sparse at the higher values. If you square 46,339 it will yield the test input 2,147,302,921. Squaring 46,340 results in 2,147,395,600. There are 92,679 values between 2,147,302,921 and 2,147,395,600 that can never be tested using this approach.
This sort of situation doesn't lend itself to the inverse function testing we just discussed. You can put noisy data in the input of a filter and get clean data out, but you usually can't shove the clean data in the output of the filter to reproduce the same noisy input.
If you have an established algorithm that is believed to work, you can put faith in its results, using the first method we discussed. You just record some noisy, real-world data, process it with both algorithms, and compare the results. But if you are developing an entirely new algorithm, then you don't have an old one to tell you what the right answer is. If you do have an old one, and it gives different results than the new one, then it might test your faith. Are you really sure the old algorithm gives the right answers? There is likely to be some doubt in your mind.
In situations like these I like to use a different approach. I manufacture some realistic input data by starting with clean data and adding artificially generated, known noise. Then when I use an algorithm to separate the signal from the noise, I know exactly what the correct answer should be.
The description of a valid input signal for your algorithm depends entirely on the problem, but generally you can create a pure, valid input signal by adding several sine waves with various amplitudes, frequencies, and phase differences. You can compute the value of this input signal at the moments of interest and put it in a file called CLEAN.DAT. Use this as the input to your algorithm and store the resulting output in TRUTH.DAT.
The next step is to create some realistic corruption. Maybe the signal is corrupted by frequencies that are out of band, or periodic noise spikes. If the signal is likely to be corrupted by gaussian noise, you can use the RANDOM_NUMBERS.Noise function to create a noise waveform. If appropriate you can add a bias and/or filter the noise waveform. You should have a good idea of the characteristics of the noise in your system (if you don't, you haven't done your homework), and should be able to model it. Simulate a representative noise waveform and store it in NOISE.DAT.
The third step is to add each element in CLEAN.DAT to the corresponding element in NOISE.DAT, and store the result in INPUT.DAT. You now have a realistic, noisy input signal for testing your algorithm.
Finally, process the noisy input and store the result in OUTPUT.DAT. You can compare OUTPUT.DAT to TRUTH.DAT to see how well the algorithm performed in the presence of noise. You can do this as often as you want with different kinds of input signals and different kinds of noise.