The Simple Algebraic Math Library (SAML) is a low-level library to handle the classic objects of symbolic calculus: arbitrary big integers, reals of arbitrary precision, rationals, polynomials, tensors, et cetera. This document describes the internals of the library and other implementation details.
The SAML package is free; this means that everyone can use it and redistribute it on a free basis. Some files in the SAML archive were not written by me, and have their own copyright notice -- they are included only for the convenience of the installer. Some files are so trivial that they've been explicitly released in the Public Domain. All other files, including this documentation, are copyrighted by me (Thierry Bousch) and redistributable according to the GNU Public Licence. Read the file `COPYING' for the details of this license.
The following program creates two variables, stores 123 and 456 into them, adds the second variable to the first one and displays the result.
/* File example.c */
#include <stdio.h>
#include "saml.h"
int main (void)
{
mref_t v1, v2;
/* Create two mrefs */
v1 = mref_new();
v2 = mref_new();
/* Store the initial values */
mref_build(v1, ST_INTEGER, "123");
mref_build(v2, ST_INTEGER, "456");
/* v1 := v1 + v2 */
mref_add(v1, v1, v2);
printf("The sum is %s\n", mref_string(v1));
/* Free the mrefs */
mref_free(v1); mref_free(v2);
return 0;
}
The two mref_free() calls at the end of the program could have been
omitted, since the process will terminate and release all resources. They
would be important if this routine was part of a larger program.
Now compile and run the program:
$ gcc -o example example.c -lsaml $ ./example The sum is 579 $
You'll easily build more impressive examples, for instance by using rationals instead of integers:
mref_build(v1, ST_RATIONAL, "-777777777777777");
mref_build(v2, ST_RATIONAL, "123/456");
Note that mref_build() with type ST_RATIONAL will
automatically create rationals over the integers, whereas the
ST_RATIONAL type allows rationals over any ring.
For polynomials, the situation is not so nice; it is not clear whether
an expression like 3*xy+1 is a polynomial with integer or
rational coefficients (or worse), and it is not clear either whether
xy is a literal or the product of the two literals x and
y. Thus, the "build" operation is not supported by the
ST_POLY type, because it is not very well defined. You have to
build polynomials "by hand", first creating literals and then promoting
them to polynomials. This is obviously painful, fortunately SAML provides
a simple parser for algebraic expressions (in files `simple-lexers.c'
and `simple-parser.c') which can be used for polynomials. You still
need to create an object of the right type (a polynomial with integer
coefficients, say) before calling the parser. Our example would be:
v1 = mref_new(); v2 = mref_new();
mref_build(v2, ST_INTEGER, "0");
mref_cast(v2, ST_APOLY);
/* Now v2 contains the null polynomial */
saml_init_lexer_mem("3*xy+1", 6);
saml_parse(v1, v2);
/* Now v1 contains 3*xy+1, v2 is unchanged */
saml_init_lexer_mem("x-y", 3);
saml_parse(v2, v2);
/* And v2 contains x-y */
The parser assumes C-style identifiers (although FORM-style identifiers
like [foo-bar] are also accepted) and all operators (sum,
product) must be explicit. The string length argument is mandatory,
because this parser should also be able to work with strings which are
not zero-terminated, or could contain embedded zeros.
There are two basic data types used by SAML. The most fundamental one is the mnode (math node). Just likes inodes on the UNIX system, these objects are invisible to the user, and are only used internally.
Mnodes are described in the header file `mnode.h'. Each mnode consists of a type (represented by an integer), a reference counter, and type-specific data. These data can include pointers to other mnodes. There must be no cyclic references; an mnode should only reference mnodes which have been created before.
The type field is used for dispatching; when you want to add two
polynomials, SAML checks that both operands have the same type, and
calls the type-specific operation, in this case poly_add(). This
function will, in turn, invoke mono_add() to group similar
monomials, which will in turn call the generic function
mnode_add() to add the coefficients.
When you add two polynomials, there will be many common monomials between the two operands and their sum. This is why a reference counter has been added, to avoid unnecessary copies of objects. When a mnode is no longer needed, the library decrements the counter, and frees the object when the counter reaches zero (once again, it is very similar to hard links on a UNIX filesystem). If this mnode refers to other mnodes, they should have their reference counters decremented as well, and so on, recursively.
Error conditions (as the result of a calculation), and uninitialized
objects, are also implemented as mnodes of the special type
ST_VOID. The advantage of this scheme is that all type-specific
operations on mnodes, like poly_add(), can always assume that
their operands are valid, since the dispatcher will never call this
function if one of the operands doesn't have the type ST_POLY.
Mnodes must be allocated and freed explicitly. This is error-prone and can easily lead to crashes or memory leaks. Therefore, it is important to use mnodes only internally, and to expose a higher-level interface which takes care of allocation details and guarantees that user code can't crash the library.
This is done with mrefs (math references). If mnodes are similar to inodes, then mrefs are similar to file descriptors on UNIX. They are identified with a small integer, and must be created and deleted explicitly. But unlike mnodes, which represent a constant object, mrefs are really variables: you can assign to them, do operations on them, without any reference-counting hassle. Also, the library is able to check that the mrefs are valid, whereas it wouldn't be possible if objects were referred to by pointers. Thus, mrefs provide some simple form of garbage collecting: there can be no dangling pointers (since you don't have access to pointers) and mnodes can leak only if you don't keep track of your mrefs.
It is strongly recommended to declare mrefs as objects of type
mref_t (which is only a typedef for int), especially in
function prototypes. Otherwise you have annoying ambiguities -- see for
instance the prototype of mref_power() below. It also avoids the
horrible "hungarian notation" where one encodes the type into the
variable names, and calls some variable `mra', say, instead of the
simpler `a'.
These are math types whose mnodes don't hold references to other mnodes.
In other terms, they are completely self-contained. For now, the
following simple math types are implemented: ST_VOID,
ST_LITERAL, ST_MINT, ST_INTEGER, ST_CYCLIC
and ST_FLOAT.
ST_VOID represent either errors, or uninitialized
objects. They are couples (errno,where) where errno is
a positive error number (or 0 for an uninitialized object) and
where is a short string indicating in which routine the error
originated. Obviously, there are almost no operations possible on these
objects, and the dispatcher treats them specially anyway.
These math types are only "patterns" allowing to build complicated
types from simple ones. For instance, you can create fractions over any
ring without divisors of zero, and polynomials over any ring. For now,
the following parametrized math types are implemented:
ST_RATIONAL, ST_COMPLEX, ST_MONO, ST_POLY,
ST_UPOLY, ST_LINE, ST_MATRIX and ST_TENSOR.
mnode_gcd() to allow simplification of fractions.
ST_APOLY, the code is clumsy and doesn't
handle exceptional conditions, like exponent overflow (which is
ignored). On the other hand, for some problems this representation can
be much more efficient, so I don't want to throw it away entirely.
ST_POLY representation.
samuel(1) briefly explains them.
In previous versions of the library, you had to call saml_init()
before any other operation. Otherwise, the program would crash at the
first invocation of mref_new(). This is why you'll find this
function call at the beginning of some programs; but it's useless now,
the initialization will automagically take place at the first invocation
of mref_new().
ST_VOID, should always be initialized.
All library functions operate on mrefs (variables). New mrefs are
created with mref_new(); the return value is a small non-negative
integer which identifies the mref. When the mref is no longer needed,
you can delete it with mref_free(id) where id is the
value returned by mref_new(). New mrefs are initialized with the
special value "void", which really means "uninitialized".
mref_new(), or if it has
already been freed.
All operations taking mref identifiers as arguments will abort()
if any of these identifiers is invalid, with a beautiful error message.
We won't precise this every time.
free() these
strings yourself.
saml-mtypes.h. It would be better to call this the
"superficial type" of id because it only reflects the type of the
topmost object; thus, all polynomials have type ST_POLY, no matter
the type of their coefficients.
saml-errno.h.
ST_POLY, then id will be cast to
a polynomial with rational coefficients.
mref_cast(id,ST_POLY) would have given a polynomial with
integer coefficients.
Returns zero on success, a negative error number on failure.
INT_MAX).
Returns zero on success, a negative error number on failure.
saml_init_lexer_fd() or saml_init_lexer_mem(),
otherwise, the behaviour is undefined. The file descriptor or memory area
will be parsed and the result put into dest; it will have the same
type as the model. This function returns 0 on success and a negative
error number on failure.
The following operations are meaningful only on certain types (polynomials, tensors). Like the other ones, they return 0 on success and a negative error number on failure.
ST_MATRIX) is stored into dest.
ST_TENSOR and
ST_LITERAL. If the tensor doesn't contain the first literal, an
error is returned. This function is similar to mref_subs() but
must be different, for namespace reasons, since the same literal could
be used both as a tensorid and as a polynomial indeterminate.
We have already seen how most mref functions return a negative value when an error occurs. However, it is tedious (and error-prone) to test all return values, so other mechanisms exist.
The mnode layer, in the case of an error, returns a mnode of the special
type ST_VOID, which contains the error number and the name of the
routine where it occurred. The mref layer detects this condition, calls
the error handler with two arguments and (unless this handler caused the
program to abort) returns -errno to the caller. This
handler is declared in the `saml.h' header file as follows:
extern void (*saml_error_handler) (int reason, const char* where)
The library provides three error-handlers as examples, but the user can install his own. Here they are:
saml_errors.
saml_errors. This is useful when you
expect errors (for example, to verify user-supplied data)
and you'd like to print a more informative error message; this is what
factorint(1) does to detect malformed numbers.
Objects of type ST_VOID (uninitialized variables and run-time
errors) are also handled in a special way by the mnode dispatcher. When
a mref function takes several arguments and these arguments have
conflicting types, the dispatcher will return a "type mismatch" error
unless one of the arguments has type ST_VOID; in this
case, this argument will be the return value. In other words, errors
will be propagated without alteration (unless they meet another error,
then we'll have to drop one of them). This behaviour is similar to NaN for
floating-point numbers: errors are "sticky" and can't disappear.
NULL in
the s_mtype structure.
mnode_cast().
register_mtype().
mnode_stringify() could not be transformed into something of the
given type.
ST_APOLY polynomials can't exceed
2^32-1), and this error will be raised if the result can't fit.
This is probably better than silently ignoring the rollover.
Resizable strings are extensively used by SAML, and since user programs
may need them, I have included them in the public interface. Unlike C
strings, they're not necessarily NUL-terminated and can contain embedded
NULs; the length is part of the data structure. The gr_string
structure is defined in `saml-util.h', as follows:
typedef struct _growable_string {
size_t buflen; /* allocated length */
size_t len; /* current length */
char s[0]; /* data */
} gr_string;
It is OK to manipulate these objects directly, provided that len (the current length) never exceeds buflen (the allocated length). To create them, or increase their length, you should use the following functions, also prototyped in `saml-util.h':
gr_string, which can be different from dest if a
reallocation took place.
Important: Never forget that a gr_string can be relocated
to a new address if you add data to it. Thus, you should never ignore the
return values of the above functions. Also, don't forget to put an explicit
'\0' at its end if you want to use it as a normal C string.
To conclude, let's notice that it's very easy to shrink an existing growable
string; just reduce the len field. In particular, setting this field
to zero will clear the string. This trick is used several times in the
sample programs.
SAML provides a simple parser for the common cases where you want to read symbolic expressions from a file or an area of memory. To use it, you must first initialize the lexer, with one of the following functions:
@noindent{}and then you can call saml_parse().
The lexical conventions are simple: identifiers can be either C-style
identifiers like __foo_1 or FORM-style bracketed identifiers like
[z^2+c]. Whitespace is ignored. The complete grammar can be found
in the file simple-parser.c.
[ Some remarks about the parser ]
Mnodes are only visible inside the library. To use them, you have to include `mnode.h'. Every mnode begins with the following header, so a pointer to a mnode (which has type `s_mnode*' or `mn_ptr' for brevity) is identified with a pointer to its header:
typedef struct mnode_header {
int type; /* superficial type */
int refs; /* reference counter */
} s_mnode, *mn_ptr;
The first field contains the so-called "superficial" type, used for dispatching. The second field counts how many times the mnode is referenced -- the mnode is deleted when this counter reaches zero. But this is only a template, all mnodes will contain other data appended to the above header. The format of these data is type-specific; there can be integers, references to other mnodes, whatever. So-called standard mnodes are simply arrays of references to other mnodes; they are defined in `mnode.h' as follows:
typedef struct {
struct mnode_header hdr; /* see above */
int length; /* length of the remaining data */
mn_ptr x[0]; /* remaining data (appended) */
} std_mnode, *smn_ptr;
With these types you can use the array operations of the public interface:
mref_array(), mref_getitem() and mref_setitem(), but
be careful: SAML won't warn you if you attempt these operations on
non-standard types, and even standard types can expect class invariants,
which can be violated with these functions.
Each type, standard or not, is identified by a small non-negative
numbers; builtin types are defined in `saml-mtypes.h'. Each type has an
initialization routine, usually called from saml_init(), which is
responsible for registering the new type and conversion routines, and
maybe initialize some private data structures, like hash tables.
If you create a new type, you should insert its initialization routine
in the saml_init() routine.
s_mtype or more probably unsafe_s_mtype.
These two data structures are bitwise equivalent, but the second one is
more laxist on type verifications, which is often desirable.
The s_mtype structure, defined in `mnode.h', describes the
operations supported by this type. A NULL value means that some
operation is unsupported. Since these functions are always called via
the dispatcher, the implementor can assume that all type checking has
already been done, and the arguments have the right types. They should
return a new mnode with the same type (or a new reference to an already
existing mnode) or an error mnode.
Not everything is stored in the s_mtype structure. Type conversion
is particular because its requires a double dispatching (depending on the
initial and final type), and thus is implemented in a special way.
You can register a conversion routine with
-1 as the first argument (generic conversion routines are
supposed to take an argument of any type, and convert it to type
t2).
Assume that t1 and t2 correspond to types foo and var. The prototype of the conversion routine should be as follows:
NULL
or the address of a mnode of type bar. It should return a new
reference to an object of type bar or an error.
Generic conversion routines are usually not needed, except in rare cases where a model is needed for the conversion (matrices are an example of this, because you must know their size to convert a number to a scalar matrix).
Special functions are the constructors and destructors. File `mnode.c' provides three allocators for new mnodes:
struct mnode_header part) and one reference.
mstd_free().
ST_VOID. The
first argument is the error code, the second argument should be the name
of the function where the error occurred.
To delete or duplicate existing mnodes, it is best to use the following (inline) functions, otherwise it gets messy:
Most routines (addition, for instance) return a new mnode, more
precisely, they return a new reference to a mnode. Keep this in mind
when you write operations for parametrized types: count how many
references you create (most operations will create one reference) and
how many you delete (with unlink_mnode()). If you unlink too many
references, you will sooner or later get a crash because of a dangling
pointer. On the other hand, if you forget to remove references to
temporary objects, they will never be freed and it is a memory leak. As
a debugging aid to track this kind of problems, the library provides
three global variables, called:
extern long nb_mnodes_allocated; extern long nb_mnodes_reserved; extern long nb_mnodes_freed;
These three counters are initially zero. The first counter is
incremented every time __mnalloc() is called (perhaps via
mstd_alloc()), and the last counter is incremented every time
destroy_mnode() is called (via unlink_mnode() normally).
The second counter should be incremented each time you allocate a mnode
you don't plan to ever release (for good or bad reasons); for instance,
the module `Integer.c' preallocates two mnodes for 0 and 1 in its
initialization routine, and so increments nb_mnodes_reserved by
two units. Reserved mnodes should be exceptional.
If all the reference-counting is correct, the counters will verify the following relation, after all mrefs have been freed:
nb_mnodes_allocated == nb_mnodes_reserved + nb_mnodes_freed
Yes, I have heard about automatic garbage collecting. It is nice for the
programmer but not so nice for the computer; in my experiments, the
garbage-collected version of samuel used three times as much
memory as the normal reference-counting version, and it was slightly
slower. You can't have your cake and eat it, too.
Most SAML routines expect arguments of the same type. You cannot directly
add an integer and a fraction; you must convert the integer to a fraction
first, and then call mref_add() or mnode_add(). At the mnode
layer, two functions allow you to convert objects:
It should be noted that mnode_promote() can be useful even if arg
and model have the same superficial type; you might use it, for
instance, to convert a polynomial with coefficients in some ring K to a
polynomial with coefficients in L, if there exists a natural map
form K to L.
The only problem with mnode_promote() is that it requires a "model",
i.e. a mnode which has already the right type. It's a chicken-and-egg
problem: how do you create the first object with this type?
This is why mnode_cast() is useful, it allows you to create new types
from simpler ones.
cast and promote work
This is not pretty nor easy to describe, but I'll try anyway; in any case
you should read `promote.c' for more detail. Conversion routines
between types are registered through register_CV_routine() and
stored into a dictionary whose keys are the (type1,type2) pairs.
Let type1 and type2 the original type and the destination
type, respectively. If there is a conversion routine from type1 to
type2, then we execute it (and if it fails we return the error
to the caller). Otherwise, if type1=type2 or type1 is
ST_VOID, we return an unmodified copy of the initial object.
Otherwise, we look for a "generic" conversion routine to type2;
if there is one, we call it, and if it succeeds we return the result.
Finally (i.e., if there was no generic conversion routine or if it failed)
then we call mnode_make() to "cast" the object to type2,
and report either success or failure.
This is hairy and probably more complicated then needed. More importantly,
it's basically wrong, as it supposes that cast and promote
are essentially the same operation. But in fact, the role of cast
is to create new types: for instance, I have an object of type A and
want to create a polynomial whose coefficients have type A.
The current scheme, for instance, doesn't allow you to create nested
"Unary polynomials" using cast().
Type conversion is a big mess and should be redone entirely. Also, some very important things are missing from the API. There's no way to get the numerator and denominator of a fraction, for instance, because it isn't a "standard" type (i.e., an array) like unary monomial. More generally there's no way to break a complex object like a polynomial or a tensor into simpler objects: there are no iterators. Sure, it would expose some implementation details, but I think it's unavoidable, in some cases it is really necessary.
I think the right way to fix these two problems would be to add to each "opaque" type (like "Rational"), a "broken" type which would be an array-type -- here, an array of length two, containing the numerator and the denominator. The broken type would allow no operations, apart from accessing the pieces, and conversion routines to and from the opaque type. This would be much cleaner than the current ad-hockery and would give us a safe way to peek into objects.
There should be some kind of "call by name" as in Objective-C to call arbitrary methods, but then we should choose a convention for passing arguments and the retrieving the result(s).
We should really implement mnode_info() someday. It's supposed to
get "hints" about the mnode from an upper layer of the library, for
example, what's its estimated complexity (this is important for the Bareiss
algorithm), are multiplications really expensive, should we avoid
denominators, etc.
For C programmers, mrefs are a great advantage over mnodes because they
handle the gory details of reference counting, and they're mutable
objects like most other C types. In the case of Python, these
considerations are not relevant, because the language does the
reference-counting automagically and has no problem with immutable
objects. Therefore, the class Mathnode defined by the
saml1 module is only a thin layer around mnodes, and the
semantics are pretty much the same; in particular, the objects are
immutable. The main difference is that "error" mnodes are reported as
exceptions instead of class instances.
saml1
This module defines the class Mathnode, and the functions
MnAllocStats(), MnArray() and resultant().
It also defines the symbolic constants ST_* corresponding
to the builtin SAML types, and the Saml*Error exceptions.
saml1.MnAllocStats()
saml1.MnArray(mtype, tuple)
saml1.resultant(poly1, poly2, lit)
self have been the first, second, or third argument?)
Mathnode class
Most of the interesting operations on mnodes are accessible as methods
of the class Mathnode. Usually the self argument will be
the first argument of the corresponding operation on mnodes. This class
also offers attributes to get the math-type and the reference counter of
the underlying mnode. Here is the complete list of attributes and methods
(note that all attributes are read-only, since the objects are immutable):
self.refs
self.type
self.stype
self.__members__
self.__methods__
Mathnode(mtype, string)
mnode_build(), and is supported only for "simple" math-types like
integers, rationals or literals. More complicated types (like arrays,
polynomials and tensors) can't be obtained this way; you must use
either MnArray() to create an array type, or cast() or
promote() a simpler object.
self.cast(new@_type)
mnode_cast().
self.promote(model)
mnode_promote().
self.parse(string)
saml_parse(), which is not a function
on mnodes but still a very useful one.
self.zero()
mnode_zero().
self.one()
mnode_one().
self.negate()
self, this corresponds to
mnode_negate().
self.sqrt()
mnode_sqrt()
self.det()
mnode_det()
self.ne0()
mnode_notzero().
self.lt0()
mnode_isneg().
self.differ(arg)
mnode_differ().
self.lessthan(arg)
mnode_lessthan().
self.info(what)
mnode_info()
(which is not yet implemented currently, but the stubs are here).
self.gcd(arg)
mnode_gcd().
self.pow(exponent)
mnode_pow().
self.mvl(lit1, lit2)
mnode_move_lit().
self.len()
self.item(index)
The usual operations (addition, substraction, multiplication, division and remainder) are accessible with the normal binary operators. The binary comparison operators are also usable; the only problem is that they can't raise an exception if an error occurs. Finally, integers can be coerced to Mathnodes, so it's possible to write `z+1' or `if 3<z<4' and have the conversions performed automatically.
Mnodes representing errors (that is, mnodes for which
saml_errno() would return a non-zero value) are not considered as
Mathnode objects, but are reported as exceptions. There is a one-to-one
correspondance between the saml1 exceptions and the possible
return values of saml_errno(), except out-of-range values which
are all reported as SamlUnknownError (this is consistent with the
behaviour of saml_strerror(), btw). Have a look at saml1.c
for the complete list of exceptions.
Calling SAML from Python is much simpler than the same job in C, because all objects are created and destroyed automatically. For instance, our first sample program (adding 123 and 456) would be simply
from saml1 import * print Mathnode(ST_INTEGER,'123') + Mathnode(ST_INTEGER,'456')
Note that the second Mathnode() wasn't absolutely necessary, we
could have relied on automatic coercition, and written
from saml1 import * print Mathnode(ST_INTEGER,'123') + 456
In any case, we need at least one object of the "right" type,
then we can rely on promote() and parse() and automatic
conversion of integers. For example, to compute
(1+X+X^2/2)^100, we need polynomials with rational coefficients.
A solution is
from saml1 import *
model = Mathnode(ST_RATIONAL,'0').cast(ST_APOLY)
print model.parse('(1+X+X^2/2)^100')
In this particular example we could have omitted the assignment to
model and done the whole thing in one line, but usually we'll need
the model for lots of other things, so we'll keep it preciously.
Other examples can be found in the `python' subdirectory.
Coercition should work with Python long integers as well, and there should be methods to convert a Mathnode to an integer when possible. Standard mathnodes (i.e., array types) should behave like Python arrays.
Perhaps we should ensure that different Mathnodes always point to different mnodes; thus, Mathnode identity would be equivalent to mnode identity. It would allow us to use some Mathnodes (for instance, literals) as dictionary keys.
Integers are implemented as custom mnodes (but put zero in the
length field so that they can use the standard destructor), with
the following data structure:
typedef struct {
int blocks;
__u32 d[0];
} big_int;
The blocks field contains the number of limbs, and the sign of
the number. Zero has zero blocks. Numbers are written in base 10^9 (the
highest power of 10 which fits in a 32-bit word) in little-endian
format. For space-efficiency reasons, two mnodes are preallocated for
"0" and "1".
Most algorithms are taken from Knuth, Seminumerical Algorithms. The square root algorithm is the simplest I could think of, and is obviously very inefficient for large numbers. I only needed it to implement Brillhart-Morrison's algorithm, so performance wasn't crucial.
This is also a custom type with the standard destructor. Reals are stored as follows:
typedef struct {
int blocks;
int exponent;
__u16 d[0];
} big_float;
Here blocks is always positive, and a power of two; it is the
size of the d[?] array. The number is stored under the form
m.2^e, where m is between 1/2 (included) and 1 (excluded)
and e is some integer. The exponent field has the same sign
as the number, and its absolute value is 2^30+e. The mantissa is
expressed in base 2^16, and stored in big-endian format, i.e.,
d[0] is the most significant limb, in particular its most
significant bit is 1, unless the number is zero.
Here again, most algorithms are taken from Seminumerical Algorithms. Someday we should rewrite this type to use 32-bit limbs, implement the square root, and all usual transcendental functions. Later...
All cyclic integers (ST_CYCLIC) are stored in a hash table, which
is resized when appropriate. When a new cyclic integer must be returned,
we examine if it already exists in the table; if it does, we just
copy_mnode() it, otherwise we create the new object, and insert
it into the hash table.
It is exactly similar to what we do for literals, but the goal here is mainly to save memory, and to avoid extra malloc/free operations which are much more expensive than a hash table lookup.
Important: The cyclic_invert() routine assumes that the
modulus n is prime. Otherwise it will just raise your number to
the power n-2. Maybe I should fix this someday.
Polynomials (ST_POLY) are a standard type. A polynomial with
n terms (with the convention that n=0 for the null
polynomial) is represented as an array of n+1 monomials, where the
first one is zero. Why, you might ask? This is because we need to keep
type information about the coefficients of the polynomial, even if the
polynomial is zero. On the other hand, it complicates the algorithms if
we allow some terms to be zero; so this term is not really part of the
polynomial, it is only used to store type information.
Monomials are a custom type with a custom destructor. Their representation consists of the coefficient, followed by (literal, exponent) pairs. This representation is particularly suited to very sparse polynomials, but requires many allocations/deallocations and reference-count modifications.
Dense unary polynomials with an unnamed variable (ST_UPOLY) are
also a standard type. A polynomial of degree d is represented as
an array of d+1 coefficients; the null polynomial is represented
like other constant polynomials. They are used by other parts of the
library to implement polynomial substitution, or polynomial gcd (using
Collins' method).
The first implementation, while perfectly suited to very sparse polynomials, is not particularly time- or space-efficient when the number of literals is small. All monomials are allowed to use their own set of literals, which is usually redundant, and it complicates all operations.
In this alternative implementation (ST_APOLY), all literals are
grouped at the beginning of the object -- so, they don't have to be
repeated for each monomial; monomials are simply a (pointer to a)
coefficient and an array of exponents. For common problems, in appears
to be much faster and less memory-hungry than the previous approach. On
the other hand, it will be very inefficient in some cases, for instance
imagine a polynomial with is the sum of thousand different literals.
After some work, it seems to work nicely now. It is not a standard type. The set of coefficients is sparsely encoded, but the exponents are densely encoded.
[Description of how data are represented]
Tensors are the solution to a serious limitation of SAML, the absence of an external multiplication. Since one cannot multiply objects of different types, the only solution to implement vectors, matrices, numbers, and still be able to combine them together, is to consider them as special cases of one type: the tensor.
Tensors are represented by a header containing its order (the number of literals involved), then (literal, start, len) triples, and finally the data. Special functions are provided to iterate over the indices. The representation is always dense.
Rewrite all the multiprecision stuff to use base 2^32 exclusively, and unify the routines for integer and real arithmetic.
The "new" implementation of polynomials (type ST_APOLY) uses
too much memory and CPU time when there are too many literals. We should
devise something else, allowing for moderately sparse exponent vectors,
but still time- and space-efficient when the number of literals is small.
It would allow us to throw away ST_POLY entirely.
samuel(1) is a user-friendly front-end to SAML, and an unvaluable
tool to debug the library. It can be used either in interactive mode or
in shell scripts. It doesn't allow rational functions, yet. There are
other limitations; look at its manpage. User-friendliness has been
improved with readline(3) support, and the possibility to filter
all the output through an external program. As an illustration of
this, colorize(1) formats and colorizes polynomials; to use it,
just pass the option `--post=colorize' to samuel and
watch the fun.
The syntax was designed to avoid keywords because they might
clash with literal names. The option -b (or --batch) is
useful in non-interactive mode, to avoid ambiguities between variable
and literal names.
The current version of samuel does no longer use the stack provided
by yacc to store temporary expressions, but handles the stack
itself. This has at least two advantages. The most visible one is that it
avoids memory leaks when a parse error occurs. Also, it simplifies the
handling of variable-length argument lists (like determinants).
Security issue: samuel can execute arbitrary code. Thus, you
should not feed it with unverified data. Reading and writing to pipes is
especially dangerous.
induce(1) is heavily inspired from make(1), but it handles
algebraic expressions instead of files. It is a purely declarative
language; there are no variables (once a node has been computed, its
value can no longer be modified) and no loops. This has the advantage
that intermediary results can be memoized, even across multiple
invocations of induce. The implementation consists of several
more-or-less independant parts.
induce computes the smallest directed acyclic graph (DAG)
containing all the "dependencies" of the roots (the nodes whose values
have been requested by the user). If it detects cyclic dependencies, a
fatal error is raised. If this graph exists but is infinite or simply
too big, induce will die of memory exhaustion. Each vertex of the
DAG requires between 60 and 70 bytes approximately; it contains
references to its parents (immediate dependencies), a pointer to a
bytecode string, and a reference counter, which counts how many children
of this vertex haven't been evaluated yet.
induce evaluates all vertices of the DAG, starting from the
leaves (vertices which don't depend on other nodes), and going to the
roots. There is usually more than a way to evaluate vertices in an order
compatible with the partial ordering of the DAG. We try to order
calculations in order to minimize the number of active mrefs at a time,
using some simple heuristics. Since this problem is equivalent to
optimal register allocation, which is known to be NP-complete in the
general case, we cannot expect that these heuristics will work in all
situations. And indeed, they are particularly inefficient when the DAG
is a tree (i.e., when there are no common subexpressions).
Anyway. This order guarantees that when a vertex is evaluated, all its
parents, if any, have been evaluated first. So we execute the bytecode
pointed to by the vertex, store the resulting mref in the vertex, and
decrement the reference counters of the parents; if the counters reaches
zero, we can free these vertices -- and, more importantly, the mrefs
holding their values.
When the roots are read from standard input, they are processed one line
at a time (so induce can be used as a bidirectional filter, in
line mode). The above procedure is repeated for each line. A consequence
of this is that induce doesn't remember the dependencies and the
results of the previous lines; this is intentional, and ensures that the
process size remains bounded. If you want memoizing between lines, mark
as `precious' the appropriate nodes.
This program is rather straightforward. Factoring integers is a very
specialized task, thus it deserves a separate program (at least if one
believes in the Unix philosophy of "small sharp tools"). I've taken the
same notations as Knuth in Seminumerical Algorithms. The coolest
feature is the option -m (or --mail) which puts
factorint in the background, and sends a mail to you when the
factorization is complete.
See the manpage for details; factorint uses the Rabin-Miller test
(aka "strong Fermat test") for pseudo-primality testing. Factorization uses
Pollard's rho method, and then the Multiple-Polynomial Quadratic sieve.
Factorization can be extremely long. On my machine (a 486dx2/66 running Linux) it takes two or three days to factorize a seventy-digit number; Moreover, the program eats up to fifteen megabytes of memory and the same amount of disk space.
Put tests, loops and control structures in the samuel interpreter?
I'm not sure it's really a good idea; we'd have to redesign the interpreter
entirely, whereas such functionality can be provided by external tools
(either using pipes, or by interfacing with scripting languages). Let's
keep it simple.
On the other hand induce lacks many commodities: it would be nice
to have (internal) functions and loops and sums over a range of indices,
but I doubt it'll ever be implemented, because these things (and much much
more) are trivial to do in Python. Maybe we could rewrite induce
in python.
The program factorint still has a few quirks with numbers whose
factorization contains big primes raised to exponents greater than one.
Sure, it recognizes perfect squares but MPQS will choke on numbers like
p^2.q^3 where p and q are two big primes.
Moreover, it might be useful to implement Lenstra's algorithm which
uses elliptic curves.
This document was generated on 16 November 1998 using the texi2html translator version 1.54.