Other SETL native and utility packages

One of the most essential features of very high level languages is the ability to interface easily and efficiently with lower level languages. This is important for two reasons. (I) Such interfaces make it possible to reuse the vast treasure of public domain code and libraries available on the Net, (i) particularly critical routines can be brought to high efficiency by coding them in lower level languages.

When SETL2 was revised in 1990, this need was recognized by Kirk Snyder, who added an interim Call-out and Call-back capability to the language. The present system provides a considerably more mature and powerful version of this important capability. Using it, the SETL system has been supplied with a collection of native packages providing useful libraries of high-efficiency functions in various areas. Besides the extensive Tk package described in the preceding chapter, these now include operations for string manipulation and matching, analysis of strings using regular expressions, image processing, numerical computation, computation with modular integers, and others. Still other native packages make important internal SETL system functions like parsing, compilation, and library manipulation available for external use within the SETL language. This includes the interface to Unix decribed below.

An interface to the Python language is also provided, making it possible to access all of Python's facilities, including its Web access operations.

The high-efficiency native packages are supplemented by a family of utility packages written in SETL itself. These provide a variety of string, graph, I/O, and integer procedures.

In this report we will first describe this collection of libraries in a manner helpful to their would-be user. Subsequently we descend to the 'C' level, to give the layer of detail needed by would-be implementors of additional **native package**s. This will make it necessary to describe parts of the internal structure of the SETL2 system and the way in which **native package**s are implemented.

Every **native package** must have a SETL specification, which
must be compiled into a SETL library to make the **native package** available
to the system.

An example of a SETL specification for a **native package** is:

native packagestring_utility_pak;procedurebreakup(obj,chars); -- break string into tuple at any occurrence of a c in chars. -- works recursively on tuples of strings, etc.procedurejoin(tuple,chars); -- joins tuple of strings into string by inserting 'chars' betweenproceduresingle_out(obj,chars); -- break string into tuple by separating occurrence of a c in chars -- into singletons. works recursively on tuples of strings, etc.proceduresegregate(obj,chars); -- break string into tuple by separating runs of c's in chars -- from nonruns. works recursively on tuples of strings, etc.procedurekeep_chars(string,chars); -- suppress chars of string that are not in 'chars'proceduresuppress_chars(string,chars); -- suppress chars of string that are 'chars'procedurecase_change(string,mode); -- capitalize if mode is "lu"; decapitalize if mode is "ul"; -- and more ...endstring_utility_pak;

Although SETL performs no type checking at compilation time (and no type declarations appear in a native package specification), the specification does document the functions provided by the compiled package.

As the above example illustrates, a native package specification is identical
to the specification part of a standard SETL package, with the exception
that the keyword **native** appears in the package header. Any variable name can be used as a placeholder for a parameter
in the **native package** Specification, but it is important the number of
arguments match the number expected by the package implementation.

Once a **native package** is included by some SETL code that **use**s it, the functions provided by the will
**native package** will become available to the user exactly as if they were built-in functions.

For example, a user of the *breakup* operation supplied by the **native package** whose specifier is seen above, would write:

program test_breakup; use string_utility_pak; print(breakup({"a,b,c,d",["d,e","f"]},",")); end;to use it. When the above program is compiled, the system will verify statically that the breakup function is called correctly, i.e. with two parameters, but of course will not check the types of the arguments statically, since SETL treats these dynamically.

When program execution begins, the SETL runtime system will realize that the program seen above
requires a **native package**, in this case 'string_utility_pak', and will
search the appropriate directories for the presence of the DLL containing the
code. If no such DLL can be found, the program will abort with an error message announcing that
certain Native functions cannot be resolved.
Otherwise, the DLL is loaded, and the system checks that all the functions
declared in the **native package** specification are defined in the DLL.

string_utility_pak | Provides various commonly used string manipulation primitives, in high-performance form. These include string split and join operations, character suppression operations, and upper/lowercasing. A library of high-performance operations on strings represented as flat byte-arrays is also provided. |

file_pak | Provides various file and directory manipulation operations, including creation, deletion, listing, and navigation. |

rx_pak | High-performance regular-expressions package derived from the Gnu library. |

stringm_pak | High-performance exact and approximate string matching operations, including Knuth-Morris-Pratt and Boyer-Moore matching algorithms, Karp-Rabin probabilistic match, Aho-Corasik match of list of patterns to a string, suffix-tree match, edit-distance based approximate string matching and alignment, and detection of longest common substring of a pair of strings. |

tk_pak | Interface to the Tk widget set. This is the basis for SETL's interactive Graphical Interface, described in Chapter 10. |

parser_pak | Interface to SETL's parse and compile functions. 'Parse' analyses a SETL source file and returns its parse tree as a nested collection of tuples, making the standard SETL error messages available if there are parse errors. A 'compile' function, which takes a string as representing a SETL source file as input and compiles it into the currently active library is also provided. |

ide_pak | This package gives access to the string contents of the currently active SETL edit window and makes it possible to create the interactive source-code analysis and manipulation utilities described below. Its functions also allow the currently selected area of an edit window to be detected and manipulated. A second group of functions in this package allow the library list currently in effect to be read and modified, packages to be reloaded after compilation, and programs to be rerun. |

setldb_pak | Extensive collection of SETL database operations, described in Chapter 12. |

libnr_pak | Large library of high-performance numerical routines, derived from the collection described in the well-known treatise 'Numerical Recipes in C.' |

libntl_pak | Package of integer, modular arithmetic, and modular polynomials package derived for Victor Shoup's well-known 'NTL' library. |

unix_pak | Makes the facilities of Unix available to SETL programs on systems allowing this. |

python_pak | Interface to the Python language and all of its many libraries. In particular,this gives access to Python's internet routines, including its URL and Socket routines. |

grlib_pak | Extensive library of image-manipulation functions. Multiplane images are handled in both byte and floating formats. |

zip_pak | A Lempel-Ziv based package of compression and Zip-File access utilities. |

sound_pak | Procedures for manipulation and analysis of sound in AIFF and related formats. |

midi_pak | Procedures for manipulation of music in MIDI format. |

movie_pak | Procedures for manipulation and display of video image sequences. |

native packagestring_utility_pak; -- ************ general string utilities ************procedurebreakup(obj,chars); -- break string into tuple at any occurrence of a c in chars. -- works recursively on tuples of strings, etc.procedurejoin(tuple,chars); -- joins tuple of strings into string by inserting 'chars' betweenproceduresingle_out(obj,chars); -- break string into tuple by separating occurrence of a c in chars -- into singletons. works recursively on tuples of strings, etc.proceduresegregate(obj,chars); -- break string into tuple by separating runs of c's in chars -- from nonruns. works recursively on tuples of strings, etc.procedurekeep_chars(string,chars); -- suppress chars of string that are not in 'chars'proceduresuppress_chars(string,chars); -- suppress chars of string that are 'chars'procedurecase_change(string,mode); -- capitalize if mode is "lu"; decapitalize if mode is "ul";procedurehexify(string); -- converts string to its hex representation -- ************ operations on flat strings ************procedureflat_create(length); -- creates a flat-string object, as a flat C-type buffer of specified lengthprocedureflat_len(obj); -- returns the length (in characters) of a flat-string objectprocedureflat_to_setl(obj); -- converts a flat-string object to the SETL string with the same charactersprocedureflat_from_setl(string); -- converts a SETL string to the flat-string object with the same charactersprocedureflat_clone(obj); -- returns copy of a flat-string object (as flat-string object)procedureflat_add(obj1,obj2); -- returns concatenation of two flat-string objects (as flat-string object)procedureflat_reverse(obj); -- returns reverse of a flat-string object (as flat-string object)procedureflat_slice(obj,i,j); -- returns indicated slice of a flat-string object -- (as flat-string object)procedureflat_slice_end(obj,i); -- returns indicated tail slice of a flat-string object -- (as flat-string object)procedureflat_set_slice(obj,i,obj2); -- sets indicated slice of a flat-string object -- (from another flat-string object, in place)procedureflat_translate(obj,map); -- translates characters of flat-string object, returns resultprocedureflat_reverse_translate(obj,map); -- translates and reverses characters of flat-string object, returns resultprocedureflat_set_char(obj,i,c); -- sets indicated character of a flat-string object -- (from SETL char, in place)procedureflat_get_char(obj,i); -- gets indicated character of a flat-string object (as SETL char)procedureflat_translate_all(obj,off,minlen,codestr,rar); --procedureflat_file_get(filename,start_char,end_char); -- reads flat string from specified section of a fileprocedureflat_file_put(filename,start_char,obj); -- write flat string to specified location in a fileprocedureflat_match_scores(obj1,obj2,off,repeats); --procedureflat_toto_prepare(flats); --procedureflat_toto_print(obj); --procedureflat_toto_clear(obj); --procedureflat_toto_match(obj,stg,maxnumber); --endstring_utility_pak;

**General string utilities**. The first parameter of 'breakup' can be a string, tuple or set of strings, tuple of tuples of strings, etc. Its second parameter must be a string of characters. The same is true of 'single_out' and 'segregate'. If the first parameter of 'breakup' is simply a string 'stg', 'breakup' breaks 'stg' into a tuple by cutting it at any occurrence of any c in its second parameter, 'chars'. (It two characters of 'chars' follow in immediate succession, two cuts are made, and an empty string between them results.) The characters belonging to 'chars' are removed.

'single_out' behaves in almost exactly the same way, but retains the characters belonging to 'chars', which appear as single characters in the result tuple. 'single_out' behaves similarly, but makes only one cut for each run of characters belonging to 'chars'; these runs are left in the result tuple, in which they alternate with run of characters not belonging to 'chars'.

The small example

programtest; -- string_utility testusestring_utility_pak;endtest;

illustrates these rules. Its output is

["a", "b", "c", "", "d"] ["a", ",", "b", ".", "c", ",", ".", "d"] ["a", ",", "b", ".", "c", ",.", "d"]

The first parameter of 'breakup', 'single_out', or 'segregate' can also be a tuple or set of strings, tuple of tuples of strings, etc., in which case 'breakup', 'single_out', or 'segregate' is applied recursively to the elements or components of the set or tuple until strings to which it can be applied directly are reached. This is illustrated by the example

programtest; -- string_utility test 2usestring_utility_pak;endtest;

which produces the output

{["A", "B", "C", "", "D"], ["a", "b", "c", "", "d"]} [["a", ",", "b", ".", "c", ",", ".", "d"], ["A", ",", "B", ".", "C", ",", ".", "D"]] {["A", ",", "B", ".", "C", ",.", "D"], ["a", ",", "b", ".", "c", ",.", "d"]}

The recursive action of these three operations makes it easy to deal with arrays of objects represented by punctuated strings with multiple layers of punctuation marks. For example, if a 3-dimensional tensor of integers is represented in the string form

stg := "1,2,3;11,12,13;21,22,23:101,102,103;111,112,113;" + "121,122,123:201,202,203;211,212,213;221,222,223";

the following lines of code will read it and produce the output seen:

procedureunstr_rec(obj);return ifis_set(obj)then{unstr_rec(x):xinobj}elseif is_tuple(obj)then[unstr_rec(x):xinobj]else unstr(obj)end if;endunstr_rec;

[[[1, 2, 3], [11, 12, 13], [21, 22, 23]], [[101, 102, 103], [111, 112, 113], [121, 122, 123]], [[201, 202, 203], [211, 212, 213], [221, 222, 223]]]

The 'join' operation is a kind of inverse of 'breakup'. join(tuple,chars); takes a tuple of strings as its first argument and joins this into a string by inserting the characters of 'chars' between its components. An example is

join(breakup("a,bb,ccc",","),"@@@")

which produces the output

a@@@bb@@@ccc

Another example: to replace all the runs of (possibly multiple) whitespace characters in a string by single spaces, one can simply use

join([x: x in segregate(stg," \t") | x(1) notin " \t"]," ")

'case_change' is a simple capitalization/decapitalization utility which operates in the two modes "lu" and "ul". It capitalizes its first string argument if its mode parameter (second parameter) is "lu"; decapitalizes if its mode is parameter is "ul",The following line of code

which produces the output

A.B.C a.b.c

shows this.

'keep_chars' and 'suppress_chars' are simple character-filtering utilities. 'keep_chars' suppresses all characters of its first string argument that are not found in its second argument 'chars'; 'keep_chars' suppresses all characters of its first string argument that are found in its second argument 'chars'. This is shown by the example

which produces the output

abc abcABC

forkin[1..50]loopstg := 1000 * " ";forjin[1..#stg]loopstg(j..j) := "ab";end loop;end loop;

performs its 50,000 insertions in strings of length averaging 1200 in 4 seconds (450 MHz Mac processor), but the code

forkin[1..5]loopstg := 10000 * " ";forjin[1..#stg]loopstg(j..j) := "ab";end loop;end loop;

requires 43 seconds to perform the same number of insertions in strings of length averaging 12000. The reason is that the internal SETL string representation is as a linked list of short string sections. For work with longer strings, and provided that the operations involved do not modify the length of the strings involved, a C-like 'flat' string representation, which allows high-efficiency machine-level addressing of individual, will be much more efficient. A library of operations on such objects, which we shall call 'flat strings', is provided by string_utility_pak. The mutually inverse operations

flat_from_setl(stg) and flat_to_setl(obj)convert standard SETL strings to and from 'flat string' objects. This is shown in the following example,

which produces the output

Hello World11

This example also shows the use of the 'flat_len' operation, which returns the length (in characters) of its flat-string argument.

flat_set_char(obj,i,c) and flat_get_char(obj,i)

provide access to individual characters of flat strings, the first of these setting a designated character and the second retrieving a single character. Examples are

obj := obj2 := flat_from_setl("Hello World");

producing the result

H Jello World Jello World

Note that this operation (characteristically for native objects, but uncharacteristically for SETL) has 'pointer' semantics; it changes the object on which it acts 'in place'. The same is true for the 'flat_set_slice' operation, as seen in the example

obj := obj2 := flat_from_setl("Hello World");

This produces the output

ello ello World HXXXX World HXXXX World

The operation flat_slice(obj,i,j) returns the indicated slice of a flat-string object (as a flat-string object). The operation flat_slice_end(obj,i) returns indicated tail slice of a flat-string object, from the indicated character to the end, (as a flat-string object). The operation flat_set_slice(obj,i,obj2) sets the indicated slice of a flat-string object, starting at the indicated position, and for as many characters as ore in obj2 , in place (s the preceding example shows).

flat_clone(obj) creates a copy of the indicated flat string. flat_add(obj1,obj2) returns the concatenation of two flat-string objects (as a flat-string object). These appear in the following small example, whose output is 'Hello World!Hello World!'.

Hello World!Hello World!programtest; -- flat-string cloning and addingusestring_utility_pak; flat_stg := flat_from_SETL("Hello World!");endtest;

flat_translate(obj) and flat_reverse_translate(obj)

The output produced isprogramtest; -- flat_translate and flat_reverse_translateusestring_utility_pak; flat_stg := flat_from_SETL("Hello World"); map_stg := "" +/ [char(j): jin[0..255]]; -- all 256 characters map to themselves map_stg(abs("H") + 1) := "J"; -- change two of the character mappings map_stg(abs("o") + 1) := "A"; map := flat_from_SETL(map_stg); -- convert to flats stringendtest;

JellA WArld dlrAW AlleJ

In the next few paragraphs, we give a brief syntactic and semantic summary of regular expressions and their use. Regular expressions (which are the chocolate truffles of programming: small but amazingly rich) are patterns which can be matched efficiently to strings. Each such pattern is itself represented by a string, which must conform to the simple, very condensed syntax described in the next few paragraphs. This syntax treats a small number of characters, namelynative packagerx; -- derived from the GNU regular expression packageprocedureregcomp(rwrx,pattern,flags); -- sets rx to the compiled form of the regular expression 'pattern' -- and returns 0 or an error code. -- flags is an integer whose bits allow for special options,including -- extended = (treat as extended expression, see GNU documentation) -- case = (ignore case); nosub = (don't return matching parts info) -- line = (match list of substrings delimited by newlines separately) -- for additional options, see the detailed documentation on the -- Free Software Foundation websiteprocedureregexec(rx,s,flags,rwpos); -- runs compiled rx against string, returning 0 or an error code -- if the match succeeds, pos is returned as a tuple [] defining -- the way in which the whole expression -- and its subparts matched -- flags is an integer whose two low order bits allow for -- the following special options: -- NOT_BEGLINE (value = 1): don't regard start of string as line start -- NOT_ENDLINE (value = 2): don't regard end of string as line endprocedureregerror(err,rx); -- returns detailed error string if rx failed and -- returned 'err' as its error codeendrx;

in a special way, as syntactic/semantic markers, designating such things as subpattern repetition, and alternation. (Each of these characters can then be used 'as itself' within other strings by being preceded with the backslash character '\'). As a 'gotcha' for the unwary, a few other important characters must be preceded with a backslash to function as syntactic/semantic markers; if unescaped they represent themselves. These are the characters

Yet a further nuisance arises from the fact that to be written within SETL quotes the backslash character '\' must be doubled to '\\'.

The regular expression matching process works left-to-right through the 'target' string being matched, and tries to find the leftmost-longest substring of the target which matches the pattern. More specifically, the normal matching process begins first at the first character of the target string, but then if no match can be found, tries again starting at the second, third, .. characters of the target, until a match can be found; the match found is then the longest substring, starting at the given position, which matches the pattern. For simple matches,the match process generates and return a pair [i,j], where i is the starting position of the match, and j is the position of the first subsequent character which cannot be matched, i.e. 1 past the position of the last matched character. (When matches involving patterns and sub-patterns are performed, a more complex tuple is returned, as explained below.) if no match can be found, an empty tuple is returned.

The 'flags' parameters appearing in 'regcomp' and in 'regexec' allow regular expression compilation and matching to proceed in a variety of modes, e.g. case-sensitive and case-insensitive matching. Details concerning these flags and their use are given at the end of the present section.

For the rest one must understand what strings are matched by each possible pattern string. The rules determining this are as follows.

- Any non-special string matches itself. For example, the code
**program**test; -- regular expressions example 1**use**rx; regcomp(regx,"Bc",0); -- compiles "Bc" into 'regx' regexec(regx,"abcdaBcd",0,pos); -- matches compiled "Bc" to "abcdaBcd"**print**(pos);**end**test;produces the output

[[6, 8]]

which indicates that the matched section of the target runs from its 6th to its 7th character.

- the character . matches any single character. Thus the code
**program**test; -- regular expressions example 2**use**rx; regcomp(regx,"b..a",0); -- compiles "b..a" into 'regx' regexec(regx,"abcdaBcd",0,pos); -- matches compiled "b..a" to "abcdaBcd"**print**(pos);**end**test;matches the substring 'bcda' of the target and so produces the output

[[2, 6]]

- bracketed constructs like [a-z0-9$&] consisting of a concatenation of character ranges and single characters match any character that is either explicitly listed or belongs to one of the ranges listed. For example,
**program**test; -- regular expressions example 3**use**rx; regcomp(regx,"ab[a-df-z!][a-df-z!]d",0); -- compiles "ab[a-df-z!][a-df-z!]d" into 'regx' regexec(regx,"Aab!cdzBcd",0,pos); -- matches compiled "ab[a-df-z!][a-df-z!]d" to "Aab!cdzBcd"**print**(pos);**end**test;produces the output

[[2, 7]]

since it matches the substring 'ab!cd' of the target.

- [^chars_spec] matches the set of characters complementary to the characters explicitly listed in'chars_spec' or belonging to one of the ranges listed. For example,
**program**test; -- regular expressions example 4**use**rx; regcomp(regx,"[^a-d][a-df-z!]d",0); -- compiles "[^a-d][a-df-z!]d" into 'regx' regexec(regx,"aab!cdzBcd",0,pos); -- matches compiled "[^a-d][a-df-z!]d" to "aab!cdzBcd""**print**(pos);**end**test;produces the output

[[4, 7]]

since it matches the substring '!cd' of the target.

- [[:builtin:]] matches the set of characters designated by one of a list of keywords designating important character sets. The sets provided are designated by the keywords alnum (all alphanumeric characters), alpha (all alphabetic characters), lower (all lowercase alphabetic characters), upper (all uppercase alphabetic characters), digit (all digits), xdigit (all hexadecimal digits), punct (all punctuation marks),
blank (a blank character), space (all whitespace), cntrl (control characters), graph (all printable characters except space). For example,
**program**test; -- regular expressions example 5**use**rx; regcomp(regx,"[[:punct:]]c[[:digit:][:punct:]][[:digit:][:punct:]]",0); regexec(regx,"aab!c9,Bcd",0,pos);**print**(pos);**end**test;since it matches the substring '!c9,' of the target,and so produces [[4, 8]].

- To this point we have listed constructs matching single characters and character runs of prespecified length. Now we begin to list regular expression forms designating repetitions, and so capable of matching character runs of arbitrary length. These are all postfixes. The first of these is '\+'. When postfixed to an arbitrary regular subexpression, this designates an arbitrary positive number of repetitions (at least one) of the subexpression. An example is
**program**test; -- regular expressions example 6**use**rx; regcomp(regx,"[[:digit:][:punct:]]\\+",0); regexec(regx,"aab!9,99!!888Bcd",0,pos);**print**(pos);**end**test;which matches the substring '!9,99!!888' of the target, and so produces [[4, 14]].

- when postfixed to a subitem of a regular expression, the construct \{m,n\}, where m and n are both integers, calls for at least m and at most n repetitions of the preceding item. This is shown by the example
**program**test; -- regular expressions example 7**use**rx; regcomp(regx,"[[:digit:][:punct:]]\\{2,9\\}",0); -- compiles "Bc" into 'regx' regexec(regx,"aab!9,99!!888Bcd",0,pos); -- matches compiled "Bc" to "abcdaBcd"**print**(pos);**end**test;which produces the output [[4, 13]] since it matches the 9 characters '!9,99!!88' of the target.

- The suffix '\?' is exactly equivalent to '{0,1}', e.g.
**program**test; -- regular expressions matching example 8**use**rx;**print**(regcomp(regx,"\\(fe\\)\\?d",0)); regexec(regx,"fed",0,pos);**print**(pos); regexec(regx,"d",0,pos);**print**(pos);**end**test;is seen to match successfully in both cases since the 'optional' group '\(fe\)\?' matches both the 'fe' in'fed' and the empty string in 'd'.

- The suffix '*' is like '\+', but allows 0 repetitions of the preceding item. and so can match a null string.
- The next two regular expression constructs serve to group subexpressions and allow one of several such subexpressions to be matched. The first of these is simply \(subexp\), which groups a subexpression. The more general form \(subexp1|subexp2|subexp3|...\) lists alternative patters separated by the sign '|'. These subpatterns are tried in sequence; the one leading to the longest match is then selected for the overall match.
A first example of this is

**program**test; -- regular expressions example 9a**use**rx; regcomp(regx,"c\\(a\\+\\)",0); -- compiles 'c\(a\+\)' into 'regx' regexec(regx,"caaaabb",0,pos); -- matches compiled 'c\(a\+\)' to "caaaabb"**print**(pos);**end**test;which we shall compare to

**program**test; -- regular expressions example 9b**use**rx; regcomp(regx,"ca\\+",0); -- compiles 'ca\+' into 'regx' regexec(regx,"caaaabb",0,pos); -- matches compiled 'c\(a\+\)' to "caaaabb"**print**(pos);**end**test;The difference between these two examples is that in the second form the regular expression used is simply 'ca\+', which matches 'c' followed by indefinitely many a', while in the first version the regular expression is the more sophisticated 'c\(a\+\)', in which the same repeated a's are matched, but as a group. The second example generates the expected match tuple [[1, 6]], but the first example generates the less expected [[1, 6], [2, 6]]. This is because besides the first, widest pair representing the entire range of characters matched, each subgroup on a regular expression generates an additional pair showing the range matched by the subgroup. In the example shown, this is the run of characters 'aaaa' in the second through fifth positions. However, if the group as a whole is repeated in the regular expression, only the last of its matches generates such a pair. The effect of this can be seen in the example

**program**test; -- regular expressions example 9c**use**rx; regcomp(regx,"c\\(a\\)\\+",0); -- compiles 'c(a)\+' into 'regx' regexec(regx,"caaaabb",0,pos); -- matches compiled 'c(a)\+' to "caaaabb"**print**(pos);**end**test;Here, since we match the repeated group '(a)\+', only its final match generates a pair into the resulting tuple, which is [[1, 6], [5, 6]], indicating that the last match of the repeated group is to the character 'a' appearing 5th in the target string. (Compare this to the [[1, 6], [2, 6]] returned when 'c\(a\+\)' is matched instead.

The example

**program**test; -- regular expressions example 9b**use**rx; regcomp(regx,"c\\(a\\+\\)\\(b\\+\\)",0); -- compiles 'c(a)\+(b)\+' into 'regx' regexec(regx,"caaaabb",0,pos); -- matches compiled 'c(a)\+' to "caaaabb"**print**(pos);**end**test;which returns [[1, 8], [5, 6], [6, 8]], makes another rule apparent. This is the fact that a series of successively matched subexpressions generates a series of successive pairs into the result returned. Each subexpression generates a pair indicating the range of characters that it has matched, or, if it is repeated, the range of characters matched by the last of its repetitions.

The general rule for construction of the tuple returned upon successful match of a regular expression involving (possibly nested) subexpressions and repetitions can now be stated as follows. Successful match generates a series of pairs corresponding to character ranges. The first such pair defines the total range covered by the match. This is followed by a sequence of pairs, grouped into successive subsequences corresponding to the series of top-level subgroups of the regular expression. For each such subgroup, the run of characters matched by the (final repetition of) the subgroup (if repeated) is indicated. If the subgroup contains nested sub-subgroups, the same rule is applied recursively, the pairs generated by any single subgroup always following successively in the overall list of pairs returned.

All this is shown in the example

**program**test; -- regular expressions example 10**use**rx; regcomp(regx,"\\(\\(a\\+b\\+\\)\\+\\(c\\+\\)\\)\\(d\\+\\)",0); regexec(regx,"aaaabbaabbbcccddd",0,pos);**print**(pos);**end**test;which returns the list of pairs

[[1, 18], [1, 15], [7, 12], [12, 15], [15, 18]]

and repays careful study. The value returned reflects the following facts:- the full regular expression '\(\(a\+b\+\)\+\(c\+\)\)\(d\+\)' matches the whole target string "aaaabbaabbbcccddd";
- its regular subexpression '\(\(a\+b\+\)\+\(c\+\)\)' is not repeated, and matches the first 14 characters of the target string;
- within this subexpression, the last repetition of the repeated sub-subexpression \(a\+b\+\) matches characters 7-11, and is followed by the sub-subexpression \(c\+\), which matches characters 12-14;
- the final part '\(d\+\)' of the regular expression matches characters 15-17.
- Within a parenthesised subexpression of a regular expression, any number of alternatives can be separated by the sign '\|', designating 'or'. For example,
**program**test; -- regular expressions matching example 11**use**rx; regcomp(regx,"\\(\\(ab\\)\\+\\|\\(cd\\)\\+\\)\\+",0); regexec(regx,"ababcdcdcd",0,pos);**print**(pos);**end**test;returns [[1, 11], [5, 11], [5, 5], [9, 11]] since the alternative '\(\(ab\)\+\\(cd\)\+\)\+' matches either runs of 'ab' or runs of 'cd',and so the whole target string "ababcdcdcd" is matched.

- As matching proceeds, the matcher keeps track of the subexpressions it encounters and the substrings they match. This makes it possible to use the construct '/n', where n can be any digit from 1 to 9, to refer to 'the substring matched by the n-th subexpression' (but of course this construct can only be used in positions following the n-th subexpression.) This gives a way of making later matches depend on the outcome of previously matched regular subexpressions. A simple example is '\(.\+)\1', which matches any string that is composed of two identical halves, since `\(.*\)' matches the first half, which can be anything, but then '\1' must match exactly the same string. Another example is
**program**test; -- regular expressions matching example 12**use**rx; regcomp(regx,"\\(cd\\+\\)\\+\\(ab\\+\\)\\+z\\1\\(\\2\\)",0); regexec(regx,"cdcddababbzcddabb",0,pos);**print**(pos);**end**test;which returns [[1, 18], [3, 6], [8, 11], [15, 18]] since its first part '\(cd\+\)\+\\(ab\+\)\+' matches 'cdcddababb' and associates the two subexpressions it contains with 'cdd' and 'abb' (the last character ranges these repeated regular subexpressions match); then '\1' and '\(\2\)' match the second 'cdd' and 'abb' respectively. The latter is a subgroup match, and so generates the pair [15, 18] into the result tuple. This example also shows that the construct '/n' can be included in grouped subexpressions, etc/

- The special characters '^' and '$' are used to pin matches to the start and end of a target string respectively. Both match null strings, but '^' matches only at the very start, and '$' only at the very end, of the target string. This is shown in the following example:
**program**test; -- regular expressions example 13**use**rx; regcomp(regx,"car$",0); -- compiles 'car$' into 'regx' regexec(regx,"carcar",0,pos); -- matches compiled 'car$' to "carcar"**print**(pos);**end**test;returns [[4, 7]] since it can only match the second occurrence of 'car' in 'carcar'. With the '$' removed it would instead match the first occurrence, and so returns [[1, 4]]. Similarly

**program**test; -- regular expressions example 14**use**rx; regcomp(regx,"^car",0); regexec(regx,"acar",0,pos);**print**(pos);**end**test;fails to match, since it can only match at the start of the target string. Without the '^' it would return [[2, 5]], since a match starting at the second character would be found.

- A few other special constructs which pin matches to specific position at the start, end, or middle of words are provided. These are `\>' (end of word), `\<' (beginning of word) ,`\b' (beginning or end of word), and `\B' (not beginning or end of word). All of these constructs match null strings in their indicated positions. This is shown in the composite example
**program**test; -- regular expressions example 15**use**rx; regcomp(regx,"cat\\>",0); regexec(regx,"thecatcat catcatcat",0,pos);**print**(pos); regcomp(regx,"\\<cat",0); regexec(regx,"thecatcat catcatcat",0,pos);**print**(pos); regcomp(regx,"\\Bcat\\B",0); regexec(regx,"the catcat catcatcat",0,pos);**print**(pos); regcomp(regx,"\\bcat\\B",0); regexec(regx,"thecatcat catcatcat",0,pos);**print**(pos);**end**test;which returns

[[7, 10]] [[11, 14]] [[15, 18]] [[11, 14]]

reflecting the fact that

- 'cat\>' can only match at the end of a word, hence the second occurence of 'cat';
- '\<cat' can only match at the start of a word, hence the third occurence of 'cat';
- '\Bcat\B' can only match in the middle of a word, hence the fourth occurence of 'cat';
- '\bcat\B' can only match in the start of a word of more than 3 letters, hence the third occurence of 'cat'.

The 'flags' parameters appearing in 'regcomp' and in 'regexec' allow regular expression compilation and matching to proceed in a variaety of modes, e.g. case-sensitive and case-insensitive maching. To set these flags, one passes appropriate integer 'flags' parameter values to 'regcomp' and 'regexec' respectively. Each flag corresponds to a bit position in this integer, which must be 1 for the corresponding flag to be set. The following table lists the 'regcomp' flags.

EXTENDED (value = 1) | Treat the pattern as an extended regular expression, rather than as a basic regular expression. Basic regular expressions are regular expressions having the sytax detailed in the preceding paragrphs. Extended regular expressions differ from these in that instances of the characters
do not need to be preceded by the backslash character '/'. That is, we can write '(', ')', '+', '?', '|', '{', and '}' instead of '\(', '\)', '\+', '\?', '\|', '\{', and '\}' respectively, a welcome relief. In extended regular expressions, we can also use {m} as an abbreviation for {m,m}, and can use {m,} to mean 'at least m occcurences'. Examples are given below. |

ICASE (value = 2) | Ignore case when matching letters. |

NEWLINE (value = 4) | Don't treat the end-of-line character as an ordinary character; instead, treat each line as a target to be matched. This allows '^' to match at the start of any line, and '$' to match at the end of any line, but prevents '.' from matching the end-of-line character. Examples are given below. |

Here are examples of the use of these flags.

programtest; -- regular expressions example 9a.1userx; regcomp(regx,"c(a+)",1); -- compiles 'c(a+)' into 'regx' regexec(regx,"caaaabb",0,pos); -- matches compiled 'c(a+)' to "caaaabb"endtest;/

is the extended regular expression version of Example 9a, and returns the same result. To use and ordinary regular expression we would have to write "c\\(a\\+\\)", as in Example 9a, instead. Similarly

programtest; -- regular expressions matching example 11.1userx; regcomp(regx,"((ab)+|(cd)+)+",1); regexec(regx,"ababcdcdcd",0,pos);endtest;

is the extended regular expression version of Example 9a, and returns the same result. Plainly "((ab)+|(cd)+)+" is much more readable than "\\(\\(ab\\)\\+\\|\\(cd\\)\\+\\)\\+".

The example

programtest; -- regular expressions example 16userx; regcomp(regx,"\\(abc\\)\\+",2); regexec(regx,"ABCabc",0,pos);endtest;

which returns [[1, 7], [4, 7]], or still better its extended regular expression version

programtest; -- regular expressions example 16buserx; regcomp(regx,"(abc)+",3); regexec(regx,"ABCabc",0,pos);endtest;

shows the use of case-insensitive matching.

The example

programtest; -- regular expressions matching example 17userx; regcomp(regx,"e$",4); regexec(regx,"Oh, Jack be nimble\nJack be quick",0,pos);endtest;

which returns

[[18, 19]] [[20, 21]] []

shows the alternative treatment of end-of-line characters, and the failure of the third match in this example shows also that if this option is used matches can never extend over more than one line.

It is often useful to simplify the use of a **native package** by representing the raw functions it provides in some more convenient object syntax. In this section we present such a 'wrapper class', called 'rxp', for the regular expressions **native package** described in the preceding section. This allows us to write the regular expression and matching operation of the preceding section as 'pattern_object * target', where the 'pattern_object' is created from a regular expression string by writing

Here 'pattern' is a regular expression string of the kind described in the preceding section, and 'flags' is an integer composed of bitflags, as described there. If the matching operation pattern_object * target succeeds it returns exactly the kind of position tuple described in the preceding section.

This object class also allows 'pattern_object ** target' to be used to match a regular expression repeatedly to string

Note that the class seen below interprets the error codes returned by 'regcomp' and so explains their meaning.

The class code is as follows.

classrxp; -- SETL regular SETL regular-expression classprocedurecreate(pattern,flags); -- compile pattern into regular expressionendrxp;class bodyrxp; -- SETL regular SETL regular-expression classuserx; var native_rx; -- the native compiled form of the regular expressionprocedurecreate(pattern,flags); -- compile pattern into regular expressionconstREG_NOSUB := 2; flags ?:= REG_NOSUB; -- use no flags as defaultif(errc := regcomp(native_rx,pattern,flags)) = 0then return;end if; message_string :=caseerrc -- otherwise there is an errorwhen1 => "Bad brackets \\{...\\}"when2 => "Bad regular expression syntax"when3 => "Baduseof repetition operator, e.g. *,?,+ \\{...\\}"when4 => "Bad collating element"when5 => "Bad characterclassname"when6 => "Terminal occurrence of \\"when7 => "Bad digitindigit construct"when8 => "Unbalanced square brackets"when9 => "Unbalanced parentheses \\(...\\)"when10 => "Unbalanced braces \\{...\\}"when11 => "Bad endpointinarangeexpression"when12 => "Ran out of memory during compilation"end case;abort("**** REGULAR EXPRESSION COMPILATION ERROR: " + message_string);endcreate;procedureself * stg; -- match regular expression to string, returning [] or match-tupleif(errc := regexec(native_rx,stg,0,pos_tuple)) > 1then-- ran out of memoryabort("**** REGULAR EXPRESSION EXECUTION ERROR: " + "Ran out of memory");end if;return iferrc = 0thenpos_tupleelse[]end if;end;procedureself ** stg; -- match regular expression repeatedly to string, -- returning [] or match-tuple matches := []; -- these will be collected match_offset := 0; -- offset to apply to match tuplewhilestg /= ""loopif(errc := regexec(native_rx,stg,0,pos_tuple)) > 1then-- ran out of memoryabort("**** REGULAR EXPRESSION EXECUTION ERROR: " + "Ran out of memory");end if;iferrc > 0or(nexc := pos_tuple(1)(2)) = 1then returnmatches;end if; -- nothing matched matcheswith:= [[x + match_offset,y + match_offset]: [x,y]inpos_tuple]; -- apply the match offset to the tuple returned match_offset +:= (nexc - 1); -- advance the match offset to the last matched character stg := stg(nexc..); -- continue with the remaining part of the stringend loop;returnmatches;end;endrxp;

The following small test program illustrates the use of this object class.

programtest; -- regular expressions class exampleuserxp;endtest;

The following small test program illustrates the use of this object class. It returns

[[2, 5]] [[[2, 5]], [[5, 8]], [[8, 12]]],

the first pair corresponding to the three characters 'abb' matched by the first, uniterated, match operation, and the list of three pairs corresponding to the threefold match iteration triggered by the second match operation.

The full regular expression for this entire lexical analyzer would be a totally incomprehensible coded string of several hundred characters. To have a chance of writing this correctly (which will certainly involve some amount of debugging) we must find some way of decomposing it into the string-sum of meaningful submodules. A technique generally applicable to large regular expressions is to think of it as a (perhaps recursive) alternation of concatenated subparts which can be written separately, debugged, and then combined. To make an alternation out of its separate alternatives we have then only to write

"("alternativeTo make a concatenation out of its separate pieces we write_{1}+ "|" + alternative_{2}+ "|" + ... + alternative_{k}")"

"pieceThis is done to construct the following lexical analyzer. Note that the analyzer is written as an extended regular expression to reduce clutter. At the top level, it sees the string to be analyzed as a concatenation of whitespace, variable names, quoted strings, based numerals, ordinary numerals, operation words, comments, and end-line characters. If a line ends with a comment, the following end-line character is included in the comment. Otherwise the end-line character which terminates a line counts as whitespace._{1}+ piece_{2}+ ... + piece_{k}

The lexical analyzer appears as the routine 'get_tokens'in the following package, which is supplied with the SETL system.

packageextended_SETL_lexer;constextra_opsigns := "~!@$%^&"procedureget_tokens(stg); -- decompose a string into generalized SETL tokensendextended_SETL_lexer;package bodyextended_SETL_lexer;userx,string_utility_pak; var lexer := OM; -- the 'escaped' characters allowed within SETL strings -- are \\, \0, \n, \r, \f, \t, \",\xddconstD_hex := "\\\\x[0-9a-f][0-9a-f]"; -- regular expression for hex stringsconstD_escaped := "\\\\f|\\\\0|\\\\t|\\\\n|\\\\r|\\\\\\\\|\\\\\\\""; -- regular expression for escaped special charactersconstD_opsigns_and_parens "`~!@#$%^&*()+={}[|':/?><,.[`;"; -- regular expression for generalized opsigns/parensconstD_opsigns := "([~!@#$%^&*+=:(){}[/?><.]|])"; -- regular expression for generalized opsignsconstquote := "\""; -- quote characterconstquote_eol := "[\"\n\r]"; -- quote and endline charactersconstnot_quote_eol_char := "[^\"\n\r]"; -- characters other than quotes and endlinesconstescaped_quote := "\\\\\\\""; -- regular expression for quote mark preceded by backslashconststart_quoted_string := quote + "(" + not_quote_eol_char + "|" + escaped_quote + ")+"; -- regular expression for start of complete or incomplete -- quoted string, without terminating quote or endline characterconstD_quoted := start_quoted_string + quote_eol; -- regular expression for complete and incomplete quoted stringsconstD_name := "[a-zA-Z]([[:alnum:]]|_)*"; -- regular expression for SETL variable namesconstD_white := "( |\t|\r|\n)*"; -- regular expression for whitespaceconstD_num := "[0-9][_0-9]*(.[0-9][_0-9]*((E|e|e-|e-)[0-9]+)?)?"; -- regular expression for real and floating numeralsconstD_based_num := "[0-9]+#[0-9a-z][_0-9a-z]*(.[0-9a-z]" + "[_0-9a-z]*#((E|e|e-|e-)[0-9]+)?)?"; -- regular expression for real and floating numerals -- with non-decimal basesconstD_not_monadic_end := "([~!@#$%^&*+=-?><.:/]*[*+=:?><./])";constD_opword := "([~!@#$%^&*+=(){},;[:/?><.]|]|" + D_not_monadic_end + "|[a-zA-Z][[:alnum:]]*)";constppe := {"procedure","program","class","package","use","var","const","end"}; -- declaration headers and tailsconsticl:= {"if","case","loop"}; -- control structure headers (omitting headers for iterative loops, -- which, as tokens, do not differ from variable names)procedureget_tokens(stg); -- decompose a string into generalized SETL tokensif#stg = 0then return[];end if; -- null string caseifstg(#stg)notin"\n\r"thenstg +:= "\n";end if; -- add a terminating endline if needediflexer = OMthen-- we must compile the regular expression for the lexer regcomp(lexer,"(" + D_white + "|" + D_name + "|" + D_quoted + "|" + D_based_num + "|" + D_num + "|" + D_opword + ")",1); -- use extended regular expression syntaxend if; tokens := [ ]; -- will collectwhilestg /= ""loopif#stg > 1andstg(1..2) = "--"then-- we have a comment comment :=break(stg,"\n\r"); -- this runs to the line end comment := comment +span(stg,"\n\r \t"); -- and includes the line-end character tokenswith:= comment; -- collect the comment string as a tokenelseifregexec(lexer,stg,1,pos) = 0then-- the regular expression finds a token [strt,endd] := pos(1); -- get token start and end tokenswith:= stg(strt..(endd - 1)min#stg); -- collect the tokenifendd = 1then endd := 2;;end if; -- if the token is a null string then take 1 char stg := stg(enddmin(#stg + 1)..); -- continue with remainder of input stringelse-- if no token can be found,then exitexit;end if;end loop;returntokens; -- return the list of tokensendget_tokens;endextended_SETL_lexer;

The package specifier is

native packagestringm; -- string matching package ------------------------------------ -- Exact String Matching ------------------------------------procedurekmp(string,pattern); -- Knuth-Morris-Pratt match of a pattern to a string -- returns tuple of all starting points of matchesprocedurekmp_compile(pattern); -- optimized Knuth-Morris-Pratt match; returns compiled patternprocedurekmp_exec(string,pat); -- optimized Knuth-Morris-Pratt match using pre-compiled pattern -- returns tuple of all starting points of matchesprocedurebm(string,pattern); -- Boyer-Moore match of a pattern to a string -- returns tuple of all starting points of matchesprocedurebm_compile(pattern); -- optimized Boyer-Moore match; returns compiled patternprocedurebm_exec(string,pat); -- optimized Boyer-Moore match using pre-compiled pattern -- returns tuple of all starting points of matchesprocedurekr(string,pattern); -- Karp-Rabin probabilistic match of a pattern to a string -- returns tuple of all starting points of matches ------------------------------------ -- Exact Matching to a Set of Strings ------------------------------------procedureac_compile(tuple); -- Aho-Corasik match of list of patterns to a string; compiles patternprocedureac_init(pattern,text); -- Aho-Corasik match of list of patterns to a string; binds pattern to textprocedureac_next_match(pattern); -- Aho-Corasik match of list of patterns to a string; -- returns next match -- (as a triple [tuple_index_of_pattern_matched, lic,len]), or OM ------------------------------------ -- Exact Matching using Suffix Trees ------------------------------------procedurest_compile(tuple); -- suffix-tree match of list of patterns to a string; compiles patternprocedurest_match(pattern,text); -- suffix-tree match of list of patterns to a string -- returns tuple of all match points ------------------------------------ -- Edit Distance ------------------------------------procedureedist(s1,s2); -- edit distance between two stringsprocedureexedist(s1,s2,w); -- weighted edit distance between two strings -- the weights w should be supplied as an integer quadruple: -- [insert_weight,delete_weight,substitute_weight,match_weight]procedureetrans(s1,s2); -- best edit sequence between two strings -- returns a descriptor consisting of characters M (matching character) -- S (replace character) D (delete), I (insert), followed by integer distanceprocedureexetrans(s1,s2,w); -- best edit sequence between two strings, using weights -- returns a descriptor consisting of characters M (matching character) -- S (replace character) D (delete), I (insert), followed by integer distance -- the weights w should be supplied as an integer quadruple: -- [insert_weight,delete_weight,substitute_weight,match_weight] ------------------------------------ -- Longest matching substring ------------------------------------procedurelcseq(s1,s2); -- finds longest matching subsequence of two strings; returns a pair -- [common_string,length]. Note that this is an increasing sequence, not -- necessarily contiguous ------------------------------------ -- Approximate matching operations ------------------------------------procedurepwscores(a); -- compile a list of character-score triples into the form -- expected by simil and simltproceduresimil(s1,s2,pws); -- finds the 'optimal alignment distance' between two stringsproceduresimilt(s1,s2,pws); -- finds the 'optimal alignment' of two strings, representing this -- as a descriptor of the same kind that 'exetrans' returnsendstringm;

**Exact String Matching**. Three exact matching procedures for matching to single strings, namely the Knuth-Morris-Pratt, the roughly equivalent Boyer-Moore procedure, and the potentially more efficient Karp-Rabin probabilistic algorithm constitute the first group of operations provided. The first two of these are made available in two forms, one of which matches a pattern string directly to a target, while the second begins by compiling the pattern into a set of tables and then matches this compiled form of the pattern to the target string. This second form is intended for use in situations in which a single pattern is to be matched repeatedly to the members of a collection of target strings.

A first example of the use of these routines is

programtest; -- string matching example 1usestringm;endtest;

which returns

[1, 3, 6] [1, 3, 6] [1, 3, 6]

in each case the starting characters of the three match positions. Note that the position '3' is found by all these procedures, even though it lies within the zone covered by the match to its left.

The 'precompiled' form of this same example is

programtest; -- string matching example 1usestringm;endtest;

which returns the same result.

**Exact Matching to a Set of Strings**. Two methods, the Aho-Corasik method and a method using suffix Trees, are provided to match a list of patterns to a target string. As seen in the following example, the Aho-Corasik scheme is set up as a triple of procedures. The first of these compiles a list of strings into a form suitable for efficient matching, the second binds the resulting pattern object to the target string to be matched, and the third acts something like a SET iterator, retuning a new match each time it is called, until a final call returns OM, indicating that no more matches exist.

programtest; -- exact string matching to a list of patternsusestringm; compiled_pat := ac_compile(["abc","aa","cde"]); -- compile list of patterns ac_init(compiled_pat,"abaabcde"); -- bind to target stringendtest;

Since each successful call to 'ac_next_match' returns a triple of the form [index_of_matching_pattern,match_start,match_length] this example generates the output

[2, 3, 2] [1, 4, 3] [3, 6, 3] <om>

Note that, even in the presence of overlapping matches, the algorithm reports every point at which one of the strings in the pattern list can be matched to the target string.

**Approximate Matching**.

This group of functions can be used to compute edit distances, longest common subsequences and string similarity.

**Edit Distance**.
The *edit distance* between two strings is defined as the minimum number
of edit operations needed to transform the first string into the
second.

The edit operations considered In our implementation are insertions, deletions and substitutions. Two versions 'edist' and 'exedist' of the edit distance primitive are provided. The first of these computes edit distance using standard weights assigned to the edit operation. The second, extended version allows explicit assignment of a weight to each edit operation. These weights must be passed as a tuple of four integer weights having the form:

The default weights used in the standard version are simply [1,1,1,0].

The two routines 'etrans' and 'exetrans' work with the same edit distance as 'edist' and 'exedist', but
also generate an *edit transcript*. This is a sequence of one-character edit operation codes in which D, I, S, and M designate
Deletion, Insertion, Replacement, and Match respectively. The following example illustrates the use of 'etrans' and the way in which the edit transcript that it returns can be used to set up an alignment between two similar but not identical strings.

programtest; -- approximate string matchingusestringm;forcintranscriptloopcasecwhen"M" => revised_s1 +:= s1(loc_in_s1 +:= 1); revised_s2 +:= s2(loc_in_s2 +:= 1); indicator_stg +:= ".";when"S" => revised_s1 +:= s1(loc_in_s1 +:= 1); revised_s2 +:= s2(loc_in_s2 +:= 1); indicator_stg +:= "|";when"D" => revised_s1 +:= s1(loc_in_s1 +:= 1); revised_s2 +:= "-"; indicator_stg +:= "|";when"I" => revised_s1 +:= "-"; revised_s2 +:= s2(loc_in_s2 +:= 1); indicator_stg +:= "|";end case;end loop;endtest;

This produces the output

["MMIMMMSMD", 3] aa-ababab ..|...|.| aababaaa-

The first line is the standard pair [edit_transcript,edit_dist] returned by 'etrans'. The three remaining lines are (i) the first string, with insertions marked by dashes; (ii) an indicator string, showing where the aligned first and second strings differ; (iii) the second string, with insertions marked by dashes.

By increasing the cost of insertions and deletions, as in the following variant, we can of course force a different alignment.

programtest; -- approximate string matchingusestringm;forcintranscriptloopcasecwhen"M" => revised_s1 +:= s1(loc_in_s1 +:= 1); revised_s2 +:= s2(loc_in_s2 +:= 1); indicator_stg +:= ".";when"S" => revised_s1 +:= s1(loc_in_s1 +:= 1); revised_s2 +:= s2(loc_in_s2 +:= 1); indicator_stg +:= "|";when"D" => revised_s1 +:= s1(loc_in_s1 +:= 1); revised_s2 +:= "-"; indicator_stg +:= "|";when"I" => revised_s1 +:= "-"; revised_s2 +:= s2(loc_in_s2 +:= 1); indicator_stg +:= "|";end case;end loop;endtest;

This generates the alignment

["MMDMMMMSS", 5] aaabababb ..|....|| aa-babaaa

This approach is that embodied in the group of three procedures 'pwscores', 'simil', and 'similt'. pwscores(score_mat) accepts a scoring matrix as parameter and returns a compiled form of this matrix as a opaque 'score matrix object'. The scoring matrix should be supplied as a list of triples [char_1,char_2,similarity_value], where scores for the nominal character '0' indicating a blank are also required. As seen in the next example, the procedures 'simil', and 'similt' then use this 'score matrix object' to find the optimal alignment of two strings containing the characters for which similarity values have been assigned.

programtest; -- approximate string matching using a similarity matrixusestringm; pws_obj := pwscores([["a","b",1],["a","a",2], ["b","b",2],[0,"a",0],[0,"b",0]]); -- compile character similarity matrixforcintranscriptloopcasecwhen"A" => revised_s1 +:= (c1 := s1(loc_in_s1 +:= 1)); revised_s2 +:= (c2 := s2(loc_in_s2 +:= 1)); indicator_stg +:=ifc1 = c2then"."else"|"end if;when"D" => revised_s1 +:= s1(loc_in_s1 +:= 1); revised_s2 +:= "-"; indicator_stg +:= "|";when"I" => revised_s1 +:= "-"; revised_s2 +:= s2(loc_in_s2 +:= 1); indicator_stg +:= "|";end case;end loop;endtest;

The code seen here generates the alignment

["AADAAAAAA", 14] aaabababb ..|....|| aa-babaaa

Note that 'simil' returns only an optimal similarity score; 'similt', like 'exetrans', returns a pair [transcript,score].

The 'similt' alignment technique is commonly used in computational genomics to align sequences of characters representing the amino acids out of which proteins are built. The standard characters for the 23 amino acids occuring in ordinary proteins are "ABCDEFGHIKLMNPQRSTVWXYZ". A commonly used protein-alignment scoring matrix, the so-called 'Blosum-62' matrix derived by collecting statistics from large numbers of related proteins, is (in the same order of amino acids, and with very low scores understood for pairs containing blanks)

4, -2, 4, 0, -3, 9, -2, 4, -3, 6, -1, 1, -4, 2, 5, -2, -3, -2, -3, -3, 6, 0, -1, -3, -1, -2, -3, 6, -2, 0, -3, -1, 0, -1, -2, 8, -1, -3, -1, -3, -3, 0, -4, -3, 4, -1, 0, -3, -1, 1, -3, -2, -1, -3, 5, -1, -4, -1, -4, -3, 0, -4, -3, 2, -2, 4, -1, -3, -1, -3, -2, 0, -3, -2, 1, -1, 2, 5, -2, 3, -3, 1, 0, -3, 0, 1, -3, 0, -3, -2, 6, -1, -2, -3, -1, -1, -4, -2, -2, -3, -1, -3, -2, -2, 7, -1, 0, -3, 0, 2, -3, -2, 0, -3, 1, -2, 0, 0, -1, 5, -1, -1, -3, -2, 0, -3, -2, 0, -3, 2, -2, -1, 0, -2, 1, 5, 1, 0, -1, 0, 0, -2, 0, -1, -2, 0, -2, -1, 1, -1, 0, -1, 4, 0, -1, -1, -1, -1, -2, -2, -2, -1, -1, -1, -1, 0, -1, -1, -1, 1, 5, 0, -3, -1, -3, -2, -1, -3, -3, 3, -2, 1, 1, -3, -2, -2, -3, -2, 0, 4, -3, -4, -2, -4, -3, 1, -2, -2, -3, -3, -2, -1, -4, -4, -2, -3, -3, -2, -3, 11, 0, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, -1, -1, 0, 0, -1, -2, -1, -2, -3, -2, -3, -2, 3, -3, 2, -1, -2, -1, -1, -2, -3, -1, -2, -2, -2, -1, 2, -1, 7, -1, 1, -3, 1, 4, -3, -2, 0, -3, 1, -3, -1, 0, -1, 3, 0, 0, -1, -2, -3, -1, -2, 4

**Longest common subsequence**. The final native procedure supplied by STRINGM_PAK is lcseq(s1,s2), which finds longest common subsequence of two strings, returning a pair [common_string,length]. Note that the sequence found is an ordered sequence of characters, not necessarily contiguous in either of the two strings being compared. That is, as the following example shows, the procedure finds a longest common subsequence, not the longest common run.

programtest; -- longest common subsequenceusestringm;time forall good men","What hath God wrought?"));endtest;

The pair ["ht a od ", 8] is returned.

**String alignment using 'similt'**.

programtest; -- string alignment using 'similt'usestringm; scores_list := [["a","b",1],["a","c",1],["a","c",1],["a",0,-2], ["b","a",1],["b","c",1],["b","c",1],["b",0,-2], ["c","a",1],["c","b",1],["c","c",1],["c",0,-2], ["d","a",1],["d","b",1],["d","c",1],["d",0,-2]]; scores_obj := pwscores(scores_list);endtest;

PARSER_PAK provides a small but important collection of functions which give access to SETL's internal parse operations. Its specifier is

The 'compile' function can be used to compile any SETL string dynamically that could otherwise be compiled as a SETL program, package, or class.If compilation is successful, the byte-code produced by the compilation is written to the current library in the standard manner, as explained in Section XXX.native packageparser; -- package of SETL parse functionsprocedureparse(str); -- parse a string, returning its parse tree as a nested set of tuples -- this handles all SETL 'compilation units'procedureparse_expr(str); -- parse a string, returning its parse tree as a nested set of tuples -- this handles any string that is a valid 'program body'procedurecompile(str); -- compile a string, writing result to current libraryproceduresetl_num_errors(); -- read number of errors generated during last previous -- 'parse' or 'compile' operationproceduresetl_err_string(err_num); -- get text of indicated error message generated during -- last previous 'parse' or 'compile' operationendparser;

**The 'parse' and 'parse_expr' functions**. The 'parse' and 'parse_expr' functions give more intimate access to the internal structure of the parse trees generated by the SETL system as a compilation intermediate, and so can be used to support variant forms of SETL and other quite different systems which can make use of the SETL syntax. When a compilation succeeds, the 'parse' function returns the syntax tree of the code it has compiled. This is returned as a nested set of tuples, structured in the manner explained below. 'parse' can handle any SETL 'compilation unit', i.e. anything that the ordinary SETL compiler can handle. 'parse_expr' is a variant form of 'parse', designed to make it a bit easier to work with isolated expressions. 'parse_expr' can handle sequences of expressions and statements, which must end in semicolons. It works as follows: the string submitted to as argument to 'parse_expr' is 'wrapped' by prefixing it with an initial line 'program test;' and appending a final line 'end;'. The resulting string is then parsed as a SETL program. If this operation succeeds, the relevant part of the resulting parse tree is returned. Otherwise **OM** is returned,in which case the number of error messages which are generated is available via the procedure 'setl_num_errors', and the individual error messages via
'setl_err_string'. The following small example shows the kind of constructs that can be parsed successfully:

The output produced by the preceding code is seen below. Details concerning the manner in which this output represents parse trees are explained a bit later.programtest; -- examples of successful parsesuseparser;forjin[1..3]loopx + 1; x := {w * w: winy};end loop;")); -- loops can be parsedproceduref(z);return{z};endf;")); -- any 'program body' (hence package body and class body) can be parsedendtest;

["ast_list", ["ast_add", "X", "1"], ["ast_add", "Y", "1"]] ["ast_list", ["ast_add", "X", "1"], ["ast_assign", "X", ["ast_add", "Y", "1"]]] ["ast_list", ["ast_for", ["ast_iter_list", ["ast_in", "J", ["ast_arith_tup", ["ast_list", "1"], "3"]]], ["ast_null"], ["ast_list", ["ast_add", "X", "1"], ["ast_assign", "X", ["ast_genset", ["ast_mult", "W", "W"], ["ast_iter_list", ["ast_in", "W", "Y"]], ["ast_null"]]]]]] ["ast_list", ["ast_add", "X", "1"], ["ast_add", "Y", "1"]]'parse' handles error messages in the same way as 'parse_expr' and 'compile'. Here is an example showing how error messages generated by failed parses can be recovered.

The output produced by this example isprogramtest; -- example of failed parse, and error message recoveryuseparser;forjin[1..3]loop" + " x +, 1; x := {w * w: w inn y}; endloop;"));forjin[0..ne]loop print(setl_err_string(j));end loop;endtest;

3 [2:26] *ERROR* => syntax error at , [2:29] *ERROR* => syntax error at ; [2:46] *ERROR* => syntax error at INN

**Parse tree structure.** The parse trees returned by the 'parse' function reflect the nested structure of the inputs analyzed.

The overall list of statements is represented by a top-level list node of the form

Constructs involving built-in SETL operators and syntactic forms are represented by special nodes.

Variables and calls to functions other than built-in infix and prefix operators are represented using nodes of the form "_of, followed by the name of the function and its list of arguments. Examples are:

b(x,y) is represented by ["_of", "B", ["_list", "X", "Y"]] b() is represented by ["_of", "B"] b; is represented by the list element "B"

As these examples show, in parse trees all names are capitalized; this gives SETL its case independence.

The node-types used are:

[x,y] is represented by ["_enum_tup", "X", "Y"] {x,y} is represented by ["_enum_set", "X", "Y"] [x] is represented by ["_enum_tup", "X", "Y"], etc. [ ] is represented by ["_nulltup"], etc. [x..y] is represented by ["_arith_tup", ["_list", "X"], "Y"], etc. [x,y,z..w] is represented by ["_arith_tup", ["_list", "X", "Y", "Z"], "W"], etc. -- (note that the parser handles this general case) b(x..y) is represented by ["_slice", "B", "X", "Y"] b(x..) is represented by ["_end", "B", "X", "Y"] b := y is represented by ["_assign", "B", "Y"] --- all the other assignment forms are derived from this. -- For example, b(x..) := y is represented by ["_assign", ["_end", "B", "X"], "Y"] b(x..y) := z is represented by ["_assign", ["_slice", "B", "X", "Y"], "Z"] b(x) := z is represented by ["_assign", ["_of", "B", ["_list", "X"]], "Z"] b(x,y) is represented by ["_of", "B", ["_list", "X", "Y"]] b(x,y) := z is represented by ["_assign", ["_of", "B", ["_list", "X", "Y"]], "Z"] -- superfluous parentheses are elided in the parse tree, -- simplifying it but creating some problems -- when parse trees must be back-correlated with -- the token tuples from which they originate.The setl binary operators ** * /

"_EXPON" "_MULT" "_DIV" "_MOD" "_QUESTION" "_ADD" "_SUB" "_MAX" "_MIN" "_EQ" "_NE" "_LT" "_LE" "_GT" "_GE" "_IN" "_NOTIN" "_SUBSET" "_INCS" "_AND" "_OR" "_FROM" "_FROMB" "_FROME" "_NPOW"The SETL monadic operators '

"_UMINUS" "_NELT" "_ARB" "_DOMAIN" "_RANGE" "_POW"Set and tuple formers like [x: y in s | condition] are represented by structures like

["_gentup", "X", ["_iter_list", ["_in", "Y", "S"], "CONDITION"], etc.If the CONDITION is omitted in the source text it is omitted in the parse tree also. If the leading expression is omitted, "_gentup" and "_genset" are replaced by "_gentup_noexp" and "_genset_noexp" respectively.

All such iterative constructs involve an '"_iter_list' node, whose elements are the successive iterators. These are all expressions, e.g. [y: y = f(x) | condition] is represented by

["_gentup", "Y", ["_iter_list", ["_eq", "Y", ["_of", "F", ["_list", "X"]]]], "CONDITION"].Conversely the parser will accept arbitrary expressions in this position, e.g. [y: u + v = f/x | condition] is represented by

["_list", ["_gentup", "Y", ["_iter_list", ["_eq", ["_add", "U", "V"], ["_div", "F", "X"]]], "CONDITION"]]In the full SETL compiler such constructs are rejected during a post-parse phase; but programs which use the parser directly may find use for them.

Similar remarks apply to other allowed iterators, e.g. 'exists y = f(x) | condition' is represented by

["_exists", ["_ex_iter", ["_eq", "Y", ["_of", "F", ["_list", "X"]]]], "CONDITION"]and likewise for 'forall y = f(x) | condition'. The 'for' iterator construct

for y = f(x) | condition loop a; b; c; end loop;is represented by

["_for", ["_iter_list", ["_eq", "Y", ["_of", "F", ["_list", "X"]]]], "CONDITION", ["_list", "A", "B", "C"]],the body of the loop being represented by a final component of "_list" type. 'while' and 'until' iterators like

while condition loop a; b; c; end loop;are similarly represented by

["_while", "CONDITION", ["_list", "A", "B", "C"]], etc.Conditionals like

if c1 then a; b; c; elseif c2 then d; else e; end if;are represented by

["_if_stmt", "C1", ["_list", "A", "B", "C"], ["_if_stmt", "C2", ["_list", "D"], ["_list", "E"]]]and the corresponding conditional expressions, e.g.

if c1 then a elseif c2 then d else e end if;by

["_if_expr", "C1", "A", ["_if_expr", "C2", "D", "E"]].A 'case' statement like

case a when b,b2 => c; otherwise => d; end case;is represented by

["_case_stmt", "A", ["_list", ["_when", ["_list", "B", "B2"], ["_list", "C"]]], ["_list", "D"]]If the 'otherwise' clause is omitted so is the final component of the subtree representing the case statement, e.g

case a when b,b2 => c; when b3 => d; end case;is represented by

["_case_stmt", "A", ["_list", ["_when", ["_list", "B", "B2"], ["_list", "C"]], ["_when", ["_list", "B3"], ["_list", "D"]]]]similarly case expressions like

case a when b,b2 => c otherwise => d end caseare represented by trees like

["_case_expr", "A", ["_list", ["_when", ["_list", "B", "B2"], "C"]], "D"]Alternative case statement forms like case when b => c; otherwise => d; end case; are represented by trees like

["_guard_stmt", ["_list", ["_when", "B", ["_list", "C"]]], ["_list", "D"]];Again, omission of the 'otherwise' clause simply causes omission of the final component of the corresponding subtree.

Alternative case expression forms like

case when b => c otherwise => d end case;are represented by trees like

["_guard_expr", ["_list", ["_when", "B", "C"]], "D"]The 'return' statement is treated like a monadic operator, so

return x is represented by ["_return", "X"] return is represented by ["_return"]Similarly

exit is represented by ["_exit"] continue is represented by ["_continue"] stop is represented by ["_stop"]Assignment targets like [a,-,b,-,-] containing placeholders are represented by nodes like

["A", _placeholder, "B", _placeholder, _placeholder]Since no semantic checks are performed by the parser, constructs of this kind are allowed in any position.

Quoted strings within Setl Expressions are represented by themselves, with their enclosing double quotes.

Since no semantic checks are performed by the parser, it can handle a somewhat wider range of constructs than the SETL language proper allows, including various expression forms useful in a logic system.

This includes Iterators without restrictions in set formers and quantifiers, e.g {x + y: x, y in z | x > 0}, unrestricted existentals and universals like 'exists x, y in z | x > 0' and 'forall x, y in z | x > 0', and even mover general constructs like 'forall OM | expn(x,y,z,u,v)', which can be used to represent universal quantification all the variables in a general expression expn(x,y,z,u,v). These are seen in the following example

The parse-trees produced areprogramtest; -- parsing of generalized SETL-like expressionsuseparser;inz | x > 0};")); -- set former with unrestricted 'iterator'inz | x > 0;")); -- existential with unrestricted 'iterator'forallx, yinz | x > 0;")); -- universal with unrestricted 'iterator'forallOM | expn(x,y,z,u,v);")); -- universal over 'all free variables'endtest;

["ast_list", ["ast_genset", ["ast_add", "X", "Y"], ["ast_iter_list", "X", ["ast_in", "Y", "Z"]], ["ast_gt", "X", "0"]]] ["ast_list", ["ast_exists", ["ast_ex_iter", "X", ["ast_in", "Y", "Z"]], ["ast_gt", "X", "0"]]] ["ast_list", ["ast_forall", ["ast_iter_list", "X", ["ast_in", "Y", "Z"]], ["ast_gt", "X", "0"]]] ["ast_list", ["ast_forall", ["ast_iter_list", "OM"], ["ast_of", "EXPN", ["ast_list", "X", "Y", "Z", "U", "V"]]]]

If compilation fails, error messages can be recovered using 'setl_num_errors' and 'setl_err_string' in the manner already explained in connection with our discussion of the 'parse' operation.

Here is an example of the use of 'compile' to generate and execute a series of entirely different procedures, all called 'changes_dynamically'.

The three entirely different results produced by the three quite different routines ultimately called are:packagetest_pak;procedurechanges_dynamically(x); -- declare a procedure which will change dynamicallyendtest_pak;programtest; -- parsing of generalized SETL-like expressionsuseparser,ide_pak; -- define header and trailer lines for -- packages and procedures to be compiled var changes_dynamically; -- declare global name for the procedure to be created dynamically header_line := "package bodytest_pak;\nprocedure changes_dynamically(x);\n"; trailer_line := "endchanges_dynamically;\nend test_pak;"; proc_bodies := ["returnx + x;\n", "return str(x) +str(x);\n", "return cos(float(x));\n"]; -- these three strings will be used as procedure bodiesforpbodyinproc_bodiesloopres := compile(header_line + pbody + trailer_line); -- wrap the string containing the procedure body, -- and compile it pak := reload_library_package("TEST_PAK"); -- force load of the newly compiled package body changes_dynamically := pak("CHANGES_DYNAMICALLY"); -- extract the 'changes_dynamically' routineend loop;proceduremy_procedure(x); -- auxiliary static routinereturnchanges_dynamically(x);endmy_procedure;endtest;

4 22 -0.41614683655Note the following facts about the preceding example: (i) Both the name "TEST_PAK" of the package with which we work and the name "CHANGES_DYNAMICALLY" of the procedure retrieved from it must be given as uppercase strings, since (a) this is what 'reload_library_package' and (b) 'reload_library_package' produces a map from capitalized procedure names to the corresponding procedure object. (ii) in order for 'library_package' or 'reload_library_package' to find a package, the package's specifier must have been compiled previously.

These operations allow the IDE edit window to be used as an input console and are crucial to the programming of SETL Help Scripts of the kind described in Section XXX.

The IDE_PAK specifier is

nativepackageide_pak; -- procedures useful only in programs invoked by Help scriptsprocedureget_buffer(); -- get contents of edit windowprocedureget_length(); -- get length of string in edit windowprocedureset_buffer(from_pos,to_pos,string); -- insert string into edit windowprocedureget_selection(); -- get boundaries of current selection in edit windowprocedureset_selection(from_pos,to_pos); -- set boundaries of current selection in edit windowprocedureget_position();procedureget_parameter(); -- Help Scripts, especially those of the form <SCRIPT=Personal4>, -- can be preceded by a text line which serves as a personalizing -- parameter. 'get_parameter()' returns this line -- procedures of general utilityprocedureide_rerun(commandline);procedureget_library_list(); -- get current library list, as a comma-separated string -- of file namesprocedureset_library_list(liblist); -- set current library list, from a comma-separated string -- of file namesprocedurereload_library_package(pak); -- force reload of the named package,which may have been -- newly compiledendide_pak;

We now give a variety of examples illustrating the use of these operations. Since the operations of the first group all refer implicitly to the contents of a SETL IDE edit window, they must be executed as SETL Help Scripts, by dragging a Help item into such a window, in the manner explained in Section XXX. For this it is convenient to use one of the 4 special help items available under the Help Panel tools tab as 'Personal1'-'Personal4'. Each of these is set up so that when dropped on a SETL IDE edit window it executes a program of the same name (which must have been precompiled). For example, the 'Personal1' Help Script is set up as

<Personal1> ANY PARAMETER STRING YOU LIKE <SCRIPT=Personal1>

so that when dropped on an edit window it executes the program named 'Personal1' (which must have been precompiled). Suppose that this has been supplied as

programPersonal1; -- 13:59:58useide_pak; stg := get_buffer(); -- read edit window contentstime()); -- write time to edit windowendPersonal1;

and that the 'Personal1' Help item is dropped on an Edit window containing the text of this small program when characters 87-97 of the Edit window have been selected. Then the output printed will be

461 [87, 97] 462 ANY PARAMETER STRING YOU LIKE

and the current time will replace the original '13:59:58' in the edit window.

As a second example we give the script of the 'capitalize' Help tool item, which is simply

programtest;useide_pak,string_utility_pak; [ix1,ix2] := get_selection(); -- read edit window selection boundariesif abs(ix2 - ix1) /= (ix2 - ix1)then[ix1,ix2] := [1,get_length()];end if; -- if no selection, take whole buffer buffer_part := get_buffer()(ix1..ix2); -- read selected part of edit window set_buffer(ix1,ix2,case_change(buffer_part,"lu")); -- replace by capitalized versionendtest;

The following example shows the use of 'set_library_list' and 'get_library_list'. In order to minimize the risk of including a nonexistent library file in the library list (which would disable SETL) we set the library list to the appropriately comma-separated but redundant 'setl2.lib,setl2.lib'. The output printed by '**print**(get_library_list());' is then of course 'setl2.lib,setl2.lib'.

An example of the use of reload_library_package is given in the preceding section.programtest; -- library-list manipulationuseide_pak; set_library_list("setl2.lib,setl2.lib");endtest;

This small package typifies the small interfaces that can be used to give SETL access to other fully dynamic languages, i.e. languages which support dynamic execution of statements in them. Unless very high efficiency is needed, the interface can consist of just one command which must be added to the external language, and two SETL routines. The command C added to the external language can simply pass a string back to SETL. One of the two SETL interface procedures should then which pass commands from SETL to the external language. The second must capture the string passed by the command C added to the external language.

Note that to use this scheme, one must prepare a slightly modified version of the external language, embodying the added command C mentioned above.

python_pak is just such an interface, for the increasingly popular Python language. The modified version of the Python interpreter that is needed is supplied with the SETL system as the binary module 'PythonCore'. This must be placed in the same folder as the SETL interpreter.

Note that when running in default mode, Python opens an annoying console window that SETL does not use and that the SETL-associated version of Python will only use when it needs to write error messages. To prevent this window from appearing unnecessarily you should invoke the 'EditPythonPrefs' utility that comes with Python and check the 'Delay Console Window Until Needed' option offered under the 'Default Startup Options' section of the Python preferences dialog.

The python_pak specifier is simply

native packagepython; -- SETL interface to Python languageprocedurepython_exec(command); -- pass command string to Python for executionprocedurepython_getstring(); -- gets string returned to SETL by setl.pass_string(stg)endpython;

The command added to the modified 'PythonCore' Python interpreter supplied with SETL is

setl.pass_string(stg)Here, 'stg' can be any string_valued expression in the Python language. This command converts it to a SETL string, which is returned to SETL. Note that this command implicitly uses Python's object syntax; 'setl.pass_string' is a normal Python method call, addressed to the object 'setl' which is built-in to 'PythonCore'.

Using these primitives it is easy to set up SETL procedures which make it easy to use Python, and in particular its extensive built-in libraries. Two such procedures are

procedurepython_calc(expn); -- evaluate a Python expression, returning the value as a string -- request Python to execute 'result = expn', passing 'result' to SETL python_exec("result = " + expn + "\nsetl.pass_string(str(result))");returnpython_getstring(); -- return value of Python expression 'str(result)' to SETLendpython_calc;procedurepython_val(vname); -- get the value of the indicated Python variable python_exec("setl.pass_string(str(" + vname + "))");returnpython_getstring(); -- return value of Python expression 'str(vname)' to SETLendpython_val;

An extensive collection of Internet access an file/folder manipulation utilities, built on the corresponding Python operations, is then available in the **package** net_utilities supplied with the SETL system, The specifier for this package is

(The commented-out lines in this list document a few procedures available in the body ofpackagenet_utilities; -- package of Python-derived internet utilities var telnet_prompt := "%"; -- the telnet prompt characterprocedurenet_init(); -- obligatory initializationprocedurepython_calc(expn); -- evaluate a Python expression, returning the value as a stringprocedurepython_val(vname); -- get the value of the indicated Python variableprocedurehttp_get(url); -- transmit http request in 'get' formprocedurehttp_post(url,stg); -- transmit http request in 'post' form ('old' version)procedurehttp_post_multi(url,map); -- transmit http request in 'post' form ('multi' version)procedureftp_read(url); -- transmit request for ftp file read by blocksprocedureftp_read_lines(url); -- transmit request for ftp file read by linesprocedureftp_list(url); -- transmit request for ftp file-listing of directoryproceduremail_send(host,who_from,whoall_to,subj,msg); -- transmit request for mail sendproceduremail_summary(mail_host_user_pwd,nlines); -- transmit request for summary of mailboxprocedureread_gzip(file_name); -- read a Gzipped file -- ******** telnet routines ********proceduretelnet_open(host); -- open connection to telnet host; name can be qualified by portproceduretelnet_close(host); -- close connection to telnet host;proceduretelnet_read_until(host,stg); -- read telnet stream up to indicated markerproceduretelnet_expect(host,regex_list,secs); -- read telnet stream up to regexp, or timeoutproceduretelnet_read_timeout(host,stg,secs); -- read telnet stream up to marker, or timeoutproceduretelnet_write(host,stg); -- write telnet commandproceduretelnet_response(host,stg); -- write telnet command and get responseproceduretelnet_put_file(host,file_name,dest_file_name); -- send a text file via telnetproceduretelnet_put_stg(host,dest_file_name,stg); -- send a string to start a file via telnetproceduretelnet_add_stg(host,dest_file_name,stg); -- send a string to start adding to a file via telnetproceduretelnet_end_stg(host); -- end a string sent via telnetproceduretelnet_continue_stg(host,dest_file_name,stg); -- add a string to the end of a file via telnetproceduretelnet_get_file(host,source_file_name,dest_file_name); -- get a text file via telnetprocedurelogin(host,who,pwd); -- telnet login procedure --proceduresanitize(stg); -- sanitize a string containing quotes and backslashes -- ******** file utilities ********procedurefile_chdir(path); -- change current working directory to "path"procedurefile_getcwd(); -- get current working directoryprocedurefile_listdir(path); -- get contents of indicated directoryprocedurefile_mkdir(path); -- make a new directoryprocedurefile_remove(path); -- remove a fileprocedurefile_rmdir(path); -- remove a directoryprocedurefile_rename(src,dst) ; -- rename a file or directoryprocedurefile_stat(path); -- file status tuple, in order st_mode, ino, dev, --nlink, uid, gid, size, atime, mtime, ctimeprocedurefile_kind(path); -- returns 1 of directory, 0 for pathprocedurefile_set_times(path,creation,modif); -- set file creation and modification timesprocedurefile_set_creator_type(path,creator,ftype); --set file creator and type 4-char Macintosh fieldsprocedurepython_ctime(floating_secs); -- convert python floating point time to standard rependnet_utilities;

**The 'http' routines.** It is generally not hard to examine HTML pages visually and unravel their structure. Done successfully, this makes many Web information services available for use as a kind of subprocedure for programs you may wish to construct. For example, the well-known 'CNNfn' web site makes frequently updated stock- and bond-market quotes available in 4 forms: Dow-Jones average, NASDAQ average, Standard&Poor average, and Bond average. Suppose we want to access this data for continual use within a program (perhaps a Web application) that we are writing. As of Feb. 2001, CNN is putting this data into an HTML table with the HTML comment of the form '' near the start of the table (as can be seen by inspecting the HTML source text of their page). Hence we can begin to locate our desired data by taking the page text between the first occurrence of '' and the next following occurrence of ''. These occurrences are located by the code lines

stg := http_get("www.cnnfn.com"); -- download cnnfn home page dowlocs := kmp(stg,""); -- get occurrence locations of '' etlocs := kmp(stg,""); -- get occurrence locations of ''If we print the variables 'dowlocs' and 'etlocs' the result is

[14844, 19200] [2043, 3186, 3270, 3913, 12619, 14589, 20521, 26362, 30704, 30722, 30950, 33395, 33417, 34608, 35112, 35181, 42182, 42206, 51780, 53415, 54100, 56451, 56474, 56513, 56535, 57498]This makes it clear that '' has only two occurrences, both within the HTML table containing the data we are interested in. Hence

exists locwill locate the end of this table. Breaking the range of text from the first occurrence of '' to the end of the table into lines and examining the result, we easily see that the desired data is found in lines 4,8,12,16,27,28,30,31,33,34,36 and 37 of the table. The following program assembles these observations into code that extracts the stock-market data we want from the 'CNNfn' home page.inetlocs | loc > dowlocs(2)

Note that the data tuple as first constructed has the formprogramtest; -- getting market quotes from cnnfn websiteusepython;usenet_utilities,stringm,string_utility_pak;constmarket_names:= ["DOW","NASDAQ","S&P","BOND"]; stg := http_get("www.cnnfn.com"); -- download cnnfn home page dowlocs := kmp(stg,""); -- get occurrence locations of '' etlocs := kmp(stg,""); -- get occurrence locations of ''if#dowlocs /= 2or not(exists locinetlocs | loc > dowlocs(2))then-- assumed cnnfn page format has changedstop; -- give warning and stopend if; tup := breakup(stg(dowlocs(1)..loc),"\n\r"); --breakpage section between '' and '' into lines tup := [tup(j): jin[4,8,12,16,27,28,30,31,33,34,36,37]]; -- get lines containing desired data data := [back3(line): lineintup]; -- extract numerical data from right end of linesforjin[1..4] | data(2 * j + 4)(1) = "-"loopdata(j) := "-" + data(j); -- capture signs for first 4 data itemsend loop;formname = market_names(j)loopj2 := 2 * j + 3;end loop;procedureback3(stg); -- extract numerical data from right end of lines; -- remove other HTML clutterforjin[1..3]looprbreak(stg,"<"); -- back over 3 HTML tagsrmatch(stg,"<");end loop; data :=rbreak(stg,">"); -- extract numerical data, back to next tag on right data := suppress_chars(data,"nbsp;"); -- suppress 'nbsp;'return"" +/ [ifc = "&"then" "elsecend if: cindata]; -- change remaining '&' signs to blanksendback3;endtest;

[" 43.45", " 61.94", " 11.51", " 1/8", "10903.32", "-0.40%", "2427.72", "-2.49%", "1318.80", "-0.87%", "99 9/32", " 5.42%"]The first 4 items in this tuple lack the '-' sign that they need if the market is falling. (This is because these signs appear in cnnfn's page as icons). But the necessary signs are readily extracted from the corresponding percentage changes in data items 6, 8, 10,and 12. Doing this, our data tuple reads

["- 43.45", "- 61.94", "- 11.51", " 1/8", "10903.32", "-0.40%", "2427.72", "-2.49%", "1318.80", "-0.87%", "99 9/32", " 5.42%"]As of the early morning of Feb. 14, 2001 the output of the above program is

DOW 10903.32 - 43.45 -0.40% NASDAQ 2427.72 - 61.94 -2.49% S&P 1318.80 - 11.51 -0.87% BOND 99 7/32 3/16 5.42%As of 10:13 AM of Feb. 14 it produces

DOW 10855.29 - 48.03 -0.44% NASDAQ 2421.29 - 6.43 -0.26% S&P 1312.06 - 6.74 -0.51% BOND 99 5/16 3/32 5.42%(This illustrates the remark allegedly made by J.P. Morgan when asked to predict the future course of the stock market; "The market will fluctuate.").

Here is another example, The 'cnnfn' page we have been discussing also provides a stock quotation service. One can enter comma-separated 3-letter codes for several stocks (e.g. IBM,CNN) into a small text box and click on a 'Go' button. An HTML page containing the requested data is then generated and returned. If this is done in a browser it will be seen that the URL being requested is simply

That is, the parameter for the request (the part of this URL that follows the separating '?' character) is simply 'symbols=ibm%2Ccnn', essentially the two stock symbols separated by the 'Internet sanitized' form '%2C' of the comma character. This easy visual inspection makes it plain that the 'http_get' call

should retrieve this same data, as indeed it does. It is generally easy to read 'get'-generated URLs in just this way and so to write small routines which make equivalent Web services available as SETL procedures built using 'http_get'.

**http_post** The http_post operation has the form

It is generally used when a relatively large 'stg' parameter needs to be transmitted along with a Web URL request. In interactive pages, this data is commonly collected from For example, the cnnfn pare we have been discussing

**The 'ftp' routines.** 'ftp', or 'file transfer protocol' is a standard communication convention for the download of substantial files. Many free file collections accessible by ftp are available on the Web. The ftp_read, ftp_read_lines, and ftp_list procedures make these files available to SETL. The first two of these routines have the form ftp_read(url) and ftp_read_lines(url), and make the requested file available as a SETL string. For example, many files relevant to the Python system (including full source files) are available via ftp at www.python.org/pub/python/, and either of the expressions

ftp_read("www.python.org/pub/python/README")or

ftp_read_lines("www.python.org/pub/python/README")will retrieve the small README file available at that location. The two procedures differ only in the manner that the file is fetched over the Web, and so both produce the same result. ftp_read_lines requests that the file be read by lines, and ftp_read that it be read by blocks. The second method should be substantially more efficient for large files and more appropriate for binary files.

The ftp_list procedure has the form ftp_list (url), but here the url should reference a directory of other files rather than a bottom-level file. For example, a directory of python binaries is available at ftp://www.python.org/bin/ftp-exec/, so

ftp_list("www.python.org/bin/ftp-exec"can be used to retrieve a UNIX-form listing of this directory, which appears as

drwxr-xr-x 2 root bin 512 May 19 1998 . drwxr-xr-x 3 root wheel 512 May 19 1998 .. -r-xr-xr-x 1 root bin 17500 Jun 11 1997 ls -rwxr-xr-x 1 root bin 163 Jun 11 1997 rmake -rwxr-xr-x 1 root bin 1022 Jun 11 1997 sigjob -rwxr-xr-x 1 root bin 851 Jun 11 1997 sigjoblog -rwxr-xr-x 1 root bin 19 Jun 11 1997 tme -r-xr-xr-x 1 root bin 6816 Jun 11 1997 touch

As noted in in Section 10.19 of the preceding Chapter, all higher-level Internet communication function sare built upon the 'socket' mechanism, which provides the basic capability of sending byte-streams reliably between two computers, with certain (but not detailed) notification of communication failures.

Once socket communication has been opened between two computers, communication proceeds by the exchange of byte streams, which can be thought of as strings.

To understand the issues which any such conversation must involve we can think of the computer-to-computer socket communication as an exchange between two persons, initiated when a phone rings and one of them picks it up. The person called is normally expected to say something, as an indication that communication has been initiated and that further responses will be forthcoming. Whatever is said is expected to be in an understandable language. If after calling someone what you hear on the telephone is a mysterious series of clicks or beeps, or something in a completely unknown language, for example, 'Nov Schmoz kapop', you will normally hang up rather than trying to proceed with the conversation. Likewise, if the person called says 'Hello, who is this?', and the caller responds only with things like 'Nov Schmoz Kapop', the person called will normally terminate the conversation. Only when the conversation proceeds in conformity with mutually-understood conventions can it successfully continue.

A second issue has to do with the termination of messages coming from one or the other side. Suppose somebody telephones you and that you respond by saying 'Hello, who is there?' The caller then begins to speak in perfectly comprehensible English, but then goes on endlessly without ever stopping, or still worse, pauses occasionally, but without ever clearly finishing and then resumes. In this case, sensible conversation is impossible since one never knows when a response is called for. This shows that productive bilateral communication must depend on clear, mutually understood signals defining the termination of messages.

A final consideration relates to communication delays. Suppose a phone conversation goes like this:

- You telephone somebody.
- They say "Hello, who's calling?"
- You identify yourself successfully.
- They then respond "Hold on a minute, the doorbell is ringing."
- You wait.
- A hour goes by.
- At this point you don't know what has happened at the other end and whether the communication will ever continue.

To recover from such situations the most sensible procedure is to wait for some reasonable length of time and then hang up if no response is received.

The telnet routines lie closer to the raw process of socket communication than do the other Web-related routines provided by NET-UTILITIES and therefore must handle all of these conversation-coordination problems. Telnet communication proceeds by the exchange of punctuated strings. Since the punctuation characters used play a special role they cannot be used in quite the ordinary way and so must be specially 'sanitized' when they appear in ordinary strings. A subroutine 'sanitize' is provided for this purpose, (but is used internally by all the other routines, and so need not be applied explicitly.)

'telnet_open' initiates communication with a specified host; 'telnet_close' ends this communication. 'telnet_write' sends a string to a designated host and 'telnet_read_until' reads the response from this host. However, since we must be able to recognize the punctuation mark which ends the host's response 'telnet_read_until' has a second parameter 'stg' defining the sequence of characters in the host's response which signal the end of this response. This operation is provided in three variant forms: 'telnet_read_until' waits for the occurrence of a single specified string in the response stream and is prepared to wait forever. (Waiting forever is only a sensible thing to do in a multi-threaded environment. In a single threaded SETL environment waiting forever will hang SETL and possible also the system under which it is running.) The 'telnet_read_timeout' version of 'telnet_read_until' has a third time-limit parameter, given in seconds which indicates how long one should wait until hanging up if a complete response has not yet been received. The third variant 'telnet_expect' of the same primitive allows one to wait, not just for the occurrence of a specific string in the response stream, but for any string matching a specified regular expression. This gives one a way of waiting for the occurrence of any one of a collection of strings, any one of which might signal the end of a response. For example, if one is communicating with a UNIX system and uses the UNIX 'more' operation, the response generated will end with one kind of prompt character if the whole of a requested file is returned, but with a different prompt character if only part of the requested file has been sent. So the resulting communication can end with either of these characters; one must be prepared to resume one's own end of the conversation in either case.

One additional procedure in the group of telnet routines, 'telnet_put_file' transmits a text file via telnet, simply by sending its lines one after another.

It is instructive to trace through the series of communications which occur when one accesses a remote computer via telnet and logs in the ordinary way. One begins by signaling the remote computer via 'telnet_open'. It is then the remote computer's responsibility to send back a message, which can consist of any amount of announcement and advertising material that it wants to send, but which must end with the word 'login:'. This word is understood to mark the end of the remote computer's initial response, and so cannot occur elsewhere in this response. Thus to get this initial response one executes 'telnet_read_until' with 'login:' as its second parameter. Once the whole of this login request has been received, it is the calling client's computer's responsibility to send an account name. The remote computer will respond to this account name by sending a second message, which again can be anything, but which must end in the word 'Password:', which signals that a password is expected. The 'client' computer must then reply to this response by sending an acceptable password. Once this is done the remote computer will check the password, and respond with one of two types of messages. The login request may be rejected because the account does not exist or the password is wrong, but if the account exists and the password is correct it is the remote computer's responsibility to reply with a message ending in a mutually-understood prompt character or sequence, typically a '$' or '%'. At this point the client will understand that it has been put in communication with whatever type of 'shell program' the remote computer uses to deal with remote logins. Typically this is a UNIX shell. From this point on communication must proceed in conformity with the command line conventions defined by this UNIX shell. That is, the client computer is expected to send valid UNIX commands or other character sequences defined by UNIX, and the remote computer is expected to respond in the standard UNIX fashion using messages that always end with standard prompt characters or strings.

The remaining routines of the 'telnet' group are utility composites which could be programmed using the routines we have described already, but which are provided to facilitate common operations. 'telnet_put_stg' opens a file on the remote computer and writes a string to it; 'telnet_add_stg' adds a specified string to a file that has already been created. 'telnet_continue_stg' . 'telnet_end_stg' signals the end of a string being sent via telnet. Examples illustrating the use of these operations are given below. 'telnet_get_file' reads an entire file via telnet. 'login' is a utility procedure which handles all of the steps of a login operation, given the Web address of the remote host, plus an account name and password.

The following example shows the use of some of the main telnet primitives. Note that to execute it successfully, you will need to change the values shown for 'host', 'account_name', and 'password' to whatever values are appropriate for somehost, account, and password to which you have access. The telnet_prompt character will also have to be changed if the host system you are accessing uses a different prompt.

programtest; -- communication via telnetusenet_utilities,get_lines_pak; host := "falce.mrl.nyu.edu"; account_name := "setl"; -- define host to be used, and account name password := "nine98"; telnet_prompt := "$"; -- define password, and expected prompt telnet_open(host); -- open communication with the hostformyfile!!!"); -- create a Unix file and write to it telnet_add_stg(host,"myfile"," more contentsformyfile!!!"); -- add to the Unix file telnet_end_stg(host); -- close the Unix fileendtest;

The code shown above logs in explicitly. Telnet sessions more typically begin with a call to the 'login' procedure, whose calling sequence is

login(host,who,pwd);This does much the same thing as the first part of the code shown above. It is simply

Rewritten to use this utility, the code above is justprocedurelogin(host,who,pwd); -- telnet login procedure telnet_open(host); -- open connection to telnet host; telnet_read_until(host,"login: "); -- read telnet stream up to indicated marker telnet_write(host,who); telnet_read_until(host,"Password: ");returntelnet_response(host,pwd); -- write telnet pwd and get responseendlogin;

programtest; -- communication via telnetusenet_utilities; login(host,who,pwd); -- loginformyfile!!!"); -- create a Unix file and write to it telnet_add_stg(host,"myfile"," more contentsformyfile!!!"); -- add to the Unix file telnet_end_stg(host); -- close the Unix fileendtest;

file_mkdir(path); file_rmdir(path); and file_listdir(path);respectively create, delete, and list the contents of a directory specified in whatever way the operating system expects. (For example, on the Macintosh a file path would have a colon-separated form something like 'my_disk:my_directory:my_subdirectory:my_file', while under Unix it would have a '/'-separated form like 'my_master_directory/my_directory/my_subdirectory/my_file'. When a directory is listed using 'file_listdir' the listing generated is a tuple of tuples tf, each corresponding to one of the files in the directory. The tuples tf contain the file's name, followed by whatever additional information the operating system maintains for files. For example, on the Macintosh this extra information is a pair [creator,file_type] of 4-letter codes, the first designating the application that has (nominally) created the file, and the second designating the type of the file. Thus a typical Macintosh file-listing triple would be something like

["junkfile", "SETL", "TEXT"]If 'path' identifies a file then the expression

file_kind(path)retrieves all components of the kind of tuple (e.g. triple) described above (but not the file's name). If the file identified by 'path' is a directory or does not exist, then an empty tuple is returned.

file_remove(path);deletes a file (but file_rmdir(path); should be used instead if 'path' identifies a directory rather than an ordinary file). file_remove(path); will only work as expected if the file identified by 'path' is not being held open by any currently active application (i.e. has been properly closed). Similarly, file_rmdir(path); will only work if the directory identified by 'path' is empty.

Many operating systems make use of a notion of 'current working directory'. This is the directory to which an application (in our case, SETL) will write files whose name are not prefixed by any information defining the folder or drive on which they are found. The path name of this directory is retrieved by the expression

file_getcwd()returns the string name of this file.

All of these considerations are reflected in the following short example (written for the Macintosh).

programtest; -- file utilities demousenet_utilities;close(open(testfile := test_directory + ":junkfile","TEXT-OUT")); -- create a file in the test directoryendtest;

Notice that (i) after creating a new directory called 'junkdir' we create a file named 'junkfile' within it by opening such a file for writing, but then close this file immediately so that it can be removed; (ii) we attempt to delete 'junkdir' but this does not work initially since the directory is not empty. However, it does work once 'junkfile' is removed by executing file_rmdir(test_directory). (iii) files created using the SETL 'open' operation have all-zero Macintosh 'creator' and 'file_type' codes (as distinct from files created using the edit facilities of the SETL Interactive development environment, which have the creator code "SETl" and the file_type code "TEXT").

All this is seen in the output of the preceding sample program, whose output is

Diana:FromGiuseppeLATEST [] ["R*ch", "TEXT"] [["junkfile", "R*ch", "TEXT"]] [["junkfile", "R*ch", "TEXT"]] []The use of the operation

file_set_creator_type(file,creator,ftype)is illustrated by the following short code (taken from the SETL IDE help tool 'standardize_types'), which converts all the files in a list of files represented by separate lines in a SETL IDE edit window to the creator type and file type of first file in the list.

The operationprogramtest; -- converts all files in list to creator type and file type -- of first file in the listusenet_utilities,ide_pak,string_utility_pak; -- get lines in SETL edit window and Macintosh kind-codes -- for first file [creator,ftype] := file_kind((files := breakup(get_buffer(),"\n\r"))(1)); -- set same kind set same kind-codes for every other file in listforfileinfiles(2..)loopfile_set_creator_type(file,creator,ftype);end loop;endtest;

file_rename(src,dst)changes the name of the file named 'src' to 'dst'. Note however that both 'src' and 'dst' must be simple file names, neither qualified by an indication of the folder to which it belongs. Hence the 'file_rename' operation cannot be used to change the folder to which a file belongs. Indeed, to apply this operation to a file one must first change the 'current working directory' to be the directory containing this file. This can be done using the operation

file_chdir(path); -- change current working directory to "path"

This group of procedures manipulates opaque objects representing vectors, matrices, and 3-index tensors with floating, double-precision, integer, and long integer values.native packagelibnr;procedurenr_get_tensor(tuple_of_tuple_of_tuple); -- create tensor from tuple_of_tuple_of_tupleprocedurenr_get_matrix(tuple_of_tuple); -- create matrix from tuple_of_tuple of floatsprocedurenr_get_dmatrix(tuple_of_tuple); -- create double precision matrix from tuple_of_tuple of floatsprocedurenr_get_imatrix(tuple_of_tuple); -- create integer matrix from tuple_of_tuple of intsprocedurenr_get_vector(tuple); -- create vector from tuple of floatsprocedurenr_get_dvector(tuple); -- create double precision vector from tuple of floatsprocedurenr_get_ivector(tuple); -- create integer vector from tuple of intsprocedurenr_get_lvector(tuple); -- create long integer vector from tuple of intsprocedurenr_get_cvector(tuple); -- create unsigned byte vector from tuple of unsigned ints -- no greater than 255procedurenr_get_zero_vector(dim,stg); -- create a zero vector of a specified type; "d", "f", or "i"procedurenr_vect_max_index(vect); -- index of maximum element of a vectorprocedurenr_vect_reverse(vect); -- reverse a vectorprocedurenr_get_element(nrobject,index_tuple); -- get component of vector, matrix, etc using tuple of indicesprocedurenr_get_object(nrobject); -- return setl version of a matrix, vector or tensorprocedurenr_ptr_contents(opaque,typ);procedurenr_ptr_update(opaque,typ,element);procedurenr_get_refcount(obj); -- return the setl reference count of an nr-objectprocedurenr_print(opaque); -- print a native matrix, vector or tensorprocedurenr_mat_sum(matA, matB); -- matrix sum of two matrices of like shapeprocedurenr_mat_diff(matA, matB); -- matrix difference of two matrices of like shapeprocedurenr_mat_prod(matA, matB); -- matrix product of two matrices of compatible shapeprocedurenr_matconst_sum(matA, c); -- matrix sum of matrix with constantprocedurenr_matconst_diff(matA, c); -- matrix difference of matrix and constantprocedurenr_const_diffmat(matA, c); -- matrix difference of constant and matrixprocedurenr_matconst_prod(matA, c); -- matrix product of matrix and constantprocedurenr_vect_prod(matA, matB); -- componentwise vector productprocedurenr_cvect_conj(vect); -- complex conjugate of a complex vector represented -- in r,i,r,i,... formprocedurenr_cvect_prod(v1,v2); -- complex product of two complex vectors represented -- in r,i,r,i,... formprocedurenr_mat_clone(mat); -- clone of matrix or vectorprocedurenr_mat_fix(fmat); -- convert a matrix or vector to fixprocedurenr_mat_float(fmat); -- convert a matrix or vector to floatprocedurenr_mat_double(fmat); -- convert a matrix or vector to doubleprocedurenr_mat_transpose(fmat); -- transpose a matrix or vectorprocedurenr_rows_and_cols(mat); -- return number of rows and columns of any matrixprocedurenr_get_slice(matOrVect, tuple); -- return a rectangular portion or one dimensional slice of -- a vector or matrix. The tuple argument is [i,j] or [l,t,r,b].procedurenr_set_vect_slice(rwVect, index, replacement); -- set specified slice of a vector, starting at index, -- containing replacement --#################################################################### --# --# CHAPTER 2 - Solution of Linear Algebraic Equations --# --#################################################################### -- -- ******* Matrix inversion and decomposition routines ******* --proceduregaussj(rwfmatA,rwfmatB); -- Gauss-Jordan solution of Ax = b for square A; -- replaces A by its inverse, and B by the solutionprocedureludcmp(rwfmatA,rwivectB,rwfloatD); -- LUP decomposition. The matrix A is replaced by -- its LU decomposition, with (implicits) 1's in the diagonal of -- L, and the integer vector B (which must be suppied as an input) -- returned as the permutation needed -- in the factorization; D (which must be suppied) -- is the sign of this permutation.proceduremprove(fmatA,fmatALUD,ivectINDX,fvectB,rwfvectX); -- improves the accuracy of the solution of a linear -- equation obtained by the LU methodprocedurelubksb(fmatA,ivectB,rwfvectD); -- solution-by backsubstitution routine. Here fmatA and ivectB -- should be the matrix and vector returned -- by the ludcmp routine. The vector fvectD is converted -- into the solution of Ax = D, where A is the -- original square matrix.proceduretridag(fvectA,fvectB,fvectC,fvectR,rwfvectU); -- solution of tridiagonal system. A, B, and C are the subdiagonal, -- diagonal, and superdiagonal respectively. -- r is the right-hand side; U, which must be supplied as an input, -- is the output. The first component of -- A and the last component of C must be suppied, but are ignored.procedurebanmul(fmatA,integerM1,integerM2,fvectX,rwfvectB); -- multiplies a pair of matrices stored in band diagonal form.procedurebandec(rwfmatA,integerM1,integerM2,rwfmatAL,rwulvectINDX,rwfloatD); -- calculates lup factorization of a band-diagonal matrix.procedurebanbks(fmatA,integerM1,integerM2,fmatAL,ulvectINDX,rwfvectB); -- backsubstitution routine for matrices stored in band diagonal form.proceduresvdcmp(rwfmatA,rwfvectW,rwfmatV); -- singular value decomposition of the matrix A, which must be -- an m by n matrix with a larger number m of rows -- than n of columns. The matrix is factored into the product -- of an orthogonal matrix U of the same dimensions, -- which replacesproceduredsvdcmp(rwdmatA,rwdvectW,rwdmatV); -- double precision singular value decomposition of a matrix. -- See the comment to 'svdcmp'proceduresvbksb(fmatU,fvectW,fmatV,fvectB,rwfvectX); -- solution-by backsubstitution routine which solves Ax = B -- using the matrices U, V, and W procuced by svdcmp. -- the vector X, which must beproceduredsvbksb(dmatU,dvectW,dmatV,dvectB,rwdvectX); -- double precision backsubstitution routine. -- See the comment to 'svbksb'procedurecyclic(fvectA,fvectB,fvectC,floatALPHA, floatBETA,fvectR,rwfvectX); -- -- ******* Sparse matrix routines ******* --proceduresprsin(fmatA,floatTHRESH,rwfvectSA,rwulvectIJA); -- converts a square matrix A into row-indexed sparse form. -- THRESH is the cutoff value for nonzero elements. -- SA is the packed array of retained elements; -- IJA is an integer array that indicates the matrix position -- of these elements. -- The (rather complex) storage rule used is as follows:proceduresprsax(fvectSA,ulvectIJA,fvectX,rwfvectB); -- multiplies standard vector by sparse matrixproceduredsprsax(dvectSA,ulvectIJA,dvectX,rwdvectB); -- multiplies standard vector by sparse matrix, -- double precisionproceduresprstx(fvectSA,ulvectIJA,fvectX,rwfvectB); -- multiplies standard vector by transpose of sparse matrixproceduredsprstx(dvectSA,ulvectIJA,dvectX,rwdvectB); -- multiples standard vector by transpose of sparse matrix, -- double precisionproceduresprstp(fvectSA,ulvectIJA,rwfvectSB,rwulvectIJB); -- construct the transpose of a sparse matrixproceduresprspm(fvectSA,ulvectIJA,fvectSB,ulvectIJB, rw fvectSC,rwulvectIJC); -- forms product of two sparse matricesproceduresprstm(fvectSA,ulvectIJA,fvectSB,ulvectIJB, floatTHRESH,rwfvectSC,rwulvectIJC); -- multiply sparse matrix SA by the transpose of sparse matrix SB, -- returning the result as SC -- -- ******* Routines for specialized matrices ******* --procedurenrlinbcg(dvectB,rwdvectX,integerITOL,doubleTOL, integerITMAX,rwintegerITER,rwdoubleERR,vfuncASOLVE,vfuncATIMES); -- solution of Ax = b by conjugate gradient method. -- subroutine for linbcg, which mimics multiplication by -- a sparse matrix, not explicitly given -- subroutine for linbcg, which mimics multiplication by -- inverse of a preconditioner matrix, -- not explicitly givenproceduresnrm(dvectSX,integerITOL); -- compute l2 or max morm of vector SX, signallled by ITOL <= 3procedurevander(dvectX,rwdvectW,dvectQ); -- solution of Vandermonde linear systemproceduretoeplz(fvectR,rwfvectX,fvectY); -- solution of Toeplitz linear systemprocedurecholdc(rwfmatA,rwfvectP); -- Cholesky decomposition of matrixprocedurecholsl(fmatA,fvectP,fvectB,rwfvectX); -- solution of linear system using Cholesky decompositionprocedureqrdcmp(rwfmatA,rwfvectC,rwfvectD,rwintegerSIGN); -- QR decomposition of matrixprocedurersolv(fmatA,fvectD,rwfvectB); -- solution of upper triangular linear systemprocedureqrsolv(fmatA,fvectC,fvectD,rwfvectB); -- solution of linear system using QR decompositionprocedureqrupdt(fmatR,rwfmatQT,rwfvectU,fvectV); -- update QR decomposition after additon of 1-D vectorprocedurerotate(rwfmatR,rwfmatQT,integerI,floatA,floatB); -- subroutine for qrupdt: apply specified Jacobi rotation -- to a pair of matrices --#################################################################### --# --# CHAPTER 3 - Interpolation and Extrapolation --# --####################################################################procedurepolint(fvectXA,fvectYA,floatX,rwfloatY,rwfloatDY); -- polynomial interpolation of vector of valuesprocedureratint(fvectXA,fvectYA,floatX,rwfloatY,rwfloatDY); -- rational interpolation of vector of valuesprocedurespline(fvectX,fvectY,floatYP1,floatYPN,rwfvectY2); -- compute second derivatives of cubic spline thru given pointsproceduresplint(fvectXA,fvectYA,fvectY2A,floatX,rwfloatY); -- spline interpolation of vector of values, using output of 'spline'procedurelocate(fvectXX,floatX,rwulongJ); -- binary search in real arrayprocedurehunt(fvectXX,floatX,rwulongJLO); -- binary search in real array, using estimated target positionprocedurepolcoe(fvectX,fvectY,rwfvectCOF); -- solves Vandermonde equation for coefficients -- of interpolating polynomial given values at enough pointsprocedurepolcof(fvectXA,fvectYA,rwfvectCOF); -- variant, slightly mopre stable, solution of -- Vandermonde equation for coefficients -- of interpolating polynomial given values at enough pointsprocedurepolin2(fvectX1A,fvectX2A,fmatYA,floatX1,floatX2,rwfloatY,rwfloatDY); -- intepolates 2-variable polynomial given values -- on rectangular array of pointsprocedurebcucof(fvectY,fvectY1,fvectY2,fvectY12,floatD1,floatD2,rwfmatC); -- (n must be 4) calculates coefficients of bicubic spline -- interpolation on a single grid square (subroutine for splie2 and splin2)procedurebcuint(fvectY,fvectY1,fvectY2,fvectY12,floatX1L,floatX1U,floatX2L, floatX2U,floatX1,floatX2,rwfloatANSY,rwfloatANSY1,rwfloatANSY2); -- interpolates bicubic spline on a single grid square -- (subroutine for splie2 and splin2)proceduresplie2(fvectX1A,fvectX2A,fmatYA,rwfmatY2A); -- computes second-derivative table for 2-D bicubic splineproceduresplin2(fvectX1A,fvectX2A,fmatYA,fmatY2A,floatX1,floatX2,rwfloatY); -- computes 2-D bicubic spline interpolation using -- second-derivative table produced by 'splie2' --#################################################################### --# --# CHAPTER 4 - Integration of Functions --# --####################################################################proceduretrapzd(ffuncF,floatA,floatB,integerN); -- computes successive refinements for extended trapezoidal rule; -- subroutine for 'qtrap'procedureqtrap(ffuncF,floatA,floatB); -- numerical 1-D function integration by trapezoidal ruleprocedureqsimp(ffuncF,floatA,floatB); -- numerical 1-D function integration by Simpson's ruleprocedureqromb(ffuncF,floatA,floatB); -- numerical 1-D function Romberg integrationproceduremidpnt(ffuncF,floatA,floatB,integerN); -- computes sucessive stages of refinement for extended midpoint rule -- open interval integration (subroutine for 'qromo') -- Romberg integration adapted to use with improper integrals -- over open intervalproceduremidinf(ffuncF,floatA,floatB,integerN); -- computes sucessive stages of refinement for extended midpoint -- rule open interval integration in 1/x -- (subroutine for 'qromo')proceduremidsql(ffuncF,floatA,floatB,integerN); -- variant of 'midpt' allowing for sqrt singlarity at lower limit -- of integration range (subroutine for 'qromo')proceduremidsqu(ffuncF,floatA,floatB,integerN); -- variant of 'midpt' allowing for sqrt singlarity at upper limit -- of integration range (subroutine for 'qromo')proceduremidexp(ffuncF,floatA,floatB,integerN); -- variant of 'midpt' allowing for infinite upper limit of -- integration (subroutine for 'qromo')procedureqgaus(ffuncF,floatA,floatB); -- numerical 1-D function integration by Gauss points-and-weights methodproceduregauher(rwfvectX,rwfvectW,n); -- get Hermite polynomial points-and-weightsproceduregaucof(rwfvectA,fvectB,floatAMU0,rwfvectX,rwfvectW); -- computes recurrence coeffients for polynomials orthogonal -- in a weight with given moment, given coeffients defining recurrence -- coeffients for an auxiliary set of orthonal polynomialsproceduregaujac(rwfvectX,rwfvectW,n,floatALF,floatBET); -- get Jacobi polynomial points-and-weightsproceduregaulag(rwfvectX,rwfvectW,n,floatALF); -- get Laguerre polynomial points-and-weightsproceduregauleg(floatX1,floatX2,rwfvectX,rwfvectW,n); -- get Legendre polynomial points-and-weightsprocedureorthog(fvectANU,fvectALPHA,fvectBETA,rwfvectA,rwfvectB); -- calculates recurrence coefficients for general orthogonal -- polynomials, given moments -- integration of function over 3-D region with curved boundaries -- by recursive use of Gauss methodprocedurepythag(floatA,floatB); -- sqrt(a**2 + b**2)proceduredpythag(doubleA,doubleB); -- sqrt(a**2 + b**2), double precision --#################################################################### --# --# CHAPTER 5 - Evaluation of Functions --# --####################################################################procedurechebev(floatA,floatB,fvectC,floatX); -- evaluates Chebyshev polynomial representation, given its coefficientsprocedurechebft(floatA,floatB,rwfvectC,ffuncFUNC); -- fits Chebyshev approximation with specified number of terms -- to specified functionprocedurechder(floatA,floatB,fvectC,rwfvectCDER); -- computes coefficients of differentiated Chebyshev approximationprocedurechint(floatA,floatB,fvectC,rwfvectCINT); -- computes coefficients of integrated Chebyshev approximationprocedurechebpc(fvectC,rwfvectD); -- convert Chebyshev polynomial sum to ordinary polynomialprocedurepcshft(floatA,floatB,rwfvectD); -- map polynomial to a translated intervalprocedurepccheb(fvectD,rwfvectC); -- convert polynomial to Chebyshev polynomial sumprocedurepade(rwdvectCOF,rwfloatRESID); -- converts polynomial approximation to function into coeffients -- of Pade rational approximationprocedureratlsq(dfuncFUNC,doubleA,doubleB,integerKK,dvectCOF,rwdoubleDEV); -- find Chebyshev polynomia-basedl rational approximation -- to given functionprocedureeulsum(rwfloatSUM,floatTERM,integerJTERM,rwfvectWKSP); -- accelerates series convergence by Euler differencing technique. -- Should be called repeatedly with successive terms of series.procedureddpoly(fvectC,doubleX,rwfvectPD); -- evaluates poolynomial, with specifed number of derivatives, -- given coeffientsprocedurepoldiv(fvectU,fvectV,rwfvectQ,rwfvectR); -- polynomial division routine for polys with real coefficientsproceduredfridr(ffuncFUNC,floatX,floatH,rwfloatERR); -- numerical approximation of derivative by Ridder extrapoloation method --#################################################################### --# --# CHAPTER 6 - Special Functions --# --####################################################################procedurebeschb(doubleX,rwdoubleGAM1,rwdoubleGAM2,rwdoubleGAMpl,rwdoubleGAMmi); -- auxiliary gamma function calculation for bessjyprocedurebessi0(floatX); -- modified Bessel function I0procedurebessk0(floatX); -- modified Bessel function K0procedurebessi1(floatX); -- modified Bessel function I1procedurebessk1(floatX); -- modified Bessel function K1procedurebessi(integerN,floatX); -- modified Bessel function Inprocedurebessk(integerN,floatX); -- modified Bessel function Knprocedurebessik(floatX,floatXNU,rwfloatRI,rwfloatRK,rwfloatRIP,rwfloatRKP); -- modified Bessel functions of fractional orderprocedurebessj0(floatX); -- Bessel function J0procedurebessy0(floatX); -- Bessel function Y0procedurebessj1(floatX); -- Bessel function J1procedurebessy1(floatX); -- Bessel function Y1procedurebessj(integerN,floatX); -- Bessel function Jnprocedurebessy(integerN,floatX); -- Bessel function Ynprocedurebessjy(floatX,floatXNU,rwfloatRJ,rwfloatRY,rwfloatRJp,rwfloatRYp); -- Bessel functions of fractional orderprocedurebetai(floatA,floatB,floatX); -- incomplete beta functionprocedurebetacf(floatA,floatB,floatX); -- incomplete beta function continued fraction subroutineprocedurebico(integerN,integerK); -- binomial coefficent in floating-point formproceduredawson(floatX); --Dawson's exponential integralprocedureei(floatX); -- exponential integral for n = 1procedureellf(floatPHI,floatAK); --Legendre elliptic integral of the first kindprocedureelle(floatPHI,floatAK); --Legendre elliptic integral of the second kindprocedureellpi(floatPHI,floatEN,floatAK); --Legendre elliptic integral of the third kindprocedureerff(floatX); -- error function (integral of gaussian distribution)procedureerffc(floatX); -- complementary error functionprocedureerfcc(floatX); -- high precision complementary error functionprocedureexpint(integerNR,floatX); -- exponential integralprocedurefactrl(integerN); -- factorial functionprocedurefactln(integerN); -- log of factorial functionprocedurebeta(floatZ,floatW); -- beta functionproceduregammp(floatA,floatX); -- incomplete gamma functionproceduregammq(floatA,floatX); -- complementary incomplete gamma functionproceduregammln(floatXX); -- log of gamma functionproceduregser(rwfloatGAMSER,floatA,floatX,rwfloatGLN); -- retuens incomplete gamma function and its logarithm, -- evaluated using series representationproceduregcf(rwfloatGAMMCF,floatA,floatX,rwfloatGLN); -- incomplete gamma function from continued fraction representationprocedurerf(floatX,floatY,floatZ); -- Carlson's elliptic integral of the first kind -- Carlson's elliptic integral of the second kindprocedurerj(floatX,floatY,floatZ,floatP); -- Carlson's elliptic integral of the third kindprocedurerc(floatX,floatY); -- degenerate Carlson's elliptic integralprocedureairy(floatX,rwfloatAI,rwfloatBI,rwfloatAIP,rwfloatBIP); -- Airy integralprocedurecisi(floatX,rwfloatCI,rwfloatSI); -- sine and cosine integralsprocedureplgndr(integerL,integerM,floatX); -- Legendre polynomials Plm(x)proceduresphbes(integerN,floatX,rwfloatSJ,rwfloatSY,rwfloatSJP,rwfloatSYP); -- spherical Bessel functions Jn, Yn, and their derivativesprocedurefrenel(floatX,rwfloatS,rwfloatC); -- Fresnel integrals S aand Cproceduresncndn(floatUU,floatEMMC,rwfloatSN,rwfloatCN,rwfloatDN); --#################################################################### --# --# CHAPTER 7 - Random Numbers --# --####################################################################procedureran0(rwlongIDUM); -- first random number generator (Park-Miller method)procedureran1(rwlongIDUM); -- second random number generator (Park-Miller method with shuffle)procedureran2(rwlongIDUM); -- third random number generator (l'Ecuyer long-period method with shuffle)procedureran3(rwlongIDUM); -- fourth random number generator (Knuth subtractive method)procedureran4(rwlongIDUM); -- fifth random number generator (Park-Miller method)procedureexpdev(rwlongIDUM); -- exponentially distributed random numbers, averaging 1proceduregasdev(rwlongIDUM); -- Gaussian random numbers with zero mean, deviation 1proceduregamdev(integerIA,rwlongIDUM); -- gamma distributed random numbers, averaging 1procedurepoidev(floatXM,rwlongIDUM); -- poisson distributed random numbers (integers as floats), averaging 1procedurebnldev(floatPP,integerN,rwlongIDUM); -- binomial-distributed random numbersprocedureirbit1(rwulongISEED); -- random bit generatorprocedureirbit2(rwulongISEED); -- random bit generator, second versionprocedurepsdes(rwulongILWORD,rwulongIRWORD); -- random number generator based on DES encryption schemeprocedureranpt(rwfvectPT,fvectREGN); -- returns uniformly random point in specified rectangular regionprocedurerebin(floatRC,fvectR,rwfvectXIN,rwfvectXI); -- rebins a vector of densities into specified new set of bins -- (subroutine of 'vegas')procedurenrsobseq(rwintegerN,rwfvectX); -- generated N-dimensional Sobol random subsequenceprocedurevegas(fvectREGN,ffuncFXN,integerINIT,ulongNCALL,integerITMX, integerNPRN,rwfloatTGRAL,rwfloatSD,rwfloatCHI2A); -- Monte_Carlo integration of function of several variables -- in rectangular domain by vegas adaptive methodproceduremiser(ffuncFUNC,fvectREGN,ulongNPTS,floatDITH,rwfloatAVE,rwfloatVAR); -- Monte_Carlo integration of function of several variables -- in rectangular domain by recursive stratified sampling --#################################################################### --# --# CHAPTER 8 - Sorting --# --####################################################################procedurepiksrt(rwfvectARR); -- simple extraction sortprocedurepiksr2(rwfvectARR,rwfvectBRR); -- simple extraction sort, rearraanging second vectorprocedureshell(rwfvectA); -- shell sortproceduresort(rwfvectARR); -- quicksortproceduresort2(rwfvectARR,rwfvectBRR); -- quicksort, rearranging second vectorproceduresort3(rwfvectRA,rwfvectRB,rwfvectRC); -- rearranges tow additional sequence using sort-permutation of firstprocedurehpsort(rwfvectRA); -- heapsortprocedureindexx(rwfvectARR,rwulvectINDX); -- calculates sorted-position index for an arrayprocedurerank(ulvectINDX,rwulvectIRANK); -- -- calculates sorted-ranks for array elementsprocedurenselect(integerK,rwfvectARR); -- puts k-th largest element in poition k of an arrayprocedureselip(integerK,rwfvectARR); -- returns k-th largest element of an arrayprocedurehpsel(rwfvectARR,rwfvectHEAP); -- puts m-th largest elements of array into a subarrayprocedureeclass(rwivectNF,rwivectLISTA,rwivectLISTB); -- determines element equivalence class from list of equivalencesprocedureeclazz(rwivectNF,ifuncEQUIV); -- determines element equivalence class from pairwise -- equivalence-testing function --#################################################################### --# --# CHAPTER 9 - Root Finding and Nonlinear Sets of Equations --# --####################################################################procedurezbrac(ffuncF,rwfloatX1,rwfloatX2); -- root bracketing by progressive expansion of intervalprocedurertbis(ffuncF,floatX1,floatX2,floatXACC); -- root finding by bisection of interval known to contain rootprocedurertflsp(ffuncF,floatX1,floatX2,floatXACC); -- root finding by method of false positionprocedurertsec(ffuncF,floatX1,floatX2,floatXACC); -- root finding by secant method; secant serves as approximation -- for the true derivativeprocedurezriddr(ffuncF,floatX1,floatX2,floatXACC); -- root finding by Ridder methodprocedurezbrent(ffuncF,floatX1,floatX2,floatTOL); -- root finding by Brent high-speed methodprocedurezrhqr(fvectP,rwfvectR,rwfvectI); -- polynomial root finding by associated matrix computation. -- fvectR and fvectI are real and imaginary -- parts of rootsprocedurezbrak(ffuncF,floatX1,floatX2,fvectXB1,fvectXB2,rwintegerNB); -- preliminary to root finding; decomposes interval containing roots -- into specified number of parts, -- returns lists of stats, ends of intervas containing reversalsprocedureqroot(fvectP,rwfloatB,rwfloatC,floatEPS); -- improves accuracy of trial quadratic factor of polynomialprocedurenrmnewt(integerNTRIAL,rwfvectX,floatTOLX,floatTOLF,vfuncUSRFUN); -- multidimensional root finding/improvement by Newton method; -- requires predefined user function which supplies values and -- Jacobian of vector function to be zeroed.procedurenrtnewt(vfuncD,floatX1,floatX2,floatXACC); -- 1-D root finding by Newton method; -- initial guess is middle of inteval givenprocedurenrtsafe(vfuncD,floatX1,floatX2,floatXACC); -- 1-D root finding by safe variant of Newton method, -- using bisection technique where -- Newton techniqe gives smaller reduction of interval. -- Initial guess is middle of inteval givenprocedurelnsrch(fvectXOLD,floatFOLD,fvectG,fvectP,rwfvectX,rwfloatF,floatSTPMAX,rwintegerCHECK,ffuncFUNC); -- multidimensional root finding by modified Newton method, -- constrained to take steps that actually reduce the norm -- of the target multidimensional functionprocedurenewt(rwfvectX,rwintegerCHECK,vfuncVECFUNC); -- multidimensional root finding by globally convergent variant -- of Newton method (uses lnsrch)procedurefdjac(fvectX,fvectFVEC,rwfmatDF,vfuncVECFUNC); -- coarse numerical estimate of function Jacobian from function valuesprocedurebroydn(rwfvectX,rwintegerCHECKk,vfuncVECFUNC); -- multidimensional root finding by Broyden method; -- uses secant constrain and additional heuristic to -- approximate a multidimensional Newton-like step -- iterative convergence of initial X to a root value of A -- finds all roots of a polynomial A --#################################################################### --# --# CHAPTER 10 - Minimization and Maximization of Functions --# --####################################################################proceduremnbrak(rwfloatAX,rwfloatBX,rwfloatCX,rwfloatFA,rwfloatFB,rwfloatFC,ffuncF); -- downhill search for three points which bracket a minimumproceduregolden(rwfloatAX,rwfloatBX,rwfloatCX,ffuncF,floatTOL,rwfloatXMIN); -- golden-section subdivison search for a function minimumprocedurebrent(rwfloatAX,rwfloatBX,rwfloatCX,ffuncF,floatTOL,rwfloatXMIN); -- Brent parabolic-interpolation search for function minimumproceduredbrent(rwfloatAX,rwfloatBX,rwfloatCX,ffuncF,ffuncDF,floatTOL,rwfloatXMIN); -- Brent parabolic-interpolation search for function minimum, -- using derivativesproceduresimplx(fmatA,integerM1,integerM2,integerM3,rwintegerICASE,rwivectIZROV,rwivectIPOSV); -- simplex method for linear programming. == A is the input linear-programming tableau. -- This is a matrix with rows consisting of: -- (i) coefficients of the linear form to be minimized --(constant term always first) -- (ii) coefficients of the linear forms which must be positive -- (iii) coefficients of the linear forms which must be negative -- (iv) coefficients of the linear forms which must be zero -- The integers M1,M2, and M3 are -- (i) number of linear forms which must be positive -- (ii) number of linear forms which must be negative -- (iii) number of linear forms which must be zero -- The number of rows of A must be M + 2 = 1 + M1 + M2 + M3 + 1 -- (one spare row is used as workspace) -- The number of columns of A must be 1 + (N = number of variables) -- The ICASE integer returned is 0 if a amx was found, -- 1 if the mx is infinite, -1 if the contraints -- are unsatisfiable. -- Aside from this, the output is returned in A, -- in a permuted form represented by two integer matrices: -- TZROV contains M integers, with A(1)(1) = the max; -- A(j + 1)(1) = value at maximum of variable number -- TZROV(j), provided that TZROV(j) <= N; otherwise it is == the val of the TZROV(j) - N +1st slack variable. -- Similarly A(1)(j + 1) = value at maximum of variable -- number IZROV(j), provided that IZROV(j) <= N; -- otherwise it is the val of the IZROV(j) - N +1st slack variableproceduresimp1(fmatA,ivectLL,integerIABF,rwintegerKP,rwfloatBMAX); -- subroutine for 'simplx'. determines maximum, or maximum absolute value, -- of A matrix elements listed in LLproceduresimp2(fmatA,rwivectL2,rwintegerIP,integerKP,rwfloatQ1); -- subroutine for 'simplx'. locates a pivot element in the matrix Aproceduresimp3(fmatA,integerI1,integerK1,integerIP,integerKP); -- subroutine for 'simplx'. -- exchange of right and left-hand tableaux elementsprocedurelinmin(fvectP,fvectXI,rwfloatFRET,ffuncFUNC); -- multidimensional function minimization by downhill line searchproceduredlinmin(fvectP,fvectXI,rwfloatFRET,ffuncFUNC,vfuncDFUNC); -- multidimensional function minimization by -- downhill line search using derivativesproceduref1dim(floatX); -- subroutine for 'linmin'; represents linear section -- of function being minimizedproceduredf1dim(floatX); -- subroutine for 'dlinmin'; represents linear section -- of function being minimizedproceduremetrop(floatDE,floatT);procedureamoeba(rwfmatP,rwfvectY,floatFTOL,ffuncFUNK, --rwintegerNFUNK); -- downhill simplex method for multidimensional function minimizationprocedureamotry(rwfmatP,rwfvectY,rwfvectPSUM,ffuncFUNK,integerIHI,floatFAC); -- extrapolates function through simplex face -- (subroutine for 'amoeba')procedurepowell(fvectP,fmatXI,floatTOL,rwintegerITER,rwfloatFRET,ffuncFUNC); -- Powell quadratically convergent conjugate-direction method -- for multidimensional function minimizationprocedurefrprmn(rwfvectP,floatFTOL,rwintegerITER,rwfloatFRET,ffuncFUNC,vfuncDFUNC); -- conjugate-direction downhill search method for -- multidimensional function minimizationproceduredfpmin(rwfvectP,floatGTOL,rwintegerITER,rwfloatFRET,ffuncFUNC,vfuncDFUNC); -- variable-metric conjugate-direction downhill search method for -- multidimensional function minimizationprocedureamebsa(rwfmatP,rwfvectY,fvectPB,rwfloatYB,floatFTOL,ffuncFUNK,rwintegerITER,floatTEMPTR); -- multidimensional function minimization by simuated annealing -- combined with downhill simplex methodprocedureamotsa(rwfmatP,rwfvectY,fvectPSUM,fvectPB,rwfloatYB,ffuncFUNK,integerIHI,rwfloatYHI,floatFAC); -- function extrapolation through designated face of simplex -- (subroutine for 'amebsa')procedureanneal(fvectX,fvectY,rwivectIORDER); -- approximate solution of traveling salesman problem by simulated annealingprocedurerevcst(fvectX,fvectY,ivectIORDER,ivectN); -- returns cost function for proposed path reversal -- (subroutine for 'anneal')procedurereverse(rwivectIORDER,ivectN); -- performs path reversal (subroutine for 'anneal')proceduretrncst(fvectX,fvectY,ivectIORDER,rwivectN); -- (k is 6) calculates cost function for proposed traveling salesman -- path segment (subroutine for 'anneal')proceduretrnspt(rwivectIORDER,ivectN); -- (k is 6) determines reconfiguration acceptability in -- traveling salesman problem (subroutine for 'anneal') --#################################################################### --# --# CHAPTER 11 - Eigensystems --# --####################################################################procedurejacobi(rwfmatA,rwfvectD,rwfmatV,rwintegerNROT); -- eigenvectors and eigenvalues of real symmetric matrix -- by Jacobi methodproceduretred2(rwfmatA,rwfvectD,rwfvectE); -- reduce real symmetric matrix to tridiagonal form -- by Householder methodproceduretqli(rwfvectD,rwfvectE,rwfmatZ); -- eigenvectors and eigenvalues of real symmetric tridiagonal matrix -- by Jacobi methodprocedurebalanc(rwfmatA); -- replace A by 'balanced' matrix with the same eigenvaluesprocedureelmhes(rwfmatA); -- reduce real matrix to Hessenberg formprocedurehqr(rwfmatA,rwfvectR,rwfvectI); -- eigenvalues of real Hessenberg matrix by QR reduction methodprocedureeigsrt(rwfvectD,rwfmatV); -- sort eigenvectors and eigenvalues into descending order of eigenvalues --#################################################################### --# --# CHAPTER 12 - Fast Fourier Transform --# --####################################################################procedurefour1(rwfvectD,integerISIGN); -- vector fast Fourier transform and inverse. -- put ISIGN = 1 for direct transform; -- -1 for #fvectD/2 times inverse transform. -- dimension of fvectD must be a power of 2 -- The input (and output) of this routine is a -- real representation of a complex vector of -- length 2**n, with the real part of each complex number -- in odd position n and the corresponding imaginary part -- in the next following position n + 1.proceduredfour1(rwdvectD,integerISIGN); -- double precision vector fast Fourier transform and inverse. -- put ISIGN = 1 for direct transform; -- -1 for #fvectD/2 times inverse transform. -- dimension of fvectD must be a power of 2proceduretwofft(fvectD1,fvectD2,rwfvectF1,rwfvectF2); -- fast Fourier transform and inverse, pair of real vectorsprocedurerealft(rwfvectD,integerISIGN); -- fast Fourier transform of a real vectorproceduredrealft(rwdvectD,integerISIGN); -- double precision fast Fourier transform of a real vectorproceduresinft(rwfvectD); -- fast sine transform of real vector. -- This is self inverse, if divided by #fvectD/2 -- the first component of fvectD must be 0procedurecosft1(rwfvectD); -- fast cosine transform of real vector. -- This is self inverse, if divided by #fvectD/2procedurecosft2(rwfvectD,integerISIGN); -- second form of fast cosine transform of real vector.procedurefourn(rwfvectD,ulvectNN,integerISIGN); -- multidimensional fast Fourier transform and inverse.procedurerlft3(rwftensD,rwfmatSPEQ,integerISIGN); -- 3-dimensional fast Fourier transform for a real tensor or array. -- 1-D or multidimensional FFT of large vector stored on external media -- rewinds and renumbers the files used by 'fourfs' --#################################################################### --# --# CHAPTER 13 - Fourier And Spectral Applications --# --####################################################################procedureconvlv(fvectDATA,fvectRESPNS,ulongM,integerISIGN,rwfvectANS); -- convolution of a pair of vectors; -- the RESPNS function must be heled in wrap-around order in the -- first m elements of the n-long RESPNS vector, ANS must be twice as longprocedurecorrel(fvectDATA1,fvectDATA2,rwfvectANS); -- correlation of a pair of vectors, with result returned in ANS, -- which must be twice as long -- FFT-baased power spectrum estimationprocedurewt1(fvectA,integerISIGN,vfuncWTSTEP); -- driver for 1-D wavelet calculationprocedurewtn(fvectA,ulvectNN,integerISIGN,vfuncWTSTEP); -- driver for n-dimensional wavelet calculation; -- k is the product of all the components of NNprocedurespread(floatY,rwfvectYY,floatX,integerM); -- auxiliary routine for 'fasper'proceduredftcor(floatW,floatDELTA,floatA,floatB,fvectENDPTS,rwfloatCORRE,rwfloatCORIM,rwfloatCORFAC); -- calculate correction to Fourier estimate of integralproceduredftint(ffuncFUNC,floatA,floatB,floatW,rwfloatCOSINT,rwfloatSININT); -- sample integration routine illustrating the use of 'dftcor'procedurepwtset(integerN); -- coefficient initializing routine for pwtproceduredaub4(fvectA,integerISIGN); -- Daubechies 4-point wavelet filterprocedurepwt(fvectA,integerISIGN); -- partial wavelet subroutine for 'wt1' aand 'wtn'procedurepredic(fvectDATA,fvectD,rwfvectFUTURE); -- predicts data points using linear prediction coefficientsprocedureevlmem(floatFDT,fvectD,floatXMS); -- generates power-spectrum estimate from 'memcof' outputprocedureperiod(fvectX,fvectY,floatOFAC,floatHIFAC,rwfvectPX,rwfvectPY,rwintegerNOUT,rwintegerJMAX,rwfloatPROB); -- calculates Lomb normalized periodogramprocedurefasper(fvectX,fvectY,floatOFAC,floatHIFAC,rwfvectWK1,rwfvectWK2,rwulongNOUT,rwulongJMAX,rwfloatPROB); -- calculates Lomb normalized periodogram by fast methodprocedurememcof(fvectDATA,rwfloatXMS,rwfvectD); -- estimates linear prediction coefficients from given dataprocedurefixrts(rwfvectD); -- given linear prediction coefficients, finds all roots of the characteristic polynomial --#################################################################### --# --# CHAPTER 14 - Statistical Description of Data --# --####################################################################proceduremoment(fvectDATA,rwfloatAVE,rwfloatADEV,rwfloatSDEV,rwfloatVAR,rwfloatSKEW,rwfloatCURT); -- calculate moments of data vectorprocedurettest(fvectDATA1,fvectDATA2,rwfloatT,rwfloatPROB); -- Student test for significance of difference of meansproceduretutest(fvectDATA1,fvectDATA2,rwfloatT,rwfloatPROB); -- Student test for significance of difference of means -- if variances differprocedureavevar(fvectDATA,rwfloatAVE,rwfloatVAR); --proceduretptest(fvectDATA1,fvectDATA2,rwfloatT,rwfloatPROB); -- Student test for significance of difference of means -- in paired samplesprocedureftest(fvectDATA1,fvectDATA2,rwfloatF,rwfloatPROB); -- F test for significance of difference of variancesprocedurechsone(fvectBINS,fvectEBINS,integerKNSTRN,rwfloatDF,rwfloatCHSQ,rwfloatPROB); -- chisquare test for significance of difference of distribution, -- binned dataprocedurechstwo(fvectBINS1,fvectBINS2,integerKNSTRN,rwfloatDF,rwfloatCHSQ,rwfloatPROB); -- chisquare test for significance of difference of distribution, -- binned data, generalized caseprocedureksone(fvectDATA,ffuncFUNC,rwfloatD,rwfloatPROB); -- Kolmogrov-Smirnov test for significance of difference -- of distribution from assumed model, unbinned dataprocedurekstwo(fvectDATA1,fvectDATA2,rwfloatD,rwfloatPROB); -- Kolmogrov-Smirnov test for significance of difference -- of distribution in two data sets, unbinned dataprocedureprobks(floatALAM); -- Kolmogrov-Smirnov probability function subroutineprocedurecntab1(imatNN,rwfloatCHISQ,rwfloatDF,rwfloatPROB,rwfloatCRAMRV,rwfloatCCC); -- contingency table analysis for 2-D contingency tableprocedurecntab2(imatNN,rwfloatH,rwfloatHX,rwfloatHY,rwfloatHYGX,rwfloatHXGY,rwfloatUYGX,rwfloatUXGY,rwfloatUXY); -- entropy-based contingency table analysis -- for 2-D contingency tableprocedurepearsn(fvectX,fvectY,rwfloatR,rwfloatPROB,rwfloatZ); -- Pearson analysis for linear dependence of -- two associated data setsprocedurespear(fvectDATA1,fvectDATA2,rwfloatD,rwfloatZD,rwfloatPROBD,rwfloatRS,rwfloatPROBRS); -- Spearman rank-order correlation coefficientprocedurecrank(fvectW,rwfloatS); -- ranking subroutine for Spearman rank-order correlationprocedurekendl1(fvectDATA1,fvectDATA2,rwfloatTAU,rwfloatZ,rwfloatPROB); -- Kendall rank-order tau coefficientprocedurekendl2(fmatTAB,rwfloatTAU,rwfloatZ,rwfloatPROB); -- contingency-table based variant of Kendall rank-order -- tau coefficientprocedurenrks2d1s(fvectX1,fvectY1,vfuncQUADVL,rwfloatD1,rwfloatPROB); -- Kolmogrov-Smirnov test for significance of difference of -- distribution from assumed model, 2-D caseprocedureks2d2s(fvectX1,fvectY1,fvectX2,fvectY2,rwfloatD,rwfloatPROB); -- Kolmogrov-Smirnov test for significance of difference -- of distribution of two samples, 2-D caseprocedurequadct(floatX,floatY,fvectXX,fvectYY,rwfloatFA,rwfloatFB,rwfloatFC,rwfloatFD); -- cout of points in each of 4 quadrantsprocedurequadvl(floatX,floatY,rwfloatFA,rwfloatFB,rwfloatFC,rwfloatFD); -- sample model for 'ks2ds': uniform distribution inside a squareproceduresavgol(rwfvectY1,integerNL,integerNR,integerLD,integerM); -- calculates Savitsky-Golay filter coefficents for stated -- number of filter points to be used --#################################################################### --# --# CHAPTER 15 - Modeling of Data --# --####################################################################procedurefit(fvectX,fvectY,fvectSIG,integerMWT,rwfloatA,rwfloatB,rwfloatSIGA,rwfloatSIGB,rwfloatCHI2,rwfloatQ); -- computes maximum-likelihood fit to data -- with specified variancesprocedurefitexy(fvectX,fvectY,fvectSIGX,fvectSIGY,rwfloatA,rwfloatB,rwfloatSIGA,rwfloatSIGB,rwfloatCHI2,rwfloatQ); -- computes maximum-likelihood fit to data with -- specified variances, and variablility in both x and yprocedurechixy(doubleBANG); -- computes ch**2 coefficient as function of offset -- (subroutine of 'fitexy')procedurelfit(fvectX,fvectY,fvectSIG,rwfvectA,ivectIA,rwfmatCOVAR,rwfloatCHISQ,vfuncFUNCS); -- best least-squares fit to data with specified variancesprocedurecovsrt(rwfmatCOVAR,ivectIA,integerMFIT); -- spreads calculated covariance back into 'covar' -- (matrix (subroutine of 'lfit')proceduresvdfit(fvectX,fvectY,fvectSIG,rwfvectA,rwfmatU,rwfmatV,rwfvectW,rwfloatCHISQ,vfuncFUNCS); -- fit a specified set of nolinear functions to given data -- using the singular value decomposition techniqueproceduresvdvar(rwfmatV,rwfvectW,rwfmatCVM); -- evaluate the covarince of a fir of a model to data -- (subroutine of 'svdfit') -- procedure mrqmin(fvectX,fvectY,rwfvectSIG,rwfvectA,ivectIA, --rwfmatCOVAR,rwfmatALPHA,rwfloatCHISQ,vfuncFUNCS,rwfloatALAMDA); -- fit a specified set of nolinear functions to given data -- using the Levenberg-Marquardt method -- procedure mrqcof(fvectX,fvectY,rwfvectSIG,rwfvectA,ivectIA,rwfmatALPHA,rwfvectBETA,rwfloatCHISQ,vfuncFUNCS); -- estimate parameters for "mrqmin' -- procedure medfit(fvectX,fvectY,rwfloatA,rwfloatB,rwfloatABDEV); -- fit a line to given data by minimizing absolute devaiton -- procedure rofunc(floatB); -- estimate parameter for "medfit'procedurefgauss(floatX,fvectA,rwfloatY,rwfvectDYDA); -- calculate sum of gaussians with given centers and deviationsprocedurefleg(floatX,rwfvectPL); -- evaluate a sum of Legendre polynomials (example for svdfit)procedurefpoly(floatX,rwfvectP); -- evaluate a polynomial with given coefficients (example for svdfit) --#################################################################### --# --# CHAPTER 16 - Integration of Ordinary Differential Equations --# --####################################################################procedurerk4(fvectY,fvectDYDX,floatX,floatH,rwfvectYOUT,vfuncDERIVS); -- fourth order Runge-Kutta single step of advance; -- vfuncDERIVS returns the (vector) right-hand side -- of the differential equationprocedurerkdumb2(fvectVSTART,floatX1,floatX2,vfuncDERIVS,rwfmatY,rwfvectXX); -- fourth order Runge-Kutta advance from X1 to X2 -- with fixed stepsizeprocedurerkqs(rwfvectY,fvectDYDX,rwfloatX, floatHTRY,floatEPS,fvectYSCAL,rwfloatHDID,rwfloatHNEXT,vfuncDERIVS); -- fifth order Runge-Kutta single step of advanceprocedurerkck(fvectY,fvectDYDX,floatX,floatH,rwfvectYOUT,rwfvectYERR,vfuncDERIVS); -- fifth order Runge-Kutta advance from X1 to X2 -- with adaptive stepsizeprocedureode_solve(fvectYSTART,floatX1,floatX2,floatEPS,floatH1,floatHMIN,rwintegerNOK,rwintegerNBAD,vfuncDERIVS,rwintegerKOUNT,floatDXSAV,rwfvectXP,rwfmatYP); -- top-level ODE solution driver using rkck as individual step. -- solution is constructed from X1 to X2 with accuracy eps. H1 is -- guessed first sstepsize, HMIN is smallest acceptable stepsize. -- NOK is number of good steps taken, NBAD is number of retried steps.procedureode_solve_bs(fvectYSTART,floatX1,floatX2,floatEPS,floatH1, floatHMIN,rwintegerNOK,rwintegerNBAD,vfuncDERIVS,rwintegerKOUNT,floatDXSAV,rwfvectXP,rwfmatYP); -- top-level ODE solution driver using bsstep as individual step.procedureode_solve_stiff_bs(fvectYSTART,floatX1,floatX2,floatEPS,floatH1,rwintegerNOK,rwintegerNBAD,vfuncDERIVS,rwintegerKOUNT,floatDXSAV, floatHMIN,rwfvectXP,rwfmatYP,vfuncJACOBN); -- top-level ODE solution driver using stifbs as individual step.procedureode_solve_stiff(fvectYSTART,floatX1,floatX2,floatEPS,floatH1,floatHMIN,rwintegerNOK,rwintegerNBAD,vfuncDERIVS,rwintegerKOUNT,floatDXSAV,rwfvectXP,rwfmatYP,vfuncJACOBN); -- top-level ODE solution driver using stiff as individual step.proceduremmid(fvectY,fvectDYDX,floatXS,floatHTOT,integerNSTEP,rwfvectYOUT,vfuncDERIVS); -- single step of modified midpoint method for ODEsprocedurebsstep(rwfvectY,fvectDYDX,rwfloatXX,floatHTRY, floatEPS,fvectSCAL,rwfloatHDID,rwfloatHNEXT,vfuncDERIVS); -- single Burlisch-Stoer step for Romberg method for ODEsprocedurepzextr(integerIEST,floatXEST,fvectYEST,rwfvectYZ,rwfvectDY); -- polynomial extrapolation routine for Romberg method for ODEsprocedurerzextr(integerIEST,floatXEST,fvectYEST,rwfvectYZ,rwfvectDY); -- rational extrapolation routine for Romberg method for ODEsprocedurestoerm(fvectY,fvectD2Y,floatXS,floatHTOT,integerNSTEP,rwfvectYOUT,vfuncDERIVS); -- Stoerm method for second order ODEsprocedurestiff(rwfvectY,fvectDYDX,rwfloatX,floatHTRY,floatEPS,rwfvectYSCAL,rwfloatHDID,rwfloatHNEXT,vfuncDERIVS); -- fourth order Roesenbrock step for stiff ODEsprocedurejacobn(floatX,fvectY,fvectDYDX,rwfmatDFDY); -- sample Jacobian matrix routine for 'stiff'proceduresimpr(fvectY,fvectDYDX,fvectDFDX,fmatDFDY,floatXS,floatHTOT, integerNSTEP,rwfvectYOUT,vfuncDERIVS); -- fourth order Roesenbrock step for stiff ODEsprocedurestifbs(rwfvectY,fvectDYDX,rwfloatXX, floatHTRY,floatEPSR,fvectYSCAL,rwfloatHDID,rwfloatHNEXT,vfuncDERIVS); -- semi-implicit exrapolation for integrating stiff systems -- of differential equations --#################################################################### --# --# CHAPTER 17 - Two point boundary value problems --# --####################################################################proceduresetshoot(vfuncLOAD,vfuncSCORE,vfuncDERIVS); -- Initialization function for SHOOTproceduresetshootf(vfuncLOAD1,vfuncLOAD2,vfuncSCORE,vfuncDERIVS,integerNN2, integerNVAR,floatX1,floatX2,floatXF); -- Initialization function for SHOOTFprocedurenr_getcb_shoot();procedureshoot(fvectV,rwfvectF); -- solution of two-point boundary value problem by 'shooting' methodprocedurenr_getcb_shootf();procedureshootf(fvectV,rwfvectF); -- solve 2-point boundary value problem by the 'shooting -- to a fitting point' methodproceduresolvde(integerITMAX,floatCONV,floatSLOWC,fvectSCALV, ivectINDEXV,integerNB,rwfmatY,rwftensC,rwfmatS); -- top level routine for solution of 2-point boundary value -- problem by relaxation methodprocedurepinvs(integerIE1,integerIE2,integerJE1,integerJSF,integerJC1,integerK,rwftensC,rwfmatS); -- diagonalize square subsection of S matrix (subroutine for 'solvde')procedurered(integerIZ1,integerIZ2,integerJZ1,integerJZ2,integerJM1, integerJM2,integerJMF,integerIC1,integerJC1,integerJCF,integerKC,rwftensC,rwfmatS); -- reduce columns of the S matrix using data stored in the C matrix -- (subroutine for 'solvde') -- compute S matrix for 'solvde'procedurebksub(integerNE,integerNB,integerJF,integerK1,integerK2,rwftensC); -- backsubstitution routine for 'solvde' --#################################################################### --# --# CHAPTER 18 - Integral Equations and Inverse Theory --# --#################################################################### -- computes weights for n-point equal-integral quadrature -- for specified function whose moments are returned by 'kermom'procedurekermom(rwdvectWGHTS,floatY); -- computes moments of diagonally singular kernelprocedurequadmx(rwfmatA); -- constructs quadrature matrix for sample Fredholm equation -- of the second kind with singular kernelprocedurefred2(floatA,floatB,rwfvectT,rwfvectF,rwfvectW,ffuncG,ffuncAK); -- solution of linear Fredholm equation of the second kindprocedurefredin(floatX,floatA,floatB,fvectT,fvectF,fvectW,ffuncG,ffuncAK); -- solution value for linear Fredholm equation -- of the second kind -- using 'fred2' output and Nystrom interpolationprocedurevoltra(floatT0,floatH,rwfvectT,rwfmatF,ffuncG,ffuncAK); -- solution of system of linear Fredholm equations -- of the second kind --#################################################################### --# --# CHAPTER 19 - Partial Differential Equations --# --####################################################################proceduresor(dmatA,dmatB,dmatC,dmatD,dmatE,dmatF,rwdmatU,doubleRJAC); -- sucessive over-relaxation solution of boundary value problemproceduremglin(rwdmatU,integerNCYCLE); -- multigrid solution of linear elliptic partial -- differential equationsprocedurerstrct(dmatUC,rwdmatUF); -- transfers developing PDE solution from fine grid -- to coarse grid (subroutine for 'mglin') -- coarse-to fine solution by bilinear interpolation -- (subroutine for 'mglin')procedureaddint(rwdmatUF,dmatUC,rwdmatRES); -- coarse-to fine solution by bilinear interpolation, -- with transfer to fine grid (subroutine for 'mglin')procedureanorm2(dmatA); -- calculates Euclidean norm of matrixprocedurerelax(rwdmatU,dmatRHS); -- red-black Gauss-Seidel relaxation step -- (subroutine for 'mglin')procedureresid(rwdmatRES,dmatU,dmatRHS); -- calcuates residual for use in 'mglin' example -- (subroutine for 'mglin')procedurecopy(rwdmatAOUT,dmatAIN); -- copies matrix to matrix (subroutine for 'mglin')procedurefill0(rwdmatU); -- fills matrix with zeroes (subroutine for 'mglin')proceduremgfas(rwdmatU,integerMAXCYC); -- multigrid solution of nonlinear elliptic -- partial differential equationsprocedurerelax2(rwdmatU,dmatRHS); -- red-black nonlinear Gauss-Seidel relaxation step -- (subroutine for 'mgfas')procedurematadd(dmatA,dmatB,rwdmatC); -- adds two square matricesprocedurematsub(dmatA,dmatB,rwdmatC); -- subtracts two square matricesprocedurelop(rwdmatOUT,dmatU); -- calculates numerical approximation for left-hand side -- of nonlinear PDE (subroutine for 'mgfas')procedureslvsm2(rwdmatU,rwdmatRHS); -- solves nonlinear PDE on coarsest grid -- (n = 3; subroutine for 'mgfas) --#################################################################### --# --# CHAPTER 20 - Less-Numerical Algorithms --# --####################################################################proceduremachar(rwintegerIBETA,rwintegerIT,rwintegerIRND,rwintegerNGRD,rwintegerMACHEP,rwintegerNEGEP,rwintegerIEXP,rwintegerMINEXP,rwintegerMAXEXP,rwfloatEPS,rwfloatEPSNEG,rwfloatXMIN,rwfloatXMAX); -- determine machine float parameters: radix, mantissa size, -- smallest positive nonzero float, float precision, -- smallest negative nonzero float, etc.procedureigray(ulongN,integerIS); -- return the Gray code of the integer n -- compute 16-bit circular redundancy check -- for an array of integersprocedureicrc1(ushortCRC,ushortONEC); -- compute 16-bit circular redundancy check -- for one more item added to previously checked array -- compute and verify check character for an array of charactersprocedurempadd(rwcvectW,cvectU,cvectV); -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer add -- multiprecision integer addprocedurecaldat(ulongJULIAN,rwintegerMM,rwintegerID,rwintegerIYYY); -- convert Julian day number to calendar dateprocedureflmoon(integerN,integerNPH,rwlongJD,rwfloatFRAC); -- numerical date of n-th full moon since 1900procedurejulday(integerMM,integerID,integerIYYY); -- convert calendar date to Julian day numberendlibnr;

**Basic Procedures for Matrices and Vectors.** The first group of routines in the library, nr_get_tensor, nr_get_matrix, nr_get_imatrix, nr_get_vector, nr_get_imatrix, nr_get_vector, nr_get_dvector, nr_get_ivector, and nr_get_lvectorcrate this objects from the corresponding SETL tuples, tuples-of-tuples (for matrices), and tuples-of-tuples-of-tuples (for tensors). Two auxiliary routines, nr_get_element and nr_print serve respectively to extract specified elements from a vector, matrix, or tensor and to print vectors, matrices, and tensors. A third group, nr_mat_sum, nr_mat_diff, nr_mat_prod, nr_vect_prod, . They have analogs nr_matconst_sum, nr_matconst_diff, nr_matconst_prod, which . The related routine nr_vect_prod . nr_mat_clone. Note that these routines can be applied to .

The use of this first group of procedures is illustrated by the

Thi program produces the following outputprogramtest; -- numerical library example 1uselibnr;float(x): xin[1..10]])); -- form floating vectorin[1..10]])); -- form integer vectorendtest;

Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 Type = 11 Vector Dimensions = 10 2.000000 4.000000 6.000000 8.000000 10.000000 12.000000 14.000000 16.000000 18.000000 20.000000 Type = 11 Vector Dimensions = 10 2.000000 4.000000 6.000000 8.000000 10.000000 12.000000 14.000000 16.000000 18.000000 20.000000 Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 Type = 11 Vector Dimensions = 10 1.000000 4.000000 9.000000 16.000000 25.000000 36.000000 49.000000 64.000000 81.000000 100.000000 Type = 11 Vector Dimensions = 10 101.000000 102.000000 103.000000 104.000000 105.000000 106.000000 107.000000 108.000000 109.000000 110.000000 Type = 11 Vector Dimensions = 10 -99.000000 -98.000000 -97.000000 -96.000000 -95.000000 -94.000000 -93.000000 -92.000000 -91.000000 -90.000000 Type = 11 Vector Dimensions = 10 100.000000 200.000000 300.000000 400.000000 500.000000 600.000000 700.000000 800.000000 900.000000 1000.000000 Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 Type = 13 Vector Dimensions = 10 1 2 3 4 5 6 7 8 9 10 Type = 13 Vector Dimensions = 10 2 4 6 8 10 12 14 16 18 20 Type = 13 Vector Dimensions = 10 2 4 6 8 10 12 14 16 18 20 Type = 13 Vector Dimensions = 10 1 2 3 4 5 6 7 8 9 10 Type = 13 Vector Dimensions = 10 1 4 9 16 25 36 49 64 81 100 Type = 13 Vector Dimensions = 10 101 102 103 104 105 106 107 108 109 110 Type = 13 Vector Dimensions = 10 -99 -98 -97 -96 -95 -94 -93 -92 -91 -90 Type = 13 Vector Dimensions = 10 100 200 300 400 500 600 700 800 900 1000 Type = 13 Vector Dimensions = 10 1 2 3 4 5 6 7 8 9 10The routine nr_mat_clone makes a copy of its argument vector (or matrix). nr_mat_fix makes a copy of its argument vector (or matrix), which must be a 'libnr' floating vector (or matrix), but also converts this to a 'libnr' integer vector (or matrix). nr_mat_float makes a copy of its argument vector (or matrix), which must be a 'libnr' integer vector (or matrix), but also converts this to a 'libnr' floating vector (or matrix). nr_get_object . nr_rows_and_cols

This program produces the following output.programtest; -- numerical library example 2uselibnr;float(x): xin[1..10]]))); -- form floating vectorendtest;

[1.0000000000, 2.0000000000, 3.0000000000, 4.0000000000, 5.0000000000, 6.0000000000, 7.0000000000, 8.0000000000, 9.0000000000, 10.000000000] Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 Type = 13 Vector Dimensions = 10 1 2 3 4 5 6 7 8 9 10 Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 [10, 1]The following variant of the preceding program

The output of this last program isprogramtest; -- numerical library example 3uselibnr;float(x * y): xin[1..3]]: yin[1..3]]))); -- form 3 by 3 floating matrixendtest;

[[1.0000000000, 2.0000000000, 3.0000000000], [2.0000000000, 4.0000000000, 6.0000000000], [3.0000000000, 6.0000000000, 9.0000000000]] Type = 1 Matrix Dimensions = 3 x 3 1.000000 2.000000 3.000000 2.000000 4.000000 6.000000 3.000000 6.000000 9.000000 Type = 3 Matrix Dimensions = 3 x 3 1 2 3 2 4 6 3 6 9 Type = 1 Matrix Dimensions = 3 x 3 1.000000 2.000000 3.000000 2.000000 4.000000 6.000000 3.000000 6.000000 9.000000 [3, 3]nr_get_zero_vector forms a zero vector of specified length which can either be of floating, double precision, or integer type depending the second argument supplied to nr_get_zero_vector. nr_get_slice extracts a specified range of componenets from a vector or matrix. If the first argument of nr_get_slice is a vector, then its second argument should be a pair of integers defining the range of components to be extracted. However, if the first argument of nr_get_slice is a matrix, then its second argument should be a list of four integers i, j, k, l defining a submatrix to be extracted. These integers should be in the order left most column, top most row, right most column, bottom row. nr_set_vect_slice overwrites the components of one vector with those of another starting at a specified position. nr_vect_reverse returns the reverse of a specified vector. nr_vect_max_index returns the index of the largest component in a vector. nr_mat_transpose returns the transpose of a two by two matrix. The following example shows the use of many of these routines, both for vectors and for 2-D matrices.

This code produces the following output.programtest; -- numerical library example 4uselibnr;float(x): xin[1..10]])); -- form floating vectorfloat(x * y * y): xin[1..3]]: yin[1..3]]))); -- form 3 by 3 floating matrixendtest;

Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 Type = 11 Vector Dimensions = 3 2.000000 3.000000 4.000000 Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 2.000000 3.000000 4.000000 8.000000 9.000000 10.000000 Type = 11 Vector Dimensions = 10 10.000000 9.000000 8.000000 4.000000 3.000000 2.000000 4.000000 3.000000 2.000000 1.000000 10 [[1.0000000000, 2.0000000000, 3.0000000000], [4.0000000000, 8.0000000000, 12.000000000], [9.0000000000, 18.000000000, 27.000000000]] Type = 1 Matrix Dimensions = 2 x 2 1.000000 2.000000 4.000000 8.000000 Type = 1 Matrix Dimensions = 3 x 3 1.000000 4.000000 9.000000 2.000000 8.000000 18.000000 3.000000 12.000000 27.000000Some numerical routines need to work with complex rather than real vectors, and these are often represented by real vectors in which every odd numbered component represents the real part of a complex number and the following component represents the imaginary part of the same number. The two routines nr_cvect_conj and nr_cvect_prod are set up to work with vectors of this type. nr_cvect_conj forms the complex conjugate of such a vector by reversing all of its even numbered components, and nr_cvect_prod multiplies two complex vectors having this representation. The following short example shows the use of thse operations.

programtest; -- numerical library example 5uselibnr;float(x): xin[1..10]])); -- form floating vectorendtest;

Type = 11 Vector Dimensions = 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 Type = 11 Vector Dimensions = 10 1.000000 -2.000000 3.000000 -4.000000 5.000000 -6.000000 7.000000 -8.000000 9.000000 -10.000000 Type = 11 Vector Dimensions = 10 NaN 0.000000 0.000000 0.000000 0.000000 NaN NaN NaN NaN NaNThe very important procedure nr_mat_prod mutiplies to real matrices and returns their product. The two matrices to be multiplied must be of mutually consistent size, that is, the number of rows of the first matrix must equal the number of columns of the second matric.

The example shown produces the following output.programtest; -- numerical library example 6uselibnr;float(x * y): xin[1..3]]: yin[1..3]])); -- form 3 by 3 floating matrixin[1..3]]: yin[1..3]])); -- form 3 by 3 floating matrixendtest;

Type = 1 Matrix Dimensions = 3 x 3 1.000000 2.000000 3.000000 2.000000 4.000000 6.000000 3.000000 6.000000 9.000000 Type = 1 Matrix Dimensions = 3 x 3 14.000000 28.000000 42.000000 28.000000 56.000000 84.000000 42.000000 84.000000 126.000000 Type = 3 Matrix Dimensions = 3 x 3 1 2 3 2 4 6 3 6 9 Type = 3 Matrix Dimensions = 3 x 3 14 28 42 28 56 84 42 84 126

programtest; -- numerical library example 7: Gauss numerical library example 7: Gauss-Jordan solution of linear equationsuselibnr;cos(float(x * y)): xin[1..3]]: yin[1..3]])); -- form 3 by 3 floating nonsingular matrixfloat(y)]: yin[1..3]])); -- form floating vector as 1 by 3 matrix gaussj(m,v);float(y): xin[1..2]]: yin[1..3]])); -- form floating 2 by 3 matrix gaussj(m,v);endtest;

Type = 1 Matrix Dimensions = 3 x 3 0.540302 -0.416147 -0.989992 -0.416147 -0.653644 0.960170 -0.989992 0.960170 -0.911130 Type = 1 Matrix Dimensions = 3 x 1 1.000000 2.000000 3.000000 Type = 1 Matrix Dimensions = 3 x 3 -0.230942 -0.940919 -0.740630 -0.940919 -1.041855 -0.075572 -0.740630 -0.075572 -0.372442 Type = 1 Matrix Dimensions = 3 x 1 -4.334671 -3.251344 -2.009102 Type = 1 Matrix Dimensions = 3 x 2 1.000000 1.000000 2.000000 2.000000 3.000000 3.000000 Type = 1 Matrix Dimensions = 3 x 3 0.540302 -0.416147 -0.989992 -0.416147 -0.653644 0.960170 -0.989992 0.960170 -0.911130 Type = 1 Matrix Dimensions = 3 x 2 -3.261969 -3.261969 1.157077 1.157077 -1.803043 -1.803043

This code produces the following output.programtest; -- numerical library example 3.1: interpolation proceduresuselibnr; polint(xvo := nr_get_vector(xv := [1.0,1.6,2.0]),yvo := nr_get_vector([cos(x): xinxv]),x := 1.1,y,dy);cos(x)," ",-sin(x)); -- print interpolated values ratint(xvo,yvo,x,y,dy);cos(x)," ",-sin(x)); -- print interpolated values splint(xvo,yvo,yvo,x,y);cos(x)); -- print interpolated valuesendtest;

polint: 0.44629520178 0.00090993667254 0.45359612143 -0.89120736006 ratint: 0.44627976418 0.68634641171 0.45359612143 -0.89120736006 splint: 0.43741458654 0.45359612143Numerical procedures for recovering the coefficients of a polynomial given sufficiently many of its values are also available. These two routines polcoe and polcof have the same calling sequence but polcof uses a somewhat more accurate technique. These are seen in the following example.

This example produces the following output.programtest; -- numerical library example 3.3: polynomial coefficents from polynomial valuesuselibnr; polcoe(xvo := nr_get_vector(xv := [float(x): xin[1..5]]),pvo := nr_get_vector([polyvals(x): xinxv]),coeffs);procedurepolyvals(x); -- evaluate test polynomial with specified coefficients coeffs := [1.0,2.0,3.0,4.0,5.0]; sum := 0.0;forcincoeffs(1..(nc := #coeffs) - 1)loopsum := x * (sum + c);end loop;returnsum + coeffs(nc);endpolyvals;endtest;

recovered coeffs: 3.0000000000 Type = 11 Vector Dimensions = 5 5.000000 4.000000 3.000000 2.000000 1.000000 recovered coeffs: 3.0000000000 Type = 11 Vector Dimensions = 5 5.000000 4.000000 3.000000 2.000000 1.000000

**Special Functions.** This section of LIBNR_PAK provies means of calculating the most important 'special functions' of mathematical physics. The functions included (all of which are comprehenively documented and analyzed in the treatise by Press, Teukolsky, Vetterling, and Flannery) are as follows.

- Bessel functions: bessi0, bessk0, bessi1, bessk1, bessi, bessk, bessik, bessj0, bessy0, bessj1, bessy1, bessj, bessy, bessjy, sphbes.
- Exponential integrals, beta, and gamma functions: betai, bico, expint, factrl, factln, beta, gammp, gammq, gammln, dawson, ei.
- Elliptic integrals: ellf, elle, ellpi, rf, rj, rc, sncndn.
- Error functions: erff, erffc, erfcc.
- Legendre polynomials: plgndr.
- Miscellaneous special integrals: airy, cisi, frenel.

The following example illstrates the use of these functions.

programtest; -- numerical library example 6.1: special functionsuselibnr;endtest;

Randomization is often a very valuable numeric technique and when used in this way consumes random numbers at a great rate and can be quite sensitive to their quality, that is, the degree to which they behave in truly random fashion. For this reason libnr_pak provides numerous random number generators. It also provides several modified random number generators which generate numbers satisfying various important statistical distribution laws, for example, conform to Gauss 'normal distribution' or are exponentially or binomially distributed, etc. These procedures all appear in the following example.

The two routines 'vegas' and 'miser'programtest; -- numerical library example 7.1: random numerical library example 7.1: random-number generationuselibnr; tup := [ran := 666];forjin[1..4]loopran0(ran); tupwith:= ran;end loop;forjin[1..4]loopran1(ran); tupwith:= ran;end loop;forjin[1..4]loopran2(ran); tupwith:= ran;end loop;forjin[1..4]loopran3(ran); tupwith:= ran;end loop;forjin[1..4]loopran4(ran); tupwith:= ran;end loop;forjin[1..4]looppsdes(rana,ranb); tup +:= [rana,ranb];end loop;in[1..5]];in[1..5]];in[1..5]];in[1..5]];in[1..5]];forjin[1..4]loopranpt(pt,cube); tupwith:= nr_get_object(pt);end loop;forjin[1..4]loopnrsobseq(n,pt); tupwith:= nr_get_object(pt);end loop;str(irbit1(ran)): jin[1..40]];str(irbit2(ran)): jin[1..40]];endtest;

programtest; -- numerical library example 7.1: random numerical library example 7.1: random-number generationuselibnr; cube := nr_get_vector([1.0,1.0]); miser(cosxy,cube,50,0.1,aver,vrnc);procedurecosxy(x,y);return cos(x) *cos(y);endcosxy;endtest;

**Root Finding for Single Equations and for Sets of Equations.** The routines in this group find the roots of general smooth functions in one and several variables and of polynomials (for which ). These routines fall into several groups:

- Derivative-independent methods, of increasing sophistication, for smooth functions of one variable: rtbis, rtflsp, rtsec, zriddr, zbrent, and the auxillary root-bracketing routines zbrac.
- Newton methods for functions of one variable: nrtnewt, nrtsafe.
- Multidimensional methods: .
- Polynomial root-finding routines: .

programtest; -- numerical library example 8.1: 1 numerical library example 8.1: 1-dimensional root-findinguselibnr; lo := -0.1; hi := 0.1; zbrac(cos,lo,hi);cos,lo,hi,0.1e-8)); -- find root by bisection methodcos,lo,hi,0.1e-8)); -- find root by false position methodcos,lo,hi,0.1e-6)); -- find root by secant method methodcos,lo,hi,0.1e-8)); -- find root by Ridder methodcos,lo,hi,0.1e-8)); -- find root by Brent high find root by Brent high-speed method --cos,lo,hi,0.1e-8)); -- find root by Newton method --cos,lo,hi,0.1e-8)); -- find root by safe Newton method zrhqr(nr_get_vector(5 * [1.0]),roots_r,roots_i); -- complex roots of polynomial by matrix methodprocedurecos_with_der(x); d := -sin(x);return[cos(x),-sin(x)];endcos_with_der;endtest;

**Minimization and Maximization of Functions of one or More Variables.** The routines in this group fall into the following subgroups:

- Derivative-independent methods for smooth functions of one variable: .
- Simplex method for linear programming: .
- Multidimensional function minimization, with and without derivatives: .
- Conjugate-direction downhill search method: .

programtest; -- numerical library example 8.1: 1 numerical library example 8.1: 1-dimensional minimizationuselibnr; loval :=cos(lo := -0.1); hival :=cos(hi := 0.1); midval :=cos(mid := 0.0); mnbrak(lo,mid,hi,loval,midval,hival,cos);cos,0.001,the_min);cos,0.001,the_min);cos,lambda(x);return-sin(x);end lambda,0.001,the_min);endtest;

**Eigensystems.** These very important and much-studied routines are central to numerical treatment of all sorts of linear problems.

**Fast Fourier Transform.**

**Wavelets; Fourier And Spectral Applications.**

**Statistical Description of Data.**

**Modeling of Data.**

**Integration of Ordinary Differential Equations.**

**Two point boundary value problems.**

**Integral Equations and Inverse Theory.**

**Partial Differential Equations.**

**Miscellaneous Auxiliary Procedures.**

These same moduli also control the way in which ordinary SETL objects convert to opaque NTL objects.

Polynomials are represented by their vectors of coeficents, lowest-order coefficient first. Thus a polynomial can be printed as a tuple with leading, but not trailing, zeroes. However, a polynomial with zero leading coefficients reatis its full length, so that, for example, two polynomials which are identical except for zero leading coeficients may have fifferent reverses

The following header file lists the routines available. A few simple examples illustrating their use are then given.

Modulo 2 cases are treated in a special way for efficiency.

native packagelibntl; -- conversion routines; integer and modular integer conversionsprocedurentl_get_integer(setl_integer); -- convert a setl integer to an ntl integerprocedurentl_get_integer_modulo(setl_integer); -- convert a setl integer to an ntl 'modulo p' object -- the 'ntl_set_modulo' routine defines the current modulusprocedurentl_get_integer_poly(setl_tuple); -- convert a setl tuple to an integer polynomialprocedurentl_get_integer_poly_modulo(setl_tuple); -- convert a setl tuple to a modular polynomialprocedurentl_get_integer_matrix(setl_tuple); -- convert a setl tuple of tuples to a matrixprocedurentl_get_integer_modulo_matrix(setl_tuple); -- convert a setl tuple of tuples to a matrixprocedurentl_get_integer_vector(setl_tuple); -- convert a setl tuple to a vectorprocedurentl_get_integer_modulo_vector(setl_tuple); -- convert a setl tuple to a modular vector -- conversion routines; real conversionsprocedurentl_get_real(setl_real); -- convert a setl real to an ntl objectprocedurentl_get_real_matrix(setl_tuple); -- convert a setl real matrix to an ntl objectprocedurentl_get_real_vector(setl_tuple); -- convert a setl real tuple to an ntl object -- conversion routines; poly mod poly conversionsprocedurentl_get_modextint(setl_tuple); -- convert a setl tuple to a modulo p, mod poly quantity -- the current polynomial modulus 'poly' is defined by the -- 'ntl_set_modulo_ext' routineprocedurentl_get_modextint_vector(setl_tuple); -- convert a setl tuple of tuples to a vector with -- modulo p, mod poly components. the current polynomial -- modulus 'poly' is defined by the'ntl_set_modulo_ext' routineprocedurentl_get_modextint_poly(setl_tuple); -- convert a setl tuple of tuples to a polynomial with -- modulo p, mod poly coefficients. the current polynomial -- modulus 'poly' is defined by the'ntl_set_modulo_ext' routineprocedurentl_get_modextint_poly_vector(setl_tuple); -- convert a setl tuple of tuples of tuples to a vector whose -- elements are polynomials with modulo p, mod poly -- coefficients. the current polynomial modulus -- 'poly' is defined by the'ntl_set_modulo_ext' routine -- conversion routines; modulo 2 casesprocedurentl_get_mod2int(setl_short_integer); -- convert a setl integer to an ntl 'modulo 2' objectprocedurentl_get_mod2int_poly(setl_tuple); -- convert a setl tuple to a modulo 2 polynomialprocedurentl_get_mod2int_vector(setl_tuple); -- convert a setl tuple of tuples to a vector with -- modulo 2, mod poly elements. the current polynomial -- modulus 'poly' is defined by the'ntl_set_modulo_ext' routineprocedurentl_get_mod2int_matrix(setl_tuple); -- convert a setl tuple of tuples of tuples to a matrix with -- modulo 2 elements. the current polynomial modulus -- 'poly' is defined by the'ntl_set_modulo_ext' routineprocedurentl_get_mod2extint(setl_tuple); -- convert a setl tuple to a modulo 2, mod poly quantity -- the current polynomial modulus 'poly' is defined by the -- 'ntl_set_modulo_ext' routineprocedurentl_get_mod2extint_poly(setl_tuple); -- convert a setl tuple of tuples to a polynomial with -- modulo 2, mod poly coefficients. the current polynomial -- modulus 'poly' is defined by the'ntl_set_modulo_ext' routine -- operations which define the current integer and -- polynomial moduli, as used aboveprocedurentl_set_modulo(ntl_big_int); -- set the modulusprocedurentl_set_modulo_ext(ntl_modulo_pol); -- set polynomial modulus for modular ring operations -- reconversion to a SETL objectprocedurentl_to_setl_obj(ntl_object); -- reconvert any ntl object to a SETL object -- predicatesprocedurentl_eq(ntl_obj1, ntl_obj2); -- equality test for any two ntl objectsprocedurentl_gt(ntl_obj1, ntl_obj2); -- greater than comparison for integers and floatsprocedurentl_ge(ntl_obj1, ntl_obj2); -- greater equal than comparison for integers and floatsprocedurentl_lt(ntl_obj1, ntl_obj2); -- less than comparison for integers and floatsprocedurentl_le(ntl_obj1, ntl_obj2); -- less equal than comparison for integers and floatsprocedurentl_iszero(ntl_obj); -- test object = 0procedurentl_isone(ntl_obj); -- test object = 1 -- generic arithmetic operatorsprocedurentl_add(RR_b, RR_c); -- additionprocedurentl_sub(RR_b, RR_c); -- subtractionprocedurentl_mul(RR_b, RR_c); -- multiplicationprocedurentl_div(RR_b, RR_c); -- divisionprocedurentl_DivRem(RR_b, RR_c); -- division with remainderprocedurentl_negate(RR_b); -- negationprocedurentl_power(RR_a, n); -- power -- note that the exponent n is a SETL integer hereprocedurentl_XGCD(GF2EX_s, GF2EX_t, GF2EX_a, GF2EX_b); -- some problem with different schema for return values -- ProbIrredTest -- DetIrredTestprocedurentl_BuildIrred_mod2extint_poly(n);procedurentl_BuildIrred_modextint_poly(n);procedurentl_BuildRandomIrred(GF2EX_a);procedurentl_BuildIrred_modint_poly(n); -- constructs irreducible modular polynomial of degree nprocedurentl_BuildIrred_mod2int_poly(n); -- constructs irreducible mod 2 polynomial of degree n -- required other classes yet -- ComputeDegree -- ProbComputeDegree -- TraceMapprocedurentl_BuildFromRoots(vZZ_pE_a);procedurentl_eval(GF2EX_f, GF2E_a);procedurentl_eval2(ZZ_pX_f, ZZ_pE_a);procedurentl_eval3(ZZ_pX_f, vZZ_p_a);procedurentl_eval4(ZZ_pEX_f, vZZ_pE_a);procedurentl_interpolate(vZZ_pE_a, vZZ_pE_b); -- specifically matrix and vector operators with generic componentsprocedurentl_determinant(mRR_A); -- determinant of a matrixprocedurentl_solve(mRR_A, vRR_b); -- solution of an integer linear system, times determinantprocedurentl_inv(mRR_A); -- inverse of an integer linear system, times determinant -- specifically integer operationsprocedurentl_abs(ZZ_a); -- integer abs functionprocedurentl_DivRem(ZZ_a,ZZ_b); -- integer division with remainder -- returns pair [quotient,remainder]procedurentl_rem(ZZ_a,ZZ_b); -- integer remainder, without quotientprocedurentl_probprime(ntl_obj1,n_err); -- test for primality by probabalistic means -- with specified error probability 2 ** -n_errprocedurentl_GenPrime(length,n_err); -- generate a prime probabalistically -- with specified error probability 2 ** -n_errprocedurentl_NextPrime(ZZ_m,n_err); -- find the next prime not less than ZZ_m, -- with specified error probability 2 ** -n_err -- random integerprocedurentl_RandomBnd(ntl_big_int); -- return random integers up to stated big integer bound -- specifically polynomial operationsprocedurentl_GetCoeff(poly, n); -- extract degree n coefficient of polynomial (lowest is degree 0)procedurentl_SetCoeff(poly, n); ???????? -- set degree n coefficient of polynomial; -- return new polynomialprocedurentl_trunc(poly, n); -- truncate polynomial to degree n; -- return new polynomialprocedurentl_RightShift(poly, n); -- right shift polynomial by degree n; -- return new polynomialprocedurentl_LeftShift(poly, n); -- left shift polynomial by degree n; -- return new polynomialprocedurentl_GCD(poly_a,poly_b); -- take GCD of one polynomial modulo another; -- return new polynomialprocedurentl_reverse(polynomial); -- take reverse coefficients of polynomial; -- return new polynomial -- major polynomial factorization proceduresprocedurentl_factor(intpoly); -- factorization of an integer polynomial. -- This is 'final' routine; uses SFCanZass and -- modular factorization as subroutines.procedurentl_SquareFreeDecomp(intpoly); -- squarefree decomposition of an integer polynomial -- returns list of squarefree factors of integer polynomial, -- each showing its exponent. Note that these are all -- relatively prime. This is calculated using GCD algorithm -- modular polynomial factorizations. Note in connection with -- the following that the square-free part of a polynomial can -- always be obtained as f/GCD(f,f')procedurentl_FindRoot(ZZ_pX_a); -- returns single root of monic square free modular polynomialprocedurentl_FindRoots(ZZ_pX_a); -- returns list of roots of monic square free modular polynomialprocedurentl_SFBerlekamp(ZZ_pX_a, Clong_n); -- returns list of roots of monic square free modular polynomialprocedurentl_berlekamp(GF2EX_a, Clong_n); -- returns list of factors, with multiplicities of monic -- modular polynomial. Uses berlekamp methodprocedurentl_SFCanZass(ZZ_pEX_a); -- factorization of monic, square-free modular polynomial -- Uses "Cantor/Zassenhaus" methodprocedurentl_CanZass(GF2EX_a, Clong_n); -- factorization of monic, square-free mod 2 polynomial -- Uses "Cantor/Zassenhaus" methodprocedurentl_diff(GF2EX_a); -- polynomial derivative -- operations modulo an explicit cprocedurentl_InvMod(ntl_big_int1, ntl_big_int2); -- inverse of first integer modulo secondprocedurentl_AddMod(ZZ_a, ZZ_b, ZZ_c); -- a + b modulo cprocedurentl_SubMod(ZZ_a, ZZ_b, ZZ_c); -- a - b modulo cprocedurentl_MulMod(ZZ_a, ZZ_b, ZZ_c); -- a * b modulo cprocedurentl_SqrMod(ZZ_a, ZZ_b); -- square a modulo b -- bit manipulationprocedurentl_NumBits(ntl_big_int); -- get number of bits in an integerprocedurentl_bit(ntl_big_int, long); -- get specified bitendlibntl;

'Dense' images can exist in either byte or floating form. The significance of this is as follows. Red-green-blue images captured using a scanner or digital camera are normally quantized into 256 intensity levels in each of their three color planes. This allows each pixel of an image to be represented by three bytes of data, one for red, one for green, one for blue. This representation produces images very acceptable to the eye when displayed, since the eye's ability to discriminate between very fine shades of color is limited. However, when images in this form are manipulated for purposed of analysis, irritating technical difficulties arise. For example, if two images, each represented by byte data in the range from 0 to 255 are added, overflows can occur. So analysis code working with such images must always be written in a manner that avoids unwanted overflow. Floating images, which maintain one floating point value for the pixeled data in each image plane, eliminate this problem in a radical way. (However, before floating images can be displayed they must be converted back to byte images with data lying in the range 0 to 255).

'Sparse' images can be thought of as fragments of dense images defined only in certain small subregions of a rectangle. Such images typically arise when a dense image is analysed into regions of homogeneous or approximately homogeneous color. Since these regions may be small and scattered, it would be wasteful to represent them by a full dense image. The sparse image representation handles situtations of this kind much more efficiently. Like dense images, sparse images can exist either in byte or floating form.

Working with images involves coordination of multiple tools, since the images must first be created and then displayed. Hence, besides the routines GRLIB_PAK itself, we need to use the Tk system described in Chapter 10 to create windows in which images can be displayed and then to display them. Also, to ease use of the many routines provided by GRLIB_PAK, we prefer to wrap its facilities in an auxiliary object class called 'image'. We further provide an auxiliary package called 'image_pak' containing a varity of image-analysis utilities. We begin by outlining the image handling facilities of the first two of these auxiliary bodies of code.

**Images in Tk**.

The Tk widget system deals with images in two forms. We shall refer to the first of these as 'absolute' images and the other as 'canvas' images. 'Canvas' images can be placed directly into a Tk canvas for display. 'Absolute' images are used to create canvas images, and can themselves be created from the 'grlib image' objects discussed in this section. So the steps leading to the display of an image object are: (1) convert a 'grlib image' object to a Tk absolute image; (2) convert the Tk absolute image to a canvas image; (3) display the canvas image in a Tk canvas. The specific instructions which do this are shown in an example below.

**The 'image' class**.

Images in the SETL system are also handled in two forms. The 'grlib' native package works with a variety of native data structures, built in the C language, which are seen by SETL as opaque objects in the normal way, and which we will refer to as 'grlib images'. The data for these 'grlib images' is normally created by reading in an image file acquired using supplementary external devices and tools, e.g., a digital camera or scanner and some software package that can be used with these devices. However, instead of working with these 'grlib images' directly, we prefer to facilitate their use by wrapping them into auxiliary SETL objects of a class called 'image'. Most of the examples seen below will deal with image objects of this type rather than with the slightly more primitive 'grlib images'.

Now we turn to details. The specifier for GRLIB_PAK is

native packagegrlib; -- interface to SETL native image operations library -- basic functionsproceduregr_clone(img); -- copy an image or <om> if errorproceduregr_const(values, dt,type); -- create constant image with planes indicated by tuple of values -- size is determined by dt = [w,h]; type = 1 if float, 2 if discreteproceduregr_scale(img, t); -- scale image to fill rectangle t = [height,width]proceduregr_destroy_image(img); -- deallocate the imageproceduregr_validate(img_pointer); -- convert a pointer returned by the 'save' operation added to -- the SETL version of the Tk library to a grlib image object, -- which is returned. Note that 'save' has the Tcl syntactic form -- 'save canvas_object_name', and is available as canvas("image") -- in the SETL graphical (Tk) interface. -- string functionsproceduregr_to_string(img); -- return string encoding of image, or OM if errorproceduregr_from_string(stg); -- decode and return image from string encoding, or OM or if error -- I/O functionsproceduregr_read_tiff_image(stg); -- read tiff image from stg = source fileproceduregr_write_tiff_image(stg, img); -- write image in tiff format to stg = destination fileproceduregr_read_jpeg_image(stg); -- read jpeg image from stg = source fileproceduregr_write_jpeg_image(stg, img); -- write image in jpeg format to stg = destination file -- image informationproceduregr_height(img); -- return height of imageproceduregr_width(img); -- return width of imageproceduregr_width_and_height(img); -- return width and height of image, as a tupleproceduregr_type(img); -- return type of image; type = 1 if float, 2 if discreteproceduregr_density(img); -- return 1 if image is dense, 0 if sparseproceduregr_planes_number(img); -- return number of planes in image -- operations which replace an image by a transformed variant of itproceduregr_convert_to_float(rwimg); -- convert discrete image to floating formatproceduregr_convert_to_int(rwimg); -- convert floating image to discrete formatproceduregr_crop(rwimg, t); -- crop image to fit in rectangle t = [l,t,r,b]proceduregr_stuff(rwimg_a, img_b, t); -- stuff img_b into rectangle t = [l,t,r,b] within img_aproceduregr_invert(rwimg); -- in the float image case returns the negative of the image -- in the byte case returns 255-imageproceduregr_offset(rwimg, t); -- the tuple t has the form [x,y] both being not negative -- the operation shifts img, x units to the right and y units down -- putting it in a larger rectangle without any other changeproceduregr_set_pixel(rwimg, coo, t); -- set pixel coo=[i,j] of img to the values of the tupleproceduregr_get_pixel(img, coo); -- get a tuple containing colors at pixel coo=[i,j] of img -- binary operations on imagesproceduregr_plus(img_a, img_b); -- add two imagesproceduregr_minus(img_a, img_b); -- subtract image b from image aproceduregr_times(img_a, img_b); -- multiply two imagesproceduregr_divide(img_a, img_b); -- divide image a bt image bproceduregr_maxim(img_a, img_b); -- form maximum of two imagesproceduregr_minim(img_a, img_b); -- form minimum of two imagesproceduregr_power(img_a, img_b); -- raise image to image power -- binary operation with second argument a constant -- t is for tuple of values, k is for kind of values (1/2: float/int)proceduregr_plus_const(img_a, t, k); -- add constant to imageproceduregr_minus_const(img_a, t, k); -- subtract constant from imageproceduregr_times_const(img_a, t, k); -- multiply image by constantproceduregr_divide_const(img_a, t, k); -- divide image by constantproceduregr_maxim_const(img_a, t, k); -- maximum of image and constantproceduregr_minim_const(img_a, t, k); -- minimum of image and constantproceduregr_power_const(img_a, t, k); -- image to constant power -- unary operationsproceduregr_maxof(img); -- get maximum pixel component value of imageproceduregr_minof(img); -- get minimum pixel component value of imageproceduregr_sin(img); -- take pixelwise sin of imageproceduregr_cosin(img); -- take pixelwise cosin of imageproceduregr_tan(img); -- take pixelwise tan of imageproceduregr_exp(img); -- take pixelwise exponential of imageproceduregr_log(img); -- take pixelwise log of imageproceduregr_abs(img); -- take pixelwise absolute value of imageproceduregr_even(img); -- take pixelwise 'even' of image; pixel is 1 if orig. pixel is evenproceduregr_odd(img); -- take pixelwise 'odd' of image; pixel is 1 if orig. pixel is oddproceduregr_sqrt(img); -- take pixelwise sqrt of imageproceduregr_asin(img); -- take pixelwise arcsin of imageproceduregr_acosin(img); -- take pixelwise arccos of imageproceduregr_atan(img); -- take pixelwise arctan of imageproceduregr_fix(img); -- take pixelwise rounded byte valueproceduregr_floor(img); -- take pixelwise truncated byte valueproceduregr_ceiling(img); -- take pixelwise ceiling byte valueproceduregr_flip_h(img); -- flip image horizontally, return new copyproceduregr_flip_v(img); -- flip image vertically, return new copyproceduregr_rotate(img, n); -- return an image that is a rotation of n*pi degrees from the image img -- image analysis functionsproceduregr_sort(img); -- return a tuple containing the sorted pixels of the imageproceduregr_shrink(img); -- return the smallest rectangle containing the non-black -- pixels of the original image, or <om> if errorproceduregr_sumall(img, power); -- return the first n = power powersums of the imageproceduregr_convolve(img, i_offt, j_offt, coeff_t); -- convolve with the kernel defined by the last 3 tuples. -- these give i,j, and c respectively for the points of -- the convolution kernel.proceduregr_maxvolve(img, i_offt, j_offt, coeff_t); -- maxvolve with the kernel defined by the last 3 tuples. -- these give i,j, and c respectively for the points of -- the convolution kernel. Maxvolve works as follows: -- each pixel p of the image is offset by i,j, and then -- multiplied by c. The maximum (rather than the sum) -- of all the resulting values is formed, and becomes the -- new pixel value of p. In a multi-plane image, -- this is done in each plane.proceduregr_minvolve(img, i_offt, j_offt, coeff_t); -- like maxvolve but with min instead of maxproceduregr_self(img, pair); -- extract a specified range of planes from an image. -- the 'pair' parameter is a pair [from, to] of two zero-indexed valuesproceduregr_self_put(img, img_res, pair); -- insert a specified range of planes into an image. -- the two values [from, to] in the 'pair' = parameter are zero-indexedproceduregr_threshold(img, t); -- thresholds the image down to the values listed in the tuple tproceduregr_histogram(img, t); -- histogram of the given image, with the values listed in the tuple t -- used as bin boundariesproceduregr_gradient(height, width, kind, nplanes, fn); -- constructs a gradient image; fn must be a function of i,j and plane number, and must return a real.proceduregr_unimath(img, fn); -- applies fn to each pixel of image -- (fn is a SETL function of a real var and an integer, -- and must return a real.)proceduregr_perlin(img); -- returns first Perlin transform: sum of all pixels above and to the left of i,jproceduregr_perlin2(img); -- returns second Perlin transform: sum of all pixels to the left of i,j and above -- the rising 45 degree diagonal running left from i,jproceduregr_anaglyph(img, sc, fctr); -- converts the grayscale image 'img' into a random converts the grayscale image 'img' into a random-dot red/green -- anaglyph (viewable in 3D using red/green glasses.) The -- original image, multiplied by the given real factor, -- is used as the desired depth.proceduregr_connected_components(img); -- return a new float image with all the pixels of every connected -- component of the original image marked in a common color, -- different components being marked in different colors. -- operations applicable to sparse images, and conversion to/from sparse -- r/wproceduregr_convert_to_sparse(rwimg); -- convert a dense image to equivalent sparse formproceduregr_convert_to_dense(rwimg); -- convert a sparse image to equivalent dense formproceduregr_shrink_and_cut(rwimg); -- convert the sparse image 'img' to a version in a rectangle -- with upper left hand corner at 0,0, and return the -- smallest rectangle containing the initial image -- no r/wproceduregr_get_level_dense(img_a, t); -- generate a sparse image, consisting of all the pixels of -- the original image which at which the given tuple -- of values appearproceduregr_get_near_dense(img_a, t, dt); -- generate a sparse image, consisting of all the pixels of -- the original image which at which values within distance dt -- of the given tuple of values appear. t is a vector of color -- values, dt the corresponding vector of allowed variationsproceduregr_split(img); -- return a tuple of images, each consisting of all pixels -- at which the original image has a fixed value -- this is returned as a triple of lists image_list,offsets_list,pop_list], where -- each img in img_list is a sparse image in a rectangle with -- upper left hand corner at 0,0, each correspondingoffse in offsets_list is -- the corner of the smallest rectangle, in the original image, -- if rectangle initially containing the corresponding image portion. -- the components (these are integers) of the pop_list give the number of pixels in the -- corresponding component images -- newest additionsproceduregr_superpose(img_a, img_b); -- return the boolean super-union of image_b over image_aproceduregr_random(img, min_val, max_val); -- return the rectangle of the given image, filled with random values in the indicated (real or byte) rangeproceduregr_rectangle(width, height, kind, thick, col); -- construct a rectangle as a sparse image; kind is real or byte; col is a tuple of plane construct a rectangle as a sparse image;proceduregr_ellipse(width, height, kind, thick, col); -- construct an ellipse as a sparse image; kind is real or byte; col is a tuple of plane construct an ellipse as a sparse image;proceduregr_mask(img, mask_img); -- return the given image, masked with the sparse image 'mask_img'; result is sparseproceduregr_lut(img, lut); -- apply the lut at the given byte image -- the lut should have no more than 256 elements; -- if less than 256 elements are provided then, 0 value is used -- for all the remaining valuesproceduregr_lex_min(img); -- return a tuple with a pair [x,y] which is the position of the minimum value found -- in the image img; If the minimum is the value of the default value for a sparse image -- than the pair [ than the pair [-1, -1] is returned. If an error occursOMis returnedproceduregr_lex_max(img); -- return a tuple with a pair [x,y] which is the position of the maximum value found -- in the image img; If the maximum is the value of the default value for a sparse image -- than the pair [ than the pair [-1, -1] is returned. If an error occursOMis returnedproceduregr_stuff_in_place(rwimg_a, img_b, t); -- stuff imageb into rectangle t = [l,t,r,b] within img_aproceduregr_flatten(rwimg); -- flat a sparse image with its default colorproceduregr_tile(t, h, w); -- given a tuple of images that share the same size returns an image -- of size w by hproceduregr_check(img); -- used to do many debugging checks inside an mimageproceduregr_check_mem(); -- return the amount of memory free in the application heapproceduregr_fft(img); -- execute the fast fourier transform on the image img -- and return two imagesproceduregr_fft_inverse(img_re, img_im); -- execute the inverse fast fourier transformproceduregr_wavelet(img); -- execute the wavelet transformproceduregr_wavelet_inverse(img); -- execute the inverse wavelet transformendgrlib;

The GRLIB procedures fall into several groups:

- A first group of operations serve to copy and rescale images, and to convert Tk images into 'grlib' images.
- Two functions for converting images to and from SETL strings containing equivalent data are provided.
- Operations for reading and writing images in TIFF and JPEG format are provided.
- A fourth group of operations get information concerning images, for example, their height and width, number of planes, etc.
- A fifth group of operations converts between byte and floating images, and includes an image-cropping operation, an image reversal operation, an operation for reading and setting individual image pixels, and a few related utilities.
- Next follow a group of 'algebraic' operations which combine two images in pixel-by-pixel fashion to generate a new image. Images can be added multiplied, divided and so forth. These operations are also available in a modified form in which constant values are combined with images.
- Next follow a group of somewhat more than a dozen operations which apply specified functions to images on a pixel by pixel basis. For example, one can form the sine of an image, the logarithm of an image, etc. These nonlinear operations can be useful in image analysis. Operations for flipping and rotating of images are also available.
- Next follow a group of about a dozen additional operations which are designed for transformation and analysis of images. This includes an operation which forms the convolution of an image with a numerical kernel, allowing us to blur images in controlled fashion, for example, to form horizontal or veritcal blurs in preparation for other steps in image analysis. Other operations in this group perform image thresholding, histograming, extraction and insertion of particular image planes, construction of gradient images, etc. Another operation in this set converts a grayscale image into a so-called 'random dot anaglyph' suitable for 3-D viewing with red-green colored glasses. Yet another operation analyzes an image into regions, each of which consists of portions of the original image having identical data values.
- The next group of operations serves for interconversion of images between the two major formats supported, namely 'dense' and 'sparse' images.
- Operations which form and which invert the Fourier transform of an image and its wavelet trnasform are also provided.
- In addition to the operations listed above, miscelleous other operations are available, for example, an operation which transforms an arbitrary image by passing its pixels through a lookup table, operations which form elipses and rectangles of specified size, etc.

These operations will be reviewed in detail below, and many examples wil be given.

As already noted we normally work with SETL images in their 'image object' form. The specifier for the image object class is given below. This class simply wraps the more primitive operations of 'grlib' into a more convenient syntactic form. We will detail the operations of this class later in this section, but in order not to get involved in too many details before their significance is understood it is best to give an example first. Such an example is seen immediately below. This simply reads in a JPEG image file (the image is just a spiral) and displays the image in a Tk window. But even this must involve the preliminary steps seen in the example. We read the image; then initialize the Tk system (this opens a first Tk window; when this window is closed all windows dependent on it will disappear and the SETL program will terminate). Then we enter Tk's main event loop (simply in order to keep the display open until we are finished examining the image).

Image display in this example is handled by the standard procedure 'write_captioned_image' seen below. By following the steps of this rather basic procedure we can start to grow familiar with some of the 'grlib image' operations summarized in the preceding paragraphs, which appear in their 'image object' syntax.

'write_captioned_image' begins by examining the image passed to it: if this image is sparse it is immediately converted to a dense format (since this is the only way in which sparse images can be displayed.) Next one reads the height, width, and number of planes of the image and determines whether it is a floating or byte image. A new Tk window is then opened to hold the image to be displayed and a canvas constructed in this window. If necessary, we pad the image out to a three plane (red, green and blue) image by adding additional planes; or we by drop superfluous planes if the image has more than three. The resulting RGB image is then converted to a Tk image and displayed in the canvas that has been built.

It is worth reading the code for the 'write_captioned_image' procedure closely, since this is a good way of starting to become familiar with the most basic operations on the SETL image class. (Also, to be understood fully, examples having to do with images need to be seen rather than just considered abstractly. Hence, as with all the examples in this section, the following example should not only be read but pasted into a SETL environment and executed.)

programtest; -- test program for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral_img := image("spiral.jpg").to_float(); Tk := tkw(); Tk(OM) := "Image Example 1"; -- create the Tk interpreter write_captioned_image(spiral_img,"Spiral Image"); Tk.mainloop(); -- enter the Tk main loopprocedurewrite_captioned_image(an_image,cap); -- write image to a Tk window var canv_image; -- global used in lambda seen belowifan_image.density() = 1thenan_image := an_image.to_dense();end if; -- density is 1 for sparse images [height,width,num_planes,image_type,image_dens] := #an_image; -- image_type is 1 for floating image, 2 for byte ('grayscale') -- image_dens is 0 for dense images, 1 for sparse images toplev := Tk("toplevel","10,10"); -- create a new Tk window to hold the displayed image ca := toplev("canvas",str(height) + "," +str(width)); ca("side") := "top"; -- create a canvas in the window three_zeroes :=ifimage_type = 1then[0.0,0.0,0.0]else[0,0,0]end if; -- create a vector of 3 zeroes of appropriate type -- use the Grlib image to create a Tk imageifnum_planes = 1then-- pad out the image with additional planes extended_image := image([three_zeroes,[height,width]]); -- set up an auxiliary image of 3 planes,initially blackforjin[1..3]loopextended_image(j..j) := an_image;end loop; -- copy the black and white image into all 3 planes of the auxiliary image tk_abs_image := Tk("image",extended_image.to_int()); -- use this 3 use this 3-plane 'grlib' image to create a Tk imageelseifnum_planes < 3then-- pad out the image with additional planes extended_image := image([three_zeroes,[height,width]]); -- set up a black image of 3 planesforjin[1..num_planes]loopextended_image(j..j) := an_image(j..j);end loop; -- copy the each plane of the image into all the corresponding plane of the auxiliary image tk_abs_image := Tk("image",extended_image.to_int()); -- use the 3 use the 3-plane 'grlib' image to create a Tk imageelseifnum_planes > 3then-- drop the extra planes tk_abs_image := Tk("image",an_image(1..3).to_int()); -- use the 3 use the 3-plane 'grlib' image to create a Tk imageelsetk_abs_image := Tk("image",an_image.to_int()); -- just use the 3 just use the 3-plane 'grlib' image to create a Tk imageend if; toplev(OM) := cap; -- enter the window title canv_image := ca("image",tk_abs_image); -- create as canvas image and place it canv_image("coords,anchor") := "0,0;nw"; -- place to make the image visible toplev{"Destroy"} :=lambda(); canv_image("coords") := OM;end lambda; -- remember to destroy the canvas image when the window is closed an_image := extended_image := OM; -- release the images -- ?? should not be neededendwrite_captioned_image;endtest;

The specifier for the image class is shown below. However, many of the most important operations on objects of image class are represented not in functional syntax but by overloading SETL's infix syntax and other aspects of its built-in syntax. These operations are summarized in the following table.

img + x; | sum of two images |

img - x; | difference of two images |

-img; | invert a floating or byte image |

img with img; | superpose two images |

img * x; | product of two images, or convolution |

x * img; | convolution or function application |

img / x; | quotient of two images; or division by a constant |

img ** x; | image to an image or constant power |

img max x; | maximum of two images, or 'maxvolution' |

img min x; | minimum of two images, or 'minvolution' |

img(i..j); | subimage, consisting of designated planes |

img(i..j) := x; | set subimage consisting of designated planes |

img(x); | if x is a pair of integers, returns the indicated pixel if x is a pair of images, forms the distorted image if x is a quadruple defining a rectangle, forms the cropped image |

img(x) := y; | if y is an image and x is a pair defining an offset, does the 'stuff' operation; but if y has one more plane than img, this final plane is used as a blend plane. if y is a tuple, modifies one pixel |

img{x} := y; | in-place stuff of sparse image y into dense or sparse image img |

domain img; | flatten a sparse image with its default color |

#img; | returns tuple of length 5consisting of [height,width,number of planes,float_or_byte (= 1/2),sparse_or_dense (= 1/2)] |

As appears in the preceding table, algebraic combinations such as sums, differences, products, quotients, powers, maxima, and minima can be written in standard infix syntax, i.e., as image1 + image2, etc. The SETL slice and slice assignment syntax is used to extract designated planes of an image and to insert planes of one image into specified planes of another. The SETL indexing syntax is used to form a cropped image. The '#' operator produces a tuple of length five containing the height, width, number of planes, floating/byte type, and dense/sparse type of an image.

Examples illustrating the use of all these operations are given later in this section.

The image class represents the other operations of 'grlib_pak' in the more 'functional' syntax seen in the following specifier.

classimage; -- SETL image class using the native image libraryselwidth_comp(1),height_comp(2),nplanes_comp(3),int_or_float_comp(4),density_comp(5); -- selectors for the components of #img.constint_image_kind := 2; -- flag of byte images; flag for floating is 1procedurenative_im(); -- gets the native part of an image objectprocedurecreate(source); -- creates an image from a file or existing imageprocedurewrite(file_name); -- writes image to file_name in tiff formatprocedureto_float(); -- returns floating version of imageprocedureto_int(); -- returns byte version of imageprocedureto_random(min_val, max_val); -- returns randomized version of single_plane imageprocedureto_string(); -- convert to SETL string representationprocedurescale(t); -- scales image to t = [height,width]proceduresort(); -- returns sorted list of distinct pixel valuesprocedureshrink(); -- returns boundary of smallest subrectangle containing positive values onlyprocedurecomponents(black_val); -- returns triple: [image with connected components given constant value, -- in order of number of pixels in component; vector of topmost_leftmost -- points of components; population of component]proceduresumall(n); -- returns the vector of sums of the first n powers of the pixels of the floating image xprocedurefourier(); -- returns real fourier cosine transform of the imageprocedureperlin(); -- returns perlin transform of the imageprocedureperlin2(); -- returns alternative perlin transform of the imageprocedurethreshold(x); -- threshold of image; reduces to 0.0 or downward to thresholdprocedureflip(how); -- flips the image: flips the image: how can be "h","v", 1, 2, or 3 --procedure maxvolve(kernel); -- utility sparse convolution, returns maximum of valuesprocedureanaglyph(scale,fraction_black); -- converts grayscale image to random converts grayscale image to random-dot anaglyph; scale is depth scaleprocedurehistogram(v); -- v is an increasing list of floating values, whose length is n.proceduremax_image(y); -- returns an image representing a 'local maximum' of an image.proceduremin_image(y); -- returns an image representing a 'local minimum' of an image.procedurenearest(y); -- returns a pair of images representing the nearest positive pixel k,m -- to a given pixel i,jprocedureskeletonize(); -- Returns the skeletonized form of an image.procedureultimate_points(); -- Returns the 'ultimate_point' image of an image, after erosionprocedurewatershed(); -- Returns the 'watershed' image of an image, as defined by 'NIH Image'.procedurefilter(filter_name); -- returns image after application of designated Photoshop filter.proceduredither(dither_method_name); -- returns dithered version of imageprocedureget_level(pix); -- returns sparse image consisting of all pixels with given valueprocedureget_levels_near(pix,dpix); -- returns sparse image consisting of all pixels with value near given valueprocedureto_sparse(); -- converts dense to sparse imageprocedureto_dense(); -- converts sparse to dense imageprocedureshrink_and_cut(); -- shrink the image and return its rectangleprocedureoffset(t); -- offset a sparse image by a positive amount, enlarging its rectangle as necessaryproceduredensity(); -- image densityproceduresplit(); -- split a dense image into a list of sparse images, by pixel values. -- return a tuple of images, each consisting of all pixels -- at which the original image has a fixed value -- this is returned as a triple of lists image_list,offsets_list,pop_list], where -- each img in img_list is a sparse image in a rectangle with -- upper left hand corner at 0,0, each correspondingoffse in offsets_list is -- the corner of the smallest rectangle, in the original image, -- if rectangle initially containing the corresponding image portion. -- the components (these are integers) of the pop_list give the number of pixels in the -- corresponding component imagesproceduremask(img2); -- mask the first image by the second; if that is not sparse, just cropprocedurelut(t); -- apply the lut t to the first plane of the given byte imageproceduretile(t); -- tile the given byte image with the images in tuple t, all of which -- must have the same sizeprocedurefft(); -- direct version of fft. Takes 2**n by 2**m image, returns pair of -- 2**n by 2**(m 2**n by 2**(m - 1) + 1 images representing real and imaginary parts of -- fourier transform. In the first coordinate, the frequencies run from -- 0 to 2**n 0 to 2**n - 1; in the second coordinate, from 0 to 2**(m - 1); -- the real and imaginary parts of the full symmetrical fourier image, -- with the 0,0 frequency at ist center, are therefore obtained by -- moving the left half of each image returned to its rignt, and then -- adjoining a horizontally and vertically flipped copy to at the bottom; -- this represents the middle frequency 2**(m this represents the middle frequency 2**(m - 1) in the second coordinate -- twice.procedure inv_fft(img2); -- inverse version of fft. img.inv_fft(img2) takes a 2**n by 2**(m inverse version of fft. img.inv_fft(img2) takes a 2**n by 2**(m - 1) + 1 -- pair of images like those returned by gr_fft, and reconstructs the original imgendimage;

One group of three related operations on image objects, namely 'convolve', 'maxvolve', and 'minvolve' deserve special explanation. Of these the convolution operation written in the form img* kernel, where kernel is a SETL tuple of a form to be explained, is the best known. The 'kernel' that appears in this operation must be a list of triples, each of the form [i,j,c]. Here i and j are two integers defining offsets, and c is a numerical coefficient, which should be real if the image 'img' is of real type, but should be an integer if 'img' is a byte image. The convolution operation img* kernel acts by iterating over all the triples in 'kernel', and for each of these shifting img by the amount i,j and then multiplying it by the coefficient c. All of the translated images formed in this way are then added to give a final result.

A typical use of this operation is to blur an image. Another typical use is to detect image edges. For blurring, one uses a kernel in which the coefficients are all positive and sum to 1. A convolution operation with such a kernel represents the average value of multiple shifted copies of an original image and so accomplishes blurring. To find edges one uses a kernel whose coefficients sum to 0; hence, some of these coefficients are positive whereas others must be negative. A kernel of this sort will reduce an image to 0 in each region where the image is constant, and so tends to favor image regions in which the image's intensity values are varying rapidly, i.e., edges. A variety of such edge detection convolution kernels will be described below.

'maxvolve' and 'minvolve', non-linear variants of 'convolve', which substitute the max and the min operations respectively for addition. They are written in the form img.maxvolve(kernel), etc., and take kernels of the same form as convolve. 'maxvolve' acts by forming the same set of shifted images as does 'convolve', but then combines these images by taking their maximum rather than their sub; likewise, 'minvolve'. 'maxvolve' can therefore be used to detect all pixels which lie near some pixel of large value, and 'minvolve' can be used to detect all pixels which lie near some pixel with small values. Later in this section we will give examples of the use of 'maxvolve' and 'minvolve' to find significant geometric features in images, for example, long horizontal or vertical lines.

In addition to the SETL image class and all the operations on image objects listed in the preceding paragraphs, an extensive collection of image utilities built on the operations of the image class is provided. These operations are collected in the package 'image_pak' supplied with the SETL system. The operators in this auxiliary package are reviewed later in the present section. Many of these operations duplicate and extend the image analysis facilities available in the widely-used 'NIH image' collection.

packageimage_pak; -- package of image analysis utilities var output_image_count := 0;proceduremean(img,n); -- average a floating-point image over n by n squares (n odd)procedurelaplace(img,n); -- Laplace a floating-point image over n by n squares (n odd)procedureerode(bin_img); -- erode a binary image by removing all pixels which contact the boundaryproceduredilate(bin_img); -- dilate a binary image by adding all pixels which contact the boundaryprocedureerode_horiz(bin_img); -- erode a binary image by removing all pixels which contact the boundary in the horizontal directionproceduredilate_horiz(bin_img); -- dilate a binary image by adding all pixels which contact the boundary in the horizontal directionproceduredilate_vert(bin_img); -- dilate a binary image by adding all pixels which lie within 2 pixels of the boundary in the vertical directionprocedurerefill_vert(bin_img,vert_img,n); -- dilate a binary image by n pixels in the vertical direction to refillprocedurerefill_horiz(bin_img,horiz_img,n); -- dilate a binary image by n pixels in the horizontal direction to refillproceduremaxrank(img); -- take largest pixel value in 3 by 3 neighborhoodprocedureminrank(img); -- take smallest pixel value in 3 by 3 neighborhoodprocedureoutline(bin_img); -- outline a binary imageprocedurefind_edges(img); -- standard Sobel edge-finding algorithmprocedurefind_corners(img); -- corner-finding algorithm, based on horizontals and verticalsprocedurefind_point_up(img); -- find sharp point going upprocedurefind_point_down(img); -- find sharp point going downproceduresharpen(img,c); -- sharpen filtersproceduresharpen_limit(img,c); -- limited sharpen filtersprocedureemboss(img,kind); -- emboss filters; should be applied to grayscale imageproceduremake_gray(img); -- convert a 1-plane image to grayscaleprocedureenhance_contrast(img); -- enhanced contrast filterproceduresubtract_background(img); -- subtract the average background from an imageprocedureshadow(img,kind,c); -- shadow filtersprocedurefind_edges_horiz(img); -- horizontal edge-finding algorithmprocedurefind_edges_vert(img); -- vertical edge-finding algorithmprocedurestd_dev_image(img,rect); -- std_dev-image algorithmprocedurestandard_kernel(img,kind); -- take standard kernel transform of a floating-point imageprocedureseen_from_above(bin_img); -- a point is seen from above if it is white and the two pixels above it are blackprocedureseen_from_below(bin_img); -- if it is white and the two pixels below it are blackprocedureseen_from_left(bin_img); -- if it is white and the two pixels to its left are blackprocedureseen_from_right(bin_img); -- if it is white and the two pixels to its right are blackproceduredrop_dots_for_ab(bin_img); -- drop dots from binary imageproceduredrop_dots_for_lr(bin_img); -- drop dots from binary imageproceduredrop_dots_for_ab3(bin_img); -- drop dots from binary imageproceduredrop_dots_for_lr3(bin_img); -- drop dots from binary imageproceduremake_spherize_pair(img); -- create a distortion-image that spherizes 'img'proceduremake_distort_pair(img,f,g); -- convert a set of functions into a distortion-image pair -- routines specialized for music readingproceduremusic_edges_horiz(img); -- horizontal long thin edge-finding algorithmproceduremusic_verylong_edges_horiz(img); -- horizontal very long thin edge-finding algorithmproceduremusic_edges_vert(img); -- vertical long thin edge-finding algorithmproceduremusic_blobs(img); -- find round blobs in imageprocedurefind_notes_top(mev,img,inv_img); -- find notes at top of note stemsprocedurefind_notes_bot(mev,img,inv_img); -- find notes at the bottom of note stemsprocedurewrite_captioned_image(an_image,cap); -- write image to a Tk windowendimage_pak;

Having said all this, we begin reviewing the available operations on image objects one-by-one.
Our second example illustrates the flip, scale, **max**, **min**, and convolution operations. Note that **max** and **min**, like the other image operations written in infix form, can be applied either between a pair of image objects, or between an image object and a tuple of constants representing a constant image. Two images being combined by such an operation must be *compatible*, that is, must have the same number of planes, must both either be byte or floating images, and must both either be dense or sparse images. Also when an image object is combined with a tuple of constants using one of these operations, the tuple must be compatible with the image, i.e. its length must equal the number of planes in the image, and if the image is real (resp. byte) the components of the tuple must be integers of appropriate value.

This example also illustrates the fact that some of the infix and other SETL image operations which are written using 'overloaded' syntax can represent one or another operation depending on the structure of their non-image arguments. The **max**, **min**, and product operations illustrate this.

programtest; -- test program 2 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral_img := image("spiral.jpg").to_float(); -- read in a spiral image spiral_img := spiral_img([1,1,320,320]); spiral_mirror := spiral_img.flip("h"); -- flip the image horizontally spiral_turn := spiral_img.flip(1); -- rotate the image 90 degrees spiral_scale := spiral_img.scale([100,500]); -- rescale the image to 100 by 500 spiral_max := spiral_imgmaxspiral_turn; -- form image from brightest pixels spiral_min := spiral_imgminspiral_turn; -- form image from darkest pixels spiral_rg := spiral_imgmin[255.0,255.0,0.0]; -- suppress the blue plane spiral_turn_b := spiral_turnmin[0.0,0.0,255.0]; -- keep only the blue plane kernel := [[1,1,0.25],[-1,1,0.25],[1,-1,0.25],[-1,-1,0.25]]; spiral_blurred := spiral_img * kernel; spiral_blurred3 := spiral_blurred * kernel * kernel; Tk := tkw(); Tk(OM) := "Image Example 2"; -- create the Tk interpreter write_captioned_image(spiral_img,"Spiral Image"); write_captioned_image(spiral_mirror,"Spiral Image Mirrored"); write_captioned_image(spiral_turn,"Spiral Image Turned"); write_captioned_image(spiral_scale,"Spiral Image Rescaled"); write_captioned_image(spiral_max,"Max of 2 spirals"); write_captioned_image(spiral_min,"Min of 2 spirals"); write_captioned_image(spiral_rg,"Min of 2 spirals"); write_captioned_image(spiral_rg + spiral_turn_b,"Sum of 2 spirals"); write_captioned_image(spiral_blurred,"Spiral Image Blurred"); write_captioned_image(spiral_blurred3,"Spiral Image Blurred More"); Tk.mainloop(); -- enter the Tk main loopendtest;

The next example shows the use of sparse images. We read in the same spiral image as before, and then use the 'get_level' operation to extract the region in this image that is occupied by perfectly white pixels. This region is extracted as a sparse image. Then we use the 'shrink_and_cut' operation to cut this image back to the smallest subimage containing all the white pixels, while at the same time determining the containing rectangle of those pixels. Having done this we use the **domain** operation to reconvert the sparse image to dense form. Next we invert the pixel values of the image; since this is a byte image a black-to-white reversal results. Next various of the operations which extract and manipulate specified image planes of our three plane image are shown. One of the images formed by the code is written out to a file and then read back to verify that it has been written correctly. Finally we set up a kernel suitable for use in the maxvolve and minvolve operations (which are represented syntactically as **max** and **min**). Given the way that our kernel is set up, maxvolve sets each pixel to the maximum of the four pixels which surround it horizontally and vertically, and minvolve sets it to the minimum of these same four pixels.

Our next example shows several operations which construct rectangular and elliptical images, gradients, and 3-D random-dot anaglyphs for viewing with red-green glasses. All of these operations are represented simply as image creation operations with a variety of parameters. If the parameter given begins with the string 'oval' then the remaining parameters should be the height and width of the oval, the float/byte character of the image (float = 1, byte = 2), the thickness of the oval (which is normally an outline oval, but can be filled by setting its thickness parameter appropriately), and finally, a tuple giving the color values for as many planes as the oval should have. This tuple should be a tuple of real numbers if the image is to be real, or a tuple of byte values if a byte image is desired. The image formed will be sparse in any case since this gives a much smaller representation.programtest; -- test program 3 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral_img := image("spiral.jpg"); -- read in a spiral image, keeping it in byte form spiral_img := spiral_img([1,1,320,320]); -- crop to perfect square spiral_sparse := spiral_img.get_level([255,255,255]); -- extract region of white pixels as a sparse image [white_pixels_only,rect_containing] := spiral_sparse.shrink_and_cut(); -- cut the white cut the white-pixels-onl image to the smalest rectangle contining it, and note this rectangleforwhite_pixels_only image: ",#white_pixels_only);domainspiral_sparse; -- reconvert to an equivalent dense image spiral_reverse := -spiral_img; -- invert the pixel values spiral_copy := spiral_copy2 := spiral_img; -- make two copies of the spiral image spiral_copy(3..3) := spiral_reverse(1..1); -- set its blue plane from the reversed image -- spiral_copy2(1..2) := spiral_reverse(1..2); -- set its green and blue planes from the reversed image spiral_copy.write("spiral_reverse.jpg"); -- write the modified image to a file spiral_reread := image("spiral_reverse.jpg"); -- reread the image that has just been written kernel := [[1,1,1],[-1,1,1],[1,-1,1],[-1,-1,1]]; -- set up kernel for 'maxvolution' and 'minvolution' spiral_maxvolve := spiral_imgmaxkernel; -- form 'maxvolution' spiral_minvolve := spiral_imgminkernel; -- form 'minvolution' Tk := tkw(); Tk(OM) := "Image Example 3"; -- create the Tk interpreter write_captioned_image(spiral_reverse,"Spiral Reversed"); write_captioned_image(spiral_copy,"SpiralwithReversed Blue Plane"); -- write_captioned_image(spiral_copy2,"Spiral with Reversed Green and Blue Planes"); write_captioned_image(spiral_maxvolve,"Spiral Image 'Maxvolved"); write_captioned_image(spiral_minvolve,"Spiral Image 'Minvolved"); write_captioned_image(spiral_sparse,"Spiral White Pixels"); write_captioned_image(spiral_flat,"Spiral White Pixels 2"); write_captioned_image(white_pixels_only,"Spiral White Pixels Only"); write_captioned_image(spiral_reread,"Spiral After WriteandReread"); Tk.mainloop(); -- enter the Tk main loopendtest;

Next we show how to form a gradient image. Here again the operation appears as an image creation, but now the parameter is a tuple consisting of the image width and height, number of planes, float/byte character in the same encoding as before, and gradient function. The fifth, function-valued, component of this tuple must be a function of three parameters i, j, p, where i and j define a pixel position and p is a plane number. This function must return the color value which the indicated pixel has in the indicated plane. The value returned should be a real number if the image is floating and a byte value if the image is to be a byte image.

The final operation in this example, 'anaglyph', is a specialty routine which converts one plane floating dense images to random-dot ,red-green anaglyphic images suitable for viewing through red-green glasses. As seen in the example, the first parameter of the 'anaglyph' function is a (floating) 'scale' value which converts the grayscale levels of the input image to perceived depth values in the anaglyph. Relatively small values, generally not more than one-tenth should be used for this parameter. Values larger than this will generate anaglyphs that the eye cannot perceive. The second parameter of the anaglyph method defines the amount of green to be used in the anaglyph, versus the amount of red. For this parameter values near one-half are appropriate. The code seen below actually builds two anaglyphs, one using a grayscale oval and the other using a grayscale gradient, simply to show the appearance of a 3-D gradient.

programtest; -- test program 3 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window ellipse := image(["oval",300,200,2,40,[255,0,255]]); -- construct an elliptical outline as a sparse byte image, specifying its width rect := image(["rectangle",300,200,1,40,[255.0,0.0,255.0]]); -- construct a rectangular outline as a sparse floating image, specifying its width gradient := image([300,300,3,1,rgb_fcn]); -- gradient image from rgb_fcn e1 := image(["oval",300,200,1,40,[255.0]]).to_dense(); -- build a 1 build a 1-plane ellipse as a dense floating image e1 := e1.anaglyph(0.05,0.5); -- convert it to a random convert it to a random-dot anaglyph grad1 := gradient(1..1).anaglyph(0.05,0.5); Tk := tkw(); Tk(OM) := "Image Example 3"; -- create the Tk interpreter write_captioned_image(ellipse,"Ellipse"); write_captioned_image(rect,"Rectangle"); write_captioned_image(gradient,"Gradient"); write_captioned_image(e1,"Ellipse as Random_dot Anaglyph"); -- write an anaglyph write_captioned_image(grad1,"Gradient as Random_dot Anaglyph"); -- write an anaglyph Tk.mainloop(); -- enter the Tk main loopprocedurergb_fcn(i,j,plane); -- function for creating gradient; in this context must return triple of reals fi :=float(i); fj :=float(300 - j);return ifplane = 0thenfi * 255.0elseifplane = 1thenfj * 255.0else255.0 * 300.0 *sin((fi + fj)/50.0)end if/ 300.0;endrgb_fcn;endtest;

Our fourth example begins by exhibiting two of the operations, 'histogram' and 'sumall', provided for the statistical analysis of images. 'histogram' accepts an increasing tuple of values defining bin boundaries, and histograms an image into the corresponing bins, reporting the number of pixels with values in each of the bins. 'sumall' generates a tuple of specified length containing the first few moments of an image.

Next we demonstrate the use of two utility operations suitable for certain types of image analysis. These are called 'perlin' and 'perlin2', after Ken Perlin, who was one of the first to realize their advantages. The 'perlin' operation transforms an image by setting its pixel value at position i, j, to the sum of all pixel values lying left and above position i, j, in an original source image. This transform can be built quite rapidly by a single sweep over the pixels of an image which only needs to keep the values in the immediately preceding row while the next following row are being built. But, once this transform is available, it is easy to form averages of an image over an arbitrary rectangle simply as the sum of the 'perlin' transform values at the upper left and lower right corners of the rectangle minus the same values at the upper right and lower left corners. This gives us a very fast way of blurring an arbitrary image heavily by convolving it with a kernel containing hundreds or even thousands of identical values. (The fast Fourier transform is another, much more complex but much more flexible way of doing this efficiently.) If the same operation were attempted without the use of some such special approach it would be extremely inefficient. Moreover, the 'perlin' transform has the advantage that it allows the size of the rectangle used for averaging a source image to vary from point to point; a fast Fourier transform approach would not allow this.

programtest; -- test program 4 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral := image("spiral.jpg").to_float(); -- read the spiral imagefloat(j): jin[1..10]])); -- histogram the image into indicated binsinSecond Mode"); -- write blurred spiral Tk.mainloop(); -- enter the Tk main loopendtest;

Our next example shows the use of the 'get_levels_near', 'threshold', 'components', 'split', and 'stuff' operation. We also show a few operations which produce ordinary SETL values rather than images. The program begins by reading in the usual spiral test image, and then creates one large ellipse and two small ellipses and one small rectangle. The 'get_levels_near' operation is then used to extract all pixels of the spiral image which lie within a distance 90 in each color plane from the level 100 in that plane. The parameters of 'get_levels_near' are the central vector of color levels in the color zone to be isolated, and the distance from these central levels to be included in this zone.

Next we use the 'threshold' operation to reduce the color values in each spiral position to 25 in each color plane, after which the 'get_level' operation is used to isolate all pixels with color values at this thresholded level. Note that both 'get_level' and 'get_levels_near' return sparse rather than dense images. After this, the 'split' operation is used to decompose the thresholded image into its 'parts', where each such part consists of all pixels having a common color value. In the case of our thresholded image there is just one such part (precisely because the image has been thresholded). In other images there might be many and they might be quite small. For this reason the 'parts' are generated as sparse images. The 'split' operation returns the list of all these parts, and for each one of them its smallest enclosing rectangle and total number of pixels. These are returned as three lists, the first a list of sparse images representing all the parts of the original image, the second a list of rectangles, each being the enclosing rectangle of the corresponding part in the first list, and finally a list of integers, each being the number of pixels of the corresponding part in the first list. This is shown in the code by printing the list of parts returned by 'split'.

After this we apply the more sophisticated 'components' operation, which isolates the connected components of the image to which it is applied. If you look carefully at the thresholded spiral image you will see that it actually consists of a pair of spirals, which are intertwined but never touch each other. For this reason the thresholded spiral image has three separate connected components: the two spirals in it and the black background against which they appear. The 'components' operation returns a single image, in which all of the pixels of each separate component of the original have been given the same color value. This color value is an integer (but is represented as the corresponding real number). The components are assigned the values 0.0, 1.0, 2.0, ... in increasing order. This is shown in the code by using the 'sort' operation, which collects all the pixel values from an image and returns them in sorted order. The 'components' operation has just one parameter, which should be the vector of colors corresponding to the 'outside' of the image (this will normally be black).

After this we demonstrate the use of the 'stuff' operation, which inserts a sparse image into either a sparse or a dense image at a given position. The 'image' class requires this operation to be written as an indexed assignment.

Finally, we show the use of a very important operation, which links Tk into our image processing library by making it possible to capture the contents of any Tk canvas as a 'grlib' image. This is written as a Tk operation having the form

Here, 'canvas' can be any Tk canvas; the operation captures the current contents of this canvas and converts it into an object of 'image' type. Note that the availability of this operation allows us to use the Tk facilities interactively to create collages of images, to which text can readily be added, and then to capture these collages as 'grlib' images. Since 'grlib' images can also be displayed at specified locations in Tk canvases, we have, all in all, a powerful set of mechanisms for creating, manipulating, and analyzing images.

<programtest; -- test program 5 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral := orig_spiral := image("spiral.jpg").to_float(); -- read the spiral image and make a copy e1 := image(["oval",60,40,1,4,[255.0,0.0,0.0]]).to_dense(); -- build a 1 build a 1 build a 1 build a 1-plane ellipse as a dense floating image e2 := image(["oval",300,200,1,44,[255.0,255.0,0.0]]).to_dense(); -- build a 1 build a 1-plane ellipse as a dense floating image e3 := image(["oval",60,40,1,4,[0,0,255]]).to_dense(); -- build a 1 build a 1 build a 1 build a 1-plane blue ellipse as a dense floating image e4 := image(["rectangle",60,40,1,4,[0,255,0]]).to_dense(); -- build a 1 build a 1 build a 1 build a 1-plane green rectangle as a dense floating image tiled := e2.tile([e3,e4]); -- tile large rectangle with ellipses of various colors dark_part := spiral.get_levels_near([100.0,100.0,100.0],[90.0,90.0,90.0]); -- extract part of spiral in color range below 190.0 thesholded := (spiral.threshold([25.0,25.0,25.0]).get_level([25.0,25.0,25.0])) * [10.0,10.0,10.0]; -- bring all non bring all non-black pixels of the spiral to level [25.0,25.0,25.0]true,1,#ffffff"; -- outline is 1 pixel wide, black the_spline1("fill") := "#ff0000"; -- fill with red the_data2 := "158.00000000,10.00000000,8.00000000,310.00000000,198.00000000,114.00000000,128.00000000,28.00000000"; the_spline2 := the_canvas("polygon",the_data2); -- create a spline from the preceding data the_spline2("smooth,width,outline") := "true,1,#ffffff"; -- outline is 1 pixel wide, black the_spline2("fill") := "#0000ff"; -- fill with blue ct := the_canvas("text","Red text in the canvas"); -- add some text to the canvas ct("coords") := "30,30"; ct("anchor,font,fill") := "nw,{Times 36},red"; ct2 := the_canvas("text","Blue text in the canvas"); ct2("coords") := "130,40"; ct2("anchor,font,fill") := "center,{Times 36},blue"; Tk.update(); -- update to force draw canv_img := the_canvas("image"); -- capture the canvas contents as a grlib image write_captioned_image(canv_img,"Canvas image captured"); -- display the captured image write_captioned_image(spiral,"Spiral Stuffed"); -- write the stuffed spiral write_captioned_image(orig_spiral,"Original Spiral"); -- write the stuffed spiral write_captioned_image(e2,"Large Ellipse Stuffed"); -- write the stuffed large ellipse write_captioned_image((abs* (sin* spiral)) * 255.0,"Sin of Spiral"); -- write the stuffed large ellipseendtest;

The following example shows the 'shrink', 'max_pos_val', 'min_pos_val', 'to_spare', 'to_dense', 'mask' and the image stuff-in-place operation, which is written using curly brackets as img{offset_tuple} := sparse_img, and the available image I/O operations. We create a spiral image and then apply various operations to it. The 'shrink' operation returns the smallest rectangle containing all the non-black pixels of an image. 'max_pos_val' finds the pixel with largest value (in lexicographic order of pixel values), and the maximum numerical value occurring in any pixel value tuple. 'min_pos_val' acts similarly but finds minima rather than maxima. 'to_spare' and 'to_dense' convert an image between its spare and dense representations; the 'domain' operator can also be used to convert a spare image to the corresponding dense form. Image I/O is available for two image formats, JPEG and TIFF. The JPEG format compresses images, but at the cost of changing some of their less-noticed pixel values. The TIFF format preserves pixel values precises. To write an image in JPEG format one uses the 'write' operation supplying a file name which ends in '.jpg'. To write an image in TIFF format the name of the file written should end in '.tif'. Images are read from files by operations having the syntactic form of image creation operators, but with parameters that are file names. The image being read will be converted from JPEG format if the name of the file containing it ends in '.jpg', or if this file name ends in '.tif', then the image will be converted from TIFF format.

programtest; -- test program 6 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral := orig_spiral := image("spiral.jpg"); -- read the spiral image and make a copyandval",spiral.max_pos_val());andval",spiral.min_pos_val()); spiral_sparse := spiral.to_sparse(); -- convert the spiral to sparse, and then back to dense spiral_redensed := spiral_sparse.to_dense(); -- write the spiral in tiff and jpeg formats, and reread it spiral.write("spiral_test.jpg"); -- write image to file_name in jpeg format spiral.write("spiral_test.tiff"); -- write image to file_name in tiff format spiral2 := orig_spiral := image("spiral_test.jpg").to_float(); -- read the spiral image, jpeg format spiral3 := orig_spiral := image("spiral_test.tiff").to_float(); -- read the spiral image, tiff format e1 := (esp := image(["oval",60,40,1,4,[255.0,0.0,0.0]])).to_dense(); -- build a 3 build a 3-plane ellipse as a dense floating imagedomain(esp)," ",#esp_flat); -- do an in do an in-place stuff of small sparse ellipse into spiral orig_spiral{[10,10]} := esp; orig_spiral{[15,15]} := esp; orig_spiral{[25,25]} := esp; e2sp := image(["oval",300,200,1,40,[255.0,255.0,0.0]]); -- build a 3 build a 3-plane ellipse as a sparse floating image spiral_scaled := spiral.scale([300,200]); -- scale the spiral image to match the ellipse spiral_masked := spiral_scaled.mask(e2sp); -- mask the spiral with the sparse ellipse Tk := tkw(); Tk(OM) := "Image Example 5"; -- create the Tk interpreter write_captioned_image(spiral_redensed,"Spiral redensed"); -- write the spiral, recovered from sparse form write_captioned_image(esp_flat,"Flattened ellipse"); -- write flattened ellipse write_captioned_image(spiral2,"Spiral image jpg"); -- write the written and re write the written and re-read spiral, jpeg format write_captioned_image(spiral3,"Spiral image tiff"); -- write the written and re write the written and re-read spiral, tiff format write_captioned_image(orig_spiral,"Stuffed spiral"); -- write the spiral stuffed with the small ellipse write_captioned_image(spiral_masked,"Spiral masked"); -- write the spiral masked by the sparse ellipse Tk.mainloop(); -- enter the Tk main loopendtest;

**Fourier transform of images**. The logarithm is important because it makes the isomorphism of real addition and real multiplication of positive number explicit. Historically this was used to accelerate paper and pencil calculation of product: with a table of logarithms in hand one could multiply two numbers by taking the logarithm of each, adding these logarithms and then taking the inverse logarithm of the result. Similarly the Fourier transform is important because it makes an isomorphism between convolution of functions or of images and function or image multiplication explicit. To form the convolution of image by a kernel one can take the Fourier transform of each, multiply them, and then take the inverse Fourier transform of the result. Of course, this approach is only useful if the Fourier transform of an image can be calculated more rapidly than convolutions can be calculated directly. The time required to form the convolution of an image and a kernel is proportional to the product of the number of pixels in the image by the number of terms in the kernel. If the kernel resembles a dense image, in the sense of having about as many terms as there are pixels in an image, then the time required to calculate the convolution using a naive direct approach is proportional to the square of the number of pixels in the image, and this can be very large. For this reason it is very significant that there exists a *fast* method for calculating the Fourier transform, in time proportional to the number of pixels in the image times the logarithm of this number of pixels. The availability of this fast technique has had a revolutionary effect on all types of calculation whose crucial operation is a convolution of some kind. This now includes many image processing and restoration operations including the generation of three-dimensional from hospital CAT scans. Optical and radio astronomy have also been much effected by the use of fast Fourier transform methods.

The operations 'fft' and 'inv_fft', which appear in the following example respectively calculate and invert the Fourier transform of an image. In the form supported by our image library these operations require that the number of rows and number of columns in the image should both be powers of two (but not necessary the same power of two), since this is the case which the fast Fourier technique handles most readily. Also, as it appears in 'grlib_pak', the 'fft' operation converts an image into a pair of images, one representing the real and the other the imaginary part of the image Fourier transform. To invert the Fourier transform one supplies the real part of the Fourier image as the 'main' argument of the 'inv_fft' operation, and the imaginary part of the Fourier transform (the second image component which it returns) as auxiliary parameter of this same operation. All these details are shown in the example.

The following example also displays the two components of the Fourier transform, to give a bit of their flavor.

programtest; -- test program 7 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral := orig_spiral := image("spiral.jpg").to_float().scale([256,256]); -- read the spiral image and make a copy [reimg1fft,imimg1fft] := spiral.fft(); restored := reimg1fft.inv_fft(imimg1fft); check_black(spiral - restored); -- check that an image is black Tk := tkw(); Tk(OM) := "Image Example 5"; -- create the Tk interpreter write_captioned_image(restored,"Restored version of spiral, thru fft and inv_fft"); write_captioned_image(reimg1fft,"fft of spiral, real part"); write_captioned_image(imimg1fft,"fft of spiral, imaginary part"); Tk.mainloop(); -- enter the Tk main loopprocedurecheck_black(the_image); -- check that an image is black; otherwise write it out sorted := the_image.sort();ifsorted(1) >= -0.001andsorted(#sorted) <= 0.001thenelsenot. ******* Written out as result_image",output_image_count + 1); write_image(the_image); -- write image to an incrementing fileend if;endcheck_black;endtest;

Our next example shows the use of the Fourier transform to calculate a convolution of an image by a kernel that has about as many terms as the image has pixels.

**The Wavelet Transform**.
Like the Fourier transform, the wavelet transform of an image is an invertible, linear, image-to-image transformation. The version of this transform (and of its inverse) available in 'grlib_pak' is just as fast as the Fourier transform and, like the Fourier transform, applies to images with two **k rows and **j columns. Although it does not have anything like the Fourier transform's very special relationship to the convolution operator, the wavelet transform can be useful for image compression. To see why this is so, we need to say a bit about the mathematics of this transformation, which it is easiest to do in the 1-dimensional case rather than the 2-dimensional image case. In this case, i.e., when applied to a vector w of length 2n the wavelet transform can be thought of as the result of iterating a sequence of related operations called *preliminary wavelet transforms*. The first preliminary transform in this sequence maps a vector of length 2n into a pair of vectors w_{1}, w_{2}, of length n, and has the following special property: If the vector w being transformed is linear with a jump (that is, if its i-th component w_{i} has the form a * i + b up to a certain point and then drops to 0, the vector w_{1} has just a few non-zero coefficients, while w_{2} is a vector of half the original length, which itself is almost linear with a jump. The full wavelet transformation is obtained by mapping an initial vector w by a preliminary wavelet transform into the pair w_{1}, w_{2}, then mapping w_{2}, which is of half the length of w, by a preliminary wavelet transform for vectors of half length, and so on until the last remaining vector is of length 2. If the vector w with which this process begins is 'image-like', i.e. is well approximated by a piecewise-linear function with jumps, then the final wavelet transform can be expected to have a relatively few substantial components (in the ideal case considered in the preceding sentences, this number will be roughly logarithmic in the length #w of the original vector w.) Hence we can compress w by taking its wavelet transform, a then either dropping all of its sufficiently small coefficients or quantizing them very crudely. Note that the wavelet transform can be expected to handle jumps in a vector (i.e. edges in an image) well; Fourier compression techniques do less well in this regard, since the presence of jumps affects many coefficients in a Fourier transform, rather than just a few.

The following example shows the use of the wavelet transform and its inverse. We display the wavelet transform image, simply to show hov sparse it is compared to the Fourier transform, as displayed in the preceding example.

programtest; -- test program 8 for image classusegrlib; -- use the basic image libraryuseimage,random_pak,image_pak,string_utility_pak; -- use the image class, randomization, image_pakusetkw; -- use the main widget class var Tk; -- global holding the main Tk window spiral := orig_spiral := image("spiral.jpg").to_float().scale([256,256]); -- read the spiral image and make a copy spiral_wave := spiral.wavelet(); restored := spiral_wave.inv_wavelet(); check_black(spiral - restored); -- check that an image is black Tk := tkw(); Tk(OM) := "Image Example 5"; -- create the Tk interpreter write_captioned_image(restored,"Restored version of spiral, thru waveletand inv_wavelet"); write_captioned_image(spiral_wave,"Wavelet transform of spiral"); Tk.mainloop(); -- enter the Tk main loopendtest;

**The image_pak Library**

'image_pak' is a utility library of image analysis operations built on top of 'grlib_pak' and 'image object class'. This library includes about thirty image operations. Many (but not all of these operations) are designed for use with 'binary' images. That is, byte images each of whose pixels has the color value 0 or 255. The 'erode' operation erodes a binary image by removing all image pixels which contact the outside of the image. 'erode' is available in two addition versions, one of which removes all image pixels which touch a pixel outside the image in the horizontal direction and the other one of which does the same in the vertical direction. There is also a 'dilate' operation which can be thought of as an 'erode' operation which applies to the black rather than the white pixels in a binary image.

Another important class of utility operations available in the image_pak library are various 'feature finders'. These operations attempt to enhance image regions containing significant elementary features by raising pixel values in these regions and lowering pixel values in all other image regions. If this can be done successfully then the features can be pinpointed by subsequent thresholding operations, followed by the extraction of connected components. The aim, which has been intensively pursued for forty years (with, alas, only moderate success) is to develop computerized feature finding procedures which come close enough to the performance of the human eye to support sophisticated image-based applications. The collection of operations of this type found in 'image_pak' includes the widely used Sobel edge-finding algorithm, a corner-finding technique, specialized techniques for finding edges which run in horizontal and vertical directions, an edge finder that enhances image regions in which the image varies rapidly, and several other experimental edge finders.

A third class of utility operations included in 'image_pak' are the so-called image 'filters'. These are operations which enhance an image by removing various detectable forms of clutter from it, which enhance the contrast in an image, subtract background noise or which put the image into one or another stylized form, for example, give it an embossed appearance.

The full list of operations provided by 'image_pak' is given in its specifier, and the way that these operations work can be seen by examining the source code for this package which is distributed with the SETL system. Here we simply give a few sample routines from this collection, namely, the 'erode' procedure, the Sobel edge finder, the 'sharpen_limit' filter used to remove dust-like artifacts from an image and two feature finders.

The two basic routines 'compress' and 'uncompress' respectively apply the Lempel-Ziv method and its inverse to a string. 'open_zip' generates a special form of zip-file handle, used in the remaining operations, from a valid Zip-file name. Files opened in this way shoud be closed using the 'close_zip' routine. 'list_zip(zipobj)' produces a tuple of tuples, each of the formnative packagezip_pak;procedurecompress(string); -- compresses a string using a variant of the Lempel-Ziv methodprocedureuncompress(string); -- inverse of the 'compress' operationprocedureopen_zip(file); -- returns a 'zipobj' handle to a Zip file, used in the operations belowprocedurelist_zip(zipobj); -- return a tuple describing contents of a Zip fileprocedureget_from_zip(zipobj,item_name); -- extract item from Zip file, as stringprocedureextract_from_zip(zipobj,item_name); -- extract item from Zip file, to a file ofthe same name as item_nameproceduredebug_zip(zipobj); -- print standard Zip file contents reportprocedureclose_zip(zipobj); -- close a Zip file handleendzip_pak;

The item_names appearing here are used in the two following operations

which find and decompress the named item from the zip file to which their first parameter must be a handle produced by 'open_zip'. 'get_from_zip' returns the decompressed result as a string, while 'extract_from_zip' rites this same result to a file of the same name as 'item_name'. The auxiliary routine debug_zip(zipobj) prints a standard Zip file contents report. Note that no tool for compressing a directory into a Zip archive is provided, since it is assumed that this will be done using some externa Zip tool, e.g MacZip on the Macintosh, or WinZip under Windows.

The following short program illustrates the use of these operations.

programtest; -- zip file manipulion testusezip_pak; ziphandle := open_zip("examples_folder.zip"); -- open a zip fileendtest;

Now we begin to give the layer of detail needed by would-be implementers of additional **native package**s. This will make it necessary to describe parts of the internal structure of the SETL2 system and the way in which **native package**s are implemented.

A variety of issues arise in connection with SETL-to-C communication

- The executable C code implementing the functions provided by a
**native package**is loaded dynamically (so that large libraries of such packages can be managed comfortably) rather than being linked in statically. For this, SETL relies on the dynamic linking capabilities of the operating system under which it is running. But it must also know where these code files are found, and where the relevant entry points in them are located. - For efficiency, most all of the manipulation of objects created by
**native package**s will be handled by the C routines that these packages provide. SETL will often see these objects only as 'opaque' objects resembling atoms,to which it can apply only those few operations available for atoms. All other operations on such objects will be performed by passing them to appropriate routines of their originating**native package**s, and receiving similar values in return. But at certain points in this process values directly comprehensible to SETL, such as integers, floating values, and strings may need to be passed. And for efficiency we may even want to pass SETL composites back and forth, e.g. tuples of floats from which the C code of a**native package**will construct real arrays, or tuples of floats which**native package**C code has constructed from some internally generated real array in order that it may be passed back to SETL. - For a
**native package**to be sure that an opaque object passed to it by the SETL is correctly structured, these objects need to carry type tags extending SETL's internal object tagging system. - SETL, like many other high-level languages, manages memory automatically, but C relies on manual memory management via its 'malloc' and 'free' functions. To bridge this gap, native packages must be linked into to SETL's reference-count based memory management system. This is done by requiring native objects to carry reference counts in the positions in which SETL expects them, and by providing several macros which C routines can use when they need to increment or decrement such reference counts directly (for example, if a multilevel SETL object needs to be built-in C.) All other management of reference counts is left to SETL, which will increment and decrement these counts in its normal way as opaque objects originating with
**native package**s are assigned as the values of SETL variables or made members of sets and tuples.When the reference count of an opaque object falls to 0, all the space which it occupies must be reclaimed. Since SETL is only aware of a small part of this space and does not even know how the bulk of it is structured (for example, how is an image, video file, or Tk widget laid out?), it must rely on the

**native package**which generated the object to perform the bulk of this deallocation. For this reason,**native package**s are required to register*destructor routines*with SETL immediately after they are loaded. - SETL provides various objects of complex structure (its tuples and sets, and also strings, which SETL represents in a form devised to make insertions and concatenations efficient.) C provides essentially only one kind of composite object of indefinite length, the array. Thus a problem arises whenever C must receive a composite object from SETL or prepare such an object for return to SETL. We do not want to force the author of the necessary interfacing code to learn numerous details of the SETL implementation. To avoid this, we provide a family of C macros which represent SETL compound objects (e.g. sets and tuples) by iterations over them which C can understand, and a converse set which C can use to construct sets and tuples (also strings) by iteratively adding elements to them. Each of these families of macros is structured as (i) a
*declaration*macro, which serves to introduce the required iterator or constructor along with a C variable name for it; (ii)*block header*and*end*macros, which delimit the block of C code which will iteratively analyze or construct a compound SETL object. (iii) element*access*and*insertion*macros, for use within such iterator blocks.A few additional utility macros, which get the length or cardinality of a SETL object for use in C, or copy a SETL string to a C buffer, are also provided.

Note that the set of macros provided embody a procedural, rather than a declarative, approach to the problem of SETL to C communication.

The facilities described in the following subsections are designed to deal with the complex of issues that we have just reviewed.

Once a **native package** is loaded and all its functions resolved, the
runtime system will call a 'package initialization function'.

The name of this package-associated function is constructed by appending the suffix '__INIT' to the package name. If a function with such a name is found in the DLL, it is invoked as soon as it has been loaded.

This '__INIT' function typically contains package-specific initializations, as well as code used to register the package types with SETL. A simple example of the C code pattern used is

int32 TK__INIT( SETL_SYSTEM_PROTO_VOID) { tcl_type = register_type(SETL_SYSTEM "tcltk",internal_destructor); if (tcl_type == 0) return SETL_ERROR; return SETL_OK; }

In the code seen above, the 'register_type' routine is called by the package
initialization function. The name of the package 'tcltk' is passed
for use in error reporting by the SETL runtime system. The type
is saved, in this case in a static variable, so that the **native package**
can later check if the objects passed to its functions are of the correct
type.

A destructor routine is also passed to the SETL runtime system. Whenever the reference count of an object drops to 0, the destructor registered for its type will be called. The important conventions which apply to destructors are detailed in Section XXX below.

No function to copy the object needs to be supplied, since the SETL runtime
will only add and remove references to the native objects, but these objects
can only be created and modified by the corresponding **native package**.

struct specifier_item { unsigned sp_form; /* form code */ union { int32 sp_atom_num; /* unique atom identifier */ int32 sp_short_value; /* value of short item */ struct file_item *sp_file_ptr; /* file node pointer */ struct instruction_item *sp_label_ptr; /* label value code pointer */ struct opaque_item *sp_opaque_ptr; /* opaque item */ struct proc_item *sp_proc_ptr; /* procedure */ struct integer_h_item *sp_long_ptr; /* header of integer */ struct i_real_item *sp_real_ptr; /* real number */ struct string_h_item *sp_string_ptr; /* header of string */ struct set_h_item *sp_set_ptr; /* root of set header */ struct map_h_item *sp_map_ptr; /* root of map header */ struct tuple_h_item *sp_tuple_ptr; /* root of tuple header */ struct iter_item *sp_iter_ptr; /* iterator */ struct object_h_item *sp_object_ptr; /* object (user-defined) */ struct mailbox_h_item *sp_mailbox_ptr; /* mailbox header */ void *sp_biggest; /* dummy -- it must be as large as */ /* the largest of the above */ } sp_val; };

The way in which a **native package** must begin to access the body of a native object given a pointer to a specifier for it is therefore

specifier_ptr->sp_opaque_ptr

This much is forced by SETL. Once the body block of a native object has been accessed in this way, the rest depends on the layout of the block, which is determined by the **native package** rather than by SETL (which constrains only the first two words of such blocks, as explained just below.) An example, the way in which blocks representing matrices are handled by LIBNR_PAK, is given below.

The first field of the specifier structure seen above is the *form code* determining the type of the object. For native objects (those maintained by native packages; SETL keeps only references to such objects, which it is only able to treat as opaque atoms) the form code 'ft_opaque' is used.

Simple types do not use the reference-count used to manage compound SETL objects, since the cost of maintaining such counts it will be greater than the cost of copying small objects. But memory management of compound types, including sets, maps, and native objects, uses a reference-count scheme fundamental to the SETL run-time system. The reference counts needed for this are always held in the first word of the block to which the pointer that appears in the specifier of every compound object points. For native objects, such a block must therefore begin in the following way.

struct setl_native { /* Native Object Structure */ int32 use_count; /* Reference Count */ int32 type; /* Encodes Type and Subtype */ ...... };

Whenever a **native package** creates a native object which it will manage, it must allocate and initialize such a block, create a specifier which points to it, and pass this specifier back to the SETL runtime system. The type of the native object is
encoded in the lower 16 bits of the *type* field seen, and an optional
subtype (for use by the **native package**, notby SETL) can be stored in the higher order bits.

The use_count field must only be manipulated by two C macros called *mark_specifier* and
*unmark_specifier*.

To show how all this comes together, we give an example: the 'matrix' object type allocated and managed by the LIBNR_PAK
**native package** described above. This has the following C declaration:

struct nrmatrix { /* 'body block' for LIBNR_PAK matrices */ int32 use_count; /* required by SETL */ int32 type; /* required by SETL */ int r,c,h; /* dimensions; specific to LIBNR_PAK */ void *p; /* data pointer; specific to LIBNR_PAK */ };

The structure seen here is used to store the vectors, matrices, and tensors
of various kinds used by the numerical library. Each object's type (as seen by this library; SETL sees them only as 'opaque') is
stored in the high order 16 bits of the *type* field.
The field *p* contains the pointer to the actual vector, matrix, or tensor body area set up by the memory allocation
routines of the underlying Numerical Recipes library. We use r,c,h to hold
the dimensions of the vector, matrix, or tensor object, which SETL may need to query.

For example, when a *float-vector* (subtype *nr_fvect*) native object is allocated, the
structure seen above is initialized using the following C code:

mat = (struct nrmatrix *)(malloc(sizeof(struct nrmatrix))); /* allocate structure seen above */ mat->type = nr_fvect * 65536 + nr_type; /* set up type field */ mat->use_count = 1; /* initialize reference count */ mat->p = (void *) vector(1,rows); /* allocate C data array */ mat->r = rows; /* initialize dimensions */ mat->c = 1; /* initialize dimensions */

To avoid memory leaks which might otherwise bring the SETL system to its knees, all native objects which can no longer be referenced must be deallocated. But since their internal structure is not known to the SETL system, this must be done by the **native package** which created them in the first place. For this to happen, the **native package** must make a suitable *destructor*
available for call by SETL. When called, this destructor will deallocate the memory used by the Native object.
If a destructor is not registered during the package Initialization phase,
the SETL runtime system will only release the memory occupied by that small part of the
native object structure known to it, but not the (possibly very large amounts of) memory that may have allocated internally by the **native package** during creation of the object.

The conventions employed in writing destructors are seen in the following example, taken from LIBNR_PAK:

void internal_destructor(struct nrmatrix *body_ptr) { int subtyp; subtyp = body_ptr->type >> 16; /* find object type */ switch(subtyp) { case nr_fmat: free_matrix(body_ptr->p,1,body_ptr->r,1,body_ptr->c); /* free a single-precision matrix */ break; case nr_dmat: free_dmatrix(body_ptr->p,1,body_ptr->r,1,body_ptr->c); /* free a double-precision matrix */ break; ... } }

Here, different deallocation routines must be called, depending on the subtype of the native object being deallocated

All the functions defined in SETL **native package**s have the
prototype seen below:

void NR_GET_VECTOR( SETL_SYSTEM_PROTO int argc, /* number of arguments passed */ specifier *argv, /* argument vector */ specifier *target); /* return value */

Here, SETL_SYSTEM_PROTO is a C macro which creates an encoding of and reference to the current SETL environment (program state) when multiple SETL threads are running. Next follow the number of arguments to be passed, the array of these arguments, and a pointer to the C structure to be used to store a return value when the native function returns.

At the moment of call, the argument vector will contain pointers to SETL specifier structures (as described above). The C code can access these values using an extensive set of C macros provided with the SETL system, without knowing too many details about SETL's internal data structures. These macros are described in Section XXX below.

When ready to return, the C code should prepare its return value for receipt by SETL in the following manner.

unmark_specifier(target); /* decrement reference count */ target->sp_form = ft_opaque; /* set up specifier tag */ target->sp_val.sp_opaque_ptr = (opaque_item_ptr_type)mat; /* make specifier point to body block */ return;

Note that the C code must remove a reference to the return value before modifying it.

If a routine in a **native package** has read-write parameters (which is allowed), the (perhaps opaque) objects to be returned to SETL via these parameters must be pushed onto SETL's program stack in the manner seen in the following example:

specifier return1; ... return1.sp_form = ft_opaque; /* set up specifier tag */ return1.sp_val.sp_opaque_ptr = (opaque_item_ptr_type)fmatA; /* make specifier point to body block */ push_pstack(&return1); /* push pointer to specifier onto argument stack */

The header file *macros.h* supplied with the SETL runtime system,
provides a library of macros which can be used to build **native package**s accessing
compound SETL values without knowing (too) much about SETL's internal
data structures.

The macros in this collection fall into five functional groups: declarations, constructors, iterators, retrieval macros, and utility functions.

In this group we find macros that can be used to declare iterators and constructors for SETL sets, strings and tuples which a C routine may need to build (typically for analysis of a compound value transmitted by SETL, or for construction of a compound value to be returned to SETL).

Each of these macros has one parameter, the name of the iterator or constructor to be created. This name is be used later to operate on certain variables implicitly declared by the macro. The macros in this set are

STRING_CONSTRUCTOR(source)

TUPLE_ITERATOR(name)

TUPLE_CONSTRUCTOR(source)

SET_ITERATOR(name)

SET_CONSTRUCTOR(source)

The following C code fragment, taken from the implementation of the *join* function provided by STRING_UTILITY_PAK, illustrates the way in which these macros are ordinarily used.

void JOIN( SETL_SYSTEM_PROTO int argc, /* number of arguments passed */ specifier *argv, /* argument vector (two here) */ specifier *target) /* return value */ { TUPLE_ITERATOR(ti) STRING_ITERATOR(sa) STRING_ITERATOR(sb) STRING_CONSTRUCTOR(cs) ...Here the first iterators is used to access the string components of the first argument of JOIN (which is a tuple of strings), while the second serves for iteration over the second argument of JOIN (which is a string).

STRING_CONSTRUCTOR_ADD(name,char)

TUPLE_CONSTRUCTOR_BEGIN(name)

TUPLE_ADD_CELL(name,specifier_ptr)

TUPLE_CONSTRUCTOR_END(name)

SET_CONSTRUCTOR_BEGIN(name)

SET_ADD_CELL(name,specifier_ptr)

SET_CONSTRUCTOR_END(name)

The macros used for building sets and tuples are very similar.To build an object of either of these types, we must enclose the block of C statements used to assemble its elements between a pair of the appropriate BEGIN and END macros. For example, to building a set using a SET_CONSTRUCTOR known to our C code as 'cb', we would write code around the following skeleton:

SET_CONSTRUCTOR(cb) /* declaration of a SET_CONSTRUCTOR */ ... SET_CONSTRUCTOR_BEGIN(cb); /* open SET_CONSTRUCTOR block */ ... SET_ADD_CELL(cb,&tgt); /* add set element */ ... SET_CONSTRUCTOR_END(cb); /* close SET_CONSTRUCTOR block */

Within a block like that seen in this example (i.e. between its BEGIN and END macro calls), we can use as many SET_ADD_CELL calls as we like. Note that the macro calls must supply the constructor name, and that SET_ADD_CELL must also supply a pointer to a specifier.

Construction of SETL strings is somewhat different, since it does not use an END macro. After an opening STRING_CONSTRUCTOR_BEGIN, we simply call STRING_CONSTRUCTOR_ADD for each character that we want to append to a string. The pattern is as follows

STRING_CONSTRUCTOR(cb) /* declaration of a STRING_CONSTRUCTOR */ ... STRING_CONSTRUCTOR_BEGIN(cb); /* open STRING_CONSTRUCTOR block */ ... STRING_CONSTRUCTOR_ADD(cb,c); /* add character */ ...

ITERATE_STRING_NEXT(name)

ITERATE_STRING_CHAR(name)

ITERATE_TUPLE_BEGIN(name,specifier_ptr)

ITERATE_TUPLE_END(name)

ITERATE_TUPLE_CELL(name)

ITERATE_SET_BEGIN(name,specifier_ptr)

ITERATE_SET_END(name)

ITERATE_SET_CELL(name)

Here is an example of how we would use the macros in this group to access a tuple of integers passed in from SETL.

TUPLE_ITERATOR(ti) /* declaration of a TUPLE_ITERATOR */ ... ITERATE_TUPLE_BEGIN(ti,argv[2]) /* begin iteration over tuple */ tuple_element++; if ((tuple_element > 4) || (ITERATE_TUPLE_CELL(ti)->sp_form != ft_short)) { /* use 'current' tuple component specifier set up by iteration macro */ abend(SETL_SYSTEM "Score Tuple in EDIST must have integer elements"); } ... ITERATE_TUPLE_END(ti) /* end of tuple iteration block */

The block bounded by the ITERATE_.._BEGIN and ITERATE_.._END macros is executed once for each element of the object, and sets up a 'current specifier'. We can access this current specifier using the ITERATE_TUPLE_CELL (or , for iteration over a set, the ITERATE_SET_CELL) macro.

The analysis of strings passed in from SETL is a bit different: we use a macro to initialize the iteration (ITERATE_STRING_BEGIN) and macros to access the current character (ITERATE_STRING_CHAR) and to move to the next one (ITERATE_STRING_NEXT). But in this case there is no block ender.

Once a set, tuple, or string has been built using the macros described in Section 10.6.2, we can use the following macros to build a specifier, which can then either be passed back to the SETL, or inserted into another compound object.

TUPLE_HEADER(name)

SET_HEADER(name)

For example to return a newly constructed tuple to SETL we could flesh out the following skeleton:

TUPLE_CONSTRUCTOR(tc) /* declaration of a TUPLE_CONSTRUCTOR */ ... TUPLE_CONSTRUCTOR_BEGIN(tc); /* open TUPLE_CONSTRUCTOR block */ ... TUPLE_ADD_CELL(tc,&tgt); /* add component */ ... TUPLE_CONSTRUCTOR_END(tc); /* close TUPLE_CONSTRUCTOR block */ ... target->sp_form = ft_tuple; target->sp_val.sp_tuple_ptr = TUPLE_HEADER(tc); return;

TUPLE_LEN(specifier_ptr)

SET_LEN(specifier_ptr)

STRING_CONVERT(name,char_buffer_ptr)