1 /***********************************************************************
2 @C-file{
3 author = "Nelson H. F. Beebe",
4 version = "2.00",
5 date = "10 December 2000",
6 time = "08:44:08 MST",
7 filename = "ndiff.c",
8 copyright = "Copyright (c) 2000 Nelson H. F. Beebe. This
9 code is licensed under the GNU General Public
10 License, version 2 or later.",
11 address = "Center for Scientific Computing
12 University of Utah
13 Department of Mathematics, 322 INSCC
14 155 S 1400 E RM 233
15 Salt Lake City, UT 84112-0090
16 USA",
17 telephone = "+1 801 581 5254",
18 FAX = "+1 801 585 1640, +1 801 581 4148",
19 URL = "http://www.math.utah.edu/~beebe",
20 checksum = "31885 1715 5905 53017",
21 email = "beebe@math.utah.edu, beebe@acm.org,
22 beebe@computer.org, beebe@ieee.org
23 (Internet)",
24 codetable = "ISO/ASCII",
25 keywords = "numerical file differencing",
26 supported = "yes",
27 docstring = "This program compares two putatively similar
28 files, ignoring small numeric differences.
29 Complete documentation can be found in the
30 accompanying UNIX manual page file,
31 ndiff.man.
32
33 Usage:
34 ndiff [ -? ] [ -abserr abserr ]
35 [ -author ] [ -copyright ]
36 [ -fields n1a-n1b,n2,n3a-n3b,... ]
37 [ -help ] [ -logfile filename ]
38 [ -minwidth nnn ]
39 [ -outfile filename ]
40 [ -precision number-of-bits ]
41 [ -quick ] [ -quiet ]
42 [ -relerr relerr ]
43 [ -separators regexp ] [ -silent ]
44 [ -version ] [ -www ] infile1 infile2
45
46 The checksum field above contains a CRC-16
47 checksum as the first value, followed by the
48 equivalent of the standard UNIX wc (word
49 count) utility output of lines, words, and
50 characters. This is produced by Robert
51 Solovay's checksum utility.",
52 }
53 ***********************************************************************/
54
55 /***********************************************************************
56 Coding conventions:
57
58 (1) This program is written in strictly-conforming 1989 ANSI/ISO
59 Standard C (ANSI X3.159-1989, ISO/IEC 9899:1990). No
60 concession is made to the now-obsolete Kernighan & Ritchie style
61 of C programming, since at the time of writing, the 1989
62 language standard is over a decade old. Standard-conforming
63 C compilers have been available since at least the mid-1990s on
64 every commercially-significant desktop, server, mini, and
65 mainframe computer system.
66
67 The new ANSI/ISO C language standard, ISO/IEC 9899:1999,
68 informally known as C99, has been published, but no
69 C99-conformant compilers are yet available (in late 2000) to
70 test this program.
71
72 (2) Functions are defined in strict lexicographic order and written
73 in the 7-bit ASCII character set, which occupies the first 128
74 positions in the ISO8859-1 and Unicode character sets.
75
76 (3) Macros are spelled entirely in uppercase. Global variables are
77 spelled with an initial capital letter. Functions are spelled
78 entirely in lowercase. Embedded underscores are used to make
79 names more readable. Even if lettercase is ignored, there are
80 no clashes between local and global names, and all global names
81 are unique in the first six characters, ignoring lettercase, as
82 required by Section 3.1.2 of ANSI X3.159-1989.
83
84 (4) No preprocessor conditional may depend on any system-specific
85 symbols, except in the associated configure script and the
86 single confix.h file.
87
88 (5) No system header files may be included, except those
89 defined in ANSI/ISO Standard C. All header file #include
90 statements must be bracketed by conditionals #if HAVE_xxx
91 ... #endif. The HAVE_xxx symbols are defined in config.h,
92 and determined automatically by the autoconfigure process.
93
94 (6) The code has been prettyprinted by the Free Software
95 Foundation's GNU Project C prettyprinter, indent, using these
96 options: -bap -bcd -bl -bli0 -di4 -i4 -ncdb -nce -ncs -nfc1
97 -nip -npcs -nsc (normally stored in .indent.pro in the user's
98 home directory).
99
100 (7) In order to be usable for handling files of arbitrary precision,
101 this program uses the GNU Multiple Precision library, GMP, if it
102 is found to be available at installation time.
103
104 Otherwise, it will fall back to using extended (quadruple)
105 precision, using the data type "long double", if that is
106 supported.
107
108 Failing that, it will revert to standard double precision (C
109 data type "double").
110
111 The only operational difference between these three choices is
112 the ability to handle a wider exponent range than normal, and
113 detect numeric differences of small magnitude. The --version
114 option will display the version number, along with the
115 floating-point precision that is supported, and the machine
116 epsilon. The precise numeric value of the latter may depend on
117 an earlier --precision nnn option.
118
119 (8) In the event that the program is built using native precision,
120 rather than GMP, the code should be relatively insensitive to
121 the precise behavior of native arithmetic, including the choice
122 of base, the number of bits or digits, the exponent range, and
123 the availability of representations of infinity and NaN
124 (not-a-number).
125
126 The author has only systems implementing IEEE 754 floating-point
127 arithmetic available for testing, but there should be no
128 problems with arithmetic in IBM S/390 mainframe format,
129 Compaq/DEC VAX format, or SGI/Cray format.
130
131 Nor should the code be sensitive to the presence of longer
132 intermediate-precision hardware registers, such as on the
133 Intel x86 and Motorola 68K architectures, where registers hold
134 80-bit values in temporary real format (1-bit sign, 15-bit
135 exponent, 64-bit significand with no hidden bit), while memory
136 words hold 64-bit values in double format (1-bit sign, 11-bit
137 exponent, 53-bit significand including a leading hidden 1-bit,
138 except in the case of denormalized values near the underflow
139 limit, where there is no leading 1-bit, and normalization
140 requirements are relaxed).
141
142 Special programming care is taken when necessary to ensure that
143 longer-precision registers do not produce anomalous results,
144 such as in the computation of the machine epsilon.
145
146 ***********************************************************************/
147
148 #include "awklib.h"
149
150 #include "ndiff.h"
151
152 #if !defined(DATE)
153 #define DATE "unknown"
154 #endif
155
156 #if !defined(VERSION)
157 #define VERSION "unknown"
158 #endif
159
160 #if defined(DEBUG)
161 #define DEBUGPRINT(plist) (void)fprintf plist
162 #else
163 #define DEBUGPRINT(plist)
164 #endif
165
166 typedef struct
167 {
168 const char *name;
169 size_t min_length;
170 size_t(*optfun) (const char *next_arg);
171 }
172 option_table;
173
174 static awk_table_t Fields;
175
176 static const char *CNULL = (const char *)NULL;
177 static const char *InFile1 = (const char *)NULL;
178 static const char *InFile2 = (const char *)NULL;
179
180 /***********************************************************************
181 The precise value of the regular expression, NUMBER_PATTERN, to match
182 both an integer and a floating-point number is critical, and
183 documented in the accompanying manual page: it must match not only the
184 awk- and C-style -nnn, -n.nnn, and -n.nnne+nn, but also the Fortran
185 styles -nnn, -n.nnn, -n.D+nn, -.nnD+nn, -nD+nn, -n.nnnQ+nn,
186 -n.nnnd+nn, and -n.nnn+nnn. The Fortran forms will be converted by
187 awkfloat() to awk-form. Ada permits an nonsignificant underscore
188 between digits, so we support that too.
189 ***********************************************************************/
190
191 static const char *NUMBER_PATTERN =
192 "^[-+]?([0-9](_?[0-9])*([.]?([0-9](_?[0-9])*)*)?|[.][0-9](_?[0-9])*)([DdEeQq]?[-+]?[0-9](_?[0-9])*)?$";
193
194 static char *Program_Name = "ndiff";
195
196 static int Done = 0;
197 static int Error_Count = 0;
198 static int QUIET = 0;
199 static int SILENT = 0;
200
201 static fp_t ABSERR;
202 static fp_t Macheps;
203 static fp_t Max_Abserr;
204 static fp_t Max_Relerr;
205 static fp_t RELERR;
206 static fp_t This_Abserr;
207 static fp_t This_Relerr;
208 static fp_t Zero;
209
210 static size_t Max_Abserr_NF = 0;
211 static size_t Max_Abserr_NR = 0;
212 static size_t Max_Relerr_NF = 0;
213 static size_t Max_Relerr_NR = 0;
214 static size_t MINWIDTH = 0;
215 static size_t N_Fields = 0;
216 static size_t NRLINE = 0;
217
218 #if HAVE_GMP
219 static unsigned long int Precision = 256;
220 #endif
221
222 int main(int argc, char *argv[]);
223 static void cleanup(void);
224 static void compare_all(const char *f1line, const char *f2line,
225 const char *f1parts[], const char *f2parts[],
226 size_t n);
227 static void compare_files(const char *infile1, const char *infile2);
228 static void compare_some(const char *f1line, const char *f2line,
229 const char *f1parts[], const char *f2parts[],
230 size_t n);
231 static void compare_some(const char *f1line, const char *f2line,
232 const char *f1parts[], const char *f2parts[],
233 size_t n);
234 static int diff_field(const char *field1, const char *field2, size_t nfield);
235 static void do_args(int argc, char *argv[]);
236 static size_t do_one_arg(const char *this_arg, const char *next_arg);
237 static void error(const char *s1, const char *s2, const char *s3);
238 static void fp_abs(fp_t * result, fp_t a);
239 static void fp_assign(fp_t * result, fp_t a);
240 static void fp_clear(fp_t a);
241 static int fp_cmp(fp_t a, fp_t b);
242 static void fp_div(fp_t * result, fp_t a, fp_t b);
243 static void fp_init(fp_t * a);
244 static void fp_max(fp_t * result, fp_t a, fp_t b);
245 static void fp_min(fp_t * result, fp_t a, fp_t b);
246 static void fp_mul(fp_t * result, fp_t a, fp_t b);
247 static void fp_sub(fp_t * result, fp_t a, fp_t b);
248 static const char *fp_to_string(fp_t a, int ndigits);
249 static int ignore(const char *field);
250 static void initialize(int argc, char *argv[]);
251 static void initialize_fields(const char *fields);
252 static void machine_epsilon(fp_t * result);
253 static void maxrelerr(fp_t * result, fp_t x, fp_t y);
254 static size_t opt_abserr(const char *next_arg);
255 static size_t opt_author(const char *next_arg);
256 static size_t opt_copyright(const char *next_arg);
257 static size_t opt_fields(const char *next_arg);
258 static size_t opt_help(const char *next_arg);
259 static size_t opt_logfile(const char *next_arg);
260 static size_t opt_minwidth(const char *next_arg);
261 static size_t opt_outfile(const char *next_arg);
262 static size_t opt_precision(const char *next_arg);
263 static size_t opt_quick(const char *next_arg);
264 static size_t opt_quiet(const char *next_arg);
265 static size_t opt_relerr(const char *next_arg);
266 static size_t opt_separators(const char *next_arg);
267 static size_t opt_silent(const char *next_arg);
268 static size_t opt_version(const char *next_arg);
269 static size_t opt_www(const char *next_arg);
270 static void print_lines(FILE * file, const char *lines[]);
271 static void report_difference(const char *f1line, const char *f2line,
272 size_t nfield);
273 static void string_to_fp(fp_t * result, const char *s);
274 static void usage(void);
275 static void warning(const char *s1, const char *s2, const char *s3);
276
277 int
main(int argc,char * argv[])278 main(int argc, char *argv[])
279 {
280 initialize(argc, argv);
281 if ( (InFile1 == CNULL) || (InFile2 == CNULL) )
282 usage();
283 else
284 compare_files(InFile1, InFile2);
285 cleanup();
286
287 return ((Error_Count == 0) ? EXIT_SUCCESS : EXIT_FAILURE);
288 }
289
290
291 static void
cleanup(void)292 cleanup(void)
293 {
294 /*******************************************************************
295 Free dynamically-allocated memory prior to exiting.
296 *******************************************************************/
297
298 fp_clear(ABSERR);
299 fp_clear(Macheps);
300 fp_clear(Max_Abserr);
301 fp_clear(Max_Relerr);
302 fp_clear(RELERR);
303 fp_clear(This_Abserr);
304 fp_clear(This_Relerr);
305 fp_clear(Zero);
306
307 awk_free_table(Fields.table, N_Fields);
308 awk_terminate();
309 }
310
311
312 static void
compare_all(const char * f1line,const char * f2line,const char * f1parts[],const char * f2parts[],size_t n)313 compare_all(const char *f1line, const char *f2line, const char *f1parts[],
314 const char *f2parts[], size_t n)
315 {
316 /*******************************************************************
317 Compare all fields in f1line and f2line, assuming that they have
318 already been split into n parts in f1parts[] and f2parts[].
319
320 If any fields differ, print a diff-style report, and increment
321 global variable Error_Count,
322 *******************************************************************/
323
324 size_t k;
325
326 for (k = 1; k <= n; ++k)
327 {
328 if (diff_field(f1parts[k], f2parts[k], k) != 0)
329 {
330 report_difference(f1line, f2line, k);
331 return;
332 }
333 }
334 }
335
336
337 static void
compare_files(const char * infile1,const char * infile2)338 compare_files(const char *infile1, const char *infile2)
339 {
340 /*******************************************************************
341 Compare all lines in two files, printing a diff-style report of
342 differences. If any numeric differences have been found, print a
343 one-line report of which matching line had the largest numeric
344 difference. Finally, print a diagnostic if the files differ in
345 length.
346 *******************************************************************/
347
348 const char *f1line;
349 const char *f2line;
350 const char **f1parts;
351 const char **f2parts;
352 size_t n1;
353 size_t n2;
354
355 f1line = CNULL; /* keep optimizers happy */
356 f2line = CNULL;
357 f1parts = &CNULL;
358 f2parts = &CNULL;
359
360 NRLINE = 0;
361 while ((awk_getline(infile1, &f1line) > 0)
362 && (awk_getline(infile2, &f2line) > 0))
363 {
364 NRLINE++;
365 n1 = awk_split(f1line, &f1parts, CNULL);
366 n2 = awk_split(f2line, &f2parts, CNULL);
367 if (n1 == n2)
368 {
369 if (N_Fields == 0)
370 compare_all(f1line, f2line, f1parts, f2parts, n1);
371 else
372 compare_some(f1line, f2line, f1parts, f2parts, n1);
373 }
374 else
375 report_difference(f1line, f2line, (n1 > n2) ? n1 : n2);
376 awk_free_table(f1parts, n1);
377 awk_free_table(f2parts, n2);
378 awk_free_string(f1line);
379 awk_free_string(f2line);
380 }
381 if (QUIET == 0)
382 {
383 const char *t;
384
385 if ((fp_cmp(ABSERR, Zero) != 0) && (fp_cmp(Max_Abserr, Zero) > 0))
386 {
387 t = fp_to_string(Max_Abserr, 3);
388 (void)printf("### Maximum absolute error in matching lines = %s at line %d field %d\n",
389 t, (int)Max_Abserr_NR, (int)Max_Abserr_NF);
390 awk_free_string(t);
391
392 }
393 if ((fp_cmp(RELERR, Zero) != 0) && (fp_cmp(Max_Relerr, Zero) > 0))
394 {
395 t = fp_to_string(Max_Relerr, 3);
396 (void)printf("### Maximum relative error in matching lines = %s at line %d field %d\n",
397 t, (int)Max_Relerr_NR, (int)Max_Relerr_NF);
398 awk_free_string(t);
399 }
400 }
401 if (awk_getline(infile1, &f1line) > 0)
402 {
403 awk_free_string(f1line);
404 warning("second file is short", "", "");
405 }
406 if (awk_getline(infile2, &f2line) > 0)
407 {
408 awk_free_string(f2line);
409 warning("first file is short", "", "");
410 }
411 awk_close_infile(infile1);
412 awk_close_infile(infile2);
413 }
414
415
416 static void
compare_some(const char * f1line,const char * f2line,const char * f1parts[],const char * f2parts[],size_t n)417 compare_some(const char *f1line, const char *f2line, const char *f1parts[],
418 const char *f2parts[], size_t n)
419 {
420 /*******************************************************************
421 Compare selected fields in f1line and f2line, assuming that they
422 have already been split into n parts in f1parts[] and f2parts[].
423 The globals (N_Fields, Fields[]) define which fields are to be
424 compared.
425
426 If any fields differ, print a diff-style report, and increment
427 global variable Error_Count.
428 *******************************************************************/
429
430 size_t k;
431 size_t m;
432
433 for (k = 1; (k <= N_Fields) && (k <= n); ++k)
434 {
435 m = (size_t) awk_string_to_unsigned_long(Fields.table[k]);
436 if ((m <= n) && diff_field(f1parts[m], f2parts[m], m) != 0)
437 {
438 report_difference(f1line, f2line, m);
439 return;
440 }
441 }
442 }
443
444
445 static int
diff_field(const char * field1,const char * field2,size_t nfield)446 diff_field(const char *field1, const char *field2, size_t nfield)
447 {
448 /*******************************************************************
449 If both fields are identical as strings, return 0.
450
451 Otherwise, if both fields are numeric, return 0 if they are close
452 enough (as determined by the globals ABSERR and RELERR), or are both
453 ignorable (as determined by MINWIDTH), and otherwise return 1.
454
455 Otherwise, return 1.
456
457 The computed absolute and relative errors are saved in global
458 variables (This_Abserr and This_Relerr) for later use in diagnostic
459 reports. These values are always zero for nonnumeric fields.
460 *******************************************************************/
461
462 fp_t f1;
463 fp_t f2;
464 int retval;
465
466 fp_assign(&This_Abserr, Zero);
467 fp_assign(&This_Relerr, Zero);
468 if (strcmp(field1, field2) == 0) /* handle the commonest, and easiest, case first */
469 return (0);
470 else if (awk_match(field1, NUMBER_PATTERN) &&
471 awk_match(field2, NUMBER_PATTERN))
472 {
473 /* Handle MINWIDTH test while the fields are still strings */
474 if (ignore(field1) && ignore(field2))
475 return (0);
476
477 /* Now coerce both fields to floating-point numbers,
478 converting Fortran-style exponents, if necessary. */
479 fp_init(&f1);
480 fp_init(&f2);
481
482 string_to_fp(&f1, field1);
483 string_to_fp(&f2, field2);
484
485 fp_sub(&This_Abserr, f1, f2);
486 fp_abs(&This_Abserr, This_Abserr);
487 maxrelerr(&This_Relerr, f1, f2);
488 fp_abs(&This_Relerr, This_Relerr);
489
490 retval = 0; /* assume fields match, then try to prove they do not */
491
492 if ((fp_cmp(ABSERR, Zero) != 0) && (fp_cmp(This_Abserr, ABSERR) > 0))
493 retval = 1; /* mismatch: absolute error too large */
494 if ((fp_cmp(RELERR, Zero) != 0) && (fp_cmp(This_Relerr, RELERR) > 0))
495 retval = 1; /* mismatch: relative error too large */
496
497 if (retval == 0) /* fields match, so update the */
498 { /* maximum errors in matching lines */
499 if (fp_cmp(This_Abserr, Max_Abserr) > 0)
500 {
501 Max_Abserr_NF = nfield;
502 Max_Abserr_NR = NRLINE;
503 fp_assign(&Max_Abserr, This_Abserr);
504 }
505 if (fp_cmp(This_Relerr, Max_Relerr) > 0)
506 {
507 Max_Relerr_NF = nfield;
508 Max_Relerr_NR = NRLINE;
509 fp_assign(&Max_Relerr, This_Relerr);
510 }
511 }
512
513 DEBUGPRINT((stderr,"DEBUG 0.01: f1 = %s\n", fp_to_string(f1,10)));
514 DEBUGPRINT((stderr,"DEBUG 0.02: f2 = %s\n", fp_to_string(f2,10)));
515 DEBUGPRINT((stderr,"DEBUG 0.03: This_Abserr = %s\n", fp_to_string(This_Abserr,10)));
516 DEBUGPRINT((stderr,"DEBUG 0.04: This_Relerr = %s\n", fp_to_string(This_Relerr,10)));
517 DEBUGPRINT((stderr,"DEBUG 0.05: Max_Abserr = %s\n", fp_to_string(Max_Abserr,10)));
518 DEBUGPRINT((stderr,"DEBUG 0.06: Max_Relerr = %s\n", fp_to_string(Max_Relerr,10)));
519
520 fp_clear(f1);
521 fp_clear(f2);
522 return (retval);
523 }
524 else if (awk_is_NaN(field1) &&
525 awk_is_NaN(field2))
526 return (0);
527 else if (awk_is_negative_infinity(field1) &&
528 awk_is_negative_infinity(field2))
529 return (0);
530 else if (awk_is_positive_infinity(field1) &&
531 awk_is_positive_infinity(field2))
532 return (0);
533 else
534 return (1);
535 }
536
537
538 static void
do_args(int argc,char * argv[])539 do_args(int argc, char *argv[])
540 {
541 /*******************************************************************
542 Scan the command-line arguments and hand them each off to
543 do_one_arg() for processing. Terminate if the Done flag is
544 set, or if errors occurred.
545 *******************************************************************/
546
547
548 size_t k;
549
550 for (k = 1; k < (size_t)argc; ++k)
551 k += do_one_arg(argv[k], argv[k + 1]);
552
553 if (Done || (Error_Count > 0))
554 {
555 cleanup();
556 exit((Error_Count == 0) ? EXIT_SUCCESS : EXIT_FAILURE);
557 }
558 }
559
560
561 static size_t
do_one_arg(const char * this_arg,const char * next_arg)562 do_one_arg(const char *this_arg, const char *next_arg)
563 {
564 /*******************************************************************
565 Match a command-line argument against the option table, and invoke
566 the required processing function on it. If it is not an argument,
567 then it must be one of the command-line input files.
568 *******************************************************************/
569
570 size_t extra;
571 size_t n;
572 size_t m;
573 static option_table options[] = {
574 {"?", 1, opt_help},
575 {"abserr", 1, opt_abserr},
576 {"absolute-error", 4, opt_abserr},
577 {"author", 2, opt_author},
578 {"copyleft", 5, opt_copyright},
579 {"copyright", 1, opt_copyright},
580 {"fields", 1, opt_fields},
581 {"help", 1, opt_help},
582 {"logfile", 1, opt_logfile},
583 {"minwidth", 1, opt_minwidth},
584 {"minimum-width", 4, opt_minwidth},
585 {"outfile", 1, opt_outfile},
586 {"precision", 1, opt_precision},
587 {"quick", 4, opt_quick},
588 {"quiet", 1, opt_quiet},
589 {"relerr", 1, opt_relerr},
590 {"relative-error", 4, opt_relerr},
591 {"separators", 2, opt_separators},
592 {"silent", 1, opt_silent},
593 {"version", 1, opt_version},
594 {"www", 1, opt_www},
595 {(const char *)NULL, 0, (size_t (*)(const char *))NULL}
596 };
597
598 const char *org_arg;
599 const char *arg;
600
601 arg = org_arg = awk_tolower(this_arg);
602 if ((arg[0] == '-') && (arg[1] == '-')) /* GNU/POSIX-style --option */
603 arg += 2;
604 else if ((arg[0] == '-') && (arg[1] != '\0')) /* UNIX-style -option */
605 arg++;
606 else /* not option: expect input filename */
607 {
608 if (InFile1 == CNULL)
609 {
610 InFile1 = this_arg;
611 awk_free_string(org_arg);
612 return (0);
613 }
614 else if (InFile2 == CNULL)
615 {
616 InFile2 = this_arg;
617 awk_free_string(org_arg);
618 return (0);
619 }
620 }
621
622 extra = 0;
623 n = strlen(arg);
624 for (m = 0; options[m].name != CNULL; ++m)
625 {
626 if ((n >= options[m].min_length)
627 && (strncmp(options[m].name, arg, n) == 0))
628 {
629 extra = options[m].optfun(next_arg);
630 awk_free_string(org_arg);
631 return (extra);
632 }
633 }
634 warning("unrecognized option ", this_arg, "");
635 Error_Count++;
636
637 awk_free_string(org_arg);
638 return (0);
639 }
640
641
642 static void
error(const char * s1,const char * s2,const char * s3)643 error(const char *s1, const char *s2, const char *s3)
644 {
645 /*******************************************************************
646 Report a fatal error and terminate.
647 *******************************************************************/
648
649 (void)fprintf(stderr, "FATAL ERROR: %s%s%s\n", s1, s2, s3);
650 cleanup();
651 exit(EXIT_FAILURE);
652 }
653
654
655 static void
fp_abs(fp_t * result,fp_t a)656 fp_abs(fp_t * result, fp_t a)
657 {
658 /*******************************************************************
659 Store the absolute value of the argument in *result, which is
660 assumed to already be properly initialized.
661 *******************************************************************/
662
663 #if HAVE_GMP
664 DEBUGPRINT((stderr,"DEBUG 0.1: a = %s\n", fp_to_string(a,10)));
665 mpf_abs(*result, a);
666 DEBUGPRINT((stderr,"DEBUG 0.2: abs(a) = %s\n", fp_to_string(*result,10)));
667 #else
668 *result = ((a < 0) ? -a : a);
669 #endif
670 }
671
672
673 static void
fp_assign(fp_t * result,fp_t a)674 fp_assign(fp_t * result, fp_t a)
675 {
676 /*******************************************************************
677 Assign a to *result, where both are assumed to already be properly
678 initialized.
679 *******************************************************************/
680
681 #if HAVE_GMP
682 mpf_set(*result, a);
683 #else
684 *result = a;
685 #endif
686 }
687
688
689 static void
fp_clear(fp_t a)690 fp_clear(fp_t a)
691 {
692 /*******************************************************************
693 Clear a floating-point variable.
694 *******************************************************************/
695
696 #if HAVE_GMP
697 mpf_clear(a);
698 #endif
699 }
700
701
702 static int
fp_cmp(fp_t a,fp_t b)703 fp_cmp(fp_t a, fp_t b)
704 {
705 /*******************************************************************
706 Compare two floating-point numbers, a and b, and return -2 if
707 either is a NaN, -1 if a < b, 0 if a == b, and otherwise +1 if
708 a > b.
709 *******************************************************************/
710
711 #if HAVE_GMP
712 int n;
713
714 n = mpf_cmp(a, b);
715
716 if (n < 0)
717 return (-1);
718 else if (n == 0)
719 return (0);
720 else
721 return (+1);
722 #else
723 if (a != a) /* then a is a NaN */
724 return (-2);
725 else if (b != b) /* then b is a NaN */
726 return (-2);
727 else if (a < b)
728 return (-1);
729 else if (a == b)
730 return (0);
731 else
732 return (+1);
733 #endif
734
735 }
736
737
738 static void
fp_div(fp_t * result,fp_t a,fp_t b)739 fp_div(fp_t * result, fp_t a, fp_t b)
740 {
741 /*******************************************************************
742 Return a/b in *result, which is assumed to have been properly
743 initialized.
744 *******************************************************************/
745
746 #if HAVE_GMP
747 mpf_div(*result, a, b);
748 #else
749 *result = a / b;
750 #endif
751
752 }
753
754
755 static void
fp_init(fp_t * a)756 fp_init(fp_t *a)
757 {
758 /*******************************************************************
759 Initialize a floating-point variable.
760 *******************************************************************/
761
762 #if HAVE_GMP
763 mpf_init(*a);
764 #elif HAVE_LONG_DOUBLE
765 *a = 0.0L;
766 #else
767 *a = 0.0;
768 #endif
769 }
770
771
772 static void
fp_max(fp_t * result,fp_t a,fp_t b)773 fp_max(fp_t * result, fp_t a, fp_t b)
774 {
775 /*******************************************************************
776 Store the larger of the two arguments a and b in *result, which
777 is assumed to have been properly initialized. *result may be
778 either a or b.
779 *******************************************************************/
780
781 #if HAVE_GMP
782 if (fp_cmp(a, b) > 0)
783 mpf_set(*result, a);
784 else
785 mpf_set(*result, b);
786 #else
787 *result = (a > b) ? a : b;
788 #endif
789 }
790
791
792 static void
fp_min(fp_t * result,fp_t a,fp_t b)793 fp_min(fp_t * result, fp_t a, fp_t b)
794 {
795 /*******************************************************************
796 Store the smaller of the two arguments a and b in *result, which
797 is assumed to have been properly initialized. *result may be
798 either a or b.
799 *******************************************************************/
800
801 #if HAVE_GMP
802 DEBUGPRINT((stderr,"DEBUG 0.3: a = %s, b = %s\n", fp_to_string(a,10), fp_to_string(b,10)));
803 if (fp_cmp(a, b) < 0)
804 mpf_set(*result, a);
805 else
806 mpf_set(*result, b);
807 DEBUGPRINT((stderr,"DEBUG 0.4: min(a,b) = %s\n", fp_to_string(*result,10)));
808 #else
809 *result = (a < b) ? a : b;
810 #endif
811 }
812
813
814 static void
fp_mul(fp_t * result,fp_t a,fp_t b)815 fp_mul(fp_t * result, fp_t a, fp_t b)
816 {
817 /*******************************************************************
818 Return a*b in *result, which is assumed to have been properly
819 initialized.
820 *******************************************************************/
821
822 #if HAVE_GMP
823 mpf_mul(*result, a, b);
824 #else
825 *result = a * b;
826 #endif
827
828 }
829
830
831 static void
fp_sub(fp_t * result,fp_t a,fp_t b)832 fp_sub(fp_t * result, fp_t a, fp_t b)
833 {
834 /*******************************************************************
835 Store the difference a - b in *result, which is assumed to already
836 be properly initialized.
837 *******************************************************************/
838
839 #if HAVE_GMP
840 DEBUGPRINT((stderr,"DEBUG 0.5: a = %s, b = %s\n", fp_to_string(a,10), fp_to_string(b,10)));
841 mpf_sub(*result, a, b);
842 DEBUGPRINT((stderr,"DEBUG 0.6: fp_sub(*result,a,b) = %s\n", fp_to_string(*result,10)));
843 #else
844 *result = a - b;
845 #endif
846 }
847
848
849 static double
fp_to_double(fp_t a)850 fp_to_double(fp_t a)
851 {
852 /*******************************************************************
853 Convert a to a double precision value and return it.
854 *******************************************************************/
855
856 const char *s;
857 double d;
858
859 s = fp_to_string(a, 35);
860 d = strtod(s, (char **)NULL);
861 awk_free_string(s);
862
863 return (d);
864 }
865
866
867 static const char *
fp_to_string(fp_t a,int n_digits)868 fp_to_string(fp_t a, int n_digits)
869 {
870 /*******************************************************************
871 Convert a to a string in a standard format, and return a pointer
872 to a dynamically-allocated string containing the result. If
873 n_digits is negative or zero, then as many digits as required will
874 be used. Otherwise, the number will have 1 non-zero digit before
875 the decimal point, and n_digits - 1 digits after. The exponent is
876 always signed, and contains at least two digits, with no embedded
877 spaces.
878 *******************************************************************/
879
880 #if HAVE_GMP
881 char *s;
882 char *t;
883 mp_exp_t e;
884
885 if (mpf_sgn(a) == 0) /* mpf_get_str() returns "" for zero */
886 return (awk_dupstr("0.0e+00"));
887 else
888 {
889 t = mpf_get_str((char *)NULL, (mp_exp_t *) & e, 10, 0, a);
890 if ((n_digits > 0) && (n_digits < strlen(t)))
891 t[n_digits] = '\0';
892 s = awk_padstr(t, 1 + 40); /* extra space for e plus 128-bit integer (39 digits + sign) */
893 if (t[0] == '-')
894 (void)sprintf(s, "%c%c.%se%+.2ld",
895 (int)t[0], (int)t[1], (t[2] ? &t[2] : "0"), (long int)(e - 1));
896 else
897 (void)sprintf(s, "%c.%se%+.2ld",
898 (int)t[0], (t[1] ? &t[1] : "0"), (long int)(e - 1));
899 awk_free_string(t);
900 return ((const char *)s);
901 }
902 #else
903 char s[128]; /* big enough for double and quadruple precision */
904
905 #if HAVE_LONG_DOUBLE
906 (void)sprintf(s, "%.*Le", n_digits - 1, a);
907 #else
908 (void)sprintf(s, "%.*e", n_digits - 1, a);
909 #endif
910
911 return (awk_dupstr(s));
912 #endif
913
914 }
915
916
917 static int
ignore(const char * field)918 ignore(const char *field)
919 {
920 /*******************************************************************
921 Return 1 if field is ignorable, because it is shorter than
922 MINWIDTH and appears to be a real number. Otherwise, return 0.
923 *******************************************************************/
924
925 return ((MINWIDTH > 0) &&
926 (strlen(field) < MINWIDTH) &&
927 awk_match(field, "[.DdEeQq]"));
928 }
929
930
931 static void
initialize(int argc,char * argv[])932 initialize(int argc, char *argv[])
933 {
934 /*******************************************************************
935 Process command-line options, and initialize global variables.
936 *******************************************************************/
937
938 Program_Name = (argv[0] == (char*)NULL) ? "ndiff" : argv[0];
939 awk_initialize();
940
941 #if HAVE_GMP
942 mpf_set_default_prec(Precision);
943 #endif
944
945 fp_init(&ABSERR);
946 fp_init(&Macheps);
947 fp_init(&Max_Abserr);
948 fp_init(&Max_Relerr);
949 fp_init(&RELERR);
950 fp_init(&This_Abserr);
951 fp_init(&This_Relerr);
952 fp_init(&Zero);
953
954 fp_assign(&ABSERR, Zero);
955 fp_assign(&Max_Abserr, Zero);
956 fp_assign(&Max_Relerr, Zero);
957 fp_assign(&RELERR, Zero);
958 fp_assign(&This_Abserr, Zero);
959 fp_assign(&This_Relerr, Zero);
960 Max_Abserr_NF = 0;
961 Max_Abserr_NR = 0;
962 Max_Relerr_NF = 0;
963 Max_Relerr_NR = 0;
964
965 awk_new_table(&Fields);
966 N_Fields = 0;
967 Error_Count = 0;
968
969 machine_epsilon(&Macheps);
970
971 do_args(argc, argv);
972
973 if ((fp_cmp(ABSERR, Zero) == 0) && (fp_cmp(RELERR, Zero) == 0)) /* supply default (see man pages) */
974 {
975 #if HAVE_GMP
976 fp_t temp;
977
978 mpf_init_set_d(temp, 8.0); /* temp = 8.0 */
979 mpf_mul(temp, temp, Macheps); /* temp = 8.0 * Macheps */
980 mpf_set_d(RELERR, 1.0e-15); /* RELERR = 1.0e-15 */
981 fp_max(&RELERR, RELERR, temp);
982 mpf_clear(temp);
983 #else
984 fp_max(&RELERR, (fp_t) 1.0e-15, (fp_t) 8.0 * Macheps);
985 #endif
986 }
987 }
988
989
990 static void
initialize_fields(const char * fields)991 initialize_fields(const char *fields)
992 {
993 /*******************************************************************
994 Convert a FIELDS=n1a-n1b,n2,n3a-n3b,... specification to a list
995 of N_Fields numbers in Fields[].
996 *******************************************************************/
997
998 size_t j;
999 size_t k;
1000 size_t m;
1001 size_t n;
1002 const char **numbers;
1003 const char **parts;
1004
1005 numbers = &CNULL; /* keep optimizers happy */
1006 parts = &CNULL;
1007
1008 N_Fields = 0;
1009 n = awk_split(fields, &parts, ",");
1010 for (k = 1; k <= n; ++k)
1011 {
1012 m = awk_split(parts[k], &numbers, "-+");
1013 if (m == 1)
1014 {
1015 if (!awk_match(parts[k], "^[0-9]+$"))
1016 error("non-numeric fields value [", parts[k], "]");
1017 else if (parts[k] == 0)
1018 error("zero fields value [", parts[k],
1019 "]: fields are numbered from 1");
1020 else
1021 awk_add_element(&Fields, ++N_Fields, awk_dupstr(parts[k]));
1022 }
1023 else if (m == 2)
1024 {
1025 if ((!awk_match(numbers[1], "^[0-9]+$")) ||
1026 (!awk_match(numbers[2], "^[0-9]+$")))
1027 error("non-numeric fields range [", parts[k], "]");
1028 else if ((awk_string_to_long(numbers[1]) == 0L) ||
1029 (awk_string_to_long(numbers[2]) == 0L))
1030 error("zero value in fields range [", parts[k],
1031 "]: fields are numbered from 1");
1032 else if (awk_string_to_long(numbers[1]) >
1033 awk_string_to_long(numbers[2]))
1034 error("bad fields range [", parts[k], "]");
1035 else if ((awk_string_to_long(numbers[2]) -
1036 awk_string_to_long(numbers[1]) + 1) > 100)
1037 error("fields range [", parts[k], "] exceeds 100");
1038 else
1039 {
1040 for (j = (size_t) awk_string_to_unsigned_long(numbers[1]);
1041 j <= (size_t) awk_string_to_unsigned_long(numbers[2]);
1042 ++j)
1043 awk_add_element(&Fields, ++N_Fields,
1044 awk_unsigned_long_to_string(
1045 (unsigned
1046 long)j));
1047 }
1048 }
1049 else
1050 error("bad fields range [", parts[k], "]");
1051 awk_free_table(numbers, m);
1052 }
1053 awk_free_table(parts, n);
1054 }
1055
1056
1057 static void
machine_epsilon(fp_t * result)1058 machine_epsilon(fp_t * result)
1059 {
1060 /*******************************************************************
1061 Store the machine epsilon (the smallest number x such that (1 + x)
1062 differs from 1) in *result, which is assumed to have been properly
1063 initialized.
1064 *******************************************************************/
1065
1066 static fp_t x;
1067 static fp_t y;
1068
1069 #if HAVE_GMP
1070 fp_t half;
1071
1072 mpf_init_set_d(half, 0.5); /* half = 0.5 */
1073 mpf_init_set_d(x, 1.0); /* x = 1.0 */
1074 mpf_init(y);
1075 mpf_mul(y, half, x); /* y = 0.5*x */
1076 mpf_add_ui(y, y, 1UL); /* y = 1 + 0.5*x */
1077 DEBUGPRINT((stderr, "DEBUG: machine_epsilon(): x = %s\ty = %s\n",
1078 fp_to_string(x, 3), fp_to_string(y, 0)));
1079 while (mpf_cmp_ui(y, 1UL) != 0)
1080 {
1081 mpf_mul(x, half, x); /* x = 0.5*x */
1082 mpf_mul(y, half, x); /* y = 0.5*x */
1083 mpf_add_ui(y, y, 1UL); /* y = 1 + 0.5*x */
1084 DEBUGPRINT((stderr, "DEBUG: machine_epsilon(): x = %s\ty = %s\n",
1085 fp_to_string(x, 3), fp_to_string(y, 0)));
1086 }
1087
1088 mpf_clear(y);
1089 mpf_clear(half);
1090 mpf_set(*result, x);
1091 mpf_clear(x);
1092 #elif HAVE_LONG_DOUBLE
1093 x = 1.0L;
1094 y = 1.0L + x * 0.5L;
1095 store(&y);
1096 while (y != 1.0L)
1097 {
1098 store(&x);
1099 x *= 0.5L;
1100 y = 1.0L + x * 0.5L;
1101 store(&y);
1102 DEBUGPRINT((stderr, "DEBUG: machine_epsilon(): x = %.35Lg\ty = %.35Lg\n",
1103 (long double)x, (long double)y));
1104 }
1105 *result = x;
1106 #else
1107 x = 1.0;
1108 y = 1.0 + x * 0.5;
1109 store(&y);
1110 while (y != 1.0)
1111 {
1112 store(&x);
1113 x *= 0.5;
1114 y = 1.0 + x * 0.5;
1115 store(&y);
1116 DEBUGPRINT((stderr, "DEBUG: machine_epsilon(): x = %.35g\ty = %.35g\n",
1117 (double)x, (double)y));
1118 }
1119 *result = x;
1120 #endif
1121
1122 DEBUGPRINT((stderr, "DEBUG: Macheps = %s\n", fp_to_string(*result, 0)));
1123 }
1124
1125
1126 static void
maxrelerr(fp_t * result,fp_t x,fp_t y)1127 maxrelerr(fp_t * result, fp_t x, fp_t y)
1128 {
1129 /*******************************************************************
1130 Compute the maximum relative error of two values, x and y, and
1131 store it in *result.
1132 *******************************************************************/
1133
1134 fp_t t1;
1135 fp_t t2;
1136 fp_t t3;
1137 fp_t xabs;
1138 fp_t yabs;
1139
1140 fp_init(&t1);
1141 fp_init(&t2);
1142 fp_init(&t3);
1143 fp_init(&xabs);
1144 fp_init(&yabs);
1145
1146 DEBUGPRINT((stderr,"DEBUG 0.70: x = %s\n", fp_to_string(x,10)));
1147 DEBUGPRINT((stderr,"DEBUG 0.71: y = %s\n", fp_to_string(y,10)));
1148
1149 fp_abs(&xabs, x);
1150 fp_abs(&yabs, y);
1151
1152 DEBUGPRINT((stderr,"DEBUG 0.72: x = %s\n", fp_to_string(xabs,10)));
1153 DEBUGPRINT((stderr,"DEBUG 0.73: y = %s\n", fp_to_string(yabs,10)));
1154
1155 /* See the documentation of the -relerr option in ndiff.man for the
1156 explanation of this complex definition: */
1157
1158 if (fp_cmp(x, y) == 0)
1159 fp_assign(result, Zero); /* *result = 0 */
1160 else if ((fp_cmp(x, Zero) != 0) && (fp_cmp(y, Zero) != 0))
1161 {
1162 fp_sub(&t1, x, y); /* t1 = x - y */
1163 DEBUGPRINT((stderr,"DEBUG 1: t1 = %s = %s - %s\n", fp_to_string(t1,10), fp_to_string(x,10), fp_to_string(y,10)));
1164 fp_abs(&t1, t1); /* t1 = abs(x - y) */
1165 DEBUGPRINT((stderr,"DEBUG 2: t1 = %s\n", fp_to_string(t1,10)));
1166 fp_min(&t2, xabs, yabs); /* t2 = min(abs(x),abs(y)) */
1167 DEBUGPRINT((stderr,"DEBUG 3: t2 = %s = min(%s,%s)\n", fp_to_string(t1,10), fp_to_string(xabs,10), fp_to_string(yabs,10)));
1168 fp_div(result, t1, t2); /* *result = abs(x-y)/min(x,y) */
1169 DEBUGPRINT((stderr,"DEBUG 4: *result = %s\n", fp_to_string(*result,10)));
1170 }
1171 else if ((fp_cmp(x, Zero) == 0) && (fp_cmp(y, Zero) != 0))
1172 string_to_fp(result,"1"); /* *result = 1 */
1173 else if ((fp_cmp(y, Zero) == 0) && (fp_cmp(x, Zero) != 0))
1174 string_to_fp(result,"1"); /* *result = 1 */
1175 else
1176 fp_assign(result, Zero); /* *result = 0 */
1177
1178 fp_clear(t1);
1179 fp_clear(t2);
1180 fp_clear(t3);
1181 fp_clear(xabs);
1182 fp_clear(yabs);
1183 }
1184
1185
1186 static size_t
opt_abserr(const char * next_arg)1187 opt_abserr(const char *next_arg)
1188 {
1189 /*******************************************************************
1190 Process the -abserr option, and return the number of extra
1191 arguments consumed.
1192 *******************************************************************/
1193
1194 if (next_arg == CNULL)
1195 error("expected -abserr absolute-error", "", "");
1196 string_to_fp(&ABSERR, next_arg);
1197 return (1);
1198 }
1199
1200
1201 static size_t
opt_author(const char * next_arg)1202 opt_author(const char *next_arg)
1203 {
1204 /*******************************************************************
1205 Process the -author option, and return the number of extra
1206 arguments consumed.
1207 *******************************************************************/
1208
1209 static const char *lines[] = {
1210 "Author:",
1211 "\tNelson H. F. Beebe",
1212 "\tCenter for Scientific Computing",
1213 "\tUniversity of Utah",
1214 "\tDepartment of Mathematics, 322 INSCC",
1215 "\t155 S 1400 E RM 233",
1216 "\tSalt Lake City, UT 84112-0090",
1217 "\tUSA",
1218 "\tEmail: beebe@math.utah.edu, beebe@acm.org,",
1219 "\t beebe@computer.org, beebe@ieee.org (Internet)",
1220 "\tWWW URL: http://www.math.utah.edu/~beebe",
1221 "\tTelephone: +1 801 581 5254",
1222 "\tFAX: +1 801 585 1640, +1 801 581 4148",
1223 (const char *)NULL,
1224 };
1225
1226 print_lines(stderr, lines);
1227 Done = 1;
1228 return (0);
1229 }
1230
1231
1232 static size_t
opt_copyright(const char * next_arg)1233 opt_copyright(const char *next_arg)
1234 {
1235 /*******************************************************************
1236 Process the -copyright option, and return the number of extra
1237 arguments consumed.
1238 *******************************************************************/
1239
1240 static const char *lines[] = {
1241 "########################################################################",
1242 "########################################################################",
1243 "########################################################################",
1244 "### ###",
1245 "### ndiff: compare putatively similar files, ignoring small numeric ###",
1246 "### differences ###",
1247 "### ###",
1248 "### Copyright (C) 2000 Nelson H. F. Beebe ###",
1249 "### ###",
1250 "### This program is covered by the GNU General Public License (GPL), ###",
1251 "### version 2 or later, available as the file COPYING in the program ###",
1252 "### source distribution, and on the Internet at ###",
1253 "### ###",
1254 "### ftp://ftp.gnu.org/gnu/GPL ###",
1255 "### ###",
1256 "### http://www.gnu.org/copyleft/gpl.html ###",
1257 "### ###",
1258 "### This program is free software; you can redistribute it and/or ###",
1259 "### modify it under the terms of the GNU General Public License as ###",
1260 "### published by the Free Software Foundation; either version 2 of ###",
1261 "### the License, or (at your option) any later version. ###",
1262 "### ###",
1263 "### This program is distributed in the hope that it will be useful, ###",
1264 "### but WITHOUT ANY WARRANTY; without even the implied warranty of ###",
1265 "### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ###",
1266 "### GNU General Public License for more details. ###",
1267 "### ###",
1268 "### You should have received a copy of the GNU General Public ###",
1269 "### License along with this program; if not, write to the Free ###",
1270 "### Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, ###",
1271 "### MA 02111-1307 USA. ###",
1272 "########################################################################",
1273 "########################################################################",
1274 "########################################################################",
1275 (const char *)NULL
1276 };
1277
1278 print_lines(stderr, lines);
1279 Done = 1;
1280 return (0);
1281 }
1282
1283
1284 static size_t
opt_fields(const char * next_arg)1285 opt_fields(const char *next_arg)
1286 {
1287 /*******************************************************************
1288 Process the -fields option, and return the number of extra
1289 arguments consumed.
1290 *******************************************************************/
1291
1292 if (next_arg == CNULL)
1293 error("expected -fields fields", "", "");
1294 initialize_fields(next_arg);
1295 return (1);
1296 }
1297
1298
1299 static size_t
opt_help(const char * next_arg)1300 opt_help(const char *next_arg)
1301 {
1302 /*******************************************************************
1303 Process the -help option, and return the number of extra
1304 arguments consumed.
1305 *******************************************************************/
1306
1307 static const char *lines[] = {
1308 "Usage:",
1309 "\tndiff [ -? ] [ -abserr abserr ] [ -author ] [ -copyright ]",
1310 "\t [ -fields n1a-n1b,n2,n3a-n3b,... ] [ -help ]",
1311 "\t [ -logfile filename ] [ -minwidth nnn ] [ -outfile filename ]",
1312 "\t [ -precision number-of-bits ] [ -quick ] [ -quiet ]",
1313 "\t [ -relerr relerr ] [ -separators regexp ] [ -silent ]",
1314 "\t [ -version ] [ -www ]",
1315 "\t infile1 infile2",
1316 (const char *)NULL
1317 };
1318
1319 print_lines(stderr, lines);
1320 Done = 1;
1321 return (0);
1322 }
1323
1324
1325 static size_t
opt_logfile(const char * next_arg)1326 opt_logfile(const char *next_arg)
1327 {
1328 /*******************************************************************
1329 Process the -logfile option, and return the number of extra
1330 arguments consumed.
1331 *******************************************************************/
1332
1333 if (next_arg == CNULL)
1334 error("expected -logfile filename", "", "");
1335 if (freopen(next_arg, "w", stderr) == (FILE *) NULL)
1336 {
1337 cleanup();
1338 exit(EXIT_FAILURE); /* no stderr error left to report an error on */
1339 }
1340 return (1);
1341 }
1342
1343
1344 static size_t
opt_minwidth(const char * next_arg)1345 opt_minwidth(const char *next_arg)
1346 {
1347 /*******************************************************************
1348 Process the -minwidth option, and return the number of extra
1349 arguments consumed.
1350 *******************************************************************/
1351
1352 if (next_arg == CNULL)
1353 error("expected -minwidth width", "", "");
1354 MINWIDTH = (size_t) awk_string_to_unsigned_long(next_arg);
1355 return (1);
1356 }
1357
1358
1359 static size_t
opt_outfile(const char * next_arg)1360 opt_outfile(const char *next_arg)
1361 {
1362 /*******************************************************************
1363 Process the -outfile option, and return the number of extra
1364 arguments consumed.
1365 *******************************************************************/
1366
1367 if (next_arg == CNULL)
1368 error("expected -outfile filename", "", "");
1369 if (freopen(next_arg, "w", stdout) == (FILE *) NULL)
1370 error("could not reopen ", next_arg, " on stdout");
1371 return (1);
1372 }
1373
1374
1375 static size_t
opt_precision(const char * next_arg)1376 opt_precision(const char *next_arg)
1377 {
1378 /*******************************************************************
1379 Process the -precision option, and return the number of extra
1380 arguments consumed.
1381 *******************************************************************/
1382
1383 size_t n;
1384
1385 if (next_arg == CNULL)
1386 error("expected -precision number-of-bits", "", "");
1387 n = (size_t) awk_string_to_unsigned_long(next_arg);
1388 if (n > 0)
1389 {
1390 #if HAVE_GMP
1391 Precision = n;
1392 mpf_set_default_prec(Precision);
1393 machine_epsilon(&Macheps); /* must be updated when precision changes */
1394 #else
1395 warning
1396 ("no support for multiple-precision arithmetic has been compiled into this version",
1397 "", "");
1398 #endif
1399 }
1400 return (1);
1401 }
1402
1403
1404 static size_t
opt_quick(const char * next_arg)1405 opt_quick(const char *next_arg)
1406 {
1407 /*******************************************************************
1408 Process the -quick option, and return the number of extra
1409 arguments consumed.
1410 *******************************************************************/
1411
1412 return (0); /* TO DO: complete this complex function */
1413 }
1414
1415
1416 static size_t
opt_quiet(const char * next_arg)1417 opt_quiet(const char *next_arg)
1418 {
1419 /*******************************************************************
1420 Process the -quiet option, and return the number of extra
1421 arguments consumed.
1422 *******************************************************************/
1423
1424 QUIET = 1;
1425 return (0);
1426 }
1427
1428
1429 static size_t
opt_relerr(const char * next_arg)1430 opt_relerr(const char *next_arg)
1431 {
1432 /*******************************************************************
1433 Process the -relerr option, and return the number of extra
1434 arguments consumed.
1435 *******************************************************************/
1436 fp_t one;
1437
1438 if (next_arg == CNULL)
1439 error("expected -relerr relative-error", "", "");
1440 string_to_fp(&RELERR, next_arg);
1441
1442 fp_init(&one);
1443 string_to_fp(&one,"1");
1444
1445 if (fp_cmp(RELERR, Macheps) < 1)
1446 {
1447 char s[100];
1448 const char *t1;
1449 const char *t2;
1450
1451 t1 = fp_to_string(RELERR, 3);
1452 t2 = fp_to_string(Macheps, 6);
1453 (void)sprintf(s,"RELERR = %s is below machine epsilon %s", t1, t2);
1454 warning(s,"","");
1455 awk_free_string(t1);
1456 awk_free_string(t2);
1457 }
1458 else if (fp_cmp(RELERR, one) > 0)
1459 fp_mul(&RELERR, RELERR, Macheps);
1460
1461 fp_clear(one);
1462
1463 return (1);
1464 }
1465
1466
1467 static size_t
opt_separators(const char * next_arg)1468 opt_separators(const char *next_arg)
1469 {
1470 /*******************************************************************
1471 Process the -separators option, and return the number of extra
1472 arguments consumed.
1473 *******************************************************************/
1474
1475 if (next_arg == CNULL)
1476 error("expected -separators separators", "", "");
1477 FS = next_arg;
1478 return (1);
1479 }
1480
1481
1482 static size_t
opt_silent(const char * next_arg)1483 opt_silent(const char *next_arg)
1484 {
1485 /*******************************************************************
1486 Process the -silent option, and return the number of extra
1487 arguments consumed.
1488 *******************************************************************/
1489
1490 SILENT = 1;
1491 return (0);
1492 }
1493
1494
1495 static size_t
opt_version(const char * next_arg)1496 opt_version(const char *next_arg)
1497 {
1498 /*******************************************************************
1499 Process the -version option, and return the number of extra
1500 arguments consumed.
1501 *******************************************************************/
1502 const char *t;
1503
1504 machine_epsilon(&Macheps);
1505
1506 t = fp_to_string(Macheps, 3);
1507 (void)fprintf(stderr,
1508 "Version %s of [%s] [%s: epsilon = %s]\n",
1509 VERSION, DATE, ARITHMETIC, t);
1510 awk_free_string(t);
1511 Done = 1;
1512 return (0);
1513 }
1514
1515
1516 static size_t
opt_www(const char * next_arg)1517 opt_www(const char *next_arg)
1518 {
1519 /*******************************************************************
1520 Process the -www option, and return the number of extra arguments
1521 consumed.
1522 *******************************************************************/
1523
1524 static const char *lines[] = {
1525 "Master World-Wide Web source distribution:",
1526 "\tftp://ftp.math.utah.edu/pub/misc/ndiff-" VERSION ".tar.gz",
1527 "\thttp://www.math.utah.edu/pub/misc/ndiff-" VERSION ".tar.gz",
1528 (const char *)NULL,
1529 };
1530
1531 print_lines(stderr, lines);
1532 Done = 1;
1533 return (0);
1534 }
1535
1536
1537 static void
print_lines(FILE * outfile,const char * lines[])1538 print_lines(FILE * outfile, const char *lines[])
1539 {
1540 /*******************************************************************
1541 Print the array of lines[] on outfile, supplying a newline after
1542 each.
1543 *******************************************************************/
1544
1545 size_t k;
1546
1547 for (k = 0; lines[k] != CNULL; ++k)
1548 (void)fprintf(outfile, "%s\n", lines[k]);
1549 }
1550
1551
1552 static void
report_difference(const char * f1line,const char * f2line,size_t nfield)1553 report_difference(const char *f1line, const char *f2line, size_t nfield)
1554 {
1555 /*******************************************************************
1556 Print a diff-style difference of two lines, but also show in the
1557 separator line the field number at which they differ, and the global
1558 absolute and relative errors, if they are nonzero.
1559 *******************************************************************/
1560
1561 if (SILENT == 0)
1562 {
1563 (void)printf("%luc%lu\n", (unsigned long)NRLINE, (unsigned long)NRLINE);
1564 (void)printf("< %s\n", f1line);
1565 if ((fp_cmp(This_Abserr, Zero) != 0)
1566 || (fp_cmp(This_Relerr, Zero) != 0))
1567 {
1568 const char *abserr;
1569 const char *relerr;
1570 fp_t emult;
1571 fp_t tenk;
1572
1573 fp_init(&emult);
1574 fp_init(&tenk);
1575
1576 fp_assign(&emult, This_Relerr);
1577 fp_div(&emult, emult, Macheps); /* emult = This_Relerr / Macheps */
1578 fp_assign(&tenk, Zero);
1579 string_to_fp(&tenk, "10000");
1580 abserr = fp_to_string(This_Abserr, 3);
1581 relerr = fp_to_string(This_Relerr, 3);
1582
1583 if (fp_cmp(ABSERR, Zero) == 0) /* suppress absolute error report */
1584 {
1585 if (fp_cmp(emult, tenk) >= 0)
1586 (void)printf("--- field %d\trelative error %s\n",
1587 (int)nfield, relerr);
1588 else
1589 (void)printf("--- field %d\trelative error %s [%d*(machine epsilon)]\n",
1590 (int)nfield, relerr, (int)(fp_to_double(emult) + 0.5));
1591 }
1592 else if (fp_cmp(RELERR, Zero) == 0) /* suppress relative error report */
1593 (void)printf("--- field %d\tabsolute error %s\n", (int)nfield, abserr);
1594 else /* report both absolute and relative errors */
1595 {
1596 if (fp_cmp(emult, tenk) >= 0)
1597 (void)printf("--- field %d\tabsolute error %s\trelative error %s\n",
1598 (int)nfield, abserr, relerr);
1599 else
1600 (void)printf("--- field %d\tabsolute error %s\trelative error %s [%d*(machine epsilon)]\n",
1601 (int)nfield, abserr, relerr, (int)(fp_to_double(emult) + 0.5));
1602 }
1603
1604 awk_free_string(abserr);
1605 awk_free_string(relerr);
1606
1607 fp_clear(emult);
1608 fp_clear(tenk);
1609 }
1610 else
1611 (void)printf("--- field %d\n", (int)nfield);
1612 (void)printf("> %s\n", f2line);
1613 }
1614 Error_Count++;
1615 }
1616
1617
1618 static void
string_to_fp(fp_t * result,const char * s)1619 string_to_fp(fp_t * result, const char *s)
1620 {
1621 /*******************************************************************
1622 Convert a numeric string to an fp_t floating-point number, and
1623 return the result as a pointer to a floating-point number in an
1624 internal buffer which will be overwritten on the next call.
1625
1626 Fortran use has any of E, e, D, d, Q, or q, or even nothing at
1627 all, for the exponent letter, but awk and C only allow E and e.
1628
1629 Ada usefully permits nonsignificant underscores for
1630 readability: 3.14159265358979323846 and
1631 3.14159_26535_89793_23846 are equivalent.
1632
1633 We can safely assume that there are no leading or trailing
1634 whitespace characters, because all strings passed to this
1635 function are the result of splitting lines into
1636 whitespace-delimited fields.
1637 *******************************************************************/
1638
1639 char *t;
1640 const char *t1;
1641 const char *t2;
1642
1643 t = (char*)awk_dupstr(s);
1644 (void)awk_gsub("_", "", &t); /* remove Ada-style separators */
1645 (void)awk_gsub("[DdQq]", "e", &t); /* convert Fortran exponent letters to awk-style */
1646 if (awk_match(t, "[0-9.][-+][0-9]+$")) /* then letter-less exponent */
1647 {
1648 t1 = awk_substr(t, 1, RSTART);
1649 t2 = awk_substr(t, RSTART + 1, RLENGTH - 1);
1650 awk_free_string(t);
1651 t = awk_padstr("",(size_t)(RSTART + RLENGTH));
1652 (void)sprintf(t, "%se%s", t1, t2); /* insert exponent letter e */
1653 awk_free_string(t1);
1654 awk_free_string(t2);
1655 }
1656 if (awk_match(t, "^[ \t]*[-+]?[.]")) /* missing digit before decimal point? */
1657 { /* gmp does not recognize such numbers */
1658 char *t3;
1659
1660 t3 = awk_padstr(t,1);
1661 t3[RSTART + RLENGTH - 2] = '0'; /* supply leading zero digit */
1662 (void)strcpy(&t3[RSTART+RLENGTH-1],&t[RSTART+RLENGTH-2]);
1663 awk_free_string(t);
1664 t = t3;
1665 }
1666
1667 #if HAVE_GMP
1668 mpf_set_str(*result, t, 10);
1669 #elif HAVE_LONG_DOUBLE
1670 (void)sscanf(t, "%Lg", result);
1671 #else
1672 *result = (fp_t) awk_string_to_double(t);
1673 #endif
1674
1675 awk_free_string(t);
1676 }
1677
1678
1679 static void
usage(void)1680 usage(void)
1681 {
1682 /*******************************************************************
1683 Print a usage message on stderr, and exit with a failure code.
1684 *******************************************************************/
1685
1686 static const char *lines[] = {
1687 "\t[ -? ] [ -abserr abserr ] [ -author ] [ -copyright ]",
1688 "\t\t[ -fields n1a-n1b,n2,n3a-n3b,... ] [ -help ]",
1689 "\t\t[ -logfile filename ] [ -minwidth nnn ]",
1690 "\t\t[ -outfile filename ] [ -precision number-of-bits ]",
1691 "\t\t[ -quick ] [ -quiet ] [ -relerr relerr ]",
1692 "\t\t[ -separators regexp ] [ -silent ] [ -version ]",
1693 "\t\t[ -www ] infile1 infile2",
1694 (const char *)NULL,
1695 };
1696
1697 (void)fprintf(stderr, "Usage:\n\t%s", Program_Name);
1698 print_lines(stderr, lines);
1699 cleanup();
1700 exit(EXIT_FAILURE);
1701 }
1702
1703
1704 static void
warning(const char * s1,const char * s2,const char * s3)1705 warning(const char *s1, const char *s2, const char *s3)
1706 {
1707 /*******************************************************************
1708 Print a warning message on stderr, using emacs compile-command-style
1709 message format.
1710 *******************************************************************/
1711 if (FNR > 0)
1712 (void)fprintf(stderr, "%s:%ld:%%%%%s%s%s\n", FILENAME, (long)FNR, s1, s2, s3);
1713 else /* special case for diagnostics during initialization */
1714 (void)fprintf(stderr, "%s%s%s\n", s1, s2, s3);
1715 }
1716