<> <> AlgebraStructures Notes Dennis Arnon © Copyright 1986 Xerox Corporation. All rights reserved. 1. General Ramblings 1.1. Goals Be able to create structured algebraic structures on demand, starting from any available ground structures and using built in constructors, and perform IO and operations on elements of these structures. Motivation was to get generic procs for the domains of interest. The general idea is recursive: assume that you have a ring whose elements are some kind of REF type. You can then make new rings with three structuring operations: polynomials, matrices, and real fields. 1.2. Semantics We have a single fixed class record for all instances of structures in that class. That is, all instances of structures in the class will contain a pointer to exactly the same record in memory. A class record contains only procedures. Property lists provide a primitive subclassing mechanism. That is, for each given type of class record (group, ring field, etc.), we can instantiate variations of it with different property lists. The procedure-only principle extends to property lists - they contain only procedures. We seek to enforce the paradigm that data (e.g. matrix size, variable list, minimal polynomial) pertains to particular instances of structures, and that it goes in the instance record of each particular structure. Certain familiar structures, e.g. integers, reals, bigrats, are automatically created by the package. For some differences between algebraic structures it is clear that we need different built in kinds (classes) of structures (e.g. group, ring). Having established our builtin kinds of structures, it is clear that there can be small differences among structures of a given kind; these distinctions we make via the property list (e.g. ordered rings). Properties possessed by all structures of a given kind (e.g. identity, characteristic for rings) are in the class record for that kind, provided that (1) their values can be obtained by invoking procedures, (2) it is reasonable that the parameters of such procedures be REF ANY, or a particular structure kind, which is all we have available in the AlgebraClasses interface. Properties not meeting these requirements, and clear interface data, have to go in later interfaces. For some "properties" it is less clear where they should go; generally we put them on property list (e.g. Euclidean domain). The value of a property may be a record of operations that are available because a structure has a given property (e.g. comparison ops made possible by order property). Our view is that an algebraic structure consists of two things - a class record of operations (some are defined for all rings, and others, on the propList, are defined only for some subclass of rings), and instance data. The ring class record contains procedures for the class, where a "class" of rings is all those rings builts up using a particular constructor. We can take the definition of what proc definitions belong in the constructor interface as: those which are the same for all rings in the class. Something that allows us to get polymorphism out of our structuring interfaces, e.g. polynomials, matrices. Each particular flavor of polynomial or matrix is defined to be a just a generic polynomial or matrix, e.g. RatPolynomial: TYPE = Polynomials.Polynomial; One way that this pays benefits is that the ring instantiation procedures can be in the constructor interface, because when the class procs do their NARROWs, they get generic polynomials. The interface for a new structured ring can get a handle to its ground ring in one of two ways: either the ground ring has already been instantiated in some other interface, or we instantiate it in this interface. In either the case, the direct procs just call appropriate constructor procs, passing them the ground ring. The place where type checking really happens is when we get down to operations on the ground ring, i.e. the class procs for the ground ring will do NARROWs which will fail if the wrong ground ring is present at the bottom. The interface for a structured ring contains two things. (1) the instantiation of the ring. (2) direct procedures for operating on elements of the ring. The direct procs have inputs and outputs of the particular structured ring type. They utilize the structuring operations to do their computation. They don't have to NARROW their inputs, but they do have to NARROW the outputs they get back from structuring operations. To instantiate the ring, you define class procs which narrow their arguments and call the direct procs. Obviously there is no need to "widen" the outputs that come back from the direct procs. You might say that the NARROW done by the direct proc is wasted effort, but this is the price one pays for the convenience of having the direct procs, which anyone who is a direct client of this ring will use. The standard naming convention is that the direct procs have names which mimic the fields of a ring class, e.g. "Multiply", and the class procs simply put the word "Class" before this. If you get an op for a ring via its class record (as is done e.g. by structuring operations), you have to NARROW the result, since by definition, AlgebraClasses.UnaryOp, BinaryOp, etc. return a REF. A finished, structured ring knows what its ground ring is, and even though that information is not in its elementType declaration (e.g. we define RatPolynomial: TYPE = Polynomials.Polynomial), it will be used in ring's direct procs, and indeed, the ground ring should be clearly stated somewhere at the beginning of the interface. Thus, in order to build up a structured ring you have in mind, all the "lower order" rings have to exist, i.e. have been instantiated, and you have to know what at least the outermost of the interfaces containing those instantiations is to get your hands on the outermost ground ring you need. 1.3. Fancy operations We expect that, in general, one will want additional operations on a ring than those provided by its class record or property list. These will be implemented by a separate sequence of interfaces. Thus, if one wanted an integer polynomial factorization algorithm that did complete factorization of the integer content, it would be necessary to have an integer factorization routine. To us it is appropriate that you do not, and do not expect to, find such a procedure in the class record or property list for the integer ring. There seems little uniformity in how much genericity will be possible in such modules. Thus it seems clear that one can have generic polynomial division and real root isolation modules. Gcd seems clearly to admit genericity for prs algorithms (since pseudodivision does), perhaps for Chinese remainder lifting of variables, for multivariate polynomials over a field (use prs algorithm at univariate base), and not at all for Chinese remainder lifting of integer coefficients. There appear to be two kinds of polynomials we want to factor, rational and algebraic, possibility of genericity unclear. This point of view has an algebraic interpretation - we are saying that the only kinds of "finer grain" ring operations we are willing to build into the class structure are ordered rings and fields. Operations pertaining to euclidean domains, or unique factorization domains, have to be implemented outside the ring structure. Thus the ring structure makes it easy to get I/O and basic ring operations for algebraic polynomials, giving you a platform from which to separately do the additional work to e.g. get factorization. 1.4. Modus operandi We try to follow the rule that there is one class record for more than one domain. However the sequence constructor violates this; every time you make a new sequence structure you make a new class record. The problem here is not that we have a class that has only one domain in it, but the "lack of uniqueness" aspect: we unfailingly create the new class record every time, regardless of whether we might already have made this kind of sequences. This can be dealt with by having unique hashcodes for structures, and a registry of all created structures, and then every structure constructor checks the registry to see if a requested structure is already there. This could be handy if someone tries to recreate the basic structures, or e.g. to prevent algebraic number routines from multiply recreating the same algebraic number field. Hashcodes should probably be ropes. There is a tension between having e.g. a PolynomialOpsRec that goes on the property list, and simply having procs in the interface; currently we do both (i.e. we're redundant). The general problem is that we have a single kind of StructureClassRec defined in AlgebraClasses, that has to be made to serve all algebraic structures. A better mechanism would be to simply put REF PROC's on the property list, plus build up a list of the available ops. Only the things which absolutely every structure has, e.g. a printNameProc, would be hard fields of the class record. Algebraclasses would export a "dispatcher", which would take in the name of a desired proc and know whether to get it out of a record field or off the proplist. The general issue is - how to make all knowledge about what procs are available in a structure exported to clients by the run-time structure (ideally even including documentation, perhaps a la pop-up buttons). Then the CaminoReal user interface could be just like any other client, and would contain no hardwired knowledge of the underlying algebra system, except for knowledge of the generic mechanism by which the specifications of a structure are encoded. We then have a field in class records in which is stored a list of triples: property name, whether in class record or a proc on property list or a property value on the property list, and a documentation rope. Class record then offers procs to tell you whether a property exists at all, or to get its value on the assumption that that value is a proc, or to get its value on the assumption that that value is a bool, or to get its value on the assumption that that value is a rope, etc. Also there is a proc to do a full printout of the triple list - names, nature of values, and documentation ropes. Some other proc returns a list of (standardly modified, e.g. remove "$") names of proc-valued properties; this can then become part of the pop-up menu of ops for that structure. AlgebraClasses offers a proc that takes five args (a StructureClass variable, a PropName, a Ref that is the prop value, a specification of the nature of the value, and a documentation rope), and does two things: puts the prop on the StructureClass's proplist with the specified value, and adds the (PropName, specification of the nature of the value, documentation rope) triple to the list of triples stored in that StructureClass). The Impl module for a class then has as part of its start code a call to this proc for each op associated with the class. A basic convention is that the procs implementing structure ops should all have types defined in AlgebraClasses. Thus we insist that our structures are "run-time self-documenting". The run-time user can find out in a standard way what's available in a structure. The client is expected to look at the interfaces, or perhaps there is a manual that tries to stay up to date (one way it might stay up to date is to regenerate itself on the fly, i.e. whenever you want to display it, it goes off and actually calls the same sort of procs we have been talking about, in all the interfaces, to produce a list of what ops are available). MethodApplier: PROC [methodSelector: ATOM, argList: LIST OF Object, structure: Object _ NIL] RETURNS[value: Object]; <> ok:BOOL; movingPointer, outArgs: LIST OF Object; arg: Object; CheckArgs: PROC [methodSelector: ATOM, argList: LIST OF Object, structure: Object] RETURNS [ok: BOOL, outArgs: LIST OF Object _ NIL] ~ { desiredArgStructures: LIST OF AC.Object; outArgs: LIST OF Object _ NIL; IF NOT AC.CheckMethod[methodSelector, structure] THEN RETURN[FALSE ! "No such method"]; desiredArgStructures _ AC.DesiredArgStructures[methodSelector, structure]; WHILE desiredArgStructures # NIL AND argList # NIL DO desiredArgStructure _ desiredArgStructures.first; desiredArgStructures _ desiredArgStructures.rest; arg _ argList.first; argList _ argList.rest; IF AC.StructureEqual[desiredArgStructure, arg.structure] THEN IF outArgs = NIL THEN outArgs _ LIST[arg] ELSE outArgs.rest _ LIST[arg] ELSE { recastArg _ AC.PerformMethod[$Recast, desiredArgStructure, LIST[arg] ]; -- automatic type conversion occurs here. An alternative is to have all methods like Recast, which every structure has, in AC as concrete procs: <> IF NOT AC.StructureEqual[desiredArgStructure, recastArg.structure] THEN RETURN[FALSE ! "Type clash in arg j"]; IF outArgs = NIL THEN outArgs _ LIST[recastArg] ELSE outArgs.rest _ LIST[recastArg]; }; ENDLOOP; IF desiredArgStructures # NIL OR argList # NIL THEN RETURN[FALSE ! "Bad Number of Args"] ELSE RETURN[TRUE, outArgs]; }; <> IF structure # NIL THEN { [ok, outArgs] _ CheckArgs[methodSelector, argList, structure]; IF ok THEN RETURN[AC.PerformMethod[methodSelector, structure, outArgs] ]; }; movingPointer _ argList; WHILE movingPointer# NIL DO arg _ movingPointer.first; movingPointer _ movingPointer.rest; [ok, outArgs] _ CheckArgs[methodSelector, argList, arg.structure]; IF ok THEN RETURN[AC.PerformMethod[methodSelector, arg.structure, outArgs] ]; ENDLOOP; SIGNAL NoApplicableMethodFound; END; A basic rule: every method must return exactly one object. AC.CheckMethod, AC.DesiredArgStructures, AC.PerformMethod, etc. look in the method dictionaries of their structure arguments, and if need be, the method dictionaries of their superclasses. 6/20/87 - Consider a Structure S whose class (i.e. the value of field S.class) is C (C is an Object of flavor Class). There are two possible flavors of superclasses (i.e. flavors for C.class) in this situation: Class or Structure. The first sort is obvious - consists of methods, corresponds to inclusion of categories. The second sort is only appropriate in case S is the only Structure in its class. Then if S is a substructure of a Structure T, we can set C.class to T. The obvious example is Integers are a substructure of Rationals are a substructure of Reals. The effect is that if we are searching for methods and don't find them in C, and we see that C.class is a Structure T, then we will continue our search in T.class. It is assumed that if C.class = T, then T's recast proc can handle elements of S. At the moment only single inheritance is supported. General note on type conversion: If we ask to divide a rational number by an integer, then regardless of order in which args occur, assuming that BigRats has a recast proc that handles integers, we'll get the desired BigRat result. Consider dividing one integer by another. If we don't want it available at all, then leave Ints interface just as it is. If we want it available as truncation, then make that a method in Ints interface. If we want it to be available and yield a rational result, then have a divide method in Ints interface that imbeds args as BigRats and then invokes BigRats division method. General note on inheritance: the Smalltalk paradigm seems to assume that a class has one model: a subclass is a subset of that model consisting of all elements that have some further property. Thus the (unique) model of a subclass is a subset of the (unique) model of the class. The situation we have with algebraic structures is rather multiple models. I.e. we have a category, e.g. the real closed fields. All models in this category may have an add op, but you can't give a single particular add proc that will work for all models. 2. AlgebraStructures improvements 2.1. Integer arithmetic, quotient constructor. Combiner uses CaminoReal. Separate BigRats for CedarChest like BigCardinals, not AlgebraStructures? 2.2. New Domains Vector type: Use it or point or sequence for RatIntervals, Variables. EquationsOver, RelationsOver (similar to FormulasOver). Quantifier-free formulas over? Add existential, universal quantifiers. 2.3. Extensibility issues How should users specify new operators, or new domains? Need to specify things like canonical forms, formatting rules for legal expressions. 2.4. Recast (conversion) issues Automatic conversion where appropriate, or tell user to modify domain if he want that op? How hard is automatic domain inference for an arbitrary expression (Scott Morrison)? 2.5. Substitution facility What are the correct specifications of a general substitution facility? 3. Cad implementation 3.1. New operations Request for any AlgebraicSurfaces display op starts up the tool if not already started. Or have it running all the time, and each display request creates a new surfaceDisplay viewer. AlgebraClasses has proc to add an op to a structures property list. Or MakeSequenceStructure proc takes procs to add to prop list as args. Cad ops induced cad Alternate display showTriangulationAll (valid only for cad's of 1-, 2- and 3-space) showTriangulationSections showTriangulation0Cells showTriangulation1Cells showTriangulation2Cells rayTraceInputPolynomials (valid only for cad's of 3-space) show cells show statistics Cad's should have names, statistics Store only basis signatures with cells. Both clustering during cad construction, and clustering for disambiguation of signatures, are done with basis signatures (if basis signatures with ordinary projections don't suffice to get unambiguous clusters, then we have to go back and refine cad with respect to augmented projections). Cluster ops select containing cad rayTrace (valid only for 2-cluster in 3-space) showTriangulation (valid only for clusters in cad's of 1-, 2- and 3-space) Cell ops select containing cad rayTrace (valid only for 2-section in 3-space) showTriangulation (valid only for cells in cad's of 1-, 2- and 3-space) Polynomial sequence ops, for base coeffs that are ints or bigrats rayTraceZeros (valid only 3-variate polys) decomposeZeros (forks Vax computation; pop up menu for options (including covering set parameters); inserts cad in Active MessageSet) showZerosTriangulation (valid only for 1-, 2- and 3-variate polys) New Polynomial ops, for base coeffs that are ints or bigrats rayTrace (valid only for 3-variate poly) decomposeZeros (forks Vax computation; pop up menu for options (including covering set parameters); inserts cad in Active MessageSet) showZeroTriangulation (valid only for 1-, 2- and 3-variate poly) Sequences Basic display is one line: "sequence of n shortPNames" Alternate display: Slider bar onto popup menu popup menu entries are shortPNames of sequence elements Insert creates a blank entry in popup menu; then select item to be inserted as an operand, it is type checked, inserted. Delete modifies pop up menu Selection of Insert or Delete forces Alternate display 3.2. Adjacencies and Clusters, Graph ops, read adjacencies, compute clusters We have made a basic decision to go with basis signatures. First step is the current "cells with signature" op. This assumes that every cell in the cad has a basis signature, then this is a binary op that takes a legal (i.e. correct length) basis signature (perhaps with * elements) and returns a sequence of cells with that sig. Subr for boundary of 2-sector (i, j) in E2 I3=lelt(I,3); I2=lelt(I,2); Pass I#, I3, and I2 to a routine For each 2-sector (i, j, k) in I# Subr for boundary of 2-sector (i, j) in E2: L1 _ U1 _ 0; L1L _ L1R _ U1L _ U1R _ 0; for each 1-cell (a,b) adjacent to 2-sector (i, j) if a = i then if b = j-1 then L1 _ j-1 else if b = j+1 then U1 _ j+1 else ERROR; if L1#0 then for each 0-cell (a,b) adjacent to 1-section (i, L1) if a = i-1 then L1L _ b else if a = i+1 then L1R _ b else ERROR; if U1#0 then for each 0-cell (a,b) adjacent to 1-section (i, U1) if a = i-1 then U1L _ b else if a = i+1 then U1R _ b else ERROR; Return[ LIST(L1, LIST(L1L, L1R), U1, LIST(U1L, U1R) ) ]. ExtractClusters: PROC [cad: Cad, numComponents: CARDINAL] RETURNS [clusters: SEQ.Sequence]; <> <> << User specifies polynomials with respect to which to cluster by (*,*,@,*,*,@) tuples (@ is new FormulaOperators.Op).>> <> <> <> <> <<};>> <> <> <> <> <<};>> <<>> <> <> <<>> <> <> <<>> AdjacenciesOverAdjacency: PROC [dimension: CARDINAL, adjacency: Adjacency, cad: Cad] RETURNS [newAdjacencies: LIST OF Adjacency]; <> <> <> 3.3. Cad statistics records SmartInfo about the cad (0-cell bounding box, singularities, turning points). 3.4. New cell Ops: read covering set, display covering set 3.5. AlgebraStructures Rootfinder, new cell op: generate covering set 3.6. Compute and add adjacencies to an existing cad in Cedar Do we want to separate adjacency computation from cad computation? If so, set up a SAC-2 routine to do it; allow ComputeServer to queue up a succession of jobs, to be done in background. Step two: implement the computation in Cedar. 2D case first. 3.7. Cad databases Reading such files is slow. Need command to read some more cells of this cad and add them to the database. Cypress database? Once have distributed computations, we want to be able to do incremental updates of the database. <<3.8. RefineCad operation (Real algebraic variety Combiner)>> This is a great example for implementing cad computation related stuff in Cedar by forking off appropriate subcommands to remote SAC-2 systems. RefineCad: PROC [inCad: Cad, refiningPoly: Polynomial, cellsToRefine: Sequence] RETURNS [outCad: Cad]; <> <> <<3.9. Projection Phase>> CoarsestSqFreeBasis: PROC [inputPolynomials: SEQ.Sequence, localizationFormula: QFF.Formula _ NIL] RETURNS [contents, basis: SEQ.Sequence]; Projection: PROC [basis: SEQ.Sequence] RETURNS [basisProjection: SEQ.Sequence]; <> <<(1) we have a certain set K of polynomials in E^r which we need to make delineable on each cell and each signature-cluster in E^(r-1). >> <<(2) The collection of cells we end up with in E^r must have the boundary property: the boundary of any cell must equal the union of some collection of lower-dimensional cells.>> <> <> <> <<(1) Compute contents and basis(pp) = K of inputs, pull off initial sequences of non-constant coeffs of all basis elts, compute basis for contents union initial sequences, do a signature-invariant decomp D of E^(r-1) to get regions on which each element of K has constant degree. Note that if any element of K is identically zero at some point of E^r-1, D will have maximal regions for this property.>> <<(2) using the appropriate reducta of elements of K (ignore identically zero elements), compute projection sets for each reducta pattern that occurs on a positive-dimensional regions, compute basis. For each positive-dimensional region, use this basis to break it up to be signature-invariant with respect to its projection basis. We end up with a new decomposition D'.>> <<(3) Assume now r <= 3; maybe works for general r.>> <<(4) On each region in E^(r-1), evaluate the appropriate non-identically zero reducta of elements of K to get a set of algebraic polys; compute its basis. If any element of K is identically zero on this region, determine the appropriate algebraic polynomials whose roots we need to take account of to insure boundary property, compute their basis, merge the two bases. Isolate roots of the final basis to determine the stack over this region.>> <> <<(1) Compute P1 = contents and K = basis(pp) of inputs.>> <> <> <<(Possible precomputation for each Pijk: do a heuristic test (e.g. Groebner basis to check emptiness of complex variety), or cad, to check that real varieties stay nonempty as you take progressively larger sequences of initial non-constant coeffs; if such a variety is empty, then you don't need to continue. It could happen that appropriate tails of initial sequences of elements of K could be thrown out of P2).>> <> <> <<(2) Assume now r <= 3; maybe works for general r.>> <<(3) Extension phase: On each region in E^(r-1), evaluate the appropriate non-identically zero reducta of elements of K to get a set of algebraic polys; compute its basis. If any element of K is identically zero on this region, determine the appropriate algebraic polynomials whose roots we need to take account of to insure boundary property, compute their basis, merge the two bases. Isolate roots of the final basis to determine the stack over this region. Record K-signature of each element of the stack. >> <<3.10. Base Phase>> BasePhase: PROC [cad: Cad]; <=1 space; follow chain of induced cad's down to 1-space, contruct cells there. Fill in basisSignature's = primarySignature's of cells. If 1-space cad.localizationFormula # NIL, then discard cells which don't satisfy it.>> <<3.11. Extension Phase>> <> <> <> <> RegionExtendCadToCad: PROC [cad: Cad]; <> <> <> <<1. MaximalSignedRegions[i-1, cad, TRUE]; regions _ ExtractSignedRegions[i-1, cad, TRUE]; sort regions in order of decreasing dimension (for i-1 = 2, put degenerate 0-dimensional regions after nondegenerate).>> <<2. cad.cells _ NIL; Go down sorted list of regions, for each:>> <<2a. currentCells _ ExtendRegionToCells[i, currentRegion, cad]; newAdjacencies _ NIL;>> <> <<2b. Do this step only if i = 2 or 3. If dimension of currentRegion < i-1, then determine adjacent regions of larger dimension. Sort these into order of increasing dimension. Go down this list, and for each region on it, find all cell-cell adjacencies between currentRegion and this region. Sort these adjacencies into decreasing sum of cell dimensions. For the next such adjacency, newAdjacencies _ AdjacenciesOverAdjacency[i, adjacency, cad]. If newAdjacencies = NIL, i.e. if needed auxiliary E^(i-1) cell-cell adjacencies, or non-NIL sample points, cannot be found, then go to the next cell-cell adjacency between current region and this region. If needed auxiliary adjacencies and/or non-NIL sample points cannot be found for any cell-cell adjacency between current region and this region, give up on building adjacencies in E^i between stacks over current region and this region. >> <> <<2c. Insert cells of currentCells into cad.cells in lexicographical order. IF newAdjacencies # NIL, then enter its adjacencies into the appropriate cells of cad.cells. >> <<3. If i-space cad.localizationFormula # NIL, then delete all cells of cad.cells which don't satisfy it, and all adjacencies involving those cells .>> <<4. For i = 1, 2 or 3, if no localization, we now have a graph for the cad of E^i in which a maximal signed regions computation will yield (signed) components (in 1 and 2-space, we can have all adjacencies, in 3-space, probably only a proper subset of all adjacencies.) >> ExtendRegionToCells: PROC [dimension: CARDINAL, region: SignedRegion, cad: Cad] RETURNS [newCells: SEQ.Sequence]; <> <> <<(Let i = dimension). Method: Evaluate polys at base sample point, irreducible basis of algebraic polys, isolate roots, construct extended sample points, eval gsfd of i-variate polynomial at endpoints of isolating interval, or i-variate polynomial itself at easy primitive sector sample point, when we want its sign.>> <> <> <> <> <> <<];>> ExtendCellToStack: PROC [cell: Cells.Cell, cad: Cad] RETURNS [StacksInEi: SEQ.Sequence]; <> ExtendStacksToStacks: PROC [i: CARDINAL, StacksInEiM1: SEQ.Sequence, cad: Cad] RETURNS [StacksInEi: SEQ.Sequence]; <> <> <> 3.12. QE phase <> <> <> <> <> <> <> <<];>> MakeQEProblem: PROC [inputFormula: QFF.Formula, localizationFormula: QFF.QFFormula, inputVariables: VARS.VariableSeq, minPolyVariable: VARS.VariableSeq _ AN.defaultMinPolyVariable, fieldElementVariable: VARS.VariableSeq _ RF.defaultFieldElementVariable] RETURNS [QET.QEProblem]; EliminateQuantifiers: PROC [QET.QEProblem]; <<1. Extract polynomials of formula and localizationFormula; cad _ ProjectionPhase[polynomials, inputVariables, localizationFormula]. >> <<2. Cad.BasePhase[cad];>> <<3. For i in [2,...,k] do Cad.RegionExtendCadToCad[i, cad].>> <<4. Cad.MaximalSignedRegions[k, cad, TRUE]; regions _ Cad.ExtractSignedRegions[k, cad, TRUE]; Sort regions into order of decreasing dimension (for k = 3, put degenerate 0-dimensional regions after nondegenerate), since qe for larger dimensional regions should go faster, since they are likely to have lower degree sample points.>> <<5. For each region of regions, using representative cell = cellInEk (which has extended or primitive sample point), do:>> <<5.1 StacksInEk+1 _ Cad.ExtendCellToStacks[k+1, cellInEk, cad];>> <<5.2 For i in [k+2,...,r] do StacksInEi _ Cad.ExtendStacksToStacks[i, StacksInEiM1, cad].>> <<5.3. EliminateInnermostQuantifier[StacksInEr, formula]; >> <<5.4. For i decreasing in [k+1,...,r-1] do EliminateInnerQuantifier[i, StacksInEi, StacksInEiM1, formula]; >> <<5.5 If StacksInEk+1.baseTruthValue then insert this signed component into SignatureTable. >> <<6. DisambiguateSignatureTable, extract signatures corresponding to regions in solution set, bundle them up as a SignatureSeq>> <<7. Apply SimplifySignatureSeq; Return SignatureSeqToFormula and the final SignatureTable.>> <<>> EliminateInnermostQuantifier: PROC [StacksInEr: QET.StackSeq, formula: QFF.Formula]; <> <<>> EliminateInnerQuantifier: PROC [i: CARDINAL, StacksInEi, StacksInEiM1: QET.StackSeq, formula: QFF.Formula]; <> <<>> 3.13. Simplifier The basic strategy is to assume cells in r-space invariant with respect to bases in all dimensions. No effort to separate complete factors of input polynomials from factors of projection polynomials. If we have unambiguous signatures at this point, i.e. if for each signature, either every cluster with that signature is in solution set, or no cluster with that signature is in solution set, then we just have to simplify. If not, we have to refine the cad with respect to augmented projections and get clusters invariant with respect to these new larger bases in each dimension. Having done this we know we'll get unambiguous clusters. <> <> <> <> <> <> <<];>> <> <> <> <> <> <<];>> DisambiguateSignatureTable: PROC [inTable: QET.SignatureTable] RETURNS [outTable: QET.SignatureTable]; <> <> <> <<1. Select some polynomial P from the basis of some (proper) induced cad, such that P is not an element of outTable.polynomials (use hashing).>> <<2. Augment the secondarySignature's of all cells of all regions of this entry by sign(P).>> <<3. For each region of this entry, form the subgraph it induces, and construct signed components in this subgraph wrt secondarySignature.>> <<4. Collect the cumulative data: if the sign(s) of P on all newly formed components which are contained in the solution set differ from the sign(s) of P on all newly formed components which are external to the solution set, then we have succeeded in disambiguating this entry with P.>> <<5. Delete this entry, and create two new appropriate unambiguous entries, with appropriately augmented signatures.>> <<6. For all unambiguous entries, augment their signatures with "don't care" in the P position, and augment appropriately (i.e. with the correct P sign) the secondarySignature's of all constituent cells of all entry regions.>> <<7. For each ambiguous entry, carry out steps 2,3,4, and see if we have succeeding in disambiguating that entry with P. If so, then carry out step 5 for this entry. If not, then split this entry into up to three or more (each of which may be ambiguous or unambiguous, but at least one of which must be ambiguous if we are here) entries, corresponding to splitting each of the prior entry regions into up to three or more new regions by the sign of P.>> SimplifySignatureSeq: PROC [QET.SignatureSeq] RETURNS [QET.SignatureSeq]; <> <> SignatureSeqToFormula: PROC [QET.SignatureSeq] RETURNS [QFF.QFFormula]; <> <<3.14. CadCoodTransform>> Actually is an op on a set of polys so that we end up with a decomposition which is cylindrical with respect to two coordinate systems. 3.15. Application of Cedar cad algorithm 1. Demonstrate topologically reliable display of algebraic and parametric curves and surfaces, and display of fat plane and space curves. 2. Application of distributed cad algorithm to qe problems involving roots of quartic polynomials (full Kahan ellipse, root configuration-coefficient classification for quartics) 3. Study of parametric curves and surfaces: spline analyzer, spline combiner. 3.16. Expert system for directing problem solving by cad computation Has the knowledge base of info about root configuration from evaluation of coefficients. Tries heuristics, e.g. detecting common zeros by Grobner basis computation. Knows what computations to do to build a sufficient database for solving the qe or decision problem. 4. Semi-Algebraic Set Display 4.1. General considerations for display We are displaying semi-algebraic sets in r-space, for r = 1, 2, 3, for which we have algebraic cell decompositions (acd's). For any such semi-algebraic set, there are particular sets of r-variate, r-1-variate, ..., 1-variate polynomials involved in its description. We assume that a cylindrical algebraic decomposition (cad) that is sign-invariant for those sets of polynomials has been computed, and that the acd we have for the semi-algebraic set to be displayed is the union of certain of the cells of the cad. Display is useful for using cad algorithm in solving all types of problems, not just the cases where we are trying to display a semi-algebraic set. We will show with examples why this is so. 4.2. 2-space display Covering sets and triangulations for implicit plane curves This was what was done in [ARN83]. The idea is simple: you assume you have a acd of the curve into 0- and 1-cells, with adjacencies known. Each 1-cell is a function graph, so you just step along the appropriate x-interval, find roots of the appropriate F(x,y) at each x, then use the rest of the defining formula to put yourself on the correct "branch" of the curve. You use adjacency information to make the correct transitions from 0- to 1-cells. May be desirable to have the cad algorithm produce 1-section defining formula's which have the form "$a < x < b ~AND~ i^{th} ~real~root~of~F(x,y)$". Covering sets for semi-algebraic sets in the plane The task here is - given a particular semi-algebraic set $X$ and an $\epsilon$, generate a finite set of points so that for every point $x$ in $X$, there is some point $p$ in the covering set so that the distance from $x$ to $p$ is less than $\epsilon$.\footnote{D.T. Lee (Siam Conf. Albany, July 1985) says that this is called the "sources" problem in OR. Also mentioned Shamos thesis p. 159. Seems to be some connection with the minimum covering circle problem Easy to generalize the approach used to get covering sets for curves in the plane, by adding $y$-stepping for 2-cells. Triangulation of semi-algebraic sets in the plane Assuming a covering set has been constructed, we could use the general algorithms of Lloyd (1977) thesis, Shamos thesis, Tarjan et al, to triangulate. Probably would be smarter to triangulate as the covering set is generated, given that one knows the arrangement of points. 4.3. 2D display improvements Use adjacencies for continuous (no big 0-cells) display of curves. Display of open, closed 2-cells. Eliminate use of df's in cell display? 4.4. 3-space display Triangulation of surfaces in 3-space Algorithm: Run the cad algorithm with a bounding box. Then you get a cad of the plane in which some bounded collection of 0, 1, and 2 cells constitute the projection of the surface. Construct covering set and triangulation for this semi-algebraic set. Then use the same idea used for plane curves to triangulate the 0- and 1-cells in space, by extending over the 0- and 1-cells in the plane. For each 2-cell in space, triangulate it by simply extending over its projection in the plane, and picking the appropriate root of the appropriate $F(x,y,z)$, Then using known adjacencies of 2- and 0,1-cells in space, connect up all the triangulated objects in 3-space to get the final triangulation of the (portion of the) surface (in the bounding box). Ray tracing Roth [ROT82], Davis thesis [DAV82], Kajiya [KAJ82], Hanrahan beam tracing [HAN84]. 4.5. 3D display improvements Merge in Rauen's tool (changes to ThreeDWorld format). Need to normalize forward and up vectors? Why does Rauen have to read files? Why does he need his own internal form? Use Jules' Geometry3DApplied interface. Offer wireframes. Maureen's color solid viewing tools. Giordano's assertion of the three standard, plus a certain perspective view, as giving intuition. Rauen has alg to determine what cells a ray hits? Would be a good op to offer in interface, i.e. what cells would I see if I had xray vision? Vuillemin - draw in axes on request. Giordano's 3D interface proposal (CSL-Notebook). How do you come up with a good triangulation? How to have a good algorithm which both may, and may not, take account of adjacency information? Steps by which Rauen takes a triangulation and renders it. Eliminate use of df's in cell display? 4.6. User interface considerations Mode: Topology-priority (faster) vs. Shape-priority (slower) User selects whether he is more concerned with rapid real-time interaction in which the only constraint is to show topology accurately, or with beautiful pictures that show the true shape of objects. Mode: Smooth vs. Structural User selects whether to have the individual cells "glued together", in order to see smooth objects, or whether to see the individual cells. In the latter case, 0-cells are expanded into fat balls and 1-cells are fattened into "wires". An interesting special case of this mode is displays of entire cad's. A cad is just a particular cell decomposition of the line, plane, or all space. If the user selects geometry mode, we just see some indication that the displayed set is the entire line, plane, or space. In structure mode, we see the individual cells of the cad. A simple initial interface for surface display The user gives a trivariate polynomial, which we interpret as a surface to be displayed, and a bounding box. A cad (for the surface and bounding box), and a triangulation of each positive dimensional cell, is done on Vaxc and written to a file, which includes a specification of which cells comprise the portion of the surface in the bounding box. The display tool reads the file. Initially there are two kinds of display available: a topology-priority, structural mode (done by displaying the triangulations of the cells which comprise the surface), and a shape-priority, smooth mode (done by ray tracing the surface defining equation directly). Someone (it is unspecified who) may crawl over triangulations in an effort to "topologically optimize" them, i.e. to get by with as few polygons as possible while still guaranteeing the topological correctness of any collection of cells. At this level there are no handles on individual cells. You can only display the entire portion of the surface in the bounding box. Next stage: ray tracing 2-cells contained in surfaces Now cad files are required to contain some kind of defining formula for each 2-cell, the format of which has been negotiated. The requirement is that the ray tracer has to be able to use it in order to ray trace individual 2-cells. Ideally, the ray tracer should be able to work with the "glued-together" quantifier-free formulas of a cluster of cells. Thus we should be able to have a shape-priority, structural mode for a collection of cells lying in a surface, in which we see all the different cells and 2-cells are individually ray traced, or a shape-priority, smooth mode for a subset of surface, in which we display the union of all the cells in a smooth fashion. Now the display tool should now show the indices of all the cells currently being displayed, and allow the user to add and delete cells named by index. Future user action: define a semi-algebraic set to display As mentioned, in every display situation we have particular sets of polynomials, and a cad for those sets. The display tool assigns names to these polynomials, say A, B, ... The user is then able to write any valid quantifier-free formula using these names, e.g. "$A = 0~AND~B < 0$", or "$A = 0~\Rightarrow~B < 0$" The display tool determines the cells satisfying this formula (i.e. the cells comprising the semi-algebraic set defined by that formula), and displays them in the current display mode(s). Future user action: Pointing at cells It would be nice if the user could point, or otherwise indicate, particular cells to be deleted or added to the display. 4.7. Prior work on display of sas Contour plotting implicit curve display This is perhaps one of the oldest ideas. Rice pointed to this method [SNY83]. Sederberg (letter 10/1/84) refers to [PET84]. Note the idea arising from discussions with Dobkin of getting a "wireframe" for a surface by contour tracing. Differential equation implicit curve display This was Rice's suggestion, also found in Pavlidis' book. Pitteway's algorithm for implicit curve display Pitteway's algorithm is for implicit equations; a version of it can be given for implicit equations of any degree. Knuth (Maureen says he did this in a course some time back) asked the question - would it be competitive to plot a parametric cubic curve by implicitizing it and then applying Pitteway's algorithm? Perhaps it's worth investigating how a Pitteway algorithm applied to arbitrary algebraic curves would compete with method of [ARN83] (cf. Pratt SIGGRAPH '85 [PRA85]). Miscellaneous Farouki [FAR85]. Geisow [GEI83]. Tiller [TIL85] use of power series to generate points on a curve. Griffiths surface display [GRI78]; note the need to distinguish the problems of displaying function graphs, parametric curves and surfaces, and implicit curves and surfaces. Soiffer surface display algorithm and the references he cites. Blinn thesis and paper on generalization of algebraic surface display. 5. AlgebraicSurfaces Organization: Design and Implementation Notes 5.1. Top DF file AlgebraicSurfaces.df General DF file. Public bringover gets the load and run files, user documentation, several sample surfaces, and CedarChest bcd files needed by the load file. 5.2. Documentation AlgebraicSurfacesDoc.tioga A user-oriented description of the algebraic surfaces tool. AlgebraicSurfacesOrganization.tioga This document. CADDoc.tioga Describes the primary representation of surfaces, the Scad. Describes the file format used to load and save Scads. 5.3. Load and Run files AlgebraicSurfacesPackage.config Configuration file for the package. InstallSurfaceTool.load Load file for the package. 5.4. CAD modules Overview This is the connection to AlgebraStructures. These modules provide operations on decomposed surfaces (Scads) in the internal representation. A detailed description appears in CADDoc.tioga. CADTypes.mesa Type definitions for Scads and several related types. CADIO.mesa Provides input and output of Scads. Also a few miscellaneous operations on Scads. CADIOImpl.mesa Implementation of the CADIO interface. 5.5. Surface Tracer Overview These modules provide the ray tracer, a program which generates high-resolution images of algebraic surfaces. It is, unfortunately, very slow. SurfaceTracer.mesa Provides two public procedures: TraceSurfaces and TraceCells. TraceSurfaces produces a ray tracing of complete algebraic surfaces; TraceCells produces a tracing of specified cells in the Scads of the surfaces. SurfaceTracerImpl.mesa Implementation of the SurfaceTracer interface. LightingModels.mesa Type definitions for lighting models. ShadingModels.mesa Type definitions for surface shading models. 5.6. Surface Viewer Overview This is the connection to ThreeDWorld. These modules provide for display and manipulation of algebraic surfaces. They use polygon approximations specified by a cylindrical algebraic decomposition of the surface(s). These two-cells, along with exaggerated one-cells and zero-cells, are displayed using the ThreeDWorld and Quickviewer contexts. This is fast enough to use interactively. SurfaceViewer.mesa Includes procedures to create the color display viewer, display views of the surfaces, add and remove surfaces, mask cells, and change viewing parameters. SurfaceViewerImpl.mesa Implementation of the SurfaceViewer interface. ThreeDHacks.mesa Provides two new shape classes for ThreeDWorld. They are $FatPoint and $FatSeg, which are used to assemble 0-cells and 1-cells, respectively. A procedure, RegisterNewClasses, registers them into a 3d context. Procedures are provided to construct ThreeDWorld shapes for 0-, 1-, and 2-cells. ThreeDHacksImpl.mesa Implementation of the ThreeDHacks interface. 5.7. Surface Tool Overview These modules provide a user interface to the Algebraic Surfaces programs. SurfaceToolImpl.mesa Implementation of the user interface. Running SurfaceToolImpl adds a command, "SurfaceTool", to the Command tool. Executing this command creates a new surface tool. SurfaceIcons.mesa Interface for the icons module. Provides icons for the surface tool (user interface) window and the surface viewer (color display) window. SurfaceIconsImpl.mesa Implementation of the SurfaceIcons module. 5.8. MultiPolynomial.df Overview MultiPolynomial is a data abstraction for a polynomial in an arbitrary amount of variables. It is now included in MathPackage on CedarChest. It provides most of the internal symbolic algebra used by the Algebraic Surfaces programs. MathPackage.df DF file for the MathPackage. MultiPolynomial.mesa Interface for the MultiPolynomial module. Includes constructors, operators, and substitution and evaluation mechanisms. MultiPolynomialImpl.mesa Implementation of the MultiPolynomial module. 6. Details of AlgebraStructures representations 6.1. Distributed polynomial representation The representation is (a1, d1, a2, d2, ...), where each ai is a nonzero coefficient, each di is a degree vector, and d1 > d2 > ... in the ordering of degree vectors established by DVCompare (see below). If the variable sequence is (x1, x2, ..., xr), r >= 1, and the monomial is xi1^ei1 xi2^ie2 ... xik^eik, with each ei positive, 0 <= k <= r, and i1 < i2 < ... < ik, the monomial's degree vector is (ik, ek, ik-1, ek-1, ..., i1, e1). Note the reversed ordering. Note also that a degree vector can have any length between zero (constant term) and 2r. DVCompare: PROC [dv1, dv2: DegreeVector] RETURNS [ [-1..1] ] ~ { index1, index2, exponent1, exponent2: CARDINAL; WHILE dv1#NIL AND dv2#NIL DO index1 _ dv1.first; dv1 _ dv1.rest; exponent1 _ dv1.first; dv1 _ dv1.rest; index2 _ dv2.first; dv2 _ dv2.rest; exponent2 _ dv2.first; dv2 _ dv2.rest; IF index1 > index2 THEN RETURN[1]; IF index1 < index2 THEN RETURN[-1]; IF exponent1 > exponent2 THEN RETURN[1]; IF exponent1 < exponent2 THEN RETURN[-1]; ENDLOOP; IF dv1#NIL THEN RETURN[1]; IF dv2#NIL THEN RETURN[-1]; RETURN[0]; }; 6.2. Recursive polynomial representation exponent, reductum for nonconstant polynomials, decreasing order of exponents, last exponent may be zero. 6.3. External representation of formulas A term is either an atomic formula A = B A < B A <= B A > B A >= B A ~= B or a formula enclosed in parenthesis. A formula is either term term & formula term | formula (| has lower precedence than &) An input formula is formula $ Atomic formulas A op B are transformed internally (by the reader) into A - B op 0. 6.4. Internal representation of formulas An atomic formula is represented as (P, relation), denoting the atomic formula "P relation 0", where P is a rational polynomial of nonnegative sign, and relation is one of <, >, <=, >=, =, and #. 7. CAD Formats: Types and conversion procs for CAD's in AlgebraicSurfaces 7.1. Introduction A module providing a representation of CAD's (cylindrical algebraic decompositions), including additional display information. A file format is defined, and procedures are provided to convert between the file representation and the data representation of CAD's. 7.2. Scads Representation A Scad record contains the MultiPolynomial expression for the surface that was decomposed, and a sequence of cell records. Each cell record can represent a zero-cell, one-cell, or two-cell. This could have been done with variants, but it wasn't really worth the hassle. So, in addition to its indices, a cell record contains a degree field (= 0, 1, or 2), a sequence of vertices, and a sequence of triangles. If the degree is zero, the first vertex on the vertex sequence is used for the zero-cell's location. Anything else is ignored. If the degree is one, the vertex sequence is used and the triangle sequence is ignored. If the degree is two, both sequences are used to describe the surface for the two-cell. Defining formula Every cell in the decomposition is completely described by a quantifier-free formula. A quantifier-free formula (QFF) can be either simple or compound. A simple QFF is an equation or inequality involving a polynomial. A compound QFF is a combination, by AND or OR, of QFFs. The representation is given in CADTypes.QuantifierFreeFormula. Display information Each cell also possesses display information. For zero-cells, this information is simply the cell's coordinates. For one-cells, it is a sequence of vertices that outline a piecewise linear approximation of the cell. For two-cells, it is a sequence of vertices and a sequence of polygons that comprise a triangulation of the cell. 7.3. Surface Files - format All contents of the surface file are ignored until a BEGIN token is encountered. The surface file then described as follows: surfaceFile ::= BEGIN CAD DESCRIPTION surfaceEquation surfaceName surfaceColor surfaceCells BEGIN CELL DEFINITIONS numberOfCells instances of {cellDefinition} END CELL DEFINITIONS BEGIN CELL DISPLAY INFORMATION numberOfCells instances of {cellDisplayInfo} END CELL DISPLAY INFORMATION END CAD DESCRIPTION surfaceEquation ::= Surface polynomial (equation of surface) surfaceName ::= Name nameOfSurface | Name UNNAMED surfaceColor ::= Color r g b (REALs specifying overall color of surface) surfaceCells ::= Cells numberOfCells surfaceTriangulation ::= Triangulations numberOfTriangulations Each cell is described as follows: cellDefinitions ::= CellIndices i j k (NATs identifying cell) CellDegree degree (= 0, 1, or 2; degree of cell) DefiningFormula definingFormula A defining formula may be simple or compound. definingFormula ::= simpleDefiningFormula | compoundDefiningFormula simpleDefiningFormula ::= simpleOperation polynomial simpleOperation ::= IsPositive | IsZero | IsNegative compoundDefiningFormula ::= compoundOperation numberOfArguments definingFormula (the first argument) definingFormula (the second argument)... compoundOperation ::= And | Or numberOfArguments ::= NAT >= 1 Cell display information format. cellDisplayInfo ::= zeroD | oneD | twoD zeroD ::= Indices i j k (NATs) Degree 0 Position x y z (REALs) oneD ::= Indices i j k (NATs) Degree 1 Vertices numberOfVertices numberOfVertices instances of: {x y z (REALs)} twoD ::= Indices i j k (NATs) Degree 2 Vertices numberOfVertices numberOfVertices instances of: {x y z (REALs)} Triangles numberOfTriangles numberOfTriangles instances of: {vertex1 vertex2 vertex3 (NATs)} Each cellDisplayInfo must agree exactly in degree and indices with the corresponding cellDefinition. The degree and indices effectively act as a check that the cell display information matches up with the cell definitions. References [ARN83] D. A{\small RNON}, {\sl Topologically reliable display of algebraic curves}. Computer Graphics, 17, 1983.. [COH76] E. C{\small OHEN}, {\sl A method for plotting curves defined by implicit equations}. PROC SIGGRAPH '76 (cited in [GEI83]). [DAV82] J. D{\small AVIS}, {\sl Recursive ray tracing for the realistic display of solid models}. Report CAD-82-001 (M.Sc. thesis), Purdue CADLAB, May 1982. [FAR85] R. F{\small AROUKI}, {\sl The characterization of parametric surface sections}. MS, General Electric Company, June 1985. [GEI83] A. G{\small EISOW}, {\sl Surface interrogations}. Ph.D. thesis, University of East Anglia, July 1983. [GRI78] J. G{\small RIFFITHS}, {\sl A surface display algorithm}. Computer Aided Design, 10 (3), May 1978. [HAN83] P. H{\small ANRAHAN}, {\sl Ray tracing algebraic surfaces}. Computer Graphics, 17, 1983. [HAN84] P. H{\small ANRAHAN}, {\sl Beam Tracing}. Computer Graphics, 18, 1984. [KAJ82] J.T. K{\small AJIYA}, {\sl Ray tracing parametric patches}. Computer Graphics, 16, 1982, 245-254. [OSS83] S. O{\small CKEN}, J. S{\small CHWARTZ AND} M. S{\small HARIR}, {\sl Precise implementation of CAD primitives using rational parametrizations of standard surfaces}, Tech. Rept. No. 67, Dept. of Computer Science, Courant Institute, NYU, March 1983. [PET84] C. P{\small ETERSON}, {\sl Adaptive contouring of three-dimensional surfaces}. Computer Aided Geometric Design, 1 (1984), 61-74. [PRA85] V. P{\small RATT}, {\sl Techniques for conic splines}. Computer Graphics, 19, 1985. [ROT82] S. R{\small OTH}, {\sl Ray casting for modelling solids}. Computer Graphics and Image Processing, Vol. 18, 1982, pp. 109-144. [SNY83] W.V. S{\small NYDER}, {\sl Algorithm 531 - Contour plotting}. ACM Trans. Math. Softw., 4,3, (1978), pp. 290-294. [TIL85] W. T{\small ILLER}, {\sl Power series development about a regular point of an algebraic curve and its applications to geometric modeling}, Talk, SIAM Conf. on Geometric Modeling and Robotics, Albany, NY, July, 1985.