/* Portable Windowed Applications: Pollsters, Onlookers, and Others Example Applications Copyright 1998 by OG Software, Inc., Ann Arbor, MI 48104 This C/C++ source file contains the example applications. They are portable, of course, and must be compiled in conjunction with the appropriate version of the POO software. The application generated depends on the defined constants below. Defined Constant Application ---------------- ----------- Pollster Pollster Onlooker Onlooker Random Pseudorandom Number Generators fpMap Floating-Point Mapping PNCalc Polish Notation Calculator NICMap Convergence Map for Newton's Iteration NEigen Nonsymmetric EigenProblem */ #define Pollster /* The following is used only in Random. It is defined here because Random should be run with various values of this limit. With a limit of 12000, Random is responsive; for 12000000, it can be very unresponsive when single threaded. */ #define CHK_LIMIT 12000L /* Building an Example Application in Windows or OS/2 To build an example application, you will need this file and the appropriate POO source file. In some Windows environments, a .DEF file NAME Examples DESCRIPTION '...' EXETYPE WINDOWS CODE PRELOAD MOVEABLE DISCARDABLE DATA PRELOAD MOVEABLE MULTIPLE HEAPSIZE 1024 STACKSIZE 8192 is required, but a .RC file is not needed. For some C compilers, you must specify the 'windowed' and the 'multithreaded' options, which vary from obvious to obscure depending on your C compiler. The Windows and OS/2 versions of the software run multithreaded unless you change the source code; the multithreading issue is explored in the third example. It is not important, at least for now, that you understand what role these options play. Building an Example Application in Unix/X To build an example application, you will need this file and the POO source file for the X Window System. You will need an ANSI C compiler; a Kernighan and Ritchie C compiler will NOT compile the POO software. The POSIX I/O library is used to interface to the file system, but I am not aware of a Unix/C compiler that does not support the POSIX I/O library. When you compile, you must specify the "-lX11" option because the X library is required; this option is the equivalent of the "windowed" option. The X version of the POO software runs single threaded because the variable POO_MULTITHREADED is undef'd but includes multithreading support based on the POSIX standard, Pthreads. Since the first two examples are Pollster and Onlooker, which are interactive and run single threaded in any case, you can defer consideration of multi- threading until Random. Random, the third example, provides a safe setting for moving to the multithreaded version of the software. The "multithreading" option is likely "-lpthreads" if your system supports the POSIX threads standard. Running an Example Application The first time you run one of these examples, it will not find its memory file so that it will start with a blank slate, e.g., the windows will be created with the default colors. Both Random and fpMap run this way because they do not provide an application name to the POO software just to illustrate the consequences. The other examples are more polite so that the memory file is automatically updated when they are terminated. Window attributes and function key string assignments are saved in the memory file. You should try various colors, fonts, and geometries and some function keys and re-run the examples to verify that they are retained. Onlooker is unique in that it saves its state in the memory file. You, the user, can terminate a POO application at any time with a Ctrl+D. You can also use the EXIT or STOP commands. Both Ctrl+D and EXIT terminate the application immediately but update the memory file in the process; STOP terminates but skips the update so that the memory file remains unchanged. The examples that read input supply default values if you respond by pressing Enter (Return) whenever you are solicited for input. In X, the caret may not be drawn unless you change the 'focus' to the window. If a POO window has the focus, you can simply start typing because the POO software always accepts unsolicited input. It's usually treated as typeahead, but pollsters and onlookers interpret such input as commands. Introduction A few excerpts to illustrate the use of the API functions might have been sufficient, but complete applications are considerably more effective and interesting. The code excerpts often found in the sterile and stylized descriptions of functions and their arguments common in computing literature are rarely sufficient regardless of their sparkling clarity. They are helpful in some circumstances but belie the context essential to understanding. From my perspective, the issue was to create example applications that would illustrate the use of the API functions, test their functionality, and have interesting content. The first of these criteria, of course, was the principal one. The second has value for both of us: for me because it saved the effort of writing additional test programs, and for you because it provides simple tests of the software in the context of your system. The third criterion has value for both of us if you share my interest in, perhaps passion for, numerical computation. The examples provide numerous variations on how to use the API functions; and, if you follow the suggestions offered, you should gain familiarity with the user environment created by the software. The Pollster and Onlooker applications are the first two examples. Each consists of a single line of code because both "Pollster" and "Onlooker" are reserved as application names, as window names also. Pollsters and onlookers are built into the software: the user of any POO application can open either at any time. Although these examples provide at best lip service to the first criterion, they are quite good with respect to the other two. In particular, you should use them to gain familiarity with the user environment that is created by the POO software. Pollsters are very input oriented, while onlookers are more output oriented. The five remaining examples are numerical and good relative to the first criteria, moderate to the second, and good to the third. The first two are output only, the third is interactive, and the last two solicit input before producing their output. In the last two, if you press Enter whenever solicited for input, the example will provide default values so that you don't really have to know what to enter to see what will happen. Running them with the defaults the first time will give you a feel for what to expect if you want to explore the numerical problems involved. Indeed, this tactic is convenient in most applications. Further, after a year or so, most authors forget how to run their applications, and this tactic provides an embedded tutorial of sorts. The first numerical example, Random, is designed to illustrate why windowed applications should be run multithreaded and contains a short discussion of threads. If used as is, the POO software runs multithreaded in Windows and OS/2, but single threaded in X. If you are running in X, you should use this example to get the POO software running multithreaded. You will not be satisfied running your processor intensive applications single threaded so that you should tackle the multithreading issue with Random. The next example, fpMap, tests the floating-point output and input mappings provided by the C library; it's output only. This issue was settled years ago but seems to have revived; your C library may fail the test. The third numerical example, PNCalc, is an interactive, Polish notation calculator and is moderately useful if you're familiar with this notation, e.g., have experience using Hewlett-Packard calculators. It supports only rudimentary calculations but is quite extensible. By default, it runs with a single window, but you can request windows to see the stack and/or memory, both of which are updated automatically but can be closed and will stay closed until you again ask to see them. This example is handy for experimenting with the function keys and I/O direction. In any case, you'll clearly have to enter something if you want to have anything happen, say "2 2 +". The last two examples, NICMap and NEigen, are more in the realm of interesting numerical applications. In numerical computation, problems range from well to poorly conditioned, and methods from stable to unstable. Well-conditioned problems and stable methods are nice, particularly when embedded in large problems, but not terribly interesting. Aberrations arising from either the problem or method or both can be fascinating. NICMap employs a numerical method of dubious virtue that's used regularly in various guises for solving many different types of problems. NEigen employs an impeccable method applied to nonsymmetric eigenproblems. The purpose of NICMap is to color a rectangle in the complex plane based on whether Newton's iteration converges to a zero when it's started at each grid point in the rectangle. The convergence maps are produced in a window as a printer plot, which is a now ancient but honorable output mechanism, and their quality is poor, albeit appropriate in a monograph on text windows. This example reads the polynomial, some parameters, the rectangle, and the grid size, all of which you are expected to provide. You'll want to use a fixed size font in the map window. NEigen uses an impeccable method to compute the eigenvalues of a nonsymmetric matrix and allows you to choose the eigenvalues and their condition numbers with the caveat that it picks the last condition number. Without this caveat, either the mathematics of the eigenproblem would be wrong or the example would fail; the choice seemed clear. Given the eigenvalues and condition numbers, NEigen generates an appropriate nonsymmetric eigenproblem and then solves it and prints the computed eigenvalues. If you enter nice condition numbers, the computed eigenvalues will be close to the one's you prescribed; otherwise.... Header File A header file for the POO software is not provided because only standard data types are required and there are only five function prototypes, namely, */ extern int cappl (char *appname); extern void *copen (char *name, int size, int xsize); extern int cclose (void *tw); extern int cputs (void *tw, char *text, int len); extern int cgets (void *tw, char *text, int len); /* Placing the preceding function prototypes in a file, e.g., poo.h, should suffice if you want to use a header. */ #define BEGIN { #define END } #ifdef Pollster /*================================================================ Pollster Pollster was designed to keep a reliable log consisting of entries provided by the user. Its progenitor automated an accurate log of my computer activity. Each entry starts with a header containing the date and time and a unique integer. The header is followed by the user's input, indented two spaces for visual separation. The unique integer in the header is the displacement of the first character of the entry within the file and can be used to position the file in the window or to make internal references. Further, they provide a basis for assessing the integrity of the file since manual editing of the file invalidates them unless it is done very carefully. By default, a pollster displays the current entry at the top of the window and solicits input at the bottom. You are expected but not required to enter something in response to the unasked, and thus unbiased, question. Input can be pasted (F5) so that your entry can contain excerpts from other windows associated with the application. A Ctrl+Z (EOF) terminates the entry and closes the pollster. The output of a pollster is directed to the poll file and cannot be changed. Because no mechanism is provided for preventing the update of the poll file that occurs when the pollster is closed, the format of a poll file is examined so that you cannot open a pollster for just any file. Most files won't pass the test, but the test is not a significant impediment to a mischievous user. While running this example, you should try the various facilities for editing input. You should try the Insert, Backspace, Delete, and arrow keys. Ctrl+U deletes your input, while Ctrl+A and Ctrl+B delete the input after and before the caret, respectively. Ctrl+X, Ctrl+C, and Ctrl+V or Shift+Delete, Ctrl+Insert, and Shift+Insert, respectively, implement cut/copy/paste; both the Windows-isms and OS/2-isms are supported for user portability. You should try line re-entry by placing the pointer in a line in the window and then pressing the right mouse button; if you then move the pointer to the input rectangle and press it again, it has the same effect as pressing Enter. You might also want to assign strings to a few functions keys. You can try the scrolling keys, but you will tire of typing before you have a large enough file to make it worth the effort; wait for Onlooker. Normally, pollsters are soliciting input so that commands must be terminated with Ctrl+Enter or Shift+Ctrl+Enter. The POO software distinguishes normal input and commands by how you terminate your input rather than the more common command character. To suspend a pollster's solicitation, you can press Esc, but be careful since Esc is typamatic and closes a window if input isn't being entered. Subsequently, you can resume an entry with an NB command; until you do so, all input is interpreted as a command. You should configure the pollster just to see how the attribute commands work. Window attributes are automatically saved in the memory file so that your configuration will be used the next time you invoke Pollster. For historical reasons, I have it configured to use distinctive colors and font and a geometry of 65x10c, i.e., a small window centered in the screen. This rendition of the Pollster application uses the default file, Poll. The name of the poll file can follow "Pollster" in the string passed to cappl. */ int appl (void) { return cappl("Pollster"); } #endif #ifdef Onlooker /*================================================================ Onlooker Onlooker was designed to provide a safe and convenient mechanism for rummaging or browsing through a file system which is the first thing I want to do whenever I fire up a new system. It's aptly named because it's not allowed to change files unless you give the secret password even though onlookers are not restricted in other POO applications. An onlooker can display directories and files. For a directory, the window shows the list of file names in the directory: the "." entry is replaced with the complete path name and appears first, and the ".." entry appears as itself. By default, binary files are shown in hex dump format, but you can toggle (Ctrl+O) between text and hex dump format. The first time you run Onlooker, an onlooker is opened for the current directory. But Onlooker saves both file reference and position information in the memory file and uses it to restore its state whenever it's invoked. If position information is available for a file, an onlooker opened for the file is positioned to the saved position. The current directory of an onlooker is always set to the path of whatever file is associated with it. Usually, input in an onlooker is interpreted as a command. But, if you enter an NB, the onlooker solicits input and appends it to the VM file. It stops soliciting input after an Esc or Ctrl+Z. Unlike a pollster, input is appended unchanged, and only the VM file is changed. The corresponding file in the file system, if any, isn't affected unless you explicitly store it with a > command, which will require the secret password because you are writing into an existing file. You can paste input into an onlooker at any time, a prior NB is not required. If you direct input with a < command, the file associated with the onlooker is changed. The easiest method for opening additional onlookers is to use the right mouse button and a function key assigned the string "on ^M". If an onlooker is showing a directory, line re-entry (right mouse button) inserts the line in the directory as input and leaves the caret at the beginning of the input. If you press the function key, "on " is inserted at the beginning of the line, and the "^M" is interpreted as a Ctrl+Enter so that the software interprets the input as a command. But "on" is an abbreviation of onlooker, the next blank delimited string is a file name, and command processing ignores extraneous operands, i.e., the junk following the name. In summary, a new onlooker is opened for the file. If you prefer the manual approach, you can use either "\" or "/" as the path name separation character and environmental variables in either Unix ${...} or Windows %...% format. In this context, paths are evaluated relative to the onlooker's current directory, A CD with no operand displays the current directory, while a CD with a valid operand changes it. The file displayed in an onlooker can be changed using input direction, e.g., "<." changes it to the current directory. You'll want to try the scrolling keys and the scroll bar. The bar is at the right edge where it belongs but may or may not be drawn (see the MODE command). Whether visible or not, pressing the left button when the pointer is within the scroll bar moves the window. Shift+the left button has the same effect but doesn't require that the pointer be within the scroll bar. You'll also want to drag the text around the window with the left button. In case you get too far left or right, Ctrl+Left or Ctrl+Right reset the horizontal shift; it is also reset by Ctrl+Home and Ctrl+End but not by Home and End. Note that Left, Right, Home, and End are applied to the window or the input caret, a context dependency. You can use Ctrl+A and Ctrl+B to restrict the portion of the file associated with the onlooker, and Ctrl+U to restore the complete file. These key combinations are applied to the file or the input being entered through the window, another context dependency. But these keys do not delete data from the file, they only hide the data after the last or before the first line in the window. These file restrictions are used when cutting and pasting (F7 and F5) and when storing with > or >> commands. For example, it is quite straightforward to break a large file into pieces that will fit on a diskette and then put them back together on another system. You may want to try a FIND command or two, which treats the string operand in a somewhat unusual fashion: initial and terminal ?'s alter the search criteria and are not part of the search string. If you have two files that don't differ radically, you may want to try a CF command to compare them; the other file is specified by the number of the onlooker showing it. */ int appl (void) { return cappl("Onlooker"); } #endif #ifdef Random /*================================================================ Pseudorandom Number Generators At the beginning of the POO source code, there is either a #define or #undef for POO_MULTITHREADED, which determines whether the POO software is compiled to run single or multithreaded. You want to use the multithreaded version, and the purpose of Random is to illustrate why in a safe setting. Simply generating random integers and printing the frequencies is not much of a task for a self-respecting computer, particularly when run using a full-cycle pseudorandom number generator, but it is precisely what this example does. Random solicits no input and first prints after 1024 numbers have been generated; thereafter, the time between each output line doubles. Each line starts with the number of samples and the seconds per number followed by the frequency table. The size of the frequency table is determined by the defined constant N(4). Since Random runs for the full cycle and time is doubling, Random starts quickly but becomes processor intensive and, depending on your processor, may take up to an hour. You can terminate (Ctrl+D or Esc) or suspend (Ctrl+S or Pause) if you can get its attention, which is what makes threads pertinent. Windowing Systems and Threads The X Window System, Windows, and OS/2 Presentation Manager are conceptually the same in design and do not differ significantly when viewed in terms of their programming interface, particularly Windows and OS/2 which had more or less the same origin. These windowing systems are event driven. Since I started with X, the X terminology is more natural for me; in Windows and OS/2, events are called messages. Events may be generated by the windowing system, by the user, or by the application, e.g., maximizing a window may generate a hundred or more events for numerous windows, while pressing a key generates two or more for the focus window. In any case, these windowing systems start queuing events as soon as the first window is created and expect that the application that created the window will request and then process the events in a timely fashion. When running windowed applications, you may have noticed that the pointer often changes from the familiar arrow to an hourglass or clock and that, when this happens, some or all of the windows on your display may be unresponsive until the arrow is restored. The lack of response stems from the fact that one or more applications is ignoring its event queue, which may block even those that are looking for events. Threads is an answer, i.e., an application is busy doing something and did not start another thread to request and to process the events that the windowing system continues to queue. In short, nobody is home. Windowing systems are older than threads, and there are other ways of handling these situations, but multiple threads is by far the easiest. The POO software can be compiled to run either multithreaded or single threaded but should be run multithreaded unless your system does not support threads. When multithreaded, your application runs as one thread, and the POO software creates another thread to service the event queue so that the application remains responsive even though it may be processor intensive. Further, if you are so inclined, your application can be multithreaded if the POO software is running multithreaded, i.e., you can create additional threads. When single threaded, you are responsible for ensuring that your application calls a POO API function frequently enough to remain responsive. When single threaded, each POO API function processes all queued events, then performs its function, and then checks for events again before returning. If you run single threaded and fail to call an API function frequently enough, the underlying system may complain, perhaps bitterly, or even freeze to death. The introduction to "POO MAIN PROGRAM, API FUNCTIONS, AND THREADS" in the POO source code contains additional commentary on threads if (or when) you are interested. The Windows and OS/2 versions of the POO software use the thread support provided by these systems. Technically, Windows 3.1 does not, but C/C++ compilers for Windows 3.1 often support the Win32s interface, which includes threads but limits you, in theory, to running at most one such application at a time. These versions run multithreaded unless you change them, and the only thing you need to know about threads is how to tell your compiler that your application is multithreaded. The Windows and OS/2 versions can be run single threaded by changing the #define POO_MULTITHREADED to a #undef. Random provides a safe setting in which you can gain an appreciation for why you don't want to run single threaded. The X version of the POO software is implemented in terms of the POSIX threads standard, so-called Pthreads. The multithreading code uses only components of the POSIX threads standard. The code could be implemented using POSIX semaphores quite naturally but is instead implemented using condition variables and mutexes because they are part of the POSIX threads standard whereas semaphores are not. This is unfortunate because the semaphore implementation is much cleaner and easier to follow. If your system supports POSIX threads, you will need to specify an appropriate library option, typically "-lpthreads". The POO_..._SYSTEM variable was intended to deal with differences arising from the multithreading code so that it could be based on POSIX threads or vendor-specific threads. Eventually, this work was abandoned because I could get no access to relevant systems. If your system does not support POSIX threads, you will have to either run single threaded, which is not really feasible for any application that is processor intensive, or write appropriate code based on the available thread support. This sounds much worse than it is. Threads only account for a couple dozen lines of code and the Windows and OS/2 versions can be used as guides to what must be changed because the thread support in both Windows and OS/2 is vendor specific. Although Random becomes processor intensive, it's responsive even when run single threaded because the 'chk' nonsense ensures that an API function is called regularly, namely cputs(NULL,NULL,0), so that queued events are processed in a timely fashion. The CHK_LIMIT variable defined at the beginning of this file controls the response level when running single threaded. Although you may feel that response time is not critical, the windowing system may become cranky or worse if your application does not process the events that are queued for it, e.g., your key strokes. */ #include #include int itw (int *n, long *seed); int iad (int *n, long *seed); double dtw (long *seed); double dad (long *seed); #define N 4 /* All output is done using macros, namely BSO...ESO and BFO...EFO. Because the output lines are produced in pieces, a \v is inserted so that the window updates are skipped until the line is complete. You will probably not be able to tell the difference, but it illustrates a use of \v. */ #define BSO cputs(w, #define ESO ,0); #define BFO cputs(w,buf,sprintf(buf, #define EFO )); int appl (void) BEGIN void *w; int n, i; char buf[120+10*N]; long seed, chk; double bucket[N], cnt, t, pcnt; seed = 1; for (n = 0; n < N; ++n) bucket[n] = 0.0; cnt = 0.0; w = copen("io:Random Integers;",1,0); pcnt = 1024.0; chk = 0L; BSO "The time between output lines is doubling and will\n" ESO BSO "get lengthy. Hope you're patient; 2B to go!\n\n" ESO while (1) BEGIN ++bucket[itw(&n,&seed)]; ++cnt; if (cnt == pcnt || seed == 1) { t = (1000000.0/(double)CLOCKS_PER_SEC)/cnt; BFO "\v%11.0lf %8.3lf",cnt,t*clock() EFO for (t = 0.0, i = 0; i < N; ++i) { t += bucket[i]; BFO "\v%10.6lf",bucket[i]/cnt EFO } BSO "\n" ESO if (t != cnt) { BFO "Oops, buckets only have %.0lf != %.0lf!\n",t,cnt EFO break; } if (seed == 1) break; pcnt += pcnt; } else if (++chk == CHK_LIMIT) { chk = 0L; cputs(NULL,NULL,0); } END BSO "I quit!\n" ESO return 0; END /*---------------------------------------------------------------- Full-Cycle Pseudorandom Number Generators The pseudorandom number generators below are full-cycle generators in that the seed will assume all positive 32-bit integer values precisely once before the sequence starts again. Two complete sets are included: one for positive integers between 0 and N-1, where N is a generator argument, and one for doubles in (0,1). Numerous generators can be found in the literature, and there is a plethora of articles on their efficacy. The ones here are simply useful examples and do not include code to guard against misuse, which is not practical in functions of this size. The seed argument should be a positive integer, and nobody should tamper with it; set it and forget it. When generating random integers, the range should be a positive integer less than 32768. Failure to follow these restrictions will result in garbage. Both itw and iad return a uniformly distributed int in the range from 0 to (n-1), where n (the first argument) must be a positive integer less than 32768, and update the seed. For consistency and Fortran use, n must be passed by name but is not changed. The computation in the return is fixed-point multiplication of the integer n and the seed viewed as a 31 bit fraction. After forming the high order 31 bits of the full product, the final shift takes into account the binary point. There was a time when fixed-point computations were common. Both dtw and dad return a uniformly distributed double in the open interval (0,1) and update the seed; the actual range is 2**-32 to 1 - 2**-32, inclusive. The results produced by itw and iad are n times the results produced by the floating-point generators, but timing tests give a slight edge to the integer generators, about 25%, for Intel architectures. The *tw functions are adapted from Tausworthe and Whittlesey, CACM 11, 9, Sep 1968. Alternatively, dtw and itw can be coded using any one of the following s = s ^ (s >> 6); *seed = s = 2147483647 & (s ^ (s << (31-6))); s = s ^ (s >> 7); *seed = s = 2147483647 & (s ^ (s << (31-7))); s = s ^ (s >> 13); *seed = s = 2147483647 & (s ^ (s << (31-13))); All four variations yield full-cycle generators. The *ad functions are adapted from Ahrens and Dieter, CACM 15, 1972, and employ a multiplicative congruential method. */ int itw (int *n, long *seed) { long s = *seed; s = s ^ (s >> 3); *seed = s = 2147483647L & (s ^ (s << (31-3))); return (int)(((*n)*(s >> 15) + (((*n)*(s & 32767)) >> 15)) >> 16);} double dtw (long *seed) { long s = *seed; s = s ^ (s >> 3); *seed = s = 2147483647L & (s ^ (s << (31-3))); return ((double)s)/2147483648.0; } int iad (int *n, long *seed) { long s = *seed; s = *seed = 2147483647L & (s * 663608941L); return (int)(((*n)*(s >> 15) + (((*n)*(s & 32767)) >> 15)) >> 16);} double dad (long *seed) { long s = *seed; s = *seed = 2147483647L & (s * 663608941L); return ((double)s)/2147483648.0; } #endif #ifdef fpMap /*================================================================ Floating-Point Mapping This example, fpMap, tests the floating-point conversion mappings provided by the C library, namely the output mapping contained in sprintf and the input mapping in strtod; with minor changes, you can substitute sscanf for strtod. This is an old issue that was settled theoretically in the late sixties, and implementations of the correct mappings were in general use in the early seventies. The theory is valid for conversions between any two floating-point representations; but, in computing, one is decimal (our world) and the other radix is usually 2 (its world). If you're interested in the theory, see the following articles by D. W. Matula Base conversion mappings. Proc. AFIPS 1967 Spring Joint Comput. Conf., Vol. 30, pp. 311-318. In-and-Out Conversions. Comm. ACM 11, 2 (Jan. 1968) 47-50. A Formalization of Floating-Point Numeric Base Conversion. IEEE Trans. Computers, C-19, 8 (Aug. 1970) 681-692. Matula [1968] actually discusses the converse of the problem that is considered here. Unfortunately, theory is just knowledge, not a standard, and your C library may or may not pass this test. If the conversions are correctly rounded with adequate precision, the output mapping is 1-1 and the input mapping is onto and, as a consequence, the out-in mapping is the identity. This example tests the float and double out-in mappings to see if they are the identity, i.e., the original floating-point number is recovered. For the 32-bit and 64-bit IEEE representations, the magic numbers are 9 and 17, respectively. The single-precision tests use 7 to 10 decimal digits; and the double, 14 to 19. The tests with too few digits should fail, while the remainder should succeed (if the C library performs the conversions properly). The tests use 10000 random numbers in (0,1) so that failure is definitive but success merely suggestive unless you increase the number of trials very significantly. By increasing TRIALS by one or more orders of magnitude, you can use this example to verify that the POO is running multithreaded. The pseudorandom number generator, the first three lines of code in the inner loop, is adapted from Tausworthe and Whittlesey, CACM 11, 9, Sep 1968. Each test uses the same numbers because the seed for the generator is reset for each test. */ #include #include #define TRIALS 10000L int appl (void) BEGIN void *w; int n, k; char fmt[12], buf[80]; long seed, trial; float y; double x; w = copen("Report:Floating-Point Mapping Test;",1,1); cputs(w,"Single-precision floating-point tests.\n\n",0); for (n = 6; n < 10; ++n) BEGIN *fmt = '%'; sprintf(fmt+1,"%i.%ile",n+8,n); cputs(w,buf,sprintf(buf,"For format %s, ",fmt)); for (seed = trial = 1; trial < TRIALS; ++trial) { seed = seed ^ (seed >> 13); seed = 2147483647L & (seed ^ (seed << (31-13))); y = (float)(((double)seed)/2147483648.0); sprintf(buf,fmt,y); if (y != (float)strtod(buf,NULL)) break; } if (trial < TRIALS) k = sprintf(buf,"trial %li failed for ",trial); else k = sprintf(buf,"%li trials ending with ",trial); k += sprintf(buf+k,fmt,y); *(buf+k) = '\n'; cputs(w,buf,k+1); END cputs(w,"\nDouble-precision floating-point tests.\n\n",0); for (n = 13; n < 19; ++n) BEGIN *fmt = '%'; sprintf(fmt+1,"%i.%ile",n+8,n); cputs(w,buf,sprintf(buf,"For format %s, ",fmt)); for (seed = trial = 1; trial < TRIALS; ++trial) { seed = seed ^ (seed >> 13); seed = 2147483647L & (seed ^ (seed << (31-13))); x = ((double)seed)/2147483648.0; sprintf(buf,fmt,x); if (x != strtod(buf,NULL)) break; } if (trial < TRIALS) k = sprintf(buf,"trial %li failed for ",trial); else k = sprintf(buf,"%li trials ending with ",trial); k += sprintf(buf+k,fmt,x); *(buf+k) = '\n'; cputs(w,buf,k+1); END cputs(w,"\nMapping tests completed.\n",0); return 0; END #endif #ifdef PNCalc /*================================================================ Polish Notation Calculator This example application is a Polish notation calculator and is moderately useful if you're familiar with this notation, e.g., have experience using Hewlett-Packard calculators. It supports only rudimentary calculations but can be easily enhanced. Some of the features supported by the POO software, e.g., input direction, copy/paste, and function key strings, make it straightforward to implement extended functions and recursive computations. The Polish notation calculator is implemented with a PNC struct and a function, pnc, that processes a string in the context of a calculator defined by a PNC struct. Thus the application deals only with the I/O and runs the calculator in a single window. The user can open and close additional windows to view the stack and memory; when displayed, these windows are updated automatically. */ #include #include #include #include #include #include #include /*---------------------------------------------------------------- The PNC struct completely defines the state of a Polish notation calculator, and the function pnc processes a string relative to a PNC struct and updates it appropriately. Consequently, one could implement multiple, independent calculators. All errors are indicated by setting appropriate flags, the PNC_... constants defined after the struct. Changing the stack size has no code ramifications but impacts the user. Changing the memory size of the calculator requires code changes, i.e., names for the memory locations. */ #define PNC_MEMSIZE 8 #define PNC_STKSIZE 16 typedef struct { int flags; /* Return & flag, PNC_... */ int memsize; /* Memory size, PNC_MEMSIZE */ int stacksize; /* Stack size, PNC_STKSIZE */ int stackptr; /* Stack pointer re stack component */ char *msg; /* Cryptic error comment or help */ double stack[PNC_STKSIZE]; /* Floating-point stack */ double memory[PNC_MEMSIZE]; /* Calculator memory */ } PNC; #define PNC_NOARGS 1 #define PNC_STKOVER 2 #define PNC_ZERODIV 3 #define PNC_NUMTAX 4 #define PNC_SYNTAX 5 #define PNC_ARG 6 #define PNC_ERRORS 15 #define PNC_INIT 256 #define PNC_HELP 512 #define PNC_STACK 1024 #define PNC_MEMORY 2048 #define PNC_MEMUPD 4096 int pnc (PNC *cs, char *buf); /*---------------------------------------------------------------- int pnc (PNC *cs, char *input) processes an input string for the Polish notation calculator and updates the calculator's struct appropriately and returns the number of input characters processed successfully. If pnc returns k, the state of the calculator has been updated based on the first k characters of the input string and (input+k)... were ignored. If (input+k) is a newline or null, the complete string was presumably processed; otherwise, a cryptic error message is placed in the msg component. Errors in the C mathematical functions are handled by the default mechanism, i.e., pnc does not attempt to get involved. Polish notation (Lukasiewicz) uses parenthesis free expressions. It was originally developed for logical expressions but is clearly applicable to arithmetic expressions. Thus, for example, a+(b+c) is either a b c ++ or b c + a + while (a+b)+c is either c a b ++ or a b + c +. Formally, all operands should occur before the operators, but the mingled expressions are commonly used. Nevertheless, the notation is naturally implemented with a stack, i.e., when encountered, an operand is placed on the top of the stack, while operators are applied to the number(s) on top of the stack. Alphabetic operators, e.g., sqrt, must be followed by a character that is not alphabetic, while single character operators, e.g., an asterisk (*), may be followed by any operator without intervening white space. All whitespace characters are ignored, e.g., blanks, but serve as delimiters of course. e and pi correspond to the mathematical constants these symbols are commonly used to represent; e is exp(1) and pi 4*atan(1). eps is the machine epsilon, DBL_EPSILON, i.e., the smallest number such that 1 + eps is not 1. Unary operators remove the operand on top of the stack, perform the indicated operation on this operand, and place the result on the top of the stack, i.e., = operator . If the stack is empty, any unary operator causes an operand error. Most of the standard C mathematical functions fall in this group: abs, chs (change sign), sqrt, the six trigonometric functions, exp and ln (log), and log (log10). Binary operators remove the top two operands on the stack, perform the operation on these two operands, and place the result on the top of the stack, i.e., = operator so that the size of the stack is reduced by one. If the stack does not contain at least two operands, any binary operator causes an operand error. The four arithmetic operators and the pow function fall in this category. The eight memory locations in the PNC struct are denoted by the letters s,...,z and like the alphabetic operators must be followed by a nonalphabetic character, e.g., a blank. Values are stored in memory by entering the character immediately followed by an equal sign (=), e.g., s= stores a value, while s+ adds s to the top of the stack. This is inelegant but consistent with the notation and allows multiple stores, e.g., s=t=. The special characters >, <, ^, !, and xch function perform stack operations. The > and < operators rotate the stack with the top (bottom) of the stack moving to the bottom (top); no data is lost. The ^ operator duplicates the top of the stack; while ! eliminates the top of the stack. xch exchanges the top two stack elements. clr clears the calculator. A few special functions are implemented when a question mark (?) is the first character of the operator. If it's just a ?, the help string is returned in input and the HELP flag is set. If it's "?stack" or "?memory" (abbreviations accepted), bits are set in the flags component. These two operators are used to trigger the windows for the stack and memory, respectively. */ int pnc (PNC *cs, char *buf) BEGIN double s0, *stack; int i, k, p, size; char *e, *s, *input; char *cec[] = { " ", "OPERANDS? ", "STACK OVERFLOW ", "ZERO DIVISOR ", "ILLEGAL NUMBER ", "SYNTAX ", "ILLEGAL ARGUMENT " }; char *table = "s t u v w x y z e pi eps " "abs sqrt cos sin tan acos asin atan ln exp log pow " "chs xch clr "; char *help = "\nConstants: e pi eps\n" "Operators: + - * / > < ^ !\n" " clr chs xch ? ?stack ?memory\n" " abs sqrt exp ln log pow\n" " cos acos sin asin tan atan\n"; if (cs->flags & PNC_INIT) { for (k = 0; k < cs->stacksize ; ++k) cs->stack[k] = 0.0; for (k = 0; k < cs->memsize; ++k) cs->memory[k] = 0.0; cs->stackptr = -1; } cs->msg = NULL; cs->flags = 0; if (! buf) return 0; input = buf; stack = cs->stack; size = cs->stacksize; p = cs->stackptr; if (p >= size) { cs->flags |= PNC_STKOVER; return 0; } while (*input && *input != '\n') BEGIN if (isdigit(*input) || *input == '.') BEGIN errno = 0; s0 = strtod(input,&e); if (errno || (*e && ! isspace(*e))) { cs->flags |= PNC_NUMTAX; break; } if (p+1 >= size) { cs->flags |= PNC_STKOVER; break; } else stack[++p] = s0; input = e; END else if (*input == '+' || *input == '-') BEGIN if (isdigit(*(input+1)) || *(input+1) == '.') { errno = 0; s0 = strtod (input,&e); if (errno || (*e && ! isspace(*e))) { cs->flags |= PNC_NUMTAX; break; } if (p+1 >= size) { cs->flags |= PNC_STKOVER; break; } else stack[++p] = s0; input = e - 1; } else if (*input == '+') { if (p < 1) { cs->flags |= PNC_NOARGS; break; } --p; stack[p] = stack[p] + stack[p+1]; } else if (*input == '-') { if (p < 1) { cs->flags |= PNC_NOARGS; break; } --p; stack[p] = stack[p] - stack[p+1]; } ++input; END else if (*input == '*') { if (p < 1) { cs->flags |= PNC_NOARGS; break; } --p; stack[p] = stack[p] * stack[p+1]; ++input; } else if (*input == '/') { if (stack[p] == 0.0) { cs->flags |= PNC_ZERODIV; break; } if (p < 1) { cs->flags |= PNC_NOARGS; break; } --p; stack[p] = stack[p] / stack[p+1]; ++input; } else if (isalpha(*input)) BEGIN for (k = 1, e = input+1; *e && isalpha(*e); ++e, ++k); for (i = 0, s = table; *s; i++) { if (*(s + k) == ' ' && ! strncmp(input,s,k)) break; while (*s++ != ' '); } if (! *s) { cs->flags |= PNC_SYNTAX; break; } if (*e == '=') /* Store into memory.*/ { if (i > 7) { cs->flags |= PNC_SYNTAX; break; } else if (p < 0) { cs->flags |= PNC_NOARGS; break; } else { cs->memory[i] = stack[p]; cs->flags |= PNC_MEMUPD; ++e; }} else if (i < 11) /* Load from memory or constant.*/ { if (p+1 >= size) cs->flags |= PNC_STKOVER; else if (i < 8) stack[++p] = cs->memory[i]; /* memory */ else if (i == 8) stack[++p] = exp(1.0); /* e */ else if (i == 9) stack[++p] = 4.0*atan(1.0); /* pi */ else { s0 = 1.0 + (stack[++p] = 0.5); while (s0 != 1.0) {stack[p] = 0.5*stack[p]; s0 = 1.0 + stack[p]; } stack[p] = DBL_EPSILON; }} else if (p >= 0) switch (i-11) BEGIN case 0: /* abs */ stack[p] = fabs(stack[p]); break; case 1: /* sqrt */ if (stack[p] < 0.0) cs->flags |= PNC_ARG; else stack[p] = sqrt(stack[p]); break; case 2: /* cos */ stack[p] = cos(stack[p]); break; case 3: /* sin */ stack[p] = sin(stack[p]); break; case 4: /* tan */ stack[p] = tan(stack[p]); break; case 5: /* acos */ if (abs(stack[p]) > 1.0) cs->flags |= PNC_ARG; else stack[p] = acos(stack[p]); break; case 6: /* asin */ if (abs(stack[p]) > 1.0) cs->flags |= PNC_ARG; else stack[p] = asin(stack[p]); break; case 7: /* atan */ stack[p] = atan(stack[p]); break; case 8: /* ln */ stack[p] = log(stack[p]); break; case 9: /* exp */ stack[p] = exp(stack[p]); break; case 10: /* log */ stack[p] = log10(stack[p]); break; case 11: /* ** */ if (p < 1) cs->flags |= PNC_NOARGS; else { --p; stack[p] = pow(stack[p],stack[p+1]); } break; case 12: /* chs */ stack[p] = -stack[p]; break; case 13: /* xch */ if (p > 0) { s0 = stack[p]; stack[p] = stack[p-1]; stack[p-1] = s0; } break; case 14: /* clr */ for (k = 0; k < cs->stacksize ; ++k) cs->stack[k] = 0.0; p = cs->stackptr = -1; break; default: cs->flags |= PNC_SYNTAX; END /* switch i */ if (cs->flags & PNC_ERRORS) break; input = e; END /* if alphabetic */ else if (isspace(*input)) ++input; else if (*input == '?') BEGIN for (k = 1, e = input+1; *e && isalpha(*e); ++e, ++k); if (k == 1) { cs->flags |= PNC_HELP; cs->msg = help; } else if (k <= 6 && !strncmp(input+1,"stack",k-1)) cs->flags |= PNC_STACK; else if (k <= 7 && !strncmp(input+1,"memory",k-1)) cs->flags |= PNC_MEMORY|PNC_MEMUPD; else cs->flags |= PNC_SYNTAX; input = e; END else if (ispunct(*input)) BEGIN if (*input == '>') { s0 = stack[p]; for (k = p; k > 0; --k) stack[k] = stack[k-1]; stack[0] = s0; } else if (*input == '<') { s0 = stack[0]; for (k = 0; k < p; ++k) stack[k] = stack[k+1]; stack[p] = s0; } else if (*input == '^') { if (p+1 >= size) { cs->flags |= PNC_STKOVER; break; } else if (p < 0) stack[p = 0] = 0.0; else { ++p; stack[p] = stack[p-1]; }} else if (*input == '!') { if (p >= 0) --p; } else { cs->flags |= PNC_SYNTAX; break; } ++input; END else { cs->flags |= PNC_SYNTAX; break; } END /* while on *input */ cs->stackptr = p; if (cs->flags & PNC_ERRORS) cs->msg = cec[cs->flags & PNC_ERRORS]; return (int)(input - buf); END /*---------------------------------------------------------------- This application is quite straightforward because all it has to do is read the user's input and print the results; all the real work is done by the function pnc. The I/O processing, however, uses a couple tactics that you cannot do easily with C stream I/O. The first concerns echoing the user's input; the second, the handling of the stack and memory windows. Null lines are ignored. An EOF causes termination, but the window is not closed. Stack and Memory Windows The calculator includes stack and memory display operators, which simply set the appropriate flags for the application. Normally, once requested, the stack window is updated after each input line is processed, while the memory window is updated only if pnc sets the MEMUPD flag, which indicates memory was changed by the input. When updating these windows, the cputs return is examined to see if the window has been closed, a return of -2. If so, the open identifier is set null, and the window is not updated until again requested. In the normal course of things, this suffices. But the user may close one of these windows and immediately request it again. To deal with this possibility, if the STACK or MEMORY flag is set and the open identifier is not null, which normally means that the window already exists, cputs is used to determine whether the window really exists. The stack and memory windows are opened with minimal size output streams and truncated as necessary, i.e., the oldest output is discarded. The code uses both the \f and \0 control characters in output text passed to cputs. Their use is not essential but helps to clarify the displays of the stack and memory. If the output in these windows is directed to a file, the file will be complete. The first time you run this application, you'll want to configure the windows so that the stack and memory windows don't overlap the calculator window. The height of the stack and memory windows is best set using 16 and 8, respectively, so that you don't have to scroll them. You'll likely also want to use a fixed font. */ int appl (void) BEGIN PNC calc; void *wstk, *wmem, *wnd; int i, k, mi; char input[32*PNC_STKSIZE], result[32]; char memnames[] = { 's','t','u','v','w','x','y','z' }; cappl("PNCalc"); wstk = wmem = NULL; wnd = copen("PNC:Polish Notation Calculator",1,1); calc.memsize = PNC_MEMSIZE; calc.stacksize = PNC_STKSIZE; calc.flags = PNC_INIT; mi = 24*PNC_STKSIZE; /* The window used to solicit input and print the results is opened with the default modifiers; thus the user's input is not echoed automatically. If pnc returns k and k is positive, characters 0...(k-1) of the user's input were successfully processed by pnc and they are displayed in the window followed by the value on the top of the stack, i.e., the result. The portion of the line that could not be processed, if any, is not echoed. In general, the k-th character is the newline at the end of the user's input. */ while ((i = cgets(wnd,input,mi)) >= 0) BEGIN k = pnc(&calc,input); if (k) { cputs(wnd,input,k); cputs(wnd,"\n",1); if (calc.stackptr >= 0) cputs(wnd,result,sprintf(result, "%25.17e\n",calc.stack[calc.stackptr])); } /* When pnc encounters an error, it simply stops processing input and sets an error flag. If it returns k, the error occurred while trying to process the input starting with the k-th character. In this case, the cryptic comment followed by the input that was not processed is displayed. Using the line re-entry feature, the right mouse button, you can pick up this error output, correct it, and press Enter. */ if (calc.flags & PNC_ERRORS) { cputs(wnd,calc.msg,0); cputs(wnd,input+k,0); } else if (calc.flags & PNC_HELP) cputs(wnd,calc.msg,0); /* Stack Window If pnc finds the stack display operator in the line, it sets the STACK flag. Once requested, the stack window is updated after each input line, need it or not. When updating the stack, the complete stack is formed in a buffer prefixed with a /f so that the window update places the first line of the updated stack at the top of the window; normally, the window would be updated to show the last line at the bottom. This effectively removes the previous stack displays from the window, but the user can scroll back through them. */ if (calc.flags & PNC_STACK) { if (wstk && cputs(wstk,NULL,0) == -2) wstk = NULL; if(! wstk) wstk = copen("Stack:PNC Stack",0,-1); } if (wstk) { input[0] = '\f'; k = 1; for (i = calc.stackptr; i >= 0; --i) k += sprintf(input+k, "%2i:%25.17e\n",i,calc.stack[i]); input[k++] = '\n'; if (cputs(wstk,input,k) == -2) wstk = NULL; } /* Memory Window If pnc finds the memory display operator in the line, it sets the MEMORY and MEMUPD flags. If it changes a memory location, it sets the MEMUPD flag. If the window exists, i.e., the open identifier is not null, it's updated. When updating the memory window, the complete memory display is formed in a buffer prefixed with a /f and a \0. As with the stack update, the \f causes the window to be updated from the top of the window. The \0, which cannot be the first character in the text passed to cputs, causes the VM file used to hold the output stream to be emptied before being updated, which effectively discards the previous memory display. Thus the user cannot scroll backwards to find previous displays. */ if (calc.flags & PNC_MEMORY) { if (wmem && cputs(wmem,NULL,0) == -2) wmem = NULL; if (! wmem) wmem = copen("Memory:PNC Memory",0,-1); } if (wmem && calc.flags & PNC_MEMUPD) { input[0] = '\f'; input[1] = '\0'; k = 2; for (i = 0; i < PNC_MEMSIZE; ++i) k += sprintf(input+k, "%c:%25.17e\n",memnames[i],calc.memory[i]); if (cputs(wmem,input,k) == -2) wmem = NULL; } END /* while */ return 0; END #endif #ifdef NICMap /*================================================================ Convergence Map for Newton's Iteration This application reads/generates a complex polynomial, reads the parameters describing a rectangular region of the complex plane, and then draws a map of this rectangle based on the convergence of Newton's iteration for the zeros of the polynomial when the points in the rectangular region are used to start the iteration. The map is drawn in a window as a printer plot, which is an ancient but honorable technique for drawing pictures with a computer. The quality of the map, of course, leaves a great deal to be desired, but then this monograph is about text windows. If you press Enter often enough, the default problem is run, i.e., the polynomial x**4 - 1, the rectangle -1 <= x,y <= 1, and a horizontal grid size of .01. This problem produces a set of heart shaped regions along the lines x = y and x = -y that look fractal in nature. A quite spectacular color picture for this problem can be found in Ivars Peterson's book "The Mathematical Tourist" (New York: W. H. Freeman, 1988). Indeed, after seeing it, I wrote an application that produced pixel plots for various single-point iterations that has been lost for lack of portability; this code was extracted from its remnants. After running the default problem, you can try a few variations in the rectangle and grid. If you enter an EOF when prompted for a new rectangle, you are given a chance to enter a new problem. The default polynomial (press Enter) for any degree is interesting. The map window should be configured to use a fixed font so that the horizontal spacing of the grid 'points' is uniform. Because the map window is not opened until necessary, you will not have a chance to set the font until the first map has already been drawn. You can change the font by moving to the map window and issuing an appropriate FONT command and then use the same rectangle again by pressing Enter a few times. Subsequently, the font will be in the memory file. The maps are most easily viewed by using the left mouse button to drag the map around the window; the map size is not restricted to the window size. Horizontally, the grid can have at most 1000 points. Actually, the value you specify is adjusted to give at least 10 but not more than 1000 grid points. The vertical size of the plot depends on the font size, horizontal grid size, and the y extent of the rectangle. In this respect, the size is limited only by the virtual memory available. Thus, the plot will almost always be larger than the window. The odd structures in the maps stem principally from the iteration; other single-point iterations yield maps with entirely different types of structures. If you have a decent processor and video board, this example is not processor intensive. On the other hand, if you try to look at large rectangles or use a small grid, the process may slow and the virtual memory requirements may be large. To reduce the cost of drawing the map window after each line is added, you can press F6, which suspends the updating until you press F6 again. After the map is completed, the zeros are printed in both windows so that you'll know when the map is done. Distance, e.g., convergence, is based on the sum of the absolute values of the differences between the real and imaginary parts (the L1 norm) because it is cheaper than a complex absolute value. Convergence is based on the usual absolute/relative difference, absolute if inside the unit square and relative outside. The same technique is used to determine whether a new zero is already known. */ #include #include #include #include #include #include #define MIN_EPS DBL_EPSILON #define MAX_EPS 1.0e-6 /* Prototypes for Internal Functions */ int ReadProblem (void *w); int OldeNewt (double x, double y); /* These macros are used within the main program of the application but not in the ReadProblem function. The BSO...ESO model takes both the string and its length, but the latter can be zero if the length is not readily available. The buffer used with the BFO macro provides a page (4096 bytes). */ #define BSO(w) cputs(w, #define ESO ); #define BFO(w) cputs(w,IOB,sprintf(IOB, #define EFO )); /* This application extracts information directly from the window struct, i.e., it uses the fact that copen returns the pointer (void *) to the WND struct. From the outset, unless there was a clear reason to prohibit something, the POO software erected no barriers. There was no intent to use the open identifier in the examples, but it is clearly handy. In this case, the font size is used to adjust the vertical grid spacing to control the aspect of the rectangular map. The desired relationship is * = *, which compensates for the aspect of the font but not the screen. This abbreviated struct will not work in Windows 3.1 because one of the 32-bit things is only a 16-bit thing. But this example is very limited in Windows 3.1 anyway because the map has to fit in 64K of memory. If you are using Windows 3.1, I would suggest the Win32s interface, which avoids both problems. */ typedef struct { void *next[24]; /* A bunch of 32-bit things */ int ww, wh; /* Window dimensions */ int tw, th; /* Font sizes */ } WND; /* Global Variables The maximum degree of the polynomial is set here at 16 but is not very limiting because Newton's method usually falls apart before this limit comes into play. The map produced by this example doesn't correspond to the one that would be obtained in theory because MaxStep and MaxIter affect the results. */ #define maxN 16 #define maxZ (2*maxN+1) double X0, X1; /* x range of rectangle, X0 <= X1 */ double Y0, Y1; /* y range of rectangle, Y0 >= Y1 */ double Xdel,Ydel; /* Grid point increments */ double MaxStep; /* Max change per iteration */ double Epsilon; /* Absolute/relative error measure */ double Ar[maxN], Ai[maxN]; /* Complex coefficients, 0...N */ double Zr[maxZ], Zi[maxZ]; /* Known complex zeros, 1...N */ int N; /* Degree of the polynomial */ int Known; /* Number of known zeros */ int MaxIter; /* Max iterations before failure */ char Line[1024]; /* Next line of character 'pixels' */ char IOB[4096]; /* I/O buffer */ char *Pixel = "0+-\\/|=.:<>abcdefghijklmnopqrstuvwxyz"; /*---------------------------------------------------------------- The application is structured as two (infinite) loops. The outer loop reads the polynomial, the convergence criterion, the number of Newton iterations before declaring failure, and the maximum allowed change in the iterate per iteration. Theoretically, of course, the last three are not part of Newton's iteration but are essential in practice and significantly impact the convergence of the iteration. Varying the maximum number of iterations and change per iteration may change both the map and the time required to generate it dramatically. The inner loop reads the x and y ranges that define the rectangle and the horizontal grid spacing; the vertical grid spacing is determined from the horizontal and the current font size. Thus, you can enter a problem and draw numerous rectangles based on it, and then exit the inner loop and enter another problem, etc. To exit the inner loop, enter an EOF (Ctrl+Z). The application will terminate if you enter an EOF when asked for the degree, which is in the outer loop. When starting a new problem, you can get a good feel for what is going to happen by entering the rectangle of interest and a large grid size. The map will be quite coarse, but it only takes two Enters and a new grid value to refine the grid. The abscissa and ordinate values are printed in the format sxxxx.xxx so that you are limited to rectangles in this part of the complex plane. */ int appl (void) BEGIN void *w, *map; int k, j, h, n, m; char *s; double x, y, b, t; cappl("NICMap"); w = copen("E@IO:Newton's Iteration;",1,1); if (! w) return 0; BSO(w) "POO Application Newton\n\n" "If you simply press Enter whenever asked for\n" "input, the default problem is run, i.e.,\n" "the polynomial X**4 - 1 and the rectangle\n" " -1 <= x <= 1 and -1 <= y <= 1.\n", 0 ESO N = 4; Epsilon = MIN_EPS; MaxIter = 50; MaxStep = 1000.; map = NULL; X0 = Y0 = -1.0; X1 = Y1 = 1.0; Xdel = 0.01; /* The outer loop, which reads the problem, continues unless you close the I/O window, but ReadProblem may return an EOF. After calling ReadProblem to get the problem, it's printed in the I/O window so that you can see the problem when solicited for the rectangle. Thus, if the problem is not the one desired, you can abort and re-enter the problem by responding to the 'x range?' prompt with an EOF. If the polynomial is specified by its zeros, they are printed. Nevertheless, Known is then set to zero so that the zeros that are printed after each map consist of only those found with Newton's iteration. */ while (cputs(w,NULL,0) >= -1) BEGIN if (ReadProblem(w) < 0) break; BFO(w) "\n\nNewton's Iteration for the complex polynomial\n" "of degree %i with coefficients:\n",N EFO for (k = 0; k <= N; ++k) BFO(w) " x ^ %2i: %20.11le +i* %20.11le\n",k,Ar[k],Ai[k] EFO if (Known >= 1) { BSO(w) "with the specified zeros\n",0 ESO for (k = 1; k <= Known; ++k) BFO(w) " %20.11le +i* %20.11le\n",Zr[k],Zi[k] EFO } BFO(w) "Convergence criterion %le with at most %i iterations and\nmax " "change per iteration of %9.3lf.\n\n",Epsilon,MaxIter,MaxStep EFO Known = 0; /* With the problem now well defined, the inner infinite loop is initiated and continues until you respond to an input request with an EOF. Within this loop, this example solicits the rectangle and horizontal grid size and then produces the appropriate map. */ while (1) BEGIN /* The rectangle is defined by its horizontal (x) and vertical (y) ranges; the code handles order problems. The x grid size, Xdel, is limited to the range [.1,.001] times the x extent and is set to the appropriate end point if outside this range. The actual limit on x, b, is X1 + 0.5*Xdel to allow for the inevitable rounding errors. The vertical grid size is set later based on the font size and the horizontal grid. The prompt line that is used to solicit the ranges is replaced after the actual values are entered even though input echoing was selected when the window was opened, which is done primarily to illustrate the mechanism. Omitted values leave the pertinent variables unchanged. An EOF response to any prompt causes an exit to read a new problem. */ BFO(w) "x range [%9.3lf,%9.3lf]?",X0,X1 EFO if (cgets(w,Line,64) < 1) break; sscanf(Line,"%lf %lf",&X0,&X1); if (X0 > X1) { x = X0; X0 = X1; X1 = x; } BFO(w) "\rx range [%9.3lf,%9.3lf]\n",X0,X1 EFO BFO(w) "y range [%9.3lf,%9.3lf]?",Y0,Y1 EFO if (cgets(w,Line,64) < 1) break; sscanf(Line,"%lf %lf",&Y0,&Y1); if (Y0 > Y1) { y = Y0; Y0 = Y1; Y1 = y; } BFO(w) "\ry range [%9.3lf,%9.3lf]\n",Y0,Y1 EFO BFO(w) "x grid size[%9.3lf]?",Xdel EFO if (cgets(w,Line,64) < 1) break; sscanf(Line,"%lf",&Xdel); if ((Xdel = fabs(Xdel)) > 0.1*(X1 - X0)) Xdel = (X1 - X0)/10.0; else if (Xdel < 0.001*(X1 - X0)) Xdel = (X1 - X0)/1000.0; BFO(w) "\rx grid size %9.3lf\n",Xdel EFO b = X1 + 0.5*Xdel; /* Opening the map window is deferred until this point so that it can be integrated with the code for creating it if the user closed it. The first time through, map is NULL. Subsequently, it's not NULL but is tested to make sure that the user didn't close it. */ if (! map || cputs(map,NULL,0) < -1) { map = copen("Map:Convergence Map;",32767,4096); if (! map) { BSO(w) "Unable to open map window!",0 ESO break; }} /* The output to the map window consists of the problem statement, the map, and the zeros that were found. The following code adds the problem statement. The map window could be emptied here by issuing a cputs(map,"\f\0",2), but then you could not scroll back to previous maps. */ BFO(map) "\nPrinter Plot for Netwon's Iteration\n\n" "Polynomial of degree %i with coefficients:\n",N EFO for (k = 0; k <= N; ++k) BFO(map) " x ^ %2i: %20.11le +i* %20.11le\n",k,Ar[k],Ai[k] EFO BFO(map) "Convergence criterion %le with at most %i iterations and\n" "maximum step of %9.3lf.\n",Epsilon,MaxIter,MaxStep EFO /* The map begins and ends with the abscissa line that is generated here. The abscissa line is placed in IOB and saved so it can be put at the bottom of the map. */ m = sprintf(IOB,"\n %9.3lf|",X0); for (x = X0 + 10.0*Xdel; x < b; x += 10.0*Xdel) m += sprintf(IOB+m,"%9.3lf|",x); IOB[m++] = '\n'; BSO(map) IOB,m ESO /* The map is computed using a double loop with the outer loop on y and the inner on x so that it starts at the upper left and ends at the lower right. The value of Ydel, the vertical grid spacing, is determined from Xdel and the font size so that the aspect ratio of the map is correct. A fixed size font is clearly recommended for the map window since the horizontal spacing will be somewhat messy otherwise. Diving into the WND struct is definitely not recommended. For example, a malicious user could close the map window, which would leave map invalid and likely lead to a fatal error. Opening the window with the C modifier, which prevents the user from closing it, is only slightly safer. */ Ydel = (((WND *)map)->th)*Xdel/(((WND *)map)->tw); for (y = Y1; y > Y0 - Ydel; y -= Ydel) BEGIN /* The OldeNewt function returns the number of iterations it used if the iteration converged and a negative value otherwise. Further, if the iteration converged, OldeNewt places the zero in the next available element of (Zr,Zi). Thus, setting the character for x is straightforward: if the iteration fails, an asterisk (*) is used; otherwise, the zero WILL be found in the table, and the character is set by simply indexing into the map character set. Unfortunately, however, there is a caveat. Although it is true that a polynomial of degree N has precisely N zeros (roots) in mathematics, this is numerical computation. When OldeNewt finds a zero, it means only that the iteration converged numerically. Although one might expect fewer than N zeros, more than N can arise if the convergence criterion is poorly chosen or if the polynomial is ill conditioned. Only the first 2N zeros are saved; after that 0 is used to show that the iteration converged to yet another spurious something. Matching the new zero against the known ones is done with MAX_EPS rather than Epsilon to avoid spurious zeros, where 'spurious' implies only late relative to the other zeros and does not reflect on their quality. */ h = sprintf(Line,"%9.3lf ",y); s = Line + h; for (x = X0; x < b; x += Xdel) if (OldeNewt(x,y) < 0) *s++ = '*'; else { t = fabs(Zr[Known+1]) + fabs(Zi[Known+1]); if (t < 1.0e0) t = 1.0e0; t = t*MAX_EPS; for (k = 1; k <= Known; ++k) if (fabs(Zr[k] - Zr[Known+1]) + fabs(Zi[k] - Zi[Known+1]) < t) break; if (k > Known) k = (Known < 2*N) ? ++Known : 0; *s++ = *(Pixel+k); } *s++ = ' '; n = (int)(s - Line); /* After the inner loop is completed, any known zero in the map line is replaced by a pound sign (#). I suppose there may be cases in which a line contains a zero but none of the grid points converged to it, but this would be unusual. After the outer loop terminates, the abscissa line is appended at the bottom of the map. */ for (k = 1; k <= Known; ++k) if (y - Ydel < Zi[k] && Zi[k] <= y) { j = h + ceil((Zr[k] - X0)/Xdel); if (h <= j && j < (n-1)) Line[j] = '#'; Line[h-1] = '>'; Line[n-1] = '<'; } n += sprintf(Line+n," %9.3lf\n",y); BSO(map) Line,n ESO END BSO(map) IOB,m ESO /* Abscissa line.*/ /* Finally, the known zeros are appended to both the I/O window and the map window. If you pressed F6 to suspend window updating, the printing of these zeros in the I/O window is a signal that the map is done, i.e., it's time to press F6 in the map window. */ h = sprintf(IOB,"\nZeros previously known or located:\n"); for (k = 1; k <= Known; ++k) h += sprintf(IOB+h," %20.11le +i* %20.11le\n",Zr[k],Zi[k]); IOB[h++] = '\n'; BSO(w) IOB,h ESO BSO(map) IOB,h ESO END /* of a grid map.*/ END /* of problem.*/ return 1; END /*---------------------------------------------------------------- ReadProblem reads the next problem: the polynomial, which can be specified in a number of ways, the convergence criterion, the maximum number of of iterations before failure, and the maximum step in the iterate per iteration. This last parameter is used to avoid successive iterates that would otherwise move at warp speed, i.e., iterates close to a zero of the derivative. The polynomial can be specified by entering its coefficients, its zeros, or the index of a power series that should be truncated to obtain the polynomial. This function is coded directly in terms of the API functions and assumes that your responses to the prompting messages are echoed, i.e., the E modifier was specified when the window was opened, so that they are not explicitly echoed by this function. The main program prints a complete statement of the problem. */ int ReadProblem (void *w) BEGIN int n, k, j; char *buf, *s; double t; buf = IOB; /* The degree of the polynomial must be a positive integer between 2 and maxN. An EOF response for the degree causes termination, i.e., this function returns -1. The second value read may be { c | C | z | Z | 1-4 | anything else } and indicates whether you want to enter the coefficients or zeros or use a truncated power series or the default problem, x**N - 1. */ do { cputs(w,"\nDegree of the polynomial? ",0); if (cgets(w,buf,64) < 1) return -1; sscanf(buf,"%i",&N); } while (N < 2 || N >= maxN); Known = 0; cputs(w,"\n{Series Number|Coefficients|Zeros}? ",0); cgets(w,buf,64); for (s = buf; *s && *s == ' '; ++s); /* Coefficients The k-th complex coefficient is set to 1.0 before it's read. If you enter nothing, it remains 1.0. If you enter only one value, the coefficient is real because the imaginary part remains 0.0. If the coefficient of x**N is zero, it's set to 1.0; otherwise, the polynomial wouldn't be of degree N. The default coefficients correspond to a truncated power series that converges inside the unit disk, i.e., the comments below apply to its zeros. */ if (*s == 'c' || *s == 'C') BEGIN cputs(w,buf,sprintf(buf,"\nComplex coefficients 0...%i:\n",N)); for (k = 0; k <= N; ++k) /* Zero unless provided.*/ { cgets(w,buf,64); Ar[k] = 1.0; Ai[k] = 0.0; sscanf(buf,"%le %le",&Ar[k],&Ai[k]); } if (Ar[N] == 0.0 && Ai[N] == 0.0) Ar[N] = 1.0; END /* Zeros The k-th complex zero is set to k before it's read so that if two legitimate values are not entered, the zero has default values. Thus, if you want a real zero, you only need to enter the real part. After reading the N zeros (Zr,Zi), the coefficients (Ar,Ai) are computed. The default zeros are somewhat nasty because the coefficients grow rapidly with N and the zeros are very sensitive with respect to perturbations in the coefficients. For N small, no rounding errors are incurred while computing the coefficients so that they are exact; but, as N nears 13, the situation changes. An analysis of these polynomials can be found in "Rounding Errors in Algebraic Processes" by J.H. Wilkinson (Englewood Cliffs: N.J., Prentice- Hall, 1963). The maps produced for the default zeros are striped, but the boundaries are interesting. Even for small n, Newton's iteration fails for a lot of the grid points so that *'s abound. For small n, the most interesting rectangle is [1,n] and [-.2,.2] with small grid size, which shows interesting little structures midway between the zeros. Polynomials with multiple zeros are particularly troublesome for Newton's iteration. In some cases, multiple zeros are found as a group of almost equal zeros; in others, fewer than N zeros are located. The results depend on the parameters, e.g., the maximum step and maximum iterations, as much as the polynomial. */ else if (*s == 'z' || *s == 'Z') BEGIN cputs(w,buf,sprintf(buf,"\nEnter the %i complex zeros" " one per line:\n",N)); for (k = 1; k <= N; ++k) { cgets(w,buf,64); Zr[k] = k; Zi[k] = 0; sscanf(buf,"%le %le",&Zr[k],&Zi[k]); } Known = N; Ar[N] = 1.0e0; Ai[N] = 0.0e0; Ar[N-1] = -Zr[1]; Ai[N-1] = -Zi[1]; for (k = 2; k <= N; ++k) { Ar[N-k] = Ai[N-k] = 0.0e0; for (j = N-k; j < N; j++) { Ar[j] = Ar[j] - (Ar[j+1]*Zr[k] - Ai[j+1]*Zi[k]); Ai[j] = Ai[j] - (Ar[j+1]*Zi[k] + Ai[j+1]*Zr[k]); }} END /* Truncated Power Series Four series are provided; any value other than 1-4 elicits the default problem, x**N - 1, which is not a truncated power series. Series 1-3 are real, while 4 is complex. All converge only inside the unit disk. Truncated power series are interesting because their zeros cluster around the circle of convergence, see Hurwitz's Theorem in "The Theory of Functions" by E.C. Titchmarsh (London: Oxford Univ. Press, 1939). For truncated power series, the rectangle -2,2 and -2,2 with a grid of .05 is a good starting point. As the degree increases (n > 6), series 1, 3, and 4 tend to produce flowers, while series 2 is more individualistic. */ else BEGIN for (n = 0; isdigit(*s); ++s) n = 10*n + (*s - '0'); Ar[0] = 1.0; for (k = 0; k <= N; ++k) Ai[k] = 0.0; switch (n) BEGIN case 1: for (k = 1; k <= N; k++) Ar[k] = 1.0/k; break; case 2: Ar[0] = Ar[1] = 1.0; for (k = 2; k <= N; k++) Ar[k] = 1.0/((k-1)*(k+1)); break; case 3: for (t = 1.0, k = 0; k <= N; k++) { Ar[k] = t/((k+1)*(k+2)); t = -t; } break; case 4: for (t = 2.0*3.1415926/n, k = 1; k <= N; k++) { Ai[k] = sin(k*t); Ar[k] = cos(k*t); } break; default: n = 0; for(k = 1; k < N; ++k) Ar[k] = 0.0; Ar[N] = 1.0; Ar[0] = -1.0; break; END cputs(w,buf,sprintf(buf," Ok, series %i\n",n)); END /* It remains to get the parameters that control the iteration, which should be viewed as part of the problem statement because of their impact on the map. In each case, the prompt includes the current value, and you can respond with Enter if it's acceptable. All of these parameters are restricted. If you violate the restriction, it's changed to a suitable value without comment. An EOF response has the same effect as entering nothing because cgets returns the null string in the input buffer. */ cputs(w,buf,sprintf(buf,"\nConvergence criterion %le? ",Epsilon)); cgets(w,buf,64); sscanf(buf,"%lf",&Epsilon); if (Epsilon < MIN_EPS || Epsilon > MAX_EPS) Epsilon = MIN_EPS; cputs(w,buf,sprintf(buf,"\nMaximum iterations %i? ",MaxIter)); cgets(w,buf,64); sscanf(buf,"%i",&MaxIter); if (MaxIter < 4 || MaxIter > 500) MaxIter = 50; cputs(w,buf,sprintf(buf,"\nMaximum step %f? ",MaxStep)); cgets(w,buf,64); sscanf(buf,"%lf",&MaxStep); if (MaxStep < 1.0) MaxStep = 100.0; return N; END /*---------------------------------------------------------------- OldeNewt (x,y) applies Newton's iteration starting at (x,y) based on the complex polynomial of degree N with coefficients (Ar,Ai) and returns the number of iterations performed if it converges and minus this number if it fails to converge in MaxIter iterations. The final iterate is returned in (Zr[Known+1],Zi[Known+1]) if successful. Each iteration is restricted by requiring that ||P(z)/P'(z)|| be less than MaxStep; otherwise, P'(z) at the previous iteration is used (1.0 for the first iteration). Actually, any reasonable constant can be substituted for the derivative without doing significant damage to the iteration. */ int OldeNewt (double x, double y) BEGIN int iterate, k; double ldr, ldi, pr, pi, qr, qi, t; ldr = 1.0e0; ldi = 0.0e0; for (iterate = 0; iterate < MaxIter; ++iterate) BEGIN /* The polynomial of degree N with complex coefficients (Ar,Ai) is evaluated at (x,y) using Horner's rule to obtain both its value (pr,pi) and first derivative (qr,qi). */ pr = qr = Ar[N]; pi = qi = Ai[N]; for (k = N-1; k > 0; --k) BEGIN t = (pr*x - pi*y) + Ar[k]; pi = (pr*y + pi*x) + Ai[k]; pr = t; t = (qr*x - qi*y) + pr; qi = (qr*y + qi*x) + pi; qr = t; END /* for */ t = (pr*x - pi*y) + Ar[0]; pi = (pr*y + pi*x) + Ai[0]; pr = t; /* The increment for the next iterate is (pr,pi)/(qr,qi), but it can't be used blindly because it could be almost anywhere if (x,y) is too close to a zero of the derivative. In any case, it's wise to avoid a division by zero. If the quotient would be too large, the value of the derivative at the previous iterate, (ldr,ldi), is used. The quotient is formed by inversion and multiplication, and the inverse is saved in (ldr,ldi) in case it's needed the next time. */ if (MaxStep*(fabs(qr)+fabs(qi)) > fabs(pr)+fabs(pi)) { t = qr*qr + qi*qi; ldr = qr/t; ldi = -qi/t; } x -= (qr = pr*ldr - pi*ldi); y -= (qi = pr*ldi + pi*ldr); qr = fabs(qr) + fabs(qi); t = fabs(x) + fabs(y); if (qr < ((t > 1.0e0) ? t*Epsilon : Epsilon)) { Zr[Known+1] = x; Zi[Known+1] = y; return iterate; } END /* for ... iterate */ return -iterate; END #endif #ifdef NEigen /*================================================================ Nonsymmetric EigenProblem The previous example applied a numerical method of dubious virtue to a problem of your choice to illustrate the idiosyncrasies of the method. This example illustrates the converse: an impeccable method applied to nonsymmetric eigenproblems that you specify in terms of their eigenvalues and condition numbers except that the last condition number is computed so that the mathematics of the problem is not torn asunder. By varying the condition numbers, you can control the accuracy of the computed eigenvalues because the inevitable rounding errors provide the perturbations necessary to bring the condition numbers into play. In general, it's the order of magnitude of the condition numbers that you want to vary; different numbers of the same order of magnitude produce little difference. The main program for this application is written as an infinite loop that contains the code to read the eigenvalues and condition numbers, to generate the problem, to print the matrix, to solve the problem with orthes and hqr, and to print the specified and computed eigenvalues. As in most numerical applications, the output will look best using a fixed-size font because the output is columnar. The code that reads the eigenvalues and condition numbers, which are treated as vectors, requires some explanation because it has a number of special features, e.g., comma delimiters, correction of erroneous elements, and no echo. After all, the examples are designed to illustrate use of the API functions. In general, the I/O interface is designed to make it possible to run a modified problem easily. Default values are provided if you respond with nothing; but, the default condition numbers are random so that the problems change. First, if you respond with an EOF, the application terminates. Second, if you enter fewer values than required, all of the the remaining values, eigenvalues and condition numbers, remain unchanged. But, if you enter no condition numbers, they are defaulted to random values. Thus, pressing Enter runs problems with the same eigenvalues but different condition numbers. Third, you can use blanks or commas to delimit the elements, and consecutive commas leave the omitted element unchanged. Thus, ",,10" changes only the third value. But this is unclean since "," and " ," are not the same: the first sets one value, while the second sets one value and leaves the next one explicitly unchanged because the blank delimiters the first value and the comma delimits the second (omitted) value. Fourth, if you enter a syntactically erroneous value, two lines are appended to the output: "Correct and re-enter" and the tail of your original input starting with the erroneous value. You can use the right mouse button to re-enter the last of these lines, correct it, and press Enter. The comment is deleted after you enter your correction. Fifth, your input, with any corrections in their proper place, is echoed at the end of the original prompt. Your actual input is not echoed in the usual sense. The derivation of the algorithm for generating the eigenproblem is beyond the scope of this monograph and couldn't be produced in the form of a text file anyway. But it is quite straightforward to extract the pertinent equations from the code and to verify that the eigenvector matrix computed from the condition numbers has the desired properties. The method used to generate the eigenvector matrix treats condition numbers of 1.0 too well; specifically, it produces a zero row and column that effectively eliminates the eigenvalue from any computational effects. Thus, condition numbers must exceed 1.0. The eigenvector matrix is generated from the condition numbers as a pseudo-orthogonal Householder matrix, i.e., a matrix of the form 1 - sigma * x * y-transpose where sigma is 2 divided by the inner product of x and y. Terming such a matrix pseudo-orthogonal is reasonable because it would be orthogonal if x and y were equal; in this case, x and y are not equal, but their elements differ only in sign. These matrices are their own inverses. The algorithm used to compute x and y is such that the right (columns) and left (rows) eigenvectors provide the specified condition numbers. The computation of the eigenvector matrix also computes the last condition number. Because it must be calculated so that a sum of the condition numbers with appropriate sign changes is zero, this frequently leads to a large value for the last condition number. In effect, it inherits anything that can't be taken into account by manipulating signs. The required matrix is computed by multiplying the diagonal matrix of eigenvalues on both the left and right by the pseudo-orthogonal Householder matrix, which is a similarity transformation because the Householder matrix is its own inverse. The special forms of matrices provides considerable simplification. Finally, the generated eigenproblem is reduced to upper Hessenberg form using orthogonal Householder matrices, and the QR algorithm with shifts is applied. The actual codes are translations to C of the original Algol procedures published in Numerische Mathematik and subsequently collected in Wilkinson and Reinsch. References Smith, B.T. et al. Matrix Eigensystem Routines - EISPACK Guide. Lecture Notes in Computer Science, Vol. 6, 2-nd Edition, Berlin: Springer-Verlag, 1976. Wilkinson, J.H. The Algebraic Eigenvalue Problem. Oxford: Clarendon Press, 1965. Wilkinson, J.H. and Reinsch C. eds. Handbook for Automatic Computation, Volume II, Linear Algebra. Berlin: Springer- Verlag, 1971. */ #include #include #include #include #include #include /* PutMatrix prints a matrix one row per line, which is reasonable in a windowed environment because the user can move the window to see whatever portion of the matrix will fit in the window. tauw is a pseudorandom number generator; it produces values in (0,1) and the condition numbers are obtained by using their inverses. Both orthes and hqr are C translations of Algol procedures in the literature. */ int PutMatrix (int n, int adim, double *a, char *title); double tauw (long *seed); int orthes (int n, int adim, double *a, double *d); int hqr (int n, int adim, double *a, double *wr, double *wi); /* The BSP..ESP and BFP..EFP macros illustrate a noncommittal method for handling string and formatted output. Since this application uses only one window and one buffer, BSP and BFP take no arguments; either or both the identifier and buffer could be arguments. In many cases, BSP..ESP is used even though the length of the string is known or knowable and, in this context, is a bit inefficient. */ #define BSP cputs(Std, #define ESP ,0); #define BFP cputs(Std,Buffer,sprintf(Buffer, #define EFP )); /* The values of the defined constants N, fpf, and BufSize must be correlated. The first is the maximum allowed matrix size, fpf is the format used to print doubles, and the last is the buffer size. Since a complete row is formatted in the buffer, BufSize has to be large enough relative to N and fpf (see PutMatrix). */ #define N 32 #define fpf "%24.16le" #define BufSize 1024 void *Std; long Seed = 314159L; char Buffer[BufSize+1]; double a[N*N], d[N+1], u[N+1], v[N+1]; /* For historical reasons, numerical algorithms have been formulated based on the assumption that matrices are stored by columns, but matrices are stored by rows in C. Rather than rephrase, a macro is used to reference the matrix. Normally, I would use pointers, which would avoid the issue entirely, but there is no basis for expending any effort to optimize the orthes and hqr translations for this application. */ #define a(i,j) *(a + (j-1)*adim + (i-1)) /*---------------------------------------------------------------- The main program for this application is written as an infinite loop that contains the code to read and generate the problem, to print the matrix, to solve the problem, and to print the results. It uses only one window, which is opened with a relatively small output stream that is managed by discarding the oldest output when space is needed. */ int appl (void) BEGIN int adim, n, i, j, k; char *s, *ep; double e[N+1], cn[N+1], c[N+1], x[N+1], y[N+1], alpha, sa, t; cappl("NEigen"); Std = copen("User:Eigenproblem Application;",15000,-1); BSP "This application reads the eigenvalues and condition numbers\n" "for a nonsymmetric eigenproblem, generates a matrix with the\n" "specified eigenvalues and condition numbers > 1, solves the\n" "problem by reducing it to upper Hessenberg form and applying\n" "the QR algorithm, and prints the computed eigenvalues.\n\n" ESP adim = 0; n = 5; while (1) BEGIN /* The user input consists of the size of the matrix and its eigen- values and condition numbers for all but the last eigenvalue, which is computed appropriately when constructing the matrix of eigenvectors. If you don't specify anything or don't specify an integer, the problem size remains unchanged; it's initialized to 5. If you enter a value that is out of bounds, it's set to 5. In any case, the value is explicitly echoed because the ECHO modifier was not specified in the copen and most user's expect an echo. If the size of the problem changes, both the eigenvalues and the condition numbers are assigned default values: 1,2,... for the eigenvalues and random numbers for the condition numbers. */ BSP "\nEnter the size of the eigenproblem: " ESP if (cgets(Std,Buffer,BufSize) < 0) break; i = sscanf(Buffer,"%i",&n); if (n < 2 || n > N) n = 5; BFP "%i\n",n EFP if (adim != n) /* New problem size? */ for (adim = n, i = 1; i <= N; ++i) { e[i] = i; cn[i] = 1.0/tauw(&Seed); } /* The introduction to this example describes how the eigenvalues and condition numbers are read so that it remains only to cover some mechanics of the process. After removing leading blanks, calling the ANSI C function strtod converts the longest initial substring of the input string that looks like a floating-point number to a double and sets ep to the character following the selected substring. Thus, ep should be a blank, comma, or the newline at the end. If the ep character is acceptable and some characters were scanned (s not ep), the value is saved and s is updated. But, if nothing was scanned (s equals ep) and the ep character is a comma, the corresponding value is left unchanged. Otherwise, the value contains an invalid character (ep). In this case, the user is prompted to correct the input by putting out two lines: the first a message and the second the part of the input that must be corrected. This is convenient for the user because the remainder of the input line can be re-entered using the right mouse button and then corrected. The replacement input is read into the buffer so that the buffer ends up containing the user's corrected input as a single line. If the user enters some input, the two line error comment is deleted with "\r\r". Clearly, the explanation is longer than the code to implement it. */ BFP "Enter the %i eigenvalues:\n",n EFP if (cgets(Std,s = Buffer,BufSize) < 1) return 1; for (i = 1; i <= n; ++i) BEGIN while (*s == ' ') ++s; if (*s == '\n') break; t = strtod(s,&ep); if (*ep == ' ' || *ep == ',' || *ep == '\n') { if (s != ep) { e[i] = t; s = ep; } if (*s == ',') ++s; } else { --i; BSP "Correct and re-enter:\n" ESP BSP s ESP if (cgets(Std,s,BufSize) < 1) return 1; BSP "\r\r" ESP } END BSP "\v\b " ESP BSP Buffer ESP /* The code here is the same as that used for the eigenvalues with a couple extra elements. A condition number cannot be less than 1. But condition numbers equal to 1.0 impact the algorithm used to generate the eigenproblem, i.e., the corresponding row and column of the matrix will be zero aside from the diagonal element. Thus, values less than or equal to 1.0 are viewed as erroneous; they are treated the same as a syntax error. Second, if you respond with no input or just some blanks, the condition numbers are set to random values even though this is not likely to produce really interesting condition numbers. */ BFP "Enter the %i condition numbers:\n",n EFP if (cgets(Std,s = Buffer,BufSize) < 1) return 1; for (i = 1; i <= n; ++i) BEGIN while (*s == ' ') ++s; if (*s == '\n') break; t = strtod(s,&ep); if (! t) t = cn[i]; if (t > 1.0 && (*ep == ' ' || *ep == ',' || *ep == '\n')) { if (s != ep) { cn[i] = t; s = ep; } if (*s == ',') ++s; } else { --i; BSP "Correct and re-enter:\n" ESP BSP s ESP if (cgets(Std,s,BufSize) < 1) return 1; BSP "\r\r" ESP } END BSP "\v\b " ESP BSP Buffer ESP if (i == 1) while (i < n) cn[i++] = 1.0/tauw(&Seed); /* The computation of the pseudo-orthogonal Householder matrix is straightforward. In the initial for loop, y[i] is computed directly from the i-th condition number, and x[i] is set to +1 or -1 based on sign of the cumulative sum of the products. The last condition number is used to force the cumulative sum to zero. The final elements of x and y are computed in the second loop, and it's noteworthy that they differ only in sign. They may be zero even though it would be uninteresting. */ x[1] = 1.0; for (alpha = sa = 0.0, i = 1; i < n; ++i) { y[i] = 0.25*(cn[i] - 1.0); sa += y[i]; alpha += x[i]*y[i]; x[i+1] = (alpha > 0.0) ? -1.0 : 1.0; } y[n] = fabs(alpha); sa += y[n]; cn[n] = 1.0 + 4.0*y[n]; alpha = sqrt(1.0/(1.0 + sa)); for (i = 1; i <= n; ++i) { y[i] = sqrt(y[i]*(1.0 + x[i]*alpha)/sa); x[i] *= y[i]; } x[0] = -2.0/alpha; BFP "\nPseudo-Orthogonal Transformation: %25.18le\n",x[0] EFP for (i = 1; i <= n; ++i) BFP "%3i: " fpf fpf "\n",i,x[i],y[i] EFP /* The eigenproblem is simply the product of the Householder matrix, the diagonal matrix of eigenvalues, and the Householder matrix, which works because the Householder matrix is its own inverse. The equations used here are easily verified. */ for (t = 0.0,k = 1; k <= n; ++k) t += x[k]*y[k]*e[k]; t *= x[0]; for (j = 1; j <= n; ++j) { for (i = 1; i <= n; ++i) a(i,j) = x[0]*x[i]*y[j]*(t + e[i] + e[j]); a(j,j) += e[j]; } PutMatrix(n,adim,a,"Generated Eigenproblem"); /* Unless you really get carried away with the condition numbers, hqr should not fail; it's robust, quality software. But, when matching the computed eigenvalues to the specified values, problems can arise. For example, although the specified values are real, the eigenproblem may be so ill conditioned that some complex values are produced; in any case, the imaginary part is thrown away when matching. Condition numbers in the billions often produce complex eigenvalues. Close eigenvalues may not be matched correctly. The matching, however, affects only the presentation of the results. */ orthes(n,adim,a,c); if (hqr(n,adim,a,u,v)) { BSP "hqr failed!\n" ESP continue; } for (i = 1; i <= n; ++i) if (v[i]) { BSP "Oops, hqr found complex eigenvalues.\n" ESP BFP " " fpf fpf "\n",u[i],v[i] EFP } else v[i] = 1.0; for (i = 1; i <= n; ++i) { for (k = j = 1; j <= n; ++j) if (v[j] && fabs(e[i]-u[j]) < fabs(e[i]-u[k])) k = j; c[i] = u[k]; v[k] = 0.0; } BSP "\nSpecified Eigenvalue, " "Computed Eigenvalue, " "and Condition Number\n" ESP for (i = 1; i <= n; ++i) BFP fpf fpf fpf "\n",e[i],c[i],cn[i] EFP END /* The user will tire eventually.*/ return 0; END /*---------------------------------------------------------------- PutMatrix prints an n x n matrix with the specified title. This is very straightforward because there is no need to think about line lengths when dealing with a window; each row is simply generated as one long string of characters. From the user's perspective, any part of the matrix can be viewed by dragging the text around the window, using the left and right arrow keys, or entering an appropriate COLUMN command. As elsewhere in this example, the global string fpf is used to format floating-point quantities. Buffer is large enough to hold a complete row of the matrix. If fpf or N is changed, then compensating changes should be made in Buffer. */ int PutMatrix (int n, int adim, double *a, char *title) BEGIN int i, j; char *s; BFP "\n%ix%i %s\n",n,n,title EFP for (i = 1; i <= n; ++i) { sprintf(s = Buffer,"%2i: ",i); for (j = 1; j <= n; ++j) s += sprintf(s,fpf,a(i,j)); *s++ = '\n'; *s = '\0'; BSP Buffer ESP } return 0; END /*------------------------------------------------------------------- tauw is a full-cycle pseudorandom number generator that is used to generate default values for the condition numbers. See the example named Random for details. The number produced by tauw is in (0,1), and the condition numbers, which must be greater than 1, are obtained by inversion. */ double tauw (long *k) BEGIN *k = *k ^ (*k >> 13); *k = 2147483647L & (*k ^ (*k << (31-13))); return ((double)*k)/2147483648.0; END /*---------------------------------------------------------------- This is a translation to C of the Algol function of the same name that was published in "Similarity Reduction of a General Matrix to Hessenberg Form" by R.S. Martin and J.H. Wilkinson, Numer. Math. 12, 349-368 (1968) and in Wilkinson and Reinsch[1971]. The principal changes are the omission of the submatrix arguments k and l and inclusion of the scaling used in the EISPACK version. The submatrix is used only in connection with balanc, which is not used here so that k = 1 and l = n. In the EISPACK version, these are the igh and low arguments. All other changes are cosmetic or appropriate in C. */ int orthes (int n, int adim, double *a, double *d) BEGIN int i, j, m; double f, g, h, scale; for (m = 2; m < n; ++m) BEGIN for (scale = 0.0, i = m; i <= n; ++i) scale += fabs(a(i,m-1)); if (! scale || scale == fabs(a(m,m-1))) continue; for (h = 0.0, i = m; i <= n; ++i) { d[i] = a(i,m-1)/scale; h += d[i]*d[i]; } g = sqrt(h); if (d[m] >= 0.0) g = -g; h -= d[m]*g; d[m] -= g; for (j = m; j <= n; ++j) { for (f = 0.0, i = m; i <= n; ++i) f += d[i]*a(i,j); f /= h; for (i = m; i <= n; ++i) a(i,j) -= f*d[i]; } for (i = 1; i <= n; ++i) { for (f = 0.0, j = m; j <= n; ++j) f += d[j]*a(i,j); f /= h; for (j = m; j <= n; ++j) a(i,j) -= f*d[j]; } a(m,m-1) = scale*g; END return 1; END /*---------------------------------------------------------------- This is a translation to C of the Algol function of the same name published in "The QR Algorithm for Real Hessenberg Matrices" by R.S. Martin, G. Peters and J.H. Wilkinson, Numer. Math. 14, 219- 231 (1970) and in Wilkinson and Reinsch[1971]. Aside from the changes appropriate for C, both in expression and structure, the changes are minimal; however, the code for the exceptional shifts is not included. */ int hqr (int n, int adim, double *a, double *wr, double *wi) BEGIN int i, j, k ,l, m, na, its; double p, q, r, s, t, w, x, y, z; double macheps = DBL_EPSILON; t = 0.0; its = 0; while (n >= 1) BEGIN for (l = n; l > 1; --l) if (fabs(a(l,l-1)) <= macheps* (fabs(a(l-1,l-1)) + fabs(a(l,l)))) break; x = a(n,n); if (l == n) { wr[n] = x+t; wi[n] = 0.0; --n; its = 0; continue; } na = n-1; y = a(na,na); w = a(n,na)*a(na,n); if (l == na) { p = 0.5*(y-x); q = p*p+w; y = sqrt(fabs(q)); x += t; if (q > 0) { if (p < 0) y = -y; y += p; wr[na] = x+y; wr[n] = x-w/y; wi[na] = wi[n] = 0.0; } else { wr[na] = wr[n] = x+p; wi[na] = y; wi[n] = -y; } n -= 2; its = 0; continue; } if (its++ >= 30) return 1; for (m = n-2; m >= l; --m) { z = a(m,m); r = x-z; s = y-z; p = (r*s - w)/a(m+1,m) + a(m,m+1); q = a(m+1,m+1)-z-r-s; r = a(m+2,m+1); s = fabs(p)+fabs(q)+fabs(r); p /= s; q /= s; r /= s; if (m == l) break; if (fabs(a(m,m-1))*(fabs(q)+fabs(r)) <= macheps*fabs(p)* (fabs(a(m-1,m-1))+fabs(z)+fabs(a(m+1,m+1)))) break; } for (i = m+2; i <= n; ++i) a(i,i-2) = 0.0; for (i = m+3; i <= n; ++i) a(i,i-3) = 0.0; for (k = m; k <= na; ++k) BEGIN if (k != m) { p = a(k,k-1); q = a(k+1,k-1); r = (k != na) ? a(k+2,k-1) : 0.0; x = fabs(p)+fabs(q)+fabs(r); if (! x) continue; p /= x; q /= x; r /= x; } s = sqrt(p*p+q*q+r*r); if (p < 0) s = -s; if (k != m) a(k,k-1) = -s*x; if (l != m) a(k,k-1) = -a(k,k-1); p += s; x = p/s; y = q/s; z = r/s; q /= p; r /= p; for (j = k; j <= n; ++j) { p = a(k,j)+q*a(k+1,j); if (k != na) { p += r*a(k+2,j); a(k+2,j) -= p*z; } a(k+1,j) -= p*y; a(k,j) -= p*x; } j = (k+3 < n) ? k+3 : n; for (i = l; i <= j; ++i) { p = x*a(i,k)+y*a(i,k+1); if (k != na) { p += z*a(i,k+2); a(i,k+2) -= p*r; } a(i,k+1) -= p*q; a(i,k) -= p; } END /* for k */ END /* while n */ return 0; END #endif