/*
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