00001.tif

PDF

Download

Open document

Transcript (if available)

Content A GENERAL APPROACH TO THE DESIGN OF ROBUST ALGORITHMS FOR GEOMETRIC MODELING by Amitabh Agrawal A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Computer Science) May 1995 Copyright 1995 Amitabh Agrawal UMI Number: DP22897 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is d ep en d en t upon the quality of the copy subm itted. In the unlikely event that the author did not sen d a com plete m anuscript and th ere are m issing pag es, th e se will be noted. Also, if m aterial had to be rem oved, a note will indicate the deletion. Dissertation Publishing UMI D P22897 Published by P roQ uest LLC (2014). Copyright in th e D issertation held by the Author. Microform Edition © P roQ uest LLC. All rights reserved. This work is protected against unauthorized copying under Title 17, United S ta tes C ode P roQ uest LLC. 789 E ast Eisenhow er Parkway P.O. Box 1346 Ann Arbor, Ml 4 8 1 0 6 -1 3 4 6 UNIVERSITY O f SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES, CALIFORNIA 90007 This dissertation, written by Am t abh p Agrawal............; .............................. under the direction of hXs. Dissertation Committee, and approved by all its members, has been presented to and accepted by The Graduate School, in partial fulfillment of re quirements for the degree of DOCTOR OF PHILOSOPHY Dean o f G raduate Studies Date 1 9 .9 .5 ............ DISSERTATION COMMITTEE (\ f\ ft C hairperson D ouglas J . I e r a x d i........... Dedication This thesis is dedicated to my parents Ramesh and Pushp Lata. ] Acknowledgments : < | I would like to thank Professor Ari Requicha for bringing my attention to the j problem of floating point inaccuracies in geometric algorithms. He provided invaluable j support and guidance throughout the course of my research that led to this dissertation. I Without his guidance this thesis would have been an incomplete, hodge-podge set of j ideas. He helped temper this thesis into a more coherent piece of work. ' He taught me the value of honesty and meticulousness in research. I am indebted to him (amongst other things) for repeated proofreading of this thesis. Doug Ierardi was a constant source of ideas ranging from algebraic geometry to coffee houses. A1 Inselberg taught me how to look at things in different perspectives, and Richard Weinberg forced me to write the first page. 1 I am indebted to Dr. Ken Goldberg and Professor Yorgos Papavassilopoulos for ; their suggestions during different stages of this research. j At this stage I would also like to thank my colleagues and comrades-in-arms Tim Whalen and Antonia Spyridi who made the stay in graduate school worthwhile. Grad- school wouldn’t have been the same without them. Special mention goes to Sanjiv Patel ! I for long nights of practical advice and pin-ball games. I wish to thank all the fellow students, past and present members of the ! Programmable Automation Lab at USC and friends for their support and friendship during my stay at USC. Thank you Kamen, Jan, Miguel, Rolf, Glauco, Shahril, Mani, Bharat, Giovanni, Jung, Sung-Don, Victor, Desiree, KC, Gonj, Mansi, Alex, Joe, Bin, \ Gary, Loma, Raj, Dana and Manav. ! ( And last but not the least, my thanks go to Ms. Kusum Shori who made this place j feel more like home. i j Contents i Chapter 1 In tro d u ctio n .....................................................................................................1 1.1 Motivation . ........................................................................................................ 1 1.2 Source of E rro rs ............................................. 3 1.2.1 Nearness and Id e n tity ........................................ 5 1.3 Validity and Correctness....................................................................................... 5 1.3.1 Representation V a lid ity ............................................................................. 6 1.3.1.1 Symbolic and Numeric V alidity........................................................7 1.3.1.2 Approximate V a lid ity ........................................................................9 1.3.2 Correctness...................................................................................................10 1.3.2.1 Approximate C orrectness................................................................11 1.3.2.2 Approximate Numeric C orrectness.............................................. 12 1.4 Problem Statem ent.......................................... 14 1.5 Scope of the W o rk ............................................................................................... 15 1.6 O rganization.................................................... 15 Chapter 2 Previous W o rk ............................................................................................... 16 2.1 Correctness..................................................................................................... . 16 2.1.1 S trong Correctness......................................................................................16 2.1.2 Approximate C orrectness................... 17 2.1.3 Approximate Numeric C orrectness......................................................... 18 2.1.4 Implicit C orrectness.................................................................................. 19 2.2 Previous W o rk ...................................................................................................... 19 2.2.1 Exact Arithmetic Approaches .................................................... 20 2.2.1.1 Green and Yao [Green &Yao 1986]...............................................20 2.2.1.2 Suglhara and Iri [Sugihara & Iri 1989]............................................20 2.2.1.3 Milenkovic [Milenkovic 1988].........................................................20 2.2.1.4 Karasick, Lieber and Nackman [Karasick et al. 1990]................ 21 2.2.2 Error Analysis Approaches........................................................................ 21 2.2.2.1 Purdue / Cornell Approaches............................................................21 2.2.2.2 Milenkovic [Milenkovic 1988a, Milenkovic 1988b, Milenkovic 1989].................................................................................................. 22 2.2.2.3 Li and Milenkovic [Li & Milenkovic 1990]..................................22 2.2.2.4 Fortune [Fortune 1989]................. 22 2.2.3 Tolerance Based Approaches . . . .....................................................22 2.2.3.1 Salesin [Salesin et. al 1989, Salesin 1991].....................................23 2.2.3.2 Segal and Sequin [Segal & Sequin 1985, 1988, Segal 1990] . 24 2.2.3.3 Briiderlin and Fang [Briiderlin 1990, Briiderlin 1991, Fang 1992] ......................................................................................................24 2.2.4 General A pproaches............................................. 27 2.2.4.1 Yap [Yap 1 9 9 3 ]........................... 27 1 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 2.2A.2 Fortune and Wyk [Fortune & Wyk 1993, Fortune 1993] . . . 21 2.2.5 Related W ork.............................................................................................. 29 Chapter 3 Math Preliminaries......................................................................................... 30 3.1 Interval Arithmetic................................................................................................30 3.1.1 Interval Numbers................................... 30 3.1.2 Interval Arithmetic.................................................................................... 31 3.1.3 Floating Point Representation . . . ......................................................31 3.1.4 Floating Point A rith m etic.................. 32 3.1.5 Rounded Interval Arithmetic............... 32 3.1.6 Intervals approach exact results as precision in cre ases...................... 33 3.1.7 Implementation of Interval Arithmetic....................................................35 3.1.7.1 Standard Interval Arithmetic............................................................ 36 3.1.7.2 Range Arithmetic................................................................................36 3.1.7.3 Static Error B o u n d s .................... 36 3.2 Verifying Polynomial Identities.................... 37 3.2.1 Non-probabilistic A pproach....................................................................37 3.2.2 Probabilistic Approach. .................. 38 3.2.2.1 Verifying Rational Polynomial Identities...................................... 39 3.3 Verifying Polynomial Identities on a V a rie ty .................... 39 3.4 Gap T h eo re m .................................................. 40 Chapter 4 Approach..........................................................................................................42 4.1 Accidental Coincidences...................................................................................... 42 4.1.1 Forced Coincidences............................ 45 4.1.1.1 Coincidence by Design................. 45 4.1.1.2 Coincidence by Creation. . . .; ..................................................... 46 4.1.1.3 Coincidence By Theorem................................................................. 46 4.2 Justification............................................................................................................. 50 4.3 Approach............................................................ 51 4.3.1 Interval Arithmetic Module P t . . . .■..................................................... 55 4.3.2 Exact Computation Module PE ....................................................56 4.3.3 The Controller C r .....................................................................................57 4.4 P ro o f................. 59 4.4.1 Proof of Termination................................................................................. 59 4.4.2 Proof of Correctness............................ 64 4.4.3 Handling of Division Operator.................................................................68 4.5 Handling of Perturbations.............................. 71 4.5.1 Arbitrary Direction of M otion .................................................................72 4.5.2 Specific Direction of Motion.................................................................... 74 4.6 Solve by G eneralization................................. 74 Chapter 5 Implementation...............................................................................................75 5.1 Rationale for an Implementation................... 75 5.2 Polygonal M odeler...............................................................................................76 5.3 Correctness....................................................... 77 5.4 Profile T estin g ......................................................................................................77 5.4.1 Interval G ro w th .........................................................................................83 5.4.2 Time Spent in Interval Computation .....................................................87 5.4.3 Time Spent in the Recompute Routine.....................................................88 5.4.4 Speed of the Current Im plem entation.................................... 88 5.5 Implementation of Multi-Precision Interval Arithm etic.................................89 5.5.1 The IEEE Floating Point Standard Flags................................................. 89 5.5.2 Data S to rag e...............................................................................................90 5.5.3 Usage............................................................................................................ 91 5.5.4 Implementation............................................................................................93 Chapter 6 History Management..................................................................................... 94 6.1 Basic D esign.................................................... 95 6.1.1 Parent Pointer M anagem ent................ 95 6.1.2 Children Pointer M anagem ent..................................................................98 6.2 Ripples of Inheritance.......................................................................................100 6.2.1 Parent Pointer M anagem ent................................................................... 100 6.2.2 Children Pointer M anagem ent................................................................101 6.3 Pointer Maintenance for Subobjects................................................................101 6.3.1 Parent Pointer M anagem ent................ 102 6.3.2 Children Pointer M anagem ent................................................................102 6.4 History Graph T raversal................................ 102 6.4.1 Constructor (Ob j e c t : : Ob j e c t ( ) } ...................................................103 6.4.2 My Constructor ( v o id Ob j e c t : : M y C o n s tru c to r () ) . . . .103 6.4.3 Perturb ( v o id Ob j e c t : : P e r t u r b ( ) ) .........................................105 6.4.4 MakeConsistent (SO utcom e O b j e c t : :M a k e C o n s is te n t ()). . 106 6.4.5 RecomputeDeps (SO utcom e R o o t % :R e c o m p u te D e p s ()) . .106 6.4.6 Recompute (SO utcom e Ob j e c t : : R e co m p u te ( ) ) ..................106 6.5 Levels of Granularity....................................... 107 Chapter 7 Coincidence by Theorem: An Im plem entation..................................... 110 7.1 Exact arithmetic modulo p ................................................................................ 110 7.2 An Exact Class for an Approximate C lass......................................................110 7.3 Addition of data elements in an O b je c t......................................................... I l l 7.4 Coincidence by Theorem................................ 112 7.4.1 Randomize ( v o id Ob j e c t : : R a n d o m ize () )............................. 112 7.4.2 ComputeExact (Ob j e c t E Ob j e c t ? : C o m p u te E x a c t ()) . .113 vi 7.4.3 My Constructor ( v o id Ob j e c t : :M y C o n s tr u c to r ()) . . . . 7.5 Exam ples............................................................................................................ Chapter 8 C o n clu sio n s............................................................................................... 8.1 Contributions................................................. .................................................... 8.2 Lim itations............................................................................................... . . . 8.3 Future R e se a rc h ............................................................................................... Chapter 9 Bibliography List of Figures 1 Abstraction Process...............................................................................................1 2 Geometry Programmer’s nightm are .....................................................3 3 Nearness Implies Id e n tity .............................1 ......................................................5 4 Algorithm vs. F u n ctio n ........................................................................................6 5 Invalid Representation .................................. 7 6 Representation of a Line Segment ....................................................................7 7 Numerically Valid Representations .........; ....................................................... 9 8 Approximate Valid Representations..................................................................9 9 Algorithms as Mappings between Configuration S p a c es.............................10 10 e-Correctness....................................................................................................... 11 11 Correctness............................................................................................................13 12 Strong Correctness ........................................ 16 13 Polygon Intersection.......................................................................................... 18 14 Implicit Correctness .......................................................................................... 19 15 Epsilon G eom etry.............................................................................................. 23 16 Briiderlin ’ s Tolerances ..................................................................................... 24 17 Angle vs. Error ............................................. 27 18 How to get from P to Q ..................................................................................... 42 19 Accidental Coincidence in the Intersection of three half-spaces ............... 45 20 Coincidence by creation .. ..................... 46 21 Pappus’ theorem ............................................................................................47 22 Computation with R adicals......................... 48 23 Decision Tree for Circle Intersection P rogram .............................................. 52 24 Decision Tree for sorting points (P i)................................................................53 25 The Program as a decision tr e e ................ 54 26 Decision Tree .......... 60 27 Input Space Partition...................................... 61 28 Decision tree for classification program .........................................................65 29 Configuration Space of the branch fu nction.................................................. 72 30 Cornered .................................................................................................. 73 31 Random intersection of half-spaces ............................................................... 78 32 Connecting R o d .................................................................................................. 79 33 ANC 101 (front v ie w ).................................. 80 34 Glue G u n ..............................................................................................................81 35 Trammel ................................................ 82 36 ANC 101 interval distribution.......................................................................... 84 37 Connecting Rod interval distribution ............................................................. 84 38 Glue Gun interval distribution ........................................................................ 85 39 Trammel interval distribution.......................................................................... 85 40 Pappus Theorem interval distribution............................................................. 86 Vlll 1 41 Brianchon Theorem interval distribution 86 | 42 MyNumber (straightforward multiprecision im plem entation).................... 91 j 43 Storage of Multiprecision num bers 91 j 44 Handling of N a N s...............................................................................................92 j 45 Program as decision tree (Figure 25 reproduced).......................................... 94 j 46 Object Modification vs. Object Construction............................................... 108 ; 47 Pappus T h eo rem ...........................................• .................................................. 116 48 Theorem of Brianchon .................................................................................. 117 Abstract I 1 Geometric modeling systems are rapidly displacing traditional drafting as the I primary means for defining the geometry of mechanical parts and assemblies. 1 I I j Geometric modeling has evolved significantly over the last two decades, and modem j I j systems are very powerful. However, like their older brethren, they suffer from | reliability problems, which are caused by numerical errors in floating point i : computations. Difficulties arise primarily when there are geometric singularities or j coincidences — for example, when two planes coincide or a point lies on a surface. ! Typically, two real numbers are compared and a program branches depending on the \ results of the comparison. If the numbers are almost equal the results are unpredictable. In practice, approximate comparisons are used by introducing epsilon “fuzz factors”. I The epsilons are determined empirically, without underlying theoretical foundations. This causes inconsistencies, program crashes, and wrong results. This thesis proposes a novel approach to the robustness problem in geometric modeling. This approach brings together ideas and techniques from interval arithmetic, i constraint management, randomization, and algebraic geometry. It acknowledges that i input values have tolerances, that objects within tolerances are equivalent, and that ; certain geometric singularities must be maintained because they reflect design intent or : the laws of geometry. Unlike most other research on robustness of geometric , algorithms, this approach is systematic, general, and can be applied almost mechanically i to a great variety of problems, specifically, to any geometric algorithm using the operations +, -, * and /. The required theory and algorithms have been developed, and ' the viability of the concepts has been demonstrated by an experimental implementation i involving planar halfspaces in a 2-dimensional space. The implementation focuses on 1 I algorithms for computing the boundary of two dimensional solids defined through * I boolean operations on primitives bounded by lines. j x Chapter 1: Introduction This introductory chapter first presents the motivation for solving the problem of robustness in geometric computations. Then, floating point errors incurred in a geometric computation a*e discussed in Section 1.2. Validity of representations and correctness of algorithms are addressed in Section 1.3, leading to the problem statement in Section 1.4. The scope of the thesis follows in Section 1.5, and Section 1.6 covers the organization of the thesis. 1.1 Motivation The process of modeling real world entities can be described as shown in Figure 1 [Requicha77]. When the real world is studied, mathematical abstractions (or mathematical models) are formed leading to theories about the real world. To compute on these mathematical abstractions, representations for them are devised and manipulated through algorithms. For example, the solids in the real world correspond to r-sets1 which correspond to B-Reps (boundary representations) or CSG (constructive solid geometry) representations for r-sets. Ideally there is a one-to-one correspondence between the mathematical abstractions and the representations. Therefore, the extent to which the results computed by these representations correspond to the real world, is the extent to which the theory (mathematical abstraction; is correct. However, in the field of geometric modeling the correspondence between the mathematical abstractions and their representations breaks down. Mathematical Abstractions Real World Representations \ Figure 1: Abstraction Process L An r-set is a subset of E3 that is bounded, closed, regular, and semianalytic [Requicha 1977]. 1 The floating point computation and representation is at fault. Most mathematical abstractions (especially in the field of geometric modeling) assume the existence of a Real Euclidean space with real operations, whereas the representation space is floating point space (with a finite cardinality) with floating point operations. A simple example will illustrate the issues. In Figure 2, a polygon is intersected with a half-space defined by a-line. The polygon is defined by a list of directed edges such that the material lies on the left side of each edge. The program computes the intersection as follows. The line / is intersected with the edges of the polygon and the intersection points a, b, c and d are sorted along the line according to the corresponding values of the variable t used to parameterize the line. .The intersected polygon is then computed by starting at a\ traversing the original polygon in the direction of the arrows C a, e,f, g,d) coming to d; following I to the next sorted point c; traversing the polygon (c,j, b) to come to b; and then once again following f to come to an end at a where we terminate the cycle. This is the correct, desired intersection. However, in a floating point implementation, each arithmetic operation incurs an error. Therefore, when we are comparing the parameter values for b and e, the comparison might give an erroneous answer stating that ,, when tb > t . This simple error in our decision gets magnified to give an answer of <a, e,f, g, d, b, k, m, a>, which is grossly wrong. We could go on about how this error could get further amplified, and examples abound on how easily programs can be led astray by floating po nt computations. In practice, programmers attempt to avoid similar" errors by using epsilon “ fuzz factors” in comparisons between real variables. In the past, finding epsilons for various geometric computations was a black art. They were derived by empirical experiments and programmer’s intuitions. There exists a folklore in the PADL underground about how a major company acquired the PADL-2 solid modeler, which is written in single precision, and converted it into double precision. As soon as this was done, nothing worked, and it took several man-months to get PADL-2 working again. 2 k i m 7 7 8 Figure 2: Geometry Programmer’s nightmare In the late 1980s work began in earnest in the field of robustness in geometric algorithms in the presence of floating point errors, and since then the field has evolved steadily but slowly. Much remains to be done. 1.2 Source of Errors There are three sources of errors in geometric pro grams. 1) Real world data is imprecise. The temperature of water is not 67° but 67.2093320923...° and cannot be measured exactly. 2 Therefore, any computer representation of the real world object is an approximation and therefore, in error. 2) Some numbers cannot be represented exactly in the desired representation scheme. For example, 1 /3 does not have an exact decimal floating point representation. A The situation becomes even more complicated if we adopt a quantum- mechanical model of the world. 3 3) Floating point computation introduces errors of its own. If the input to a geometric program is considered to be exact, i.e. if the error of Type 1 above is ignored, there are two types of errors incurred during a geometric computation. (We will discuss Type 1 errors later, see Section 4.2) 1) Numerical errors when a Real value is calculated. A bound can be put on the error of the computed value, and, if the need arises, the bound can be arbitrarily tightened by redoing the computation in higher precision. This is the domain of standard numerical analysis. 2) Logical errors when a branching decision is made based on a comparison of computed values. In this case, the resulting value is not in the neighborhood of the actual value and the standard bounds are not useful. These are the errors which cause programs to crash. Figure 2 and the associated discussion provide an example. Let us look at logical errors more carefully. A geometry program can be looked upon as a decision tree, where at each node a tuple of the following type is calculated <f(x) < 0,/( * ) = 0,f(x ) > 0) What we compute, however, is not fix ) but / (x) = f i x ) ± e, because of numerical errors. Here e is a small positive number. Therefore, wrong decisions can be made in the region where \f (x) | < e . These are the errors that we should try to avoid. Note that fix ) ~ 0 situations correspond to geometric singularities or degenerate cases, e.g., parallel lines, coincident points, etc. In geometric modelers a lot of code deals with the case of fix ) = 0 . Consider a simplistic analysis where a program has been completely unfolded into a decision tree with 10 levels. If all the nodes have just two branches instead of three, then the number of leaves is reduced by a factor of about 50. This, in conjunction with the previous 4 paragraph, shows that most of the code is written for the cases where the result may not be correct. Clearly there is something wrong! 1.2.1 Nearness and Identity m Figure 3: Nearness Implies Identity Programmers typically deal with floating point comparisons by using fuzz factors, which essentially imply that if two floating point numbers are near then they are considered the same. This probably has its historical basis in trying to represent and manipulate two touching solids, where at a superficial level it makes sense. In geometry, however, it can lead to topologically inconsistent decisions, even with exact arithmetic. Consider the problem of finding the arrangements of lines on a plane in Figure 3. (A planar line arrangement is a decomposition of the whole plane into a minimal collection of non-overlapping polygons whose vertices are the intersections of the lines.) Here lines I, m and n are intersecting in three points A, B and C. If the distance between A and B is less than e , A and B will be merged into a single point whereas A and C will not be (their distance is greater than e). Therefore, we end up in a situation where two lines intersect in two points. Note that this is not a problem with the floating point computation but rather with the hypothesis that near implies identical (Nil). 1.3 Validity and Correctness i : This section presents and compares different notions of representation validity and algorithm correctness. Section 1.3.1 defines the concepts of symbolic and numeric validity (1.3.1.1), e-validity (1.3.1.2) and validity within a neighborhood (1.3.1.2). Section 1.3.2 introduces the concepts of approximate correctness (1.3.2.1) and 5 approximate numeric correctness (1.3.2.2). Approximate correctness is similar to the notions of correctness used by Hoffman et al. [Hoffman et al. 1988, Hoffman 1989], whereas approximate numeric correctness is a new definition of correctness used in this thesis. The notation used in this exposition is as follows. Sets will be denoted by capitals and instances by lower-case letters. For example, a modeling space is denoted by M and a particular model by m. Figure 4 elaborates on Figure 1 by showing the domain and range of the functions and the algorithms. When a theory is formulated, a function/ is defined that acts on an input mathematical model m ,- and maps it to an output mathematical model mQ . For example, m{ - might be an r-set, m0 a Real number a n d /th e function that associates to a solid its volume. Computationally, representation schemes Si and sQ associate representations rt and rQ to the mathematical models, so that the output rQ is produced by an algorithm a operating on rL . Math Models Representations Figure 4: Algorithm vs. Function 1.3.1 Representation Validity A representation r is valid, if there exists a mode! m belonging to M such that s (m ) = r (Figure 4). The question of validity arises because s is usually not an onto mapping and there may exist r in representation space R which do not have corresponding models. For example, in Figure 5 the boundary representation corresponds to a self intersecting object which does not belong to the modeling space of solids. 6 Figure 5: Invalid Representation We assume throughout this thesis (unless we specifically state otherwise) that representation schemes are unambiguous, i.e., given a valid r there exists only one model m that corresponds to it. 1.3.1.1 Symbolic and Numeric Validity We will illustrate the issues with a simple representation scheme for wireframes, i.e., for collections of line segments. Figure 6 shows a representation (explained below) for a wireframe consisting of two line segments. A representation consists of two types of data: numeric data and symbolic data * The numeric data is the coordinate information in a representation, typically real numbers. In the example of Figure 6, this corresponds to the three pairs of numbers (0.0,0.0), (2.0,1.0), (0.2,0.9) which are the coordinates of points A, B and C. • The symbolic data captures the connectivity information. In our example (see the information bracketed by “Sym” in the figure), 3 refers to the number of points in the representation, 2 refers to the number of edges. The first pair-of- integers entry signifies that the first edge connects point 3 (C) to point 1 (A). Figure 6: Representation of a Line Segment 7 Similarly, the second entry 1 2 denotes that the second edge connects vertices 1 (A) and 2 (B). Suppose there are m parameters in the numeric part of a representation. Then the numeric data corresponds to a point in We say that is the configuration space (C-space) associated with the representation. Perturbing the values of the numeric data yields other points of the C-space. However, • Not all points in the C-space correspond to valid representations. A similar problem arises with parameterized objects. For example if we changed the last line of numeric data in Figure 6 to (0.0,0.0) the representation would be invalid because the edge AC would collapse to a point. • Defining a representation scheme does not determine a single C-space. There are many C-spaces associated with each representation scheme. For example, a triangle and a line, segment represented in out wireframe scheme have 6 and 4 parameters respectively and therefore the C-spaces will be 9t6 and S R 4 respectively. A representation is symbolically valid if there exists some assignment of numeric data that makes the representation valid. Such a representation is imbeddable in 2-space (or the appropriate Euclic ean space). A symbolically valid representation is numerically valid if the representation is valid, i.e., if the associated numeric data results in a valid representation. Note that a numerically valid representation is always valid and hence symbolically valid, but a symbolically valid representation can be numerically invalid, and hence invalid. A symbolically valid representation r defines a mathematical relation from (a subset of) the modeling space M to the C-space %m (Figure 7). In the figure, R is the space of numerically valid representations that have the same symbolic part as r, and correspond to the hatched subset of the modeling space. 8 Figure 7: Numerically Valid Representations 1.3.1.2 Approximate Validity In this subsection, approximately valid representations are defined which are near a valid representation. These are important because in practice validity is usually not guaranteed due to numerical inaccuracies. The representation r :s e-valid if there exists a representation r ' • with the same symbolic data as r and numeric data within an e-ball of r in C- space, and • r' is valid. Note that for r to be e-valid, it has to be symbolically valid. Figure 8: Approximate Valid Representations Figure 8 extends Figure 7 by growing the R set in 9V” by e to give W. This is the set of all approximately valid representations for that symbolically valid representation. 9 We can generalize the above definition as follows. Let Xj, x2,...xm be the variables associated with a representation’s numeric data. We can associate with the representation a numeric uncertainty X by saying that (x v x 2.. .xm) e l c C-Space, where X = x / 2. ..lm and I j = [ypZj], ypzj e 31 are intervals. Attaching a low and high limit to each parameter in the representation is analogous to parametric tolerancing with + / - limits. We could generalize the notion by using arbitrary subsets X, but m- dimensional intervals will suffice for our work. We define a representation r to be valid with an uncertainty X, if there exists a valid representation r' with the same symbolic data as r and with numeric data x e X cz C- space associated with r. 1.3.2 Correctness In Figure 4, the algorithm a is correct if for all valid belonging to the representation space/?, a(rj) - rQ , such that ri and r0 are valid, unambiguous representations of models m; and mQ , and f(mi) = m0. Suppose we have a correct algorithm a, which does all its computations exactly, even on Reals (i.e., we have a Real RAM model of computation). We want to study what happens when some of these computations are done In floating point arithmetic. Because exact measurements are non-realizable in reiil life and the floating point errors are small, we must study( the neighborhoods of representations. R o 3 > R o 2 R ol Figure 9: Algorithms as Mappings between Configuration Spaces For a specific input representation there is a set of parameters Xjjc2— Xi that define the associated C-space 9v*. Let /?; be the set in 9 1* corresponding to the valid representations (see Figure 9). The algorithm a maps mutually disjoint subsets of /? ,■ into 10 different output C-spaces 9toi, 9t°2 and 9t°3 (corresponding to R ol, Ro2, ^ 03) • (For example, a polygon intersection algorithm may produce polygons with different number of edges and vertices and have different output C-spaces, even for input polygons in small /?,• regions.) Suppose, then, we take a valid input representation r with numeric data x e 91' and carry out the computations prescribed by a (correct) algorithm using infinite precision. The result is a valid output representation r' with numeric data y e 91° where 91° denotes the appropriate output C-space < 3iO J . W e call this a numerically exact computation. ' ■ 1.3.2.1 Approximate Correctness Approximate correctness is a more rigorous version of the notion of correctness proposed by Hoffman et al. [Hoffman et al. 1988, Hoffman 1989]. We start with an input representation rt that is er valid, i.e. 3 ^ such that st is valid with corresponding model mt and the distance between the numeric data of rt - and S; is less than e,. The algorithm is approximately correct if it produces an output rQ that is e0-valid and the corresponding model mQ satisfies mQ = / ( m t), w h ere/is the function that we are trying to compute. Recall that, per our definition of e-validity, rt and st must have the same symbolic part and similarly for rQ and sa. m , < £ / Figure 10: e-Correcti:ess 11 Note that rt could have an st and an s'; within an tj ball in input C-space which map to different output C-spaces (Figure 9). Therefore, two implementations that are e- correct could give answers that differ in the symbolic part of the data. Approximate correctness can be generalized using the concept of representation with an uncertainty: If the input representation rt - is valid with uncertainty X (i.e. there exists an instance x of the numeric data in X that is valid) then the algorithm is approximately correct if it produces an output rQ that is valid with an uncertainty Y and the corresponding model satisfies m 0 = / ( m 4 ) . 1.3.2.2 Approximate Numeric Correctness A slighdy different notion of correctness is the following. Consider a representation ri and an uncertainty neighborhood X; associated with its numeric data X j. Apply some inexact computation a to r,- to get a result rQ , with numeric data x Q . In addition an uncertainty X 0, for xQ is also computed. Now a is approximately numerically correct, if 3 r’-, r'Q with x't e X i and x'0 e X Q such that r’ 0 is the result of applying the exact computation to r\. Furthermore, r \ and r’ 0 have the same symbolic data as rf - and rQ respectively. Note that nothing has been said about the validity of these representations. However, if x ‘ i happens to be valid and the algorithm is correct (in the usual sense) when implemented through exact computation, then x ’ Q is also valid. So, this notion of approximate correctness in essence decouples the two issues of validity of representations and numerical behavior of algorithms. Things are easy if the input data are always valid, e.g., in solid geometry, when the input is CSG; or in problems like convex hull computation, when the input is a set of points. Then, approximate numerical correctness implies approximate correctness for algorithms which are correct when executed in exact arithmetic. Input validity for x,- is automatic and output validity of x'0 is ensured by the fact that the exact computation implements a correct algorithm. In cases where input validity is not automatic the approximate numeric correctness criterion must be augmented with a proof thatx't - is valid, to get approximate correctness. 12 The example below illustrates this definition of correctness. Suppose the problem is to compute the various intersection points of two lines in a plane. (Because this is a contrived example we actually know what the number of intersection points is, but in a more complicated example that may not be a priori known.) This corresponds to the symbolic part of the representation. I m (a) Figure 11: Correctness In Figure 11(a) / and m are the lines whose intersection point(s) need to be calculated. In accordance with the definition of approximate numeric correctness, the lines are replaced with the intervals L and M. A program would be considered approximately numerically correct if • the symbolic result (i.e., the number of intersection points) is calculated correctly for some instance of the objects within the tolerance. Here /, and mj are instances within the tolerance zones and exact computation of their intersection would produce a single point. Therefore, one is an acceptable symbolic part for the output of an approximately numerically correct algorithm. • the numerical result has an output tolerance zone that contains the exact intersection of two input lines in the tolerance zones L and M. In this example, the output p k has a zone (shown here as a circle) containing the actual intersection point of / , ■ and m y. Note that the size of this zone can be made arbitrarily small if desired (by computing with more bits). 13 Now we are ready to define the notion of approximate numeric correctness more rigorously. An exact geometric computation can be viewed as a mapping C:91n — > 91m x ZP where Z is the domain of integers and 91 is the domain of reals. Here 91m and 91" are C-spaces for the input and the output. For example, if we are calculating the boundary representation of a triangle defined by CSG on 2D linear half spaces, a half-space configuration will be defined by a point in 919. (Specifically, each line is defined by 3 parameters and therefore, 3 lines are defined by 9 parameters.) Each output vertex is defined by 2 coordinates, and therefore the three vertices are defined by a point of 916. (The symbolic output which is the number of vertices is defined in Z.) ( 9t Y * ( 3t V” p Let C ':1 2 I — >12 J x Z where both the input and output Real sets have been replaced by their power sets, but the integer output (corresponding to the symbolic result) stays the same, and C' is a floating point implementation operating on floating point numbers. Then C" is approximately numerically correct if the following condition is satisfied. C'(X1,X2,X3..Xn) = (Yj, Y2...Ym,J1,J2..Jp) withx. c 91,Y. c 9t, Jk e Z , Vi, y, k , and there exists x. < = X ., y. e Y. such that * * J J C (Vp x 2-• -xn) (yp y2- • -ymi J p ^ 2 " ^ ere ^ are input un certainties and Yj are the output uncertainties. Note that C ' maps input uncertainty zones into output uncertainty zones, whereas C maps single instances to single instances. An approximately numerically correct computation with a very large Y0 is useless in practic e. Note, however that Y0 can be reduced arbitrarily by using higher precision. j- 1.4 Problem Statement We are now ready to state precisely the main problem addressed in this thesis. In the statement below we use the definitions of C, C' and correctness given in Section 1.3.2.2. Given a C defined as a program with if statements, goto statements and algebraic operations +, x, - and - 4 - on real variables, design an imple- 14 mentation o f C' operating on floating point numbers and using floating point operations which is approximately numerically correct. 1.5 Scope of the Work The paradigm proposed in this thesis is simple, and is more general than the previous work done in the area. It applies to all programs with a certain structure rather than to a specific geometric problem. Specifically, we assume that algorithms can be modeled by decision trees, where arcs correspond to arithmetic operations and nodes to branching decisions made by comparing real numbers. The paradigm was designed to tackle robustness issues in solid modeling, but can conceivably be applied to other, non geometric problems. We assume that input data carry tolerances and that the results of computation need only have finite accuracy. These assumptions are valid in most engineering applications, but may be inappropriate in other domains, e.g. for certain mathematical computations that cannot tolerate inaccuracies. The prototype implementation deals with linear half-spaces in 2D, but the approach is applicable in any dimension. Curved objects require further development of algebraic geometry algorithms, which are needed in one of the steps of the approach. 1.6 Organization The remainder of this thesis is organized as follows. Previous work in the area is covered in Chapter 2. Chapter 3 presents the basic mathematical background necessary for the thesis. This includes various topics such as Interval Arithmetic, Polynomial Identity Detection, Polynomial Identity Detection on a Variety, and Gap theorems. The overall approach to the problem is discussed in Chapter 4. Chapter 5 presents the implementation. History management is discussed in1 Chapter 6 and theorem detection in Chapter 7. Chapter 8 identifies the contributions of this thesis and possible extensions. The Bibliography follows. 15 C hapter 2: Previous W ork This chapter covers the previous work reported in the literature. Notions of correctness of robust algorithms are discussed in Section 2.1, which presents several new concepts not discussed in section 1.3. Robust algorithms that appeared in the previous literature constitute the subject of Section 2 2. These algorithms are categorized by the approach used and discussed in subsections. The subsections are further subdivided according to specific works. > 2.1 Correctness When dealing with algorithms based on inaccurate primitives, correctness is redefined to encompass algorithms which are acceptable but not exact or correct in the usual sense. The various notions are described below in increasing order of generality. R t 1 Figure 12: Strong Correctness 2.1.1 Strong Correctness This is the notion of correctness used in the works of [Li & Milenkovic 1990] and [Salesin et al. 1989]. In both references strongly convex hulls are computed. A polygon is e-strongly convex, if the convexity property is invariant to perturbation of the vertices by e. Note that this is a stronger condition than that guaranteed by an exact convex hull algorithm. Consider the problem of convex hull evaluation of a set of points S. The correct result computed with arbitrary precision for this case is C and the computed hull (with finite precision) is C'. C' is a e-convex 8-hull of S if d(C', C) < 5 (where the distance between two polygons D and E is defined as follows: d (£>, E) = maxy (pe G £) (P> 0 ) ) ) » here D and E are considered as 16 the two dimensional sets enclosed within the boundary defined by the list of vertices), and the predicate PE (C') is true. PE means that a polygon is convex even if all the vertices are perturbed by a small amount e. This is a stronger condition than merely asking that the polygon C" be conver. More generally, let us assume in Figure 12 that the function/m aps its input ra; into an output mQ that satisfies predicate P. Now, algorithm a is strongly correct if it maps the input representation rt into an output rQ such that the model m ’ 0 that corresponds to rQ is within 8 of m0, i.e. d (m'0, m 0) < 8 , and for all m in M0 such that d (m, m'Q ) < 8 , P (m) is satisfied. Note also that mi may be non-representable and therefore may actually be a representation for some model close to mf - . This notion is important because the output from these procedures can be used as an input to other procedures and correctness can be guaranteed. For example, it is easy to write a floating point rotation routine that preserves the convexity of a strongly convex polygon. The floating point rotation produces a polygon which must be convex because it is within e of the strongly convex result of the exact, otation. Note that the output from a strongly correct program can be used as an input to an approximately correct or approximately numerically correct program. This is explained in the section on approximate correctness. 2.1.2 Approximate Correctness This notion has been covered in detail in Section 1.3.2.1. An algorithm which is e- approximately correct is called type-4 robust [Dey et al. 1992] or stable (if e is small) [Fortune 1989]. For e = this notion is called correctness [Hoffmann 1989], type-2 robustness [Dey et al. 1992] or robustness [Fortune 1989] . Note that this notion of correctness may not be used in nested algorithms [Hoffmann et al. 1988]. For example, assume the existence of a line-segment intersection algorithm a that is approximately correct (say £ is large). Now, ‘riangles ABC and DEF are intersected (see Figure 13). The procedure calls a to get the intersections of line segments. It returns null with EF and BC and G with DF and BC as arguments. Therefore, a straightforward implementation of the intersection algorithm using approximately correct primitives might give an incorrect answer (in this case, perhaps 17 DGCABG). The algorithm relies on a series of local decisions to produce a global solution, but approximately correct answers to the local problems may produce gross global errors. B Figure 13: Polygon Intersection However, an approximately correct algorithm A can use the output of a strongly correct procedure B. A finds a solution for some objec t in the neighborhood of its input. A strongly correct procedure B guarantees the satisfaction of some property pB for all points in some neighborhood of the input to A. Therefore, any point chosen from the neighborhood of the input for A will satisfy the property pB and the algorithm will work correctly. 2.1.3 Approximate Numeric Correctness The notion of approximate correctness is generali sed further to separate the notions of validity and correctness (details in Section 1.3.2) to give the notion of approximate numeric correctness. This notion also does not allow for nested algorithms. (This is expected, since approximate numeric correctness is a generalization of approximate correctness.) For the domains addressed in most of the previous work (Sections 2.2.1.1, 2.2.1.2, 2.2.2.4, 2.2.3.3), the notions of approximate numeric correctness and approximate correctness are equivalent. 18 Figure 14: Implicit Correctness 2.1.4 Implicit Correctness Another way of generalizing approximate correctness is the following. In Milenkovic’s hidden variable method for calculating line arrangements, tolerances (called channels) are defined for each line, and any ntonotonic curve within tolerance satisfies the definition of line. Here the concept of line has been generalized to a monotonic curve to simplify the problem. Therefore, rather than approximating the object in representation space we approximate it in the modeling space. The problem in generalizing such a definition of correctness is the difficulty of defining a distance function in the modeling space (in a generalized sense) as compared to a distance function in the representation space. Note that the different approximations in different applications are all problem specific. For the exemplary application implemented in this thesis (boolean operations on linear half-spaces) our paradigm and notion of correc :ness can be extended to allow for curved half-spaces within tolerance zones. But we have been unable to devise a general ! way of allowing form variations in the definitions. 2.2 Previous Work A substantial part of the previous work differs from the approach described in this thesis in being problem specific and dealing with “simpler” problems such as line arrangements or convex hull computation. The works fall in the following three broad categories: exact arithmetic approach, error analysis approach and tolerance based 19 approach. Some new work has appeared which is trying to deal with the robustness problem in a more general setting. This is classified under General Approaches discussed in section 2.2.4. Subsection 2.2.5 covers related work which is not directly applicable but deserves mention. Another comprehensive survey on the state of the art in robust geometric modeling is [Fortune 1993], 2.2.1 Exact Arithmetic Approaches These approaches discretize the problem space (usually aleph-one cardinality) in some manner resulting in finitely representable / computable objects. 2.2.1.1 Green and Yao [Green &Yao 1986] • Problem: A robust algorithm to report the intersection points of segments in a plane. The endpoints of the segments are represented exactly as points of an integer grid. The intersection points are also coerced to lie on the integer grid. • Model of Correctness: Approximate correctness. The grid size a can be used to specify the tolerance for the vertex endpoints, the result is at most a / J l away from the intersection of the model. 2.2.1.2 Sugihara and Iri [Sugihara & Iri 1989] • Problem: Polygonal modeling system. Plans to implement a 3D polyhedral system where the class of objects is trihedral polyhedra. The parameter space of the lines is discretized. The rest of the compu tation is done in exact arithmetic. • Model of Correctness: Approximate correctness. 2.2.1.3 Milenkovic [Milenkovic 1988] • Problem: Polygonal modeling system. Operates on normalized polygons. The discretization is not performed on any fixed grid; the vertices are shifted to lie on the edges, and edges are cracked to maintain certain minimum-feature separation. The main problem with this system is that the normalized vertex could be O(n) far from the original vertex where n is number of vertices in the problem. • Model of Correc tness: Approximate correctness. 20 2.2.1.4 Karasick, Licijer and Nackman [Karasick et al. 1990] • Problem: Delauney triangulation using rational arithmetic. Interval arithmetic was used to compute values. When ambiguities were detected the values were computed exactly. This could be done because of the domain of the problem considered. • Model of Correctness: Exact. Exact calculations suffer from several deficiencies: the number of bits required go up very fast; some numbers cannot be represented exactly; and the discretization of the space is problem dependent, which leads to even more bits. For example, Sugihara estimates that 80 times the precision of input bits is required for handling quadric intersections. 2.2.2 Error Analysis Approaches Forward and backwar d error analyses, along with some symbolic reasoning is done to show that the output representation satisfies certain properties. A problem with most of these approaches is the complexity of the solution even for simple problems. In the case of moderately complex problems the approach has not been proven correct. 2.2.2.1 Purdue / Cornell Approaches 2.2.2.1.1 Karasick [Karasick 1989, Hoffmann et al. 1988] • Problem: Karasick did the first pioneering work in robust geometric modelling in which he wrote a polyhedral modeler. Although not hard to crash and not provably robust the program was much more reliable than other modelers. • Model of Correctness: Approximate correctness. 2.2.2.1.2 Hopcroft and Kahn [Hopcroft & Kahn 1992] • Problem: Intersect an approximately represented convex polyhedron with a half-space. Provably correct. According to the- authors the result does not seem directly generalizable to non-convex polyhedra. • Model of Correctness: Approximate correctness. 2.2.2.1.3 Stewart [Stewart 1991] • Problem: Implement a polyhedral modeler. Provided a concept of local robustness where robustness is maintained locally and it is hoped that the 21 algorithm will be robust globally. The idea is to look at features and see what relationships are possible between them locally. These observations are used to develop rules. When there is an ambiguity about any relation in a geometric computation these rules are used to filter out the unwanted relationships. Therefore, there will be a model locally for the chosen relation. A polyhedral modeler was written which is more formal than Karasick’s but not provably robust. However, Stewart has not been able to crash it. • Model of Correctness: Approximate correctness. 2.2.2.2 Milenkovic [Milenkovic 1988a, Milenkovic 1988b, Milenkovic 1989] • Problem: Construct an arrangement of lines with P bit coefficients. The algorithm correctly computes a line arrangement in the plane which is correct for “hidden lines” that lie within 2"p of the input lines. (A hidden line is an infinite precision, non-representable line for which the arrangement of lines is correct.) • Model of Correctness: Implicit correctness. 2.2.2.3 Li and Milenkovic [Li & Milenkovic 199 0] • Problem: Find strongly convex hulls. They provide a rounded arithmetic . algorithm to compute a 8 -convex 8 -hull1 that is provably correct. • Model of Correctness: Strong correctness. 2.2.2.4 Fortune [Fortune 1989] • Problem: (a) compute the convex hull of n points; and (b) maintain a triangulation of an implicit set of n points. • Model of Correctness: Approximate correctness. 2.2.3 Tolerance Based Approaches These approaches are based on interval arithmetic; [Moore 1962, Moore 1979, Kahan 1980]. Interval arithmetic uses a model of floating point computation where a bound can be placed on each floating point error. Taking these errors into account we L An e-convex 8-hull P of a set of points S' is a polygon such that no point of S is more than 8 outside of P and it is 8 strongly convex. 22 can give provably correc. answers for any computation. For example, using interval arithmetic you can know with certainty that 2 + 2 is guaranteed to lie within [3.999, 4.001]. 2.2.3.1 Salesin [Salesin et. al 1989, Salesin 1991] • Problem : This idea has been implemented to do point location in an approximately convex2 polygon, testing approximate and strong convexity of polygons and constructing e -convex 8 -hulls of a set of points. • The approach is based on a generalization of interval arithmetic called epsilon geometry. Instead of a predicate P(x) returning TRUE or FALSE, it returns a partition of the real line into (F,U,T). If a number r lies in F, it implies that you cannot m ove* by r to make P(x) TRUE. If r lies in T, then you can move * by r to make it TRUF and if r belongs to U, we cannot say anything about it. So, instead of a predicate returning TRUE or FALSE it tells you how much you have to move the data to make it true. In Figure 15, points P and Q are a distance d apart. Assuming exact arithmetic, the predicate P Co i n c i d e n t (P , Q) will return (F = x < d, U = d < x < d, T = d < x) where * is a variable and d is the distance. W hat this means is that we cannot make the predicate true by perturbing the points by a distance less than d and we can always make the predicate true by perturbing the points more than d. U is empty in this case because everything is done in exact arithmetic;. 6-2e Figure 15: Epsilon Geom etry 2' A polygon is approximately convex if the verdces can be perturbed slightly to make it convex. Conceptually this is the opposite of strongly convex. 23 2.2.3.2 2.2.3.3 Now assume the existence of uncertainty. The location of points P and Q is known with an uncertainty e, that is, if the represented values are Pc and Qc, then d(P, Pc) < e,and d(Q, Qc) < e. Now P C o in c id e n t will return (F = x < d - 2e, U = d - 2e < x < d + 2e, T — x d + %). This means that if the points are perturbed by a value in F, the points cannot be made coincident; if the points are perturbed by a value in T, then they can be guaranteed to be coincident; and if J the points are perturbed by a value in U, then we cannot say anything about the result. Once again the approach becomes too complex, too fast and a generalization does not seem possible. Model of Correctness: strong correctness and approximate correctness. Segal and Sequin [Segal & Sequin 1985,1988, Segal 1990] Problem : Polyhedral intersection algorithm, where if two features overlap they are merged and the tolerance region grown. The approach is made more rigorous in [Segal 1990] bat it works on a case by case basis. Although it has not been proven correct, it runs without crashing on Karasick’s cube / cube test. M odel of Correctness: Approximate correctness. £ = Figure 16: B ruderlin’s Tolerances B ruderlin an d Fang [Bruderlin 1990, B riiderlin 1991, Fang 1992] Problem : Boolean polyhedral modeler similar to Segal’s [Segal 1990] in that the features that coincide are merged. However, Bruderlin and Fang’s approach maintains another tolerance to detect when the decision to merge features is inconsistent. If an inconsistency is detected, the algorithm has to be rerun with a 24 different tolerance a n d /o r higher precision. Figure 16 shows the tolerances for a point P. Every point comes with a tolerance (x) associated with it and any instance within the tolerance zone is a valid instance. If v is the maximum error that we could possibly have because of our flouting point computations, then any point within e (= % - v) of P (see Figure 16) is guaranteed to be a valid instance and any point more than 8 (= x + v) away from P is guaranteed to be distinct. In Bruderlin’s approach any two points which are within e of each other are merged, any two points which are more than 8 apart are known to be distinct, and otherwise (distance more than e but less than 8), the situation is ambiguous, in which case the computation has to be rerun with a higher tolerance (increase x). Fang [Fang 1992] points out that the above approach generates too many ambiguities for 3D modeling problems and hence too many reruns with new sets of tolerances. Also, the separated objects are not dealt with correctly in the above approach. For example, if two points p j and p; are known to be coincident (pj = p 2), then there might exist a point p 3 such that pj is distinct from p 1 but not from P2- Fang [Fang 1992] fixes the second problem by modifying the updating rules in Bruderlin’s approach. In Fang’s approach eveiy object has three tolerances (instead of two) - - e, 5 and A. If x is the tolerance of an object, i.e. any instance within x is a valid instance, and if v is the error due to floating point computation, then e = x - v and 8 = x + v initially, and A = + Conceptually, if e regions of two objects intersect and there is a valid model within the intersection, then the two objects are considered coincident. If 8 regions of two objects do not intersect then they are definite!}/ apart. A is the minimum distance of an object to a neighbor that is not coincident with it. A is used to ensure that two coincident points have the same set of points that are distinct from them. This approach, called the analytical model method is very similar to Bruderlin’s and therefore, still suffers from the first problem mentioned in the previous paragraph. This i< addressed by introducing ti e implicit model method. In the analytical model when two objects are merged, the merged object has the same 25 analytical form a r , the original object. In the implicit model method, however, the merged object only satisfies the symbolic constraints (i.e. the incidences) and the tolerances of the simple objects (i.e. vertices, edges, faces etc.). An implicit model edge need not be straight as long as it maintains the correct ' incidences. A Boolean modeler has been implemented with this approach with primitives from planes and natural quadrics (cylinder, cone and sphere) and is written in Modula-2. The modeler seems to perform well with the Karasick tests. The following problems exist with Fang’s approach. • It is proved that the boundary representation of the computed representation is homeomorphic to the boundary representation of some model in Euclidean space. This implies that a torus and a torus with a knot are considered equivalent in Fang’s approach (the implicit correctness model). • The rounding error v is chosen beforehand for the entire system. Therefore, if two lines (out of, say, two thousand lines in the problem) are intersecting at a very small angle, the entire computation for all the lines will have to be repeated with a higher v. In Figure 17, the plot of error vs. angle is shown. (The error per floating point operation is kept a constant. The plot was computed using Mathematica assuming the standard model of floating point computation; one line was kept at a small angle a and the second line was slowly rotated to generate the plot) The error can become arbitrarily large as the angle between the lines become smaller, and choosing a single v for the entire system is not a good idea. • It is not clear how this approach could be generalized to preserve other properties than incidence. The only method suggested (to ensure that two lines intersect in a single point) is at best makeshift. Therefore, the main objection to his work are that the algorithm is not provably correct and that its generalization does not seem easy. 26 2.2.4 G eneral A pproaches Recently some work has emerged which has been aimed at solving the problem in a more general setting. Although some of the work could be classified under the headings of Exact Arithmetic (Section 2.2.1) and Tolerance-Based approaches (Section 2.2.3), it is discussed separately here because it is especially relevant to our research. 500 4 0 0 300 200 100 0 . 0 1 7 4 ' 0 . 0 1 7 5 0 . 0 1 7 6 Figure 17: Angle vs. E rro r 2.2.4.1 Yap [Yap 1993] Every operation is computed in exact arithmetic. Rational numbers and algebraic numbers are used where needed. The current implementation aims at problems which are not very deep computationally (so that the computation precision can be decided on before-hand). Yap is currently working on a more extensive library. For an explanation of how our work compares with his, see the discussion under Section 2.2.4.2. 2.2.4.2 Fortune and W yk [Fortune & W yk 1993, Fortune 1993] Fortune and Wyk designed and implemented an expression precompiler. The exact precision of the input is specified. The precompiler computes the required precision of the intermediate variable., required and allocates the required storage. This takes care of 27 the storage allocation and pointer management overhead often experienced in the exact arithmetic implementations. Another feature of this implementation is a floating point filter. Each computation is done in floating point and a static error bound is calculated by the precompiler. If a decision based on this computation is ambiguous, the exact value is computed. Note that this is essentially the same idea proposed by [Benouamer et al. 1993, Kao & Knott 1990], the difference being that in the earlier works the error bound was calculated dynamically. The work reviewed in this section is more general than the work mentioned in the earlier sections but focuses on inputs of low bit complexity and non-cascaded procedures. These are some possible problems with this approach • It will not work in cascaded computations or where the objects are computed dynamically. For example, if the user is picking the points as they are being constructed, this, approach will fail. • The static bounds are too pessimistic. Note that even dynamically computed interval bounds, which are much tighter than static bounds, are considered pessimistic at times. • Floating point filtering is possible only because a symbolic value is required as an answer. If the s olution needed a number, th m in this model of correctness the exact value woulu have had to be computed. In a more general setting Fortune [Fortune 1993] suggests a possible approach which is very similar to ours. A minimal requirement fo r reliable implementation is that each geometric primitive indicate whether or not its answer is reliable. I f only reliable answers are given, the output o f the algorithm is correct. I f some unreliable answer is given, then sometimes simple tricks are appropriate. For example, if the input is a set o f points, then a second execution o f the algorithm may be successful if the points are slightly perturbed. Alternatively, in a specific situation the unreliable answer might be disambiguated through further tests, or localized use o f higher precision. 28 2.2.5 Related Work : Edelsbrunner and Mucke [Edelsbrunner & Mu eke 1988] present a technique called “Simulation of Simplicity” to deal with degenerate cases in geometry. This was a modification of the conceptual perturbation method from the linear programming literature [Dantzig 1963], In a lot of geometric algorithms degenerate cases can be avoided if the data are moved consistently in some direction by e where e is infinitely small. Using this technique, for some problems an algorithm which outputs +1, 0 or -1 can transparently answer +1 or -1 by moving the data consistently. This eliminates a lot of degenerate cases and reduces code size and programming effort considerably. The technique is not directly applicable to geometric modeling because • sometimes degenerate cases are necessary in geometric modeling • the approach uses exact arithmetic extensively which is not desirable. A similar approach was presented by Yap [Yap 1988] and more recently by Canny [Emiris & Canny 1991]. 29 Chapter 3: Math Preliminaries This chapter presents the mathematical machinery required for the thesis. Section 3.1 discusses interval arithmetic and how it can be (has been) implemented. Verifying polynomial identities is the topic of Section 3.2 and Section 3.3, and gap theorems are presented in Section 3.4. Applications to the robustness problem are discussed in later chapters. 3.1 Interval Arithmetic Interval arithmetic was launched as a discipline in 1962 with the publication of the Ph.D. thesis of R. E. Moore [Moore 1962]. Since then'a lot of work has been done in the field. For the current exposition the references followed are [Moore 1979, Knuth 1969, Kahan 1980, Sun 1989, Ely 1991, Aberth & Schaefer 1992, Tisdale 1994]. In interval arithmetic, the concept of a “real” number is extended to that of an interval. An interval, represented as [a,b], can be thought of as a new kind of number. An ordinary real number a will be considered as a degenerate interval [a, a\. Armed with the basic concept we a re : io w ready to define interval arithmetic. In practice, we cannot carry out exact real arithmetic, but we can carry out reliable bounded arithmetic, i.e., given any real computation we can give reliable arbitrarily tight lower and upperbounds on the result. These bounds define an interval. Through computation on intervals, using rounded arithmetic, we can find a resulting interval containing the exact real results. This will be explained in a greater detail in the rest of this section. Some notation follows. 3.1.1 Interval Numbers An interval is defined as a closed bounded set of real numbers [a,b] = {x\a <x b} If X is an interval then its endpoints are denoted by X/ andXM . Thus X - [X/, X J . We will not distinguish between the real number a and the degenerate interval [a,a]. Two intervals X and / are equal if the corresponding sets are the same. The intersection of X and Y is defined to be empty if X n Y = 0 . Otherwise, the intersection is the interval [m ax (X,, Y,), min (X , 1 ) ] . II u u 30 The union of X and Y is an interval iff they have a nonempty intersection, otherwise it is undefined. The interval is [min (X(, Yfi, max (Xu, Y ) ] . The width of an interval X is defined to be w(X) -- Xu - X{. The absolute value of an interval X is defined to be |X| = max (|X^|, |^ u|) • The midpoint of an interval X is defined to be m(X) = (Xu + X;)/2. 3.1.2 Interval Arithmetic We can treat the intervals as numbers, and define on them operations that produce other intervals. If X and Y are the two interval operands and 0 is any binary operator mapping reals to reals then, t X O Y = { * 0 y |jte X a y e Y} There is one problem . In the above definition the resulting set need not be a single closed segment which is our definition of an interval. Therefore, the above definition is restricted to the cases where the resulting set is indeed an interval. Otherwise the operation is undefined. For example some of the mathematical operations are defined as follows. Addition. X + Y = ^ + Yt, X u + 7 u J which is the same as {x + y | (x e X) and (y e Y) } . Negation. - X = [-X a, -X j\ = {-x\x e X} Reciprocal. 1/X = {l/x\x e X} if 0 e X, else undefined. Multiplication. XY = [min XuYp X J U ) ,max ( X /,, X(Yu, X Y t, X Y J ] 1 3.1.3 Floating Point Representation Nearly all numeric computation on a computer is done in fixed precision floating point arithmetic. The number of bits s allocated for representing a number is divided into 1 - The complexity can be reduced to computing just one multiplication for the lower and the upper bound in most cases. 31 t and u, the number of bit; for representing the mantissa2 and exponent, respectively. A real number n is represented as m x 2e . The mantissa m, is a number between 0 and 1 such that its leftmost bit in a binary representation is 1. Therefore, t bits are devoted to representing m and u for the exponent e. As only a finite number of numbers can be represented by such a scheme, an arbitrary real number often has to be approximated by another. 3.1.4 Floating Point Arithmetic When a mathematical operation is applied to two floating point numbers x andy, the result z may not belong to the domain of floating point numbers and therefore needs to be approximated. Usually, it is rounded to the nearest floating point number z, such that z = z (1 + e) where -2 r + 1 < e < 2 *+ 1. For example, assuming a decimal representation (instead of the usual binary one) let t = 3 and u = 2 measured in decimal figures. Then if we wanted to represent z = 1 / 3 the number stored would be z = 0.333. Then 8 = z /z — 1 = — 1 x 1 0 ^ . Therefore, — 10 2 < e < 10~2 as shown. Most implementations also provide an option of taking the lower or upper bound of z, denoted as zm and zM respectively such that zm < z and zm = z (1 + 8) where — 2 ? + 1 < 8 < 0 . (zM is defined similarly.) Continuing with our example, zm = 0.333 and zM = 0.334 [Moore 1979]. Note that what we are saying is that the result of any arbitrary computation provided by the computer is accurate to the last but one digit. This is a strong statement but most floating point math libraries guarantee this. In the numerical analysis jargon, a computer function is said to be ideally realized if the above property is satisfied. 3.1.5 Rounded Interval Arithmetic According to the previous subsection we can replace the exact computation r on real numbers by a floating point computation that computes an interval [rm , rM] where rm and rM are the lower and upperbounds defined in Section 3.1.4, such that the width of the interval is small and is guaranteed to contain r. If we were trying to compute 2' Knuth [Knuth 1969] prefers the term fraction because mantissa in english means “a worthless addition” and has a different meaning in mathematics. 32 z = 1 /3 , then with the existing technology we can compute zm = 0.333 and zM = 0.334 and z is guaranteed to be greater than zrn and less than z^ . Therefore, now we can modify the definition of interval arithmetic operation given in section 3.1.2 on page 31 as follows. If the exact interval arithmetic operation results in an interval [a, b], then the rounded interval arithmetic results in an interval [am, b^]. It follows that the rounded interval always include the exact one, since [am’ • An example will illustrate this better. LetX = [1 ,2 ] and Y = [2,3], then according to the exact interval arithmetic, X + 1 /Y = [1 ,2 ] + [1 /3 , 1 /2 ] (= [4 /3 , 5 /2 ] and according to the rounded interval arithmetic with a decimal base and three digits for mantissa, X + l / Y = [ ( 4 / 3 ) m, (5 /2 ) M ] = [1.33,2.5] Note that the second interval [1.33, 2.5] is guaranteed to contain the exact result [1.3333..., 2.5] and it is representable in the floating point scheme. 3.1.6 Intervals approach exact results as precision increases In this section we wiil show that the rounded intervals can be made arbitrarily close to the real intervals by increasing the precision sufficiently. Let A and B be the two intervals that are computed by an arithmetic operation which is a tree of depth N. (Depth is measured from the root.) The leaves of the tree are the individual operands and the nodes, the operators +, ■ * , - and/. W e will prove that the above is true by induction on the height of the node. An interval computed in exact arithmetic, A when computed in rounded interval arithmetic is denoted by ~A . Basis: At this level the expression is basically the individual variable. So, the intervals are the values. A = [a, a]. At this level since no computation has been done, A = A. 33 Induction Hypothesis: If the resulting interval at a node of depth n is An, then w^Anj - w (An) — > 0 as precision approaches infinity. Proof: Let the operator at depth n-1 be ® and its arguments be A and 13. Then we have to show that The inclusion part can be shown as follows: • A ® B <^A ® B . A ® B = {c\a e A, b € B, c = a ® b} . Therefore, for every c in A ® B t there exist a a and b satisfying, a e A, b e B, c = a® b. Since A c l and B c S , for all c in A < 8> B , there exist an a and b satisfying, a e A,b e B, c ■ - - - a ® b . Therefore, A ® B czA ® B . This is true by definition of rounded interval arithmetic. Having proved the inclusion, the asymptotic result can be shown in two steps. 1) w[A ® E J - w ( A ® B ) —>0, and 2) w (A ® B) - w ( A ® B ) — » 0 as precision approaches infinity. The first step is straightforward because we know that at any operation the amount of error incurred is bounded by 2 m+ 1 where m is the number of bits in the mantissa of the representation. Therefore the quantity in 1) is bounded by a quantity proportional to 2 m + 1 which becomes zero as precision = m approaches infinity. The second step will be shown by going through the various individual operations (+, -, * and /) and showing that 2) holds. Let A = (/; , AM) then A. = (Am - &Am, A m + oAm ) . Similarly B = (Bm,BM) then B = (Bm - 5B m, BM + $Bm ) . The delta quantities are the errors introduced by and that w - w (A ® B) — > 0 as precision approaches infinity. 34 rounding operations. Since A and B are the results of computations at depth n, by the inductive hypothesis all the deltas approach zero as the precision approaches infinity. Then Case I: Addition: A + B = (Am- b A m+Bm- 8 B m,AM + 8AM +BM + &BM)-. Therefore w(A+B) —w{A+B) - 5A ^ + &BM + SA^ + 85 m and the quantity on the RHS approaches zero as precision approaches infinity. Case ID: Subtraction: This is the same as addition. Case II: Multiplication: Referring to the definition of multiplication from section 3.1.2, instead of looking at all the possible sixteen cases we go through only one case here. Other cases are similar. 3x2? = HAm-S A m) ■ (Am + 84m) • (Bm + SBm)) and therefore, W<3 x J ) - w ( A x B ) = SAm (Bm + 5Bm) +SBMAM + 5Am(Bm- S B m) + 8 The above expression is snonotonically decreasing with respect to the precision and therefore approaches zero as precision approaches infinity. Case IV: Division: This will be handled as an inverse followed by multiplication. Case V: Inverse: 1/3 = ( l / ( A M + 8Aw) , l / ( A m-8 A m) ) , Therefore, w ( l / 3 ) - w ( l / A ) = l / ( A m-S A m) - l / ( A M + SAM) - l / A m + l / A M The above is a monotonically decreasing function of the precision and approaches zero as precision goes to infinity. Q.EJD. 3.1.7 Implementation of Interval Arithmetic There are three different ways of implementing interval arithmetic and each one of them provides a different speed vs. interval size trade-off. 35 3.1.7.1 Standard Interval Arithmetic The standard way of implementing interval arithmetic was demonstrated in Section 3.1.5. Another possibility within this scheme is to have a tri-valued number, such that in addition to the interval bounds, the floating point value that results from applying the required floating point operations to the floating point input numbers is also kept. This was called the Triplex arithmetic in [Tisdale 1994]. The standard representation (with two numbers) provides the best bounds and the worst time complexity. This tends to be about four to five times slower than standard precision arithmetic. 3.1.7.2 Range Arithmetic Range Arithmetic [Aberth & Schaefer 1992] represents the number in the form of € (N±r) x 10 . Note thai the storage required for r is much smaller than the storage required for JV . Also, as r increases we can reduce the amount of storage used for N. The size of the interval bounds for this representation are worse than the previous case, but the time complexity is promising. This was, however,,, not chosen for our implementation because the Aberth implementation did not use any floating point hardware. W e felt that fast floating point hardware existed (in some cases much faster than integer arithmetic 1 . trdware) and should be used if possible. 3.1.7.3 Static Error Bounds Static error bounds [Fortune & Wyk 1993] are derived by analyzing the program variables statically. This is a very efficient approach because it has very little time overhead after the preprocessing is done, but the interval bounds are very pessimistic. In general, rounded interval arithmetic is alleged to have pessimistic error bounds. Therefore, for our implementation the representation with the best size bounds, viz. standard interval arithmetic (Section 3.1.7.1), was chosen. For more details about the implementation, see section 5.5. An interesting property of the implementation is that it does multiprecision calculations with zero overhead for double precision intervals. 36 Therefore, a double precision number does not have to pay the constant overhead of multiprecision representation. 3.2 Verifying Polynomial Identities The problem may be stated as follows: Given a black box F that evaluates a polynomial given an argumentX = (XqJCi - X ^ i) in n variables, verify that the polynomial calculated by F is identically equal to zero. 2 3 For example, if / = Ox y + 0yz + Oz , then using the procedure F, we have to verify whether each of tbs coefficients is actually equal to zero. This is a trivial problem if we have an explicit representation of the polynomial, because then all we have to do is to check that each of the coefficients is indeed zerc . But we don’t have /, what we have is a procedure F that evaluates/ at any given point. Both solutions given below use exact arithmetic. 3.2.1 N on-probabilistic A pproach In case we have a univariate polynom ial/of degree d, then using the fundamental theorem o f algebra, we know that if it is not identically equal to zero then it has at most d real roots. Therefore, use F to evaluate it at more than d points and if all these points ( > d ) are zero then we know that the polynomial is identically equal to zero. 3 2 For example, l e t / = (x — 1) (x — 2) ( x - 3 ) — x - 6x + 1 lx — 6 . If we are using F to evaluate it then them are three points that we can pick at which it will evaluate to zero, i.e., 1, 2 and 3. But there is no fourth point. Any other point in the real line will evaluate to non-zero showing that the polynomial is not identically equal to zero. For a multivariate case the following theorem due to Schwartz [Schwartz 1980] is useful. Theorem: Let k be a field and let S cr k be an arbitrary subset of k. Let p (X) be a polynomial of n variables x = X j.. .xn and total degree (maximum degree in any term) d with coefficients in k. Then the equation p (x) = 0 has almost d ■ |S| solutions in Sn , where |S| is the cardinality of S . For our case the number of variables could easily be in the order of ten whereas the d will typically be low. Therefore, choosing the smallest value of |S| = 2 , we get in the 37 order of five hundred computations. So, in order to test whether a single degree multivariate polynomial in 10 variables is identically equal to zero we need at least five hundred evaluations, which is not acceptable. 3.2.2 Probabilistic Approach An examination of the previous example / = (jc - 1) • (x - 2) • (x - 3) suggests that we have to be very unlucky to pick 1,2 and 3 as our three test points. Therefore, we would know that the polynomial is not identically equal to zero long before we test it at d points. It follows from the theorem presented in Section 3.2.1 that, given a multivariate polynom ial/in n variables and maximum degree d, if the polynomial is evaluated exactly at a random value r from a finite field S, then the probability that/evaluates to zero w h e n / is not identically null is d / \S \ . Here L S I is the cardinality of the field. Note that although the computation above is exact it can be done modulo a large prime, and therefore doe j not suffer the same growth problems as other exact computations. Also, this could be computed fast by storing 2b m odp (if b bits are being used for integer computation and p is the large prime). Now whenever an overflow is encountered this value is added to the result. If p is chosen large enough (2p > 2b) then the value we compute will be zero iff the value mod p is zero. For example, let us assume a decimal representation with b - 4 . Then we can store numbers from 0 to 9999 in the 4 representation. Let p = 9973 (a prime). Then 10 m od (p ) = 27, let b = 27. Say we want to compute the addition of 9900 and 1239 modulo p, then the addition signals an overflow (result =11139), because of the overflow the leftmost bit of the number is discarded (result = 1139). We add b to the number (result = 1166) which is the right result modulo p. To see that this is indeed true let N be the result after the computation (sans overflow). Let R be the largest representable number + 1. Then b = R - p , n = N um ber A fter O verflow - N - R .Then, n+b - The number returned by our procedure = N - R + R - p = N - p = N%p. Note that this is true only for large p. 38 3.2.2.1 Verifying Rational Polynomial Identities A rational polynomial p is of the form p a/pt, where p a and pjj are polynomials in n variables. If p is the polynomial being computed by the program, and p is evaluated randomly at a value r from a finite field S, then the probability of p evaluating to zero or being undefined when p or 1 [p is not identically equs . to zero is (da + d^)/\S\. 3.3 Verifying Polynomial Identities on a Variety This section defines varieties and introduces a result by Hong [Hong 1986] that applies to varieties defined by a triangular system of equations. Variety: Given m algebraic equations in n variables, then the set of n-dimensional points satisfying these equations simultaneously is a variety. In the following exposition the base field is the complex field. Let u denote the 2s variables « j, m 2. . .u2 s, and x v x 2.>-x2l be additional 21 variables introduced by the 21 equations, then the 2s-dinensional variety in the 2s+21 dimensional space (w, Jtp Jt2---*2p defined by the following equations F l {u,xx) = 0 F2(u,x-i,x2) - 0 F2l(u,xv x2,...x2[) = 0 (A) such that the following constraints exist on the F,- s. 1) The degree of X; in each Ft is at most 2. 2) The total degree for all variables (u and x) in each is at most 4. 3) In each F; the coefficients of monomials are all integers. The sum of absolute values of all these coefficients is called the weight of F ;. The weight of Fi is not more than a constant m. Then the problem is to verify that the integer-coefficient polynomial F, 39 F {u,xv x2'-.x2l) (B) is identically equal to zero on the variety defined in (A). Equivalently, the problem is to verify that the equation F = 0 is a consequence of (A), The result by Hong, provides us with a point P e, (m}, u2, . . u2s, x v x2-- -x2l) , such that F (P ) = 0 if and only if F = 0 on (A).3 This is a powerful result because it states that we can just evaluate the polynomial at a single point and if it is zero there it is guaranteed to be zero everywhere. The drawback is that the number of bits required to represent each coordinate of P is exponential in the number of independent variables (u) and the number of equations (/). In our case, though I would be typically very small, u could be large. 3.4 G ap Theorem Algebraic numbers are the real numbers defined by the roots of an integer- 2 coefficient polynomial. For example, x - 2 = 0 , defines the algebraic but irrational number . Gap theorems provide a technique for representing exactly and operating over these algebraic numbers in a computer. A gap theorem states that roots of a polynomial c tx n + cn lxn~ l + ... + Cq, such that certain constraints exist on the coefficients, are at least separated by some positive number g. Now if the rational approximations to two roots of the polynomial are closer than g, then we know that the two roots are identical. The same result can be extended to show that if the first m bits of a root are zero, then the root is guaranteed to be zero (here m is a function of g). This gives us a good handle on representing algebraic numbers by rational approximations (binary representations). 3‘ In fact this point can be used for a specific component of the variety needed for our application. 40 G ap Theorem : [Canny 1987] Let p{d ,c) be the class of polynomials of degree d and coefficient magnitude4 c. Let ...,fn(xj .jcn) e p (d ,c) be a collection of n polynomials in n variables which has only finitely many solutions when homogenized. If (otp a ,t, ... a n) is a solution to the system, then for any j either j \ ~ n i ^ a . = 0 or a > (3 dc) j • ' For example, if we have the polynomial Sx3- 1, then d = 3, n = 1 and c = 5 and by the application of the theorem we know that any non-zero root of the polynomial has to be greater than 1.09xl0~5 . 4- An integer polynomial p is said to have coefficient magnitude c if the magnitude of all the coefficients of p is less than or equal to c. C hapter 4: A pproach This chapter contains our approach to the inaccuracy problem as specified in Section 1.4. An overview of the paradigm of coincidence-is-not-by-coincidence is the subject of Section 4.1. Section 4.2 presents a justification for the paradigm from an engineer’s point of view. A mathematical description of the approach is provided in Section 4.3 and a proof of correctness follows in Section 4.4. Discussions on handling i of perturbations are the subject of Section 4.5, and Section 4.6 investigates the implications of making the notion of correctness even weaker. 4.1 Accidental Coincidences Call on God, but row away from the rocks — A n old Indian Proverb Catastrophic errors arise when wrong branching decisions are made because of inaccurate floating point data. This happens when there are geometric singularities or coincidences. The approach proposed in this thesis is best explained by the following real world example. m Figure 18: How to get from P to Q Assume that a sailor wants to sail from P to Q (Figure 18), but the straight line path I takes him through the reefs where his boat might run aground. His engineering solution to the above problem is to avoid the reefs altogether ;md take the alternate path m, 42 because he realizes that his goal is to reach Q and not to follow the shortest path from P to Q. He relaxed his constraints and thereby made his journey much more reliable. Such decisions are taken every day in engineering to make reliable decisions in spite of inaccurate or incorrect world models. Another reason one wants to avoid a path through the reefs is that a correct path will have to be executed with a very high precision and therefore might lead to high execution costs, which is undesirable. We decided to apply the same approach of sailing around obstacles to the problem of robustness in geometric algorithms. The analogy may be explained with the help of the following code. A(C) { Compute B; 1 / / A c t u a l l y B' g e t s com puted 2 / / s u c h t h a t |B' - B |< £ 3 i f (B * < 0) 4 , • • • 5 e l s e 6 7 } • In the code fragment above, the reefs correspond to the points where things might go wrong. This is the point where we branch based on the value of a computation and the value is near zero. In line 4 of the code, if |Z ? j < e then we cannot make the decision reliably. • The path from P to Q corresponds to the execution path in the decision tree of the program. In the code, let us assume that B ' < 0 and therefore the path taken is <1, 2, 3, 4, 5>. • For the sailor to be able to sail around, he has to know the presence of trouble spots, i.e., a bound on the location of the reefi:. The tool that is used here is interval arithmetic1. If values are computed using interval arithmetic, then at step 4 if the interval B contains 0 we know that we might hit the reefs. 43 • The relaxation of constraints is done by changing a single input value to a tolerance, i.e., instead of trying to compute the exact result for a single value, the result is computed for a value within the tolerance. A similar situation arises in manufacturing, where an artifact is acceptable (correct if it is manufactured within tolerance. Therefore in a branch statement ( s t e p 4 in the code), when an ambiguity is detected (presence of reefs), the value of the input is perturbed (possible because of the relaxed constraint) within the tolerance to steer clear of the ambiguity. This paradigm, designated as ‘Coincidence is not by Coincidence’ is an exact opposite of Nearness Implies Identity (Section 1.2.1), and reflects the view that accidental coincidences can be broken for computational convenience. The underlying idea behind the paradigm of Coincidence is not by Coincidence (henceforth referred to as CNC), is the observation that in reality things are never exact, objects are never coincident, things are always manufactured to a tolerance and Casio watches keep accurate time within a few seconds per year. Therefore, if a decision depends on the value off() at x, it is a valid decision for any value of x within the tolerance. The final result is correct if the coincidence was by chance and the decision could have gone either way. However, if the coincidence was by design, it needs to be preserved (Section 4.1.1). For example, in Figure 19, three half-spaces hj, h2 and % are intersected. The arrows attached to each of the half-spaces point in the direction of the matter. Each half space comes with its own tolerance (shown by a band along each of the lines). The problem is thatp (the intersection of hj and hj) cannot be reliably classified with respect to h2. If the program decides that p is below h2, then the intersection of h2 and h3 should lie to the left of the intersection of hj and h2. This may not be the case in a program written with floating point operations and could lead to an erroneous result. However, because of the tolerances we have a solution out of this quandary. We can choose either 1" Note that interval arithmetic is not central to, the paradigm. We could have used range arithmetic or static error bounds. 44 the solution in Figure 19 (b) or the solution in Figure519 (c). If any of the half-spaces from the tolerance zone is an equally valid half-space, we can either choose a configuration such that p lies reliably under h2 or we can choose a configuration where it lies reliably above h2. Note that the answer given is valid in both cases with the difference being the presence or absence of a microfeature. (a) (b) (c) Figure 19: Acciden tal Coincidence in the Intersection of three half-spaces 4.1.1 Forced Coincidences A singularity is almost always, invariably a clue Sherlock Holmes, The Valley o f Fear When coincidence is not by coincidence, one may not be allowed to perturb the entities at will. For example, in Figure 19, the designer might insist that the three half spaces pass through the same point. In such a case even though we have the tolerances, we are not allowed to perturb the entities with respect to each other. However, if we know that the coincidence is by design and know of the truth symbolically (not based on approximate numerical computation) the approach proposed will make correct decisions as we will see later. Such forced coincidences can arise in three cases. 4.1.1.1 Coincidence by Design A designer might want entities to be constrained so as to be coincident, tangential, parallel, etc. This was the case in the example given above. Our approach caters to geometric constraints that can be defined through expressions involving parameters 45 associated with the geometric entities. For example, if we want to express the fact that the lines I and m are parallel and S( is the slope of the line I, then the line m is forced to be parallel to / by setting s < — sl . Such coincidences are taken care of automatically because the intrinsic structure of the parametric constraints is the same as the intrinsic structure of a program. Note that this does not have the full power of an arbitrary constraint management system. The constraints here are unidirectional, and do not require the solution of simultaneous equations for their satisfaction. 4.1.1.2 Coincidence by Creation In Figure 20, two lines I and m are intersected to gi ve a point P. If P is classified with respect to I, it will classify on (a coincidence). A simplistic application of CNC might try to remove this coincidence (by perturbing / or m), which will not be possible because P will always belong to I. Such coincidences can be inferred by simple history management and the more obvious ones can actually be encoded in the program. Note that this coincidence is classified separately because it can usually be handled quite efficiently. Conceptually this is a part of Coincidence by Theorem, discussed next. 4.1.1.3 Coincidence By Theorem . On two arbitrary lines / and m, three random points a; are chosen on line / and another three bj, on m. Lines are drawn connecting ar to bj (Figure 21), Then it is true that the points C, D and are collinear (Pappus’ theorem). Such coincidences need to Figure 20: Coincidence by creation 46 be detected because other wise we might try to break the collinearity of C, D and E (by perturbing I, m or the placement of the six points) and not succeed. m Figure 21: Pappus’ theorem There are two subcases to be considered. Although the scope of the problem mentioned in Section 1.4, does not include computation with radicals, we mention it here to show how the complexity of the problem changes with the addition of the radical operator. 4.1.1.3.1 Computation without Radicals These coincidences can be detected as follows. For simplicity, the example in Figure 20 is chosen. Let the equation of / be A jX + A?y + A^ = 0 and that of m be B lx + B 2y + i?3 = 0 , then the solution for the point P is ^2^3~^3^2 * “ y = Now, when we are classifying point P with respect to line /, the coordinates of P are plugged into / and we get the equation (A2B 3 - A 3B 2) A 1+ (A3B l - A 1B 3)A 2 + ( a 1b 2 ~ a 2b 1) a 3 which is identically equal to zero. 47 More formally, the detection of forced coincidences can be reduced to the problem of verifying whether a polynomial is identically equal to zero. In fact the problem itself suggests a method for detecting these. If random values are plugged in the procedure to compute the polynomial, every iteration that produces a value of zero increases the probability of the polynomial being identically equal to zero, or the coincidence being by design (see Section 4.4). 4.1.1.3.2 Computation with Radicals Consider the example in Figure 22. Here, the point P defined by the intersection of s and r is classified against t. The definitions of s, r and t follow. r Q t Figure 22: Computation with Radicals s = x2 + y2 - 1 r = x - 0.5 t = s — 2r = x2 + y2 — 2x x = 0.5 Therefore, P = 48 In the program, we evaluate the classification of P with respect to t in interval i arithmetic (as in Section 4.1.1.3.1). The resultant interval contains zero indicating the possibility of a coincidence by theorem. However, the evaluation of y contains a radical and therefore we cannot use the result developed in the previous section to detect the presence of a theorem. We need a way of guaranteeing that P will always lie on t irrespective of the perturbation that we might give to the original parameters. This problem is related to the problem of detecting whether a polynomial is identically equal to zero on a variety. Consider the variety (s,r) defined by s and r (The two points P and Q). The polynomial t is defined to be identically equal to zero on the variety (s,r) if it is identically equal to zero modulo (s,r). t = x 2 + y 2 - 2 x = x 2 + y 2 - 1 (mod r) = 0 (mod s) Therefore, t is identically equal to zero on the variety (s,r). W e are currently looking at different approaches for dealing with the first problem it zero on P) and how it relates to the second problem (t zero on (P, £>)). Hong’s approach (section 3.3) can be used to solve both of the problems but is computationally intensive. Handling of Pi Also, he made a molten sea o ften cubits from brim to brim, round in compass and five cubits the height thereof; and a line o f thirty cubits did encompass it round about Old Testament (I Kings vii.23, and 2 Chronicles iv.2) Pi in) tends to appear in geometric computations ubiquitously. Even though handling pi is not a part of our original problem definition, a geometric correctness model will not be complete if it cannot handle pi. Pi presents a problem because it is not rational and therefore, cannot be operated upon. Also it is not algebraic and therefore, does not have a defining equation we can use. However, we can replace pi by a variable, 49 through algebraic extension [Hoffman 1989] and end out with a problem formulation that we can solve as before.2 4.2 Justification The world of engineering is fraught with tolerances. The domain that is closest to ours is that of manufacturing engineering. All machines manufacture parts with a certain inaccuracy. The less the inaccuracy one is willing to tolerate, the more is the cost of manufacturing the part. In this domain the solution is to design a part with tolerances. That is, in addition to nominal geometry, tolerances are specified on geometric entities and the manufactured part instances that lie within the tolerances are acceptable. Other domains of engineering also use tolerances, albeit not necessarily on geometric entities. We adopt a similar approach. The floating point hardware (like its machining counterpart) operates with a certain inaccuracy. If all the input operands are specified with tolerances, then a solution is deemed to be correct if it is correct for some instance within the tolerance. The basic idea is that the tolerances provide us with a way of being inaccurate in our operations. Also, for most applications in CAD / CAM, which is the primary domain that we are addressing, an answer which is correct for a value close to the exact input value is acceptable. Note that if it is not, then there are serious questions about the stability of the problem being solved, and of one’s ability to convert the solution into a physically feasible action. For instance, it may represent a design of a part that cannot be manufactured, a feature that cannot be inspected or, a path that cannot be followed using machinery of finite precision. This approach is also related to using r-sets to define solids and 9t-classes to define tolerances [Requicha 1977]. Both of these concepts eliminate spurious artifacts that can arise in the modeling space but which are not meaningful in the real world. 2- W hile setting random values for this variable it is quite possible that we might actually set the value of 7 t to 3 (as mentioned in the Old Testament) 50 Also note that this paradigm also takes care of the errors of type 1 mentioned in Section 1.2 on page 3. If a piece of data in the real worM needs specification at an infinite precision, it is replaced by a toleranced value. 4.3 Approach Our approach may be summarized as follows. The original, non-robust program is modified so as to operate in interval arithmetic. The result is what we call the interval- arithmetic module. Initially, arbitrary input instances are selected within the input tolerance zones, and the interval-arithmetic module is invoked on them. When we reach a branching decision that cannot be made reliably (as indicated by interval arithmetic), we perturb the input and try again. To do this we store the history of the computation that led to the decision in question. If several re-trials continue to lead to ambiguous decisions, it is likely that the coincidence is not accidental. The overall system controller switches modes and invokes the code which handles coincidences by theorem and by design. This is called the exact computation module. It uses the history of the computation to try to prove that the coincidence results from a theorem of geometry, by using exact computation modulo p, wherep is a prime number. A more formal and detailed description follows: We consider programs with bounded loops i.e. where the number of times a loop is executed is bounded, and model them as decision trees, where the branch statements correspond to internal nodes (degree 4) and the computation corresponds to the arcs (for the sake of simplicity we assume without loss of generality that there is only one algebraic operation per arc). When we refer to the computation at a node, we refer to the computation at the arc attaching the node to its parent. Therefore, there is no computation at the root node. Each execution of the program corresponds to a path from the root (starting point of the program) to a final leaf (ending point of the program). The input to the program is i, (corresponding to in Section 1.3.2.2), the output is o (the computed value, corresponding to 9tw) and d (the decisions taken by the computation in the tree, corresponding to Zp). For example in a program to find the intersection points of two circles in 2D, the input i is the center points and the radii of the two circles. That 51 is, if the centers of the circles are (Xj, yj) and (x2i, y2} and the radii are rj and r2 respectively, then i = (xj, yj, rlt x2, y2, r2). The program is I n t e r s e c t i o n P o i n t ( x l , y l , r l , x 2 , y 2 , r2) { d i s t S q r = ( x l - x 2 ) /v2 + ( y l - y 2 ) A2; i f ( d i s t S q r > ( r l + r 2 ) A2) / / # ( i n t e r s e c t i o n p t s ) = 0 e l s e i f ( d i s t S q r = = ( r l + r2) A2) { / / # ( i n t e r s e c t i o n p t s ) = 1 ( x o l , y o l ) — Computed v a l u e ; } e l s e { / / # ( i n t e r s e c t i o n p o i n t ) = 2 ( x o l , y o l ) = com puted v a l u e ; (xo2, yo2) = com pu ted v a l u e ; } 1 The tree corresponding to this program is shown in Figure 23. dist ? rl + r2 No Intersection Compute 1 pt Compute 2 Pts Figure 23: Decision Tree for Circle Intersection Program The output o are the coordinates of the intersection points, i.e., o = {xol, yol, x o l, yol) and d is the path taken by the program i.e. d = (c) where c e { - 1, 0, 1} . 52 In the exposition below, the given program (C in Section 1.4 on page 14) is denoted by To convert jP ,- into its robust counterpart, three modules are generated: the exact computation module (P%), the interval arithmetic module (Pt) and the controller (Cr). The controller in conjunction with the other two modules provides an implementation (C' in Section 1.4 on page 14) of the program Pi that is approximately numerically correct. The proof follows in Section 4.4. In the following exposition we will illustrate all the modules by an example of a s o r t routine which sorts two points on a line. Therefore, the original program P[ is as follows. s o r t (P [2 ]) { 1 i f ( P [ l ] - > t == P [ 0 ] - > t ) { 2 P [ 0 ] = NULL; 3 } 4 e l s e i f ( P [ l ] - > t < P [ 0 ] - > t ) { 5 s w i t c h P [0 ] & P [ l ] ; 6 } 7 } 8 The equivalent decision tree will be. P[0] P[l] Figure 24: Decision Tree for soning points (Pt) 53 f Figure 25: The Program as a decision tree Some terminology follows: Terminology Pi : Original Program P t : Interval Arithmetic Module Pg : Exact Computation Module Cr : The Controller i : (i t, i2, ... in), each q is a rational number (represented in floating pi int arithmetic, perhaps with arbitrary precision) n The number of inputs (the cardinality of i) to the original program I (I;, I2 ... In), the input to the trginsformed program. Ij is the tolerance for ij. IdentO : A variable that stores the names of nodes known to be identically equal to zero. NIdentO : A variable that stores the names of nodes known to be not identically equal to zero. 54 4.3.1 Interval Arithmetic Module P t The interval arithmetic module corresponding to the sorting routine is as follows s o r t T ( P i n t [ 2], t d e n t , P ath) { 1 P a t h = NULL; 2 i f ( ( P i n t [ 1 j - > t - P i n t [ 0 ] - > t ) . c o n t a i n s (0 )) 3 i f ( t h i s _ n o d e e I d e n t 0) { 4 P a t h < —P a t h + =; 5 P i n t [0} = NULL; 6 } 7 e l s e r e t u r n ( E x c e p t i o n ) ; 8 } 9 e l s e i f ( P i n t [ 0 ] - > t > P I n t [ l ] - > t ) { 10 P a t h c—P a t h + >; 11 s w i t c h P l n t [ 0 ] and P i n t [ I ] ; 12 } 13 e l s e { P a th e - P a t h + < ; } 14 r e t u r n S u c c e s s ; 15 } 16 Therefore, s o r t T executes the sort program in interval arithmetic except that (a) it knows when a certain node is identically equal to zero (this is done to avoid recomputation), and (b) it. maintains the various decisions taken in a P a t h variable. This latter is used in case we need to do theorem-detection. While doing theorem- detection we want to take the same path irrespective of the values in the internal nodes. Note that in this example P a th has maximum length of one. When an exception is raised, P a t h is null, i.e., has length zero. More generally, • All computations are performed in interval ar thmetic. The precision of the arguments is used in deciding the precision of the individual operations. 55 • The input to the module is (a) i, PInt[2] in our example, and (b) IdentO = { Wp n2-• -nm} , where each is an internal node in the decision tree corresponding to Pt. All branches corresponding to entries in IdentO are known to be identically equal to zero. Therefore, this is where the controller uses the knowledge that some of the statements are logically true. • The control executes the decision tree with i as the input and maintains the decisions taken at the individual nodes in the variable Path. Path is a list of symbols from <, > or =. • In a branch leading to node n if the computed interval contains zero, if n e IdentO, continue with the computation, taking the branch corresponding to n = 0 , else return Path and signal Exception. • If the control reaches a leaf node, the computed interval is returned and Success is signalled. 4.3.2 Exact Computa tion Module P E This module is obtained by transforming the original program P L in the following way. s o r t E ( P a t h ) { 1 P [ 0 ] = r a n d o m ( ) ; 2 P [1] = r a n d o m (>; 3 q = l e n g t h ( P a t h ) ; 4 i f (q == 1) { 5 i f ( P a t h [0] == ">" ) { s w i t c h P [0] and P [ 1 ] ; 6 r e t u r n - 1 ; } 7 e l s e i f (P a t h [ 0 ] = _» ___ ") { P [0 ] == NULL; 8 r e t u r n - 1 ; } 9 } / / S o rtE s h o u ld n ' t be c a l l e d w i t h q = l 10 e l s e { / / q == 0 11 r e t u r n ( P [ 0 ] - > t == P [ l ] - > t ) 12 } 13 56 Here s o r t E returns 1 when the points coincide and 0 when they don’t. It returns - 1 when the program reaches a branch that it shouldn’t reach. i . Therefore, s o r t E chooses input values randomly, executes the decision tree in exact arithmetic and takes each branch based on the P a t h variable. This is done to ensure that we repeat the computation i.e. evaluate the same expression, that caused the exception in so rtT . More generally, • All computations are performed in exact integer arithmetic modulo a prime p. • The input to the module is Path, which is a list of symbols from <, > or =. Let the length of Path be q. • The control generates n random values from the field of integer modulo prime p. The decision tree is traversed and all the relevant computations are performed. • The branch decision is made based on the corresponding Path element. If the branch at depth i is reached then the Pat hi is used to determine the branch to be taken. • When a branch node at depth q is reached, the. value of the argument of the if- statement is returned. 4.3.3 The Controller Cr For our example the controller will be C o n t r o l l e r ( P T o l [ 2 ] ) { / / T h e PTol i s P w i t h 1 / / t o l e r a n c e 2 C h oose P I n t [2] g PT ol [ 2 ] ; / / any 3 * IdentO = NULL; NIdentO = NULL: 4 A: i f ( s o r t T ( P i n t [ 2 ] , Id en tO , P a th ) == E x c e p t i o n ) { 5 m = l a s t I n d e x ( P a t h ) ; 6 i f (m e N I d e n t O ) { i n c r e a s e p r e c i s i o n / c h o o s e 7 a n o t h e r v a l u e f o r P i n t [ 2 ] ; 8 g o t o A; } 9 57 f o r ( i = 0; i< k ; i + +) { 10 v a l = s o r t E ( P a t h ) ; 11 i f ( v a l = = 0 ) {// c o i n c i d e n c e 12 / / c a n be b r o k e n 13 NIdentO = NIdentO", + m; 14 ' 5 g o t o A; } 15 } 16 IdentO' = I d e n t 0+ m; 17 g o t o A; 18 r e t u r n S u c c e s s ; 19 } 20 & Therefore, when the controller is called with the specified points with tolerances (PTol [ 2 ]), the program chooses some (arbitrary) instance from the tolerances, say it chooses two points along the line from the tolerance (line 3). The interval sorting program s o r t T is called with these points. If the program is able to sort the points using interval arithmetic, no exception is returned, and the result is correct. If the sort fails, then we check whether the failure was because of a theorem (lines 10 -1 6 ), in which case we maintain it (in IdentO). Otherwise, we remember the fact in N IdentO and continue the execution b> increasing precision. The program is guaranteed to terminate (see Proof) and the computed result is returned. The exact steps for a general routine follow. The controller is the /urogram that calls both the exact computation module and the interval arithmetic module. It’s input is I. 1) Set i e I (chosen randomly), IdentO < — and NIdentO < — < j> . IdentO is a set of nodes that are known to be identically equal to zero. NIdentO is a set of nodes that are known to be identically not equal to zero. 2) Call Pt (i, IdentO, Path). If Pt returns Excci tion, let the index of the last node in Path be m 58 If m e N IdentO , increase precision and / or choose another random value £br i (see Section 4.4). Goto 2. else, Gall PE {Path), k times (see Section 4.4 for how to choose k). If for any execution, the program returns a non-zero value we know the polynomial is not identically 0; set NIdentO *— NIdentO u m Increase precision and / or choose another value for i in /; goto 2. Else, if for all executions the program returns a zero-value; se! IdentO «— IdentO u m ; goto 2. 3) Return the value computed by Pt. 4.4 Proof In the exposition below, if P is a computation, then P f refers to its floating point evaluation at a certain point, P int refers to its interval arithmetic evaluation in floating point arithmetic andPg refers to its evaluation in exact integer arithmetic modulo p. The width of an interval (m, n) is n-m in exact arithmetic. We will prove that, with probability arbitrarily close to one the program will terminate wivh a result that is approximately numerically correct. 4.4.1 Proof of Termination Here we present a 4 tep proof that our algorithm terminates with a probability that can be made as close to one as desired. Step 1: Each node computation in the program tree is a polynomial evaluation in i. Here node computation refers to the computation corresponding to the previous arc. This will be proved by induction on the depth of the node in the tree. For the sake of simplicity division is ignored as an operator for the moment. Basis: At the root node the result is ij, where ij is.a variable in i. A variable is a polynomial. Induction Hypothesis: The two operands at the kth step are both polynomials P and Q in i. 59 Induction: The Ath step is either, P + Q, P - Q or P * Q. In either of the three cases, the output is a polynomial. QED. Step 2: A t a branch b, the evaluation Pint contains a zero. Then there is a possibility thatP is identically equal to zero. Check whether P is identically equal to zero. Choose random values for i from a finite field S (in our case, integer modulo p constitutes such a field). If the multivariate polynomial/corresponding to P is non-null, the probability that the polynomial will evaluate to zero at i is d /\S \ (Section 3.2). Here d is the degree of the polynomial concerned. The probability that the polynomial will ^9 k evaluate to zero for k evaluations is (d/\S\) . Now, if a polynomial is given to us ' J c (implicitly, through a “black-box” evaluator), then the probability is (d/\S\) that we will commit an error in assuming that it is identically equal to zero on k zero- evaluations, when it is identically not equal to zero. Therefore, we can find out with a high probability^ 1 - (d/\S\) th at/5 is identically equal to zero. Also, this probability can be increased arbitrarily (by increasing k). When the probability is sufficiently high, we assume that P is identically equal to zero and proceed from there. Figure 26: Decision T ree 60 A ( x , y) { A1 c o m p u t a t i o n s ; i f (x < y) / / e l s e e l s e } : < V ) / / A 3 A3 c o m p u t a t i o n s ; r e t u r n ; i f (x > y) { / / A2 A2 c o m p u t a t i o n s ; i f (x < - y ) { / /A5 AS c o m p u t a t i o n s ; r e t u r n ; } e l s e i f (x > - y ) { / / A4 c o m p u t a t i o n s A4 c o m p u t a t i o n s ; i f (y > 0) A6 c o m p u t a t i o n s ; r e t u r n ; e l s e i f (y<0) A7 c o m p u t a t i o n s ; r e t u r n ; e l s e i f (y == 0) A10 c o m p u t a t i o n s ; r e t u r n ; } e l s e i f (x == - y ) A9 c o m p u t a t i o n s ; r e t u r n ; i f (x === y) //A 8 A8 c o m p u t a t i o n s ; r e t u r n ; -0. 0 - 0.1 0.1 Figure 27: Input Space Partition Step 3: The program eventually terminates. 61 Consider the decision tree model of a program shown in Figure 26. A sample program similar to that decision tree model appears above, and the input space partition corresponding to the program appears in Figure 27. The rectangle in Figure 27 is the tolerance zone for the two input parameters in the example. This is the space from which the input parameters are chosen. Denote by Pj, P2 ... Pm the polynomial computations required to make a decision at nodes A }, A2 ... A m. In Figure 27, the input space (assume that n— 2) has been partitioned by the zeroes of the polynomials corresponding to P 1, P2 and P4 as explained below. (This is similar to a Binary Space Partition tree decomposition of the space [Naylor 1990].) The polynomial P 2 partitions the input space into two cells according to whether P2 is positive or negative evaluated at points in each cell. (In Figure 27 the relevant portion of the line P j = 0 is labelled Pj and similarly for other P^) All instances in the cell corresponding to j 3 take.the right branch (Pj > 0), and the others the left branch at Aj. The polynomial P 2 is used to partition the remaining left-branch cell into Cj5 and C p7 (the union of Cj6, Cj7 and P4). Cj67 is then partitioned into Cj6 and Cj7 by P4. The four 2-D cells correspond to the four leaves A3, A5, A 6 and A7 in the decision tree. Thus, for example, if i lies in cell Cj3 the program will branch left at A; and left again at A2, ending at Aj. (Note that the Pi s may not be straight lines in a different case.) Let us assume that initially we choose a pointp 4 (0.05, 0) lying on P4. Suppose that the interval arithmetic results are such that P 2int and Pzint allow non-ambiguous decisions to take the left branch at A2 and the right branch at A2. Suppose further that at A4, P4mt contains zero. Therefore, we choose another instance i and reevaluate. In this reevaluation there are five possibilities: 1) All ambiguities are resolved. 2) Another point lying on P4 was chosen, e.g. (0.075,0). A singularity point is reached at the same point, e.g. P4int contains zero. 3) A point lying on P2 was chosen, e.g. (0.05, -0.05). A singularly point is reached by an ancestor node of A4, e.g. P2int contains zero. . 4) In a more complex example, a singularity could be reached at a node that is not in the path to A4. For example, if A5 was a decision node, rather than a 62 leaf, selecting as an input point a root (as in root of a polynomial) of P$ would cause an ambiguity at A5. i 5) A point lying close but not on the Pi was chosen, and caused a P int to contain zero. Note that in the above, all of these cases can interact, e.g. an execution can reach an ambiguity at A4, then a re-execution can be ambiguous at A5, the next at A 2, and so on. We can keep choosing points satisfying these properties from the input space repeatedly. We will prove that infinite looping does not happen by doing a random sampling of the input space (i.e., the space of i) such that each sample corresponds to one execution of the decision tree. Let us assume for now that none of the polynom ials/^, P2 ■ ■ ■ Pm at the internal nodes are identically equal to zero. Since the polynomial roots define a zero measure sub-space in the input parameter space, there exists at least one finitely representable point j in I (which is / j x for this example) at which none of the polynomials at the internal nodes are zero and a leaf node is reached (using Section 3.2) when we use j as input, provided we use sufficient precision. Let Pq be the accuracy required to represent j. And let be the accuracy required such that w id th ^ (P Pk (J) 5 * 0 since higher precision will narrow the computed intervals (see section 3.1.6). Let p = m ax (Pq,p 1. .-pj), where I is the number of internal nodes. Therefore at accuracy p, for all evaluations at j, none of the internal nodes interval evaluations will contain zero. Now do a random sampling of the parameter space starting at accuracy p for K /2 steps, where K is the cardinality of I at accuracy p. Then double the accuracy and repeat for K /2 steps . Continuing ad infinitum, we get the following probability of K 1 hitting j. The probability of hitting j in the m!h step is — ------, since the size of the 2 2 mK sample space at the mth step is 2mK , and the number of points chosen are K / 2 . (We ir (J )) • 1 < K ir 00 • Note that this is always possible when K in t I I K I 63 assume a uniform probability density function.) Adding the probability of hitting j for all the steps we get the expression below. p » K ; 1 1 1 , , P r o b = _ { _ + _ + _ . . . } = i Step 4: Remove from the decision tree branches that cannot be taken. In Step 3 we ignored the fact that one or several could be identically equal to zero. W e showed in Step 2 that this situation can be detected with probability as close to one as desired. A polynomial Pk at the &th node serves to divide the corresponding cell into three parts, P*<0, Pk = 0 , and P k > 0. If a polynomial is identically equal to zero, only th e middle branch will be taken and cell division is unnecessary. In effect we can remove the two other branches (and subtrees) from the tree and coalesce two nodes. Step 2 shows that we can detect when a polynomial is identically zero, and Step 4 implies that one can remove from the tree the nodes associated with such polynomials. The resulting tree satisfies the assumptions of Step 3 and therefore the program terminates in the probabilistic sense described above. 4.4.2 Proof of Correctness In this section we will prove that the approach generates a solution that is approximately numerically correct. There are two steps to the proof. According to our definition of Approximate Numeric Correctness in Section 1.3.2.2, the computation C corresponds to P^ and C to Cr. For the program to be correct, we must show that f gAn ( 9\ \ m p • the domain and range of Cr are as expected, ie ., Cr : I 2 I — > I 2 I x Z , and • there exists a value within the tolerance for which the answer is correct. Step 1: Proving that the input and output o f the program Cr correspond to C'. 64 Thft innnf t nf th < » n m m m r is the tuple (Ij, ...In) where I^ are tolerances on ik. For the purposes of the proof we encode the symbolic part of the result as the path taken by the program Pt invoked by Cr This corresponds to our notion of symbolic output of a program. For example, a program to compute the classification (symbolic output) of a point P with respect to a half-space n is as follows C l a s s i f i c a t i o n ( P, n) { I f ( n (P ) < 0) / / t h e p o i n t p l u g g e d i n t o t h e e q u a t i o n The tree corresponding to the above program appears in Figure 28. The symbolic output of the program co rresponds to the path taken by the tree. Similarly for the program to compute the intersection points of two circles (appearing on page 52), the decision tree is show on Figure 23 on page 52. Here once again the symbolic output of the program corresponds to a path in the decision tree. / / o f t h e h a l f - s p a c e r e t u r n i n ; > , e l s e i f (n(P) :== 0) / / p o i n t on t h e h a l f s p a c e r e t u r n on; . e l s e r e t u r n o u t ; Figure '28: Decision tree for classification program 65 Therefore the symbolic output can be encoded as 7) = (D^...Dp) where Di e (-1 , 0, 1) . Therefore, D e Z ^ for some k. Note that the cardinality of the tuple D would be p if Pt followed the same path as the original program Pi on the input i. This has been assumed here but will be proved in Step 2 of this proof. The output O is basically the result of computing Pj but using intervals. Therefore, if o = (o , .. .o ) is some value in the range 9? o f Pj, then the output of P t is Si f Si^\m O = where o^e. Ok and 2 . Therefore, O e \^2 J . QED Step 2: Proving that there exists a value within the tolerance fo r which the answer is correct, i.e., if Cr (7) = (0,D) , then 3 (i, o, d) such that i e 7, o e O , d = D , and Pj(i) = (o,d) . 1) The existence of an i in 7 will be shown by construction. Consider the i in Step 3 of Section 4.4.1 for which the program terminates. This i belongs to 7 by definition. Now, we will (conceptually) execute Pi on i with exact arithmetic and let the results computed be o and d respectively. 2) d = D. The result/) returned by Cr is the path taken by Pt when executed on i. Therefore, execute Pt and Pj on i using interval arithmetic and exact arithmetic respectively. Assume that repeated runs with perturbed data were done, if necessary to resolve ambiguities;, and that identically null polynomials were detected and associated branches eliminated as in steps 2 & 4 of section 4.4.1. Let us focus on the final execution that successfully computed an output result. We will show that the kth entries dk of d and Dk of D are the same for all k. At an arbitrary;step k, there are two possibilities: (a) The polynomial that determines the branching decision is not identically equal to zero. In this case after any potential ambiguities have been resolved, Pt makes the same decision as Pi because the computation Ck on which the branch is made is done with interval arithmetic. That is, if the equivalent computation in Pj is ck then > 0 =» c^ > 0 and therefore, no error is made, (The case of Ck < 0 is similar.) Ck cannot contain 0, because all the ambiguities were removed, (b) The polynomial is identically equal to 66 zero. Then ck will evaluate to zero and Pi takes the middle branch. Since the program is able to detect that the polynomial is null, the controller ensures that the middle branch is taken. Therefore, for both cases dk is the same as Dk. 3) o g O . The computation of o by P t can be looked upon as a series z . = b. 0 c. where 0 e {+, x, -} and b., c. e i u {zq...z. j} with j ' J J J J J J varying from 0 to k. This basically means that a computation at step I depends on the initial values i or the values computed until then i.e. all the ZjS with j less than I. The output o is some set of z,s (za0,...zam). Equivalently, the computation of O by Pt can be looked upon a series Z j = B j ® C j where ® e {the rounded interval arithmetic equivalent operations of +, x, -} and B-, C. e i u { ZQ... Z ^. _ 1} with j varying from 0 to k. This basically means that a computation at step I depends on the initial values i or the values computed until then i.e. all the Zjs with j less than /. The output O is the set of from (Za0...Zam), with each Zai associated with azat- W e shall prove that Zj e Zj for all j from 0 to k by induction. Therefore, o e O because o = {z .} and O = {Z .} . L a iJ * ■ a iJ Basis: To prove that Z q e ZQ. zQ = bQ0 c-q and ZQ = B Q ® CQ. By definition B q = bQ and CQ = cQ. Also, 0 C q e ® Cq , because all primitive interval operations contain the actual result. Trivially, ® cQ = B ^ ® = ZQ and therefore, Zq e Z0 • i Induction Hypothesis: z; e Z ; for all I & 0 Induction Step: To prove that Zj + 1 e Z, + 1 . 1) zj + i = cj+ i by definition. 67 2) b .^ , < = B - i because either b. , = z, for / e 0 .../ and it is true by • '7 + 1 7 + 1 7 + 1 / J J hypothesis, or bJ +1& i, and it is true by definition. Similarly, cj+ l e Cj+ r 3) If we extend the 0 operator over the intervals (implemented exactly) then bj+10cJ + i S Bj + l 0 C j ^ . 4) The ® operator is the 0 operator extended to take care of rounding, then This is so, because in calculating the 7 + 1 7 + 1 7+1 7+1 0 interval operations the results are computed and rounded outwards. [Moore 1979]. 5) B; + i ® C/ + i = zj+ i' Therefore, z7+ l “ bj+ 1 ^ cj+ I s Bj + l < > Cj+ 1 - Bj+ l ® Cj+ 1 ~ Z7+ 1 and therefore, zi + 1 e Zi+ j . QED. 4.4.3 Handling of Division Operator So far we have ignored the division operator from our proofs. In this section we convert the problem with division operator to another with no division operator and we show that the results in the two cases are the same. Conceptually, every number in the original program is replaced with two numbers viz. the numerator and the denominator. Therefore, let us assume that the original program is P and the transformed program is Q. The step in P is zk = ak © bk , where a^, bk « s /' u {Zq...z^_1} and© e {+, x , —, -5-}. The branch statement in the kth step is based on , i.e., whether ck is greater or less than zero. For every variable in P, there are two variables in Q, conceptually the denominator and the numerator. Therefore, if a is a variable in P, t hen the equivalent variables in Q 68 are denoted by x (a) and y (a) . Next we show how each statement of the program P is transformed to create Q. In our formulation y (a) is never zero. (Note that this would imply an ill defined P{). 1) zk = ak + bk in P. Then this statement is replaced by two statements. x(zk) = x(ak)y(bk) +x(bk)y(ak) , yU k) = y(ak)y(bk) 2) zk = ak- bk in P. Then this statement is replaced by x(zk) = x(ak)y(bk) - x ( b k)y(ak) y(zk) = y(ak)y(bk) 3) zk = ak x bk in P. Then the replacemens will be x(zk) = x{ak)x{bk) y y(zk) = y(ak)y(bk) 4) zk = ak/ b k in P. Then the replacement will be x(zk) = x(ak)y(bk) y (zk) = y(ak)x(bk) if (y (zk) = 0) then signal exception Note that the third statement corresponds to a divide by zero scenario for which P is undefined. Therefore, for the sake of this proof we will assume that this case doesn’t happen. % 5) The branch at the kth step depends upon ck in P. The value of sgn (x (ck) ) sgn (y (ck) ) can be used to branch in Q, the sgn function computes the sign of a variable. 69 Next we prove that P and Q produce identical results. This is done by showing that x ( c k) at an arbitrary step k, ck - ^ ^ . This basically states that conceptually at each step in Q, we have the same numbers. The proof will be by induction on the number of steps in the computation zk. Basis: Input ij, Jt (/j) = ij and y (i x) = 1. Therefore, the result is true by x { i x) i1 construction, i.e. i 1 yCij) i X (.cj) Induction Hypothesis: , for all I less than k. Proof: For different operations (+, *, - and /) the following steps can be carried out. For the sake of brevity only the case of addition is shown. ck = ak +bk x(c k) = x (a k) y (b k) + x(bk) y (a k) y{c k) = y ( a k)y{b k) In the statement above the first is the step in P and the next two steps are its replacements in Q. Using these three equations we can easily show that _ x{ck) °k ~ y (c k) ' Similarly the same steps can be undertaken for subtraction, multiplication and division. Therefore, we know that for any arbitrary computation a in P, y{a) QED 70 4.5 Handling of Perturbations Our algorithm randomly perturbs the input when an Exception is raised. In the proof given in section 4.4 we showed that random selection of a new input from a uniform distribution guarantees that the algorithm will eventually terminate with a (approximately numerically) correct answer. Our definition of approximate numeric correctness implies that any perturbation within the tolerance is acceptable. However, certain perturbation schemes are more efficient than others. Issues of efficiency are discussed in section 4.5.1 below. We assume implicitly that the direction of perturbation is irrelevant to the definition of the problem. This is the case in most problems, e.g., convex hull or Voronoi diagram computation. This has intuitive validity, because in reality an object could have a random noise affect any part of the C-space of the object. If small random perturbations could change things drastically, knowledge would have been impossible because senses themselves have some random error associated with them. For certain applications one might want to reformulate the problem because the direction of perturbation is important. For example, in computing Boolean operations, if the union of two objects in contact is to be computed and if they are moved randomly away from each other, then the resulting object might have a topology very different from the correct result. Also, for worst case analysis, if the mass of an object is desired and the application is trying to ensure that the final mass is less than a certain quantity, then in case of a singularity the geometric entities should be moved in a direction so as to increase the mass of the object; this will guarantee that the final bound on the weight is a correct bound. Note that here we are not using the definition of correctness provided in Section 1.3.3 on page 13. However, by specifying a direction of motion we are providing some additional functionality. Perturbations with specified directions are discussed in section 4.5.2, 71 4.5.1 Arbitrary Direction of Motion f = * 2y-y The zeroes are at x = l,-l y = o Figure 29: Configuration Space of the branch function Assume that the function computed in the relevant if statement ((/< 0,f= 0 ,/> 0)) 2 is f(x, y) = x y - y and the tolerances on the parameters x and y are (0,2) and (-2,2). The zeroes of this function are at x = 1,-1 and y = 0. For the tolerance region at the program’s disposal only the values of x = 1 and y = 0 are of interest. The configuration space of the parameters x and y is shown in Figure 29. The zeroes are indicated by the lines x = 1 and y = 0. The region of interest is shown ^shaded. The arguments need to be perturbed if the computed function value (interval) contains zero. This happens if the function arguments are ‘near’ one of the zero lines. Therefore, we would like to move away from the zero lines (preferably along the steepest slopes). The solution is to use grad (V) of/because the gradient of a function/ is the direction in which the function has the maximum slope. This would be the locally best direction to move. S teepest descent is a well-kno wn technique in optimization and search problems. Suppose that in our example we evaluate/ at x = 1.1 andy = 0.9. Solving for the gradient we get V /= = I 2-'* * 2- 1) * ' (1.98, 0.21) Therefore, we find that although y has a bigger tolerance, it might be a better idea to perturb the value of x. 72 Some observations: • The gradient approach will work for all analytic functions. So, we don’t have to worry about the functions with radicals in them separately (cf. Section 4.1.1.3.2). Present Position Figure 30: C ornered • Sometimes, we may need to move against the gradient. For example, in Figure 30, the present position is in a corner (the curve represents the zero of the function) and we can come out of it only by moving towards the singularity. • Perturbing by the maximum amount possible would guarantee in most cases that the singularity is resolved. However we could argue that we should try to stick to the nominal position as far as possible. (Note that once again we are trying to satisfy something more than our definition of correctness.) • How does one guarantee that one branch desingularization does not cause another branch tc become singular. Although the proof guarantees that we will eventually come out of such loops, a heuristic could be to move only one argument of a function at a time, preferably the one that does not occur in other functions. • Some sort of history management will be necessary to figure out the various variables that affect our function. • The calculation of the gradient will have to be done the hard way because the analytic, closed form o f/is usually not available. We have a black-box function evaluator and one extra evaluation per parameter will be needed to estimate the gradient. 73 • In the implementation we have so far not found it necessary to actually calculate the gradient. 4.5.2 Specific Direction of Motion This may arise in specific problems such as Boolean operations and worst case analysis, as noted earlier.3 Here the problems are not formulated according to the format recommended in Section 1.4 on page 15, and it is hard: to specify them in a more abstract setting. Also, the perturbations of the primitives are very problem specific. Two possible approaches are: • The application author provides modules for perturbation of primitives. These modules have all the knowledge accessible to the program and they have to make the decision about direction and magnitude of perturbation. • The application author provides rules which are read in at the beginning of the execution and used to decide on the motion. In this case, we have to specify the language for writing of the rules and to justify the power of this language. In the current implementation we have chosen the first approach. 4.6 Solve by Generalization A coincidence is a theorem in a certain space. The idea here is that if our line could be any connected curve inside the specified tolerance1 zone then there is no reason why we should preserve its linearity. Therefore, we can also break such coincidences by generalizing the concepts. This idea will not be further explored in this thesis. It is i mentioned here for the sake of completeness (see Implicit Correctness, Section 2.1.4). For a union operation, for example, if two objects are in contact, we want to move them towards each other. For worst-case mass analysis objects should be grown to maximize their weight. 74 Chapter 5: Implementation Asymptotic analysis keeps a student’s head in the clouds attention to implementation details keeps his feet on the ground. Jon Luis Bentley This chapter contains the description of the implementation of a proof-of-concept Boolean polygonal modeler. The purpose of an implementation of a theory is described in Section 5.1. The problem chosen for implementation is described in Section 5.2, and testing its correctness in 5.3. Section 5.4 contains the various profile tests done for this implementation. An efficient implementation of an arbitrary precision floating point representation is presented in Section 5.5. 5.1 Rationale for an Implementation The main goals of our implementation are the following. 1) Reality Check. For the author of a theory, an implementation verifies whether the abstraction process was correct or not. This is one of the best ways of corroborating correctness. For a user (critic), an implementation allows him to substitute running demos for actually reading" and understanding entire proofs.1 2) Demonstration of Unproven Conjectures. If a conjecture is backed by proof, an implementation is more of a verification (see 1 above), but if it does not have a proof 'hen an implementation is required to justify the claims. The various unproven conjectures in our case are (a) the interval growth is not too pessimistic; (b) the time spent in the recompute routines is not too high; and (c) the cost of interval arithmetic instead of floating point arithmetic is not prohibitive. The benchmarks for the system should be related to real world examples to substantiate the claims. This is not to imply that running demos is superior to understanding, only that it is faster. An implementation can be used as a filter, so that the reader only reads / undei stands some of the proofs. 75 5.2 Polygonal Modeler Boundary evaluation was chosen for the sample implementation, because it is one of the most important problems in geometric modeling. The problem was restricted to the following domain. ^ Design and implement a CSG boundary evaluator on a plane with primitive half spaces defined by straight lines. In 2D, a solid can be defined by using a CSG tree. This is a rooted tree with all interior nodes either labelled with r>* (regularized intersection) or u* (regularized union). The leaves point to a half-space in 2D. For this exposition, a half-space is Ax+By+C © 0 where ©«={<,<}. Regularized union (SJ*) and intersection (n*) are 1 1 defined as follows A v * B = c l( i( A u B )) A n * B = cl (i(A r \ B}) where cl and i are the topological closure and interior operators respectively. The boundary evaluator has the following characteristics. 1) Each half-space is at its final position. In most solid modelers a half-space is moved several times before its final instantiation. However, this motion is not the cardinal problem that we are trying to solve. 2) Coincidence by design is provided by the following features in our implementation, (a) A line can be defined by two points; (b) a point can be defined by the intersection of two lines; and (c) two lines can share the same direction vector (thereby allowing parallelism). These features can be used to force coincidences. For example, (c) can be used to make lines parallel and this truth will be maintained symbolically (instead of numerically) by our program. 3) The output is a set of edges. A set membership classification approach [Tilove 1980] to boundary evaluation was adopted. The various steps followed were 1) Intersect the lines to generate the intersection points. 2) Sort the intersection points along the line. 76 3) Classify the midpoints of the line segments (generated by the intersection points) to find the boundary points. The program was implemented in C++ with about thirty five hundred lines of code. It comprises of twenty two classes out of which one is devoted to the interval class. This does not include the multiprecision code. 5.3 Correctness •! It is hard to show correctness for an implementation. The tests that are popular for boundary evaluation algorithms are 1) A cube B is obtained by rotating a cube A by a small angle a . A and B are then intersected. The literature reports that most modelers break down and either crash or return invalid results for oq < a < 0C 2 [Laidlaw et al. 1986]. For these modelers, for a < aj, the modeler detects all the edges to be coincident and returns a cube, and for a > a.2, it detects the two cubes to be different. Our program was run on the equivalent 2D problem. For all the values of a that the program was run on, it did not crash. At a value of a = otjn (the smallest representable value), the program detected that the two squares were different. Therefore, for our modeler 0C 2 = a m and for all values of a, the program runs without crashing. 2) Another test which is commonly done is f.he intersection of several cubes which have been rotated by a random angle, resulting in a polyhedron which is an approximation of a sphere. This was done by rotating several squares and intersecting them (see Figure 31). 5.4 Profile Testing Since worst case scenarios can be constructed for this program for which the program will spend exorbitant amount of time on recomputation and detecting coincidences by theorem, it is important to have good benchmarks for profiles. The benchmarks should be sufficiently complex and be typical of the types of parts encountered in real life by a CAD / CAM modeler. The following object cross-sections were approximated by polygons and used in the tests. 77 1) Connecting Rod: Figure 32. The artifact was available in the lab. 2) ANC 101: Figure 33. A popular benchmark used in feature-finding applications (Vandenbrande 90]. 3) Glue gun: Figure 34. A popular benchmark used in fixturing applications. 4) Tram m el: figure 35. An interesting mechanism. The endpoint of the rod traverses an elliptical path. mm jU rite [Locate [Quit Figure 31: Random intersection of half-spaces 78 jlir ite I ' p T s p l a J Figure 32: Connecting Rod i jurTte jp'isplay /. A ri oeceooeGtaeoBBccoaaQc w i Figure 33: ANC 101 (front view) 80 a •Write •Locate •Show [Quit aeeeaaae Figure 34: Glue G un U nite ■Locate Figure 35: Tram m el 5.4.1 Interval Growth Interval growth is a common complaint against the usage of interval arithmetic. We attempt to show here using the given examples that this does not occur in our domain. In this section first a measure of interval badness, RelativeWidth is defined and then it is shown that the different intervals do not become too big for the computations at hand. An interval is bad if it is too big i.e. if w idth = upperbound — low erbound istoo big. However, just looking at w idth for different intervals does not make sense because the error typically is committed in the lowest bit of the mantissa and different values in the exponent can scale this error quite differently. For example, if the representation involves a base 10,3 digit mantissa and 2 digit expone nt, then the interval corresponding to 3.333 x 10 might be (^3.33 x 10 , 3.34 x 1023 j . Therefore, the uncertainty is 0.01 x 1023. Note that this uncertainty would have been quite different if the actual number was 3.333 x 109? . The measure which was chosen instead was RelativeWidth. This is defined to be ( upperbound - low erbound) / a b so lu teV a lu e . This is a more appropriate measure because it takes care of the scaling of the error by the exponent. Double precision wan used in doing all the standard arithmetic operations. In the implementation of doubles on the SPARCstation, 53 bits are devoted to the representation of the mantissa (fraction part of a floating point representation). Therefore, the smallest non-zero value of the RelativeWidth possible is 1.11 x 10-16. Histograms are drawn on a logarithmic scale showing the distribution of RelativeWidth for different interval values in the result of the various profile tests (Figure 36, Figure 37, Figure 38 and Figure 39) and the coincidehce-by-theorem tests (Figure 40 and Figure 41). (These latter tests are described in Section 7.5.) From the histogram it is apparent that most of the intervals are actually quite close to the smallest value and the - 9 worst case is achieved when the RelativeWidth is in the order of 10 for the case of Glue gun. 83 400 300 200 100 rl5 ,-09 ,-08 ,-10 ,-16 Figure 36: ANC 101 interval distribution 140 120 100 80 m 60 40 i M 20 m i-1 5 i-13 ,-10 ,-09 ,-08 i-17 i-14 ,- 1 1 Figure 37: Connecting Rod interval distribution 84 100 80 G O 40 m 20 II m ,-09 ,-16 ,-15 ,-13 ,-10 ,-08 Figure 38: Glue Gun interval( distribution S O G O 40 20 10 i-17 10 ,-16 10 15 10 ,-14 10 ,-13 10 i-12 10 -11 10 -10 10 -09 10 ■ 0 8 Figure 39: Tram m el interval distribution 85 Figure 40: Pappus Theorem interval distribution m M i f Figure 41: Brianchon Theorem interval distribution Table 1: Interval Computation: Time Profile + - % / ieee_flags Time spent ANC 101 3.9% 22.8% 22.8% 3.7% 16.6% 69.8% Connecting Rod 1.5% 31.4% 27.2% 1.5% 20.3% 81.9% Glue Gun 3.3% 20.3% 20.3% 3.3% 13.5% 60.7% Trammel 4.4% 16.5% 19.3% 4.1% 16.1% 60.4% 5.4.2 Time Spent in Interval Computation The time spent in different modules of the program is shown in Table 1 above for the different benchmarks. The time spent in the routine ieee J'lags could not be attributed to the routines individually, and therefore has been kept separately. The total time in the interval routine ranges from 60.4% all the way up to 81.9%. This was expected because all these programs involve significant amount of arithmetic computation. We measured the speed of interval computation on our machines at around seven times slower than floating point for addition and fourteen times for multiplication. The tests were done on a SPARC II, with no optimization. Professor Hall of Wisconsin Madison [Hall 1994] measured the speed at a factor c.f three for addition. This was measured on a language C-XSC that supports interval operations. Professor Hall suspects the multiplication to be only a little slower. Acting conservatively we estimate the multiplication to be about twice as slow as addition. This is deduced because our multiplication is about twice as slow as addition. Looking at the mix of computation in Table 1, we deduce that the mix is 50% addition or subtraction and 50% multiplication. Therefore, we estimate that an efficient implementation of interval arithmetic will be slower than standard floating point by a factor of 4.5 for our applications. Note that this slows us down by about a factor of four to five. Also, fast hardware implementations of interval arithmetic are possible, resulting in a hardware which is about as fast as the floating point hardware. 87 Table 2: Time spent in recompute Part Name Percent of Time Spent in Recompute ANC 101 0% Connecting Rod 0% Glue Gun 6.1% Trammel 0% Pappus Theorem 62.5% Branchion’s Theorem 73.8% 5.4.3 Time Spent in the Recompute Routine The time spent in the recompute routines was plotted for different profile examples. Also, the examples with the various theorems were chosen to see how the recompute routines behaved in their presence. This recomputation time contains both the time to perturb the arguments and compute values in interval arithmetic and also doing the coincidence by theorem in exact modulo-p arithmetic. . It is clear that the time spent in the recompute routine is negligible for the real-part benchmarks at hand. The theorem examples however, take a significant amount of time (60-70%). This is acceptable because we do not expect theorems to occur frequendy in any of the real world examples. However, if they do occur the program has suitable mechanisms for detectin g them. 5.4.4 Speed of the Current Implementation Assuming that the interval arithmetic implementation is slower than floating point by a factor of 4.5, that the. program spends eighty per cent of its time in the interval arithmetic computation, and that the time spent in the recompute routines is twenty percent, we can compute the speed of the current implementation with respect to a standard floating-point implementation. The current implementation is then slower by a factor of 4.6. This is significant because off the shelf exact integer arithmetic libraries result in a slowdown factor of 40-140 [Fortune & Van Wyk 1993], off the shelf exact 88 rational libraries, of 10,000 [Karasick et al. 1990], and Yap [Yap 1993] is aiming for a factor of ten in his quest for the ideal exact arithmetic libraries. 5.5 Im plem entation of M ulti-Precision Interval Arithm etic Arbitrary precision, rounded implementation of interval arithmetic becomes inefficient because of two factors. First, two numbers have to be stored and operated upon for each interval and, secondly, while storing each number a fixed overhead cost has to be paid even if the number is storable in a double. Range arithmetic (Section 3.1.7.2) deals with the first problem by storing the number (N ) and a range (m) (N ± m ) instead of storing lower and upperbounds (((, u), I = N — m , u = N + m ). Therefore, only one number (N) needs to be stored at the higher precision. This saves some time at the expense of some accuracy. However, the second problem still exists. Our implementation of multiprecision arithmetic uses the provisions of the IEEE floating point standard covered in Section 5.5.1. The data part of the multiprecision class is discussed in Section 5.5.2 and the method part, in 5.5.3. The results of the implementation of this class is the subject of Section 5.5.4. 5.5.1 The IE E E Floating Point Standard Flags The IEEE floating point standard provides for a possibility to flag non-representable numbers (under the usual floating point representation). Therefore infinity and negative infinity can be flagged. Also allowed is a flag which when set indicates that the value stored in the variable is not-a-number (the flag is called NaN). One of the rationales behind these flags is to let program execution continue even when a division by zero is encountered. It is observed that most programs give a reasonable answer under such circumstances. A NaN is set when a 0/0 case is encountered or if it is explicitly set in the software. If 1 and 2 below have been executed to set up the proper environment, then 3 occurs when the hardware encounters a NaN. 1) The TRAP is enabled by calling the _ fp -n o d e _ routine provided by the standard. 89 2) The interrupt handler h () (user defined routine) is defined by calling s i g v e c . 3) An operation is executed with a number for which the NaN flag is set. The hardware encounters NaN and the signal gets generated. The signal is then passed by UNIX to the interrupted process through _ s i g t r amp which calls h (). Also, note that NaN is just a flag that is set to signal an interrupt, the rest of the storage in the variable can still be used to store other information. How these flags are used in our implementation is discussed below. 5.5.2 D ata Storage The storage of multiprecision numbers is discussed in this section. The first part of the section contains a description o f class My N um ber which is a multiprecision floating point class. The usage of this class to store numbers with no overhead is covered in the latter part of the section. A floating point number has two components: mantissa and exponent. For the number 3323.23 (= 0.332323 x 104 ), the mantissa is 0.332323 and the exponent is 4 for a decimal floating point representation. The rule for deciding whether a number is a valid mantissa is that it must be smaller than one and its left-most digit must be non-zero. Usually a fixed number of bits are allocated to the mantissa and the exponent part of a floating point number. For a 32 bit storage representation, typically 24 bits are stored for the mantissa and 8 bits for the exponent. This implies that only certain numbers can be stored as floating point numbers. To alleviate that problem, a class M yNumber was written to accommodate arbitrary precision mantissa and exponent. The basic representation of these arbitrary precision numbers is in the form of a linked list. If numbers 3, 4 and 5 are represented as the successive elements of the list, then the 2 1 0 number stored is 3 x 10 + 4 x 10 + 5 x 10 - 345. Therefore, any precision numbers / can be stored resulting ii an arbitrary precision floating point implementation. 90 Mantissa CJ3-HZZEr<ZZBr<Zn Exponent I H— i R— H 1 R— H ~R—H j i Figure 42: MyNumber (straightforward multiprecision implementation) Now that M y N u m b e r has been constructed, we need to be able to access it with a zero overhead for a double precision number, paying extra overhead only if higher precision is needed. Note that most numbers encountered in our implementation are doubles, therefore, a higher cost is acceptable for multiprecision. In Figure 43, the storage of two such numbers is illustrated. The first number 0.322 is stored as a double and incurs no extra overhead. In the case of the second number, the NaN flag is set and the actual location for the number is used to store a pointer to the M y N u m b e r class. Note that in both cases the actual type o f the number (C++ type) will be d o u b l e . 5.5.3 Usage There are three possible cases when using our representation for multiprecision numbers. Let us assume that the operator is addition (+) and the two operands are a and b . d o u b l e a , b , c ; NaN d o u b l e | 0.322 | F~] d o u b l e [ M y N u m b e r Figure 43: Storage of Multiprecision numbers 91 c = a + b ; Case I This is the case whet- both the operands are normal (i.e. represented with NaN set to false). The addition takes place as usual, i.e. the standard addition operator is invoked. This case incurs zero overhead. Case II Here the operand a is normal (e.g., 2.239) and b is represented with the NaN flag set and the variable pointing to the actual multiprecision number (e.g., 2023490.223...329) (see Figure 44). When the addition operator is being executed, it signals an interrupt. The interrupt is trapped and the user written interrupt handler h () is called. The routine h ( ) , looks at the operand with the NaN flag set, finds the equivalent MyNumber and calls MyNumber: :o p e r a n d + ( d o u b l e ) . The answer is stored in c . (Note that the answer is also most probably a multi-precision number and will require similar handling.) Case III Now both a and b are multiprecision numbers. This case is similar to Case II and will be treated similarly. 2.239 a Exception Raised Nab b Exception H andling routine... 2023490.223...329 returns Figure 44: H andling of NaNs 92 5.5.4 Im plem entation The implementation bf multi-precision number classes in C++ was done by Mr. Joseph Bautista under my supervision. At this writing, the implementation of M yNum ber is complete but interrupt handlers are still under study. The current implementation of the boolean modeler does not utilize the multiprecision implementation because there is a constant overhead involved in dealing with numbers which do not require multiprecision and in all our benchmarks we have not encountered any need for multiprecision. 93 Chapter 6: History Management I f a man does not learn from his experiences, he is forced to repeat them. A straightforward implementation of the approach described in Section 4.3 is not very efficient. In Figure 45, the program is shown as a decision tree. Let us assume that a coincidence is encountered in node A4\ the approach would imply that we perturb all the arguments and re-execute the program all over a^ain. This is wasteful and the problem can be alleviated by perturbing and recomputing only the relevant objects. This is done by keeping track of dependencies amongst all the objects. Figure 45: Program as decision tree (P ;gure 25 reproduced) Every object in our implementation keeps track of its parents and children. For example, a Point object o is created by the intersection of two line objects. Then the information about the lines (parents, the objects necessary to create o) is needed to recompute the Point whenever necessary. The information about the Point (child, the object created by the lines) is needed to maintain consistency when the lines are recomputed. Sections 6.1 - 6.3 discuss the design of the history management system. The basic design is covered in Section 6.1, and the complications introduced by subclasses and subobjects are discussed in Section 6.2 and Section 6.3 respectively. The usage of 94 history to take care of the. paradigm of coincidence is mot by coincidence is presented in Section 6.4. Section 6.5 addresses the issues arising out of looking at the objects at different levels of detail. The history management is object-oriented and implemented in C++. Also, we assume that every procedure computes an object and there are no side effects, i.e., all the procedures are constructors. This is not as restrictive an assumption as it sounds. A non-constructor procedure can be used by conceptually considering it a part of the constructor method in which it is invoked (Also see Section 6.5). 6.1 Basic Design 6.1.1 Parent Pointer Management Consider the constructor for Point P which takes as arguments two lines L l and L2. (Conceptually this is the intersection point of the two lines.) Here L l and L2 are the parents of P, and we need to maintain the pointers from Point (child) to Lines (parents). c l a s s P o i n t { 1 f r i e n d P o i n t L i n e : : I n t e r s e c t (L in e S) ; 2 p u b l i c : 3 P o i n t ( L i n e s , L i n e s ) ; 4 p r i v a t e : 5 L in e * A r g u m e n t1, *A rg u m en t2 ; 6 d o u b l e x , y ; 7 } 8 P o i n t : : P o i n t (L in e SL1, L in e SL2) { 9 A r g u m e n t 1 = - SL1; A rgum ent2 = SL2; 10 P = L l . I n t e r s e c t (L2); 11 } 12 The follow ing two steps are used to accomplish t ie task. 1) S in ce ,L l andL2 are the arguments to the constructor P o i n t : : P o i n t (), storage for pointers to L l and L2 is created in the data part o f P o i n t . These are the variables A r g u m e n t 1 and A r g u m e n t 2 (line 6). For an arbitrary object A, if T a x . . . Tn a n are the arguments (where a ± is the argument 95 and T± is the type) to the constructor, then the variables Tx *A1 . . . Tn * An are created in the data part of A. 2) In the constructor P o i n t : : P o i n t , L l and L2 are stored in A rgum ent 1 and A rgum ent 2 respectively (line 10). In case of an arbitrary object, a pointer to a^is stored in A^. This design implies that there could be only one constructor per object. This is because all the arguments to the constructor (Ll, L2) are being stored in the object (Argum ent 1, A rgum ent 2) and different constructors of the same object in C++ have different arguments by definition. More than one constructor can be accommodated by using dynamic binding or keeping a u n i o n 1 of the arguments along with a flag specifying which constructor was actually used. For example, assume that P o i n t has two constructors (discussed above and the second which takes three planes as input and returns the intersection P o i n t ) , then the following code will accomplish the task using dynamic binding. The constructor for a point defined by the intersection of iwo lines will become: c l a s s L in e { 1 p r i v a t e : 2 d o u b l e A, B, C; / / r e p r e s e n t s t h e l i n e o f t h e form 3 / / Ax + By + C ,= 0 4 } 5 c l a s s P o i n t { 6 p r i v a t e : 7 P o i n t ( d o u b l e , d o u b l e ) ; / / N o t e ; c o n s t r u c t o r p r i v a t e 8 d o u b l e x , v* / / C o o r d i n a t e v a l u e s , 9 } 10 11 c l a s s P o i n t l : p u b l i c P o i n t { 12 L Union as a C++ construct 96 13 p u b l i c : 14 P o i n t l ( L i n e s , L i n e s ) ; / / T h i s c o n s t r u c t o r p u b l i c l 5 p r i v a t e : 16 L in e *A r g u m e n t1 , * A rg u m en t2; / / B ack p o i n t e r s t o 17 / / c o n s t r u c t o r a r g u m e n t s 18 } 19 20 P o i n t l : r P o i n t l ( L in e s m, L i n e s n) { 21 d o u b l e D; 22 D = n .A * m.B - n .B * m.A; 23 i f (D == 0) ( e r r o r ( ) ; } 24 x = n .B * m,C - m.B * n .C ; 25 y = m.A * n .C - n .A * m.C; 26 x = x / D; y = y / D; 27 A r g u m e n t1 =. Sm; 28 A r g u m e n t2 = Sn; 29 } 30 Similarly P o i n 1 2 could be defined. The back pointers are stored as before because both P o i n t l and P o i n t 2 are valid objects with a single constructor each and we know how to deal with single constructors. In a calling program a P o i n t will be defined as follows L i n e 1, m; P o i n t * P = new P o i n t l (1, m) ; P can be used wherever P o i n t was being used leaving dynamic binding to worry about the actual type (constructor) of the object. This approach is cleaner than the second approach where we maintain the type information ourselves, the advantage in that case being that we can do it a little more efficiently. 97 6.1.2 C hildren Pointer M anagem ent Every object maintains a list of objects dependent on it. For example, in the example shown below a L in e object is used to create more than one P o i n t . c l a s s L in e { ■ 1 d o u b l e A, B, C; 2 v o i d ** l i s t ; 3 i n t num; 4 f r i e n d c l a s s P o i n t ; 5 1 6 7 P o i n t : : P o i n t (L in e &1, L in e &m) { 8 9 / / c o d e f o r p r o d u c i n g t h e P o i n t 10 11 1 . l i s t [1 .nuit!++] = t h i s ; / / S t o r e a p o i n t e r t o t h i s l 2 / / o b j e c t i n 1 13 m .l i s t [ m .n u m + + ] = t h i s ; / / Same f o r m 14 } 15 Therefore, the following steps are taken to maintain forward pointers. 1) Variables l i s t and num (lines 3 and 4) are added to the L in e object. These will provide the storage for the forward pointers. These pointers are not typed, because in contrast to parent pointers where we always know the type and number of the parent objects, in case of children pointers neither the number nor the type is known beforehand. However, this is not a limitation because, all objects know how to compute themselves, and each object knows its own type. However, this can be done more elegandy by inheriting the storage and procedures from the R o o t class as discussed below.2 98 2) In P o i n t : : P o i n t ( ) (line 7), the forward pointers from 1 and m are stored in 1 .1 i s t and m. 1 i s t (lines 11,13). For an arbitrary object A, if a 1 . . .a n are the arguments to the constructor, then the lines a i . l i s t [ a i .num++] = t h i s ; are added for each i from 1 to n. Root class The concept of root class is introduced here to do some o f the above manipulations more elegantly. c l a s s R o o t { p u b l i c : R o o t ( i n t B l k S i z e = R o o t D e f a u l t ) ; v i r t u a l SOut.come R e c o m p u t e d ; v o i d A d d D e p e n d e n c y (R o o t *) ; p r i v a t e : Rootp& o p e r a t o r [] ( i n t ) ; i n t i n d e x ; S t o r E le m * B lk ; i n t B l k S i z e ; f r i e n d S to r E le m ; }; 1) The root clas s contains the lines 12 and 14 from the code on page 98). Basically all the forward pointer data is stored in the root class.These are stored in B l k . So, a forward dependency from L in e to P o i n t is maintained by calling A d d D ep en d en cy (L) in the constructor for P o i n t (here L is one of the arguments to the constructor). 2) The issues of memory allocation were gleaned over in the previous discussion. These are taken care of in the Root class. 2‘ Actually Root class inheritance is the only way to ensure that the correct routines get called by dynamic binding. Otherwise, we will have to remember the type of the objects in l i s t as well. 99 3) All objects need to be able to recompute themselves as needed. At some point in our program we will have a list of untyped objects (say A [ m ]) that we will have to recompute. Therefore, A i i ] -> R e c o m p u te () will be called. This routine binds dynamically to the actual object where the corresponding routine will be called. Any object that needs to use the children pointer storage facility, inherits from the Root class. 6.2 Ripples of Inheritance This section considers the modifications to the pointer management schemes described in the previous section that are needed to deal with inheritance. Consider the class Cow which inherit:' from the class A n i m a l . There are two sets of genes g l and g2. Whereas g l is enough to identify a living being as an animal, g2 is necessary to classify the animal as a cow. 6.2.1 Parent Pointer Management The code is written as follows: c l a s s A n im a l: o u b l i c R o o t{ 1 p u b l i c : 2 A n im a l ( G e n e S e t l * g l ) { 3 From shirt-sleeves to shirt-sleeves in three generations. An old american saying g, = g i ; } 4 p r o t e c t e d : G e n e S e t l *g; 5 6 7 c l a s s Cow: p u b l i c A n i m a l { Cow ( G t n e S e t l * g l , G en eS e1 ^ * g 2 ) ; p r i v a t e : G eneSet.2 *g; 10 11 8 9 100 } 12 13 Cow::Cow ( G e n e S e t l * g l , GeneSer:2 * g 2 ) : A n im a l ( g l ) 14 { 15 g — g 2 ; 16 / / C r e a t e t h e o b j e c t C o w .. . 17 } 18 Since g2 is needed to create the cow, only g2 needs to be stored within Cow (see lines 10 and 16), g l has automatically been stored in the object A n im a l. 6.2.2 Children Pointer Management To allocate the storage for children pointer a class just needs to inherit from R o o t . However, note that the Cow object does not inherit from the R o o t class directly. This design decision was made to avoid wastage of storage by storing multiple copies of a R o o t object within a single object. Therefore, if a Cow object is used as an argument of a constructor for another object M ilk , then the code is written as follows: M i l k : :M ilk (Cow *C) { C - > l i s t [ C - > n u m + + ] = t h i s ; } Therefore, the forward pointers are maintained correctly but the storage is actually maintained in the A n im a l part of the object rather than the Cow part of it. 6.3 Pointer Maintenance for Subobjects Let us assume that the Cow class has one subobject l e g , i.e., a portion of the data associated with Cow is itself an object of type L e g . In this section we will discuss the effect of subobjects on history maintenance. In this section, for the sake of simplicity ignore the fact that Cow was inherited from A n i m a l . 101 6.3.1 Parent Pointer Management As far as the storage of data is concerned, a subobject is similar to a base class and will be dealt with similarly. c l a s s Cow: p u b l i c R o o t{ p u b l i c : Cow ( G e n e S e t l * g l , G e n e S e t2 * g 2 ) ; p r o t e c t e d : G e n e S e t2 *g; l e g 1; } C o w ::Cow ( G e n e S e t l * g l , G e n e S e t2 *g2) : l ( g l ) { g = g 2 ; . . . / / c o m p u t e t h e o t h e r r e l e v a n t i n f o r m a t i o n . . . } c l a s s l e g : p u b l i c R o o t{ p r o t e c t e d : G e n e S e t l *g; } Since the l e g already contains the parent pointer information to the G e n e S e t g 1 , it was not necessary to store it again in Cow. 6.3.2 Children Pointer Management The meek shall inherit the earth. Although l e g will have a R o o t object contained in it, it is not directly accessible from Cow. Therefore, Cow will inherit from R o o t as- shown. 6.4 History Graph Traversal In the earlier part of this chapter we described the; mechanism to construct and maintain the history graph. In this section we concentrate on the traversal of the history 102 graph to perturb and recompute the objects. The standard constructor is modified to provide the following three procedures, Constructor (Section 6.4.1), MyConstructor (Section 6.4.2) and Recompute (Section 6.4.6). In addition, the procedures : Perturb (Section 6.4.3), MakeConsistent (Section 6.4.4) and RecomputeDeps (Section 6.4.5), are required to implement the perturbation / recomputation paradigm. 6.4.1 Constructor ( O b j e c t : : O b j e c t ()) The constructor performs the following functions; 1) It does the parent and children pointer assignment as shown earlier. 2) It allocates all the storage. 3) If an if-statement is present in the constructor which might cause a singularity, all the non storage-allocation statements in the constructor are moved to another procedure, M y C o n s t r u c t o r (). The example below shows how P o i n t (defined by the intersection of two lines) is constructed. Since the calculation of the intersection point involves the possibility of a singularity, the routine has been separated into a Constructor and a M y C o n s t r u c t o r (). P o i n t : : P o i n t (L in e *L1, L in e *L2) { 1 1 = L l ; 2 m = L2; / / P a r e n t p o i n t e r a s s i g n m e n t 3 L l - > l i s t [ L 3 ->num++] = t h i s ; / / C h i l d r e n p o i n t e r 4 L 2 - > l i s t [ L 2 - > n u m + + ] = t h i s ; / / a s s i g n m e n t 5 / / N o t e t h a t M y C o n s t r u c t o r () w i l l b e u s e d o n l y 6 / / i f an i f - s t a t e m e n t i s p r e s e n t i n t h e c o n s t r u c t o r ? M y C o n s t r u c t o r ( ) ; 8 } ! 9 6.4.2 MyConstructor ( v o i d Ob j e c t : :M y C o r.str u c to r ()) The basic assumption in our design is that an object knows how to construct itself. The storage allocation part of the construction has been kept in the actual constructor 103 and the object is constructed in MyConstructor. This is done to ensure that we don’t allocate the storage for the object repeatedly. The example below is the continuation of the previous example, where we wrote the constructor for P o i n t . v o i d P o i n t : :M y C o n s tr u c to r 0 { 1 s t a t i c i n t E i ; / / w i l l b e s e t e q u a l t o 2 / / z e r o by d e f a u l t 3 ’l i f (1 p a r a l l e l t o m) / / e x c e p t i o n 4 i f (E i == M axCoinc) //M a x C o in c = 3, e . g . 5 / / c h a n c e s h i g h t h a t c o i n c i d e n c e i s b y t h e o r e m 6 . . . c o i n c i d e n c e by t h e o r e m h a n d l i n g 7 e l s e {E i + + ; 8 l - > P e r t u r b ( ) ; 9 m ~ > P e r tu r b ( ) ; 10 w h i l e ( ( l - > M a k e C o n s i s t e n t ( ) = = E x c e p t io n ) 11 &&( m - > M a k e C o n s i s t e n t = = E x c e p t i o n ) ) { 12 l- > P e : c t u r b () ; 13 m - > P e r tu r b ( ) ; 14 } ' 15 M y C o n s t r u c t o r ( ) ; 16 r e t u r n ; 17 } 18 . . . c o m p u t e t h e i n t e r s e c t i o n o f l i n e s e t c . 19 } 20 When a singularity is discovered (line 4 above), we need to perturb the arguments, and recompute the objects all over again. This is done by lines 9-10 that perturb the arguments and lines 11-12 that make sure that all obj ects are recomputed consistently with no exception. If all the objects are recomputed correctly, i.e. 104 M akeC onsi s t e n t {) routines do not return E x c e p t i o n s , we can compute this object (the current P o i n t object) all over again. This can conveniently be done by calling M y C o n s t r u c t o r () recursively. If the object gets computed correctly, the program returns (line 17). In the approach presented in Section 4.3.3, when a singularity is detected (possible coincidence-by-theorem) first we check whether the eoincidence-is-by-theorem or not. This is done by executing the exact arithmetic machinery. In our implementation, however, we first try to remove the singularity by perturbing the arguments MaxCo i n c times. This happens when we recursively invoke M y C o n s t r u c t o r () and each time the static variable E i gets incremented. If E i equals M axC oinc, the exact arithmetic machinery is rolled out — this is the subject of the nex t chapter. If the exact computation determines that the coincidence is by theorem, then i| is preserved, else, we know that the coincidence can be broken (Section 4.3.3). The mechanism to break the singularity (by increasing precision) when it is known that the singularity is by coincidence, has not been shown for simplicity. 6.4.3 Perturb ( v o i d Ob j e c t : : P e r t u r b () ) P e r t u r b ( ) performs the following functions. , 1) Invalidate the object. When an object is perturbed it is not yet recomputed, therefore, it .is left in an inconsistent state. 2) Perturb all the arguments which the current object depends upon. 3) If the argument is a primitive value, choose a new value from the given tolerance. The example below, illustrates the P e r t u r b () routine for our example of % P o i n t . v o i d P o i n t : : P e r t u r b () { v a l i d = F a l s e ; l - > P e r t u r b ( ) ; m - > P e r t u r b () ; } 105 6.4.4 M akeConsistent (S Out come Ob j e c t : :M akeC onsi s t e n t ()) Consider the following scenario. P e r t u r b () has modified all the relevant arguments and their parents all the way up to the primitive elements (where we had the tolerances). New values have been chosen. Now we need to recompute all those objects again. This is done by calling M a k e C o n s i s t e n t () . For our example, L i n e : : L i n e { 1 4 p u b l i c : 2 P o i n t * P t ; / / d e f i n e d by a p o i n t P t on t h e l i n e 3 P o i n t * d i r ; / / and a d i r e c t i o n d i r 4 } 5 SOutcom e = { S u c c e s s , E x c e p t i o n } ; 6 SOutcom e L i n e : : M a k e C o n s i s t e n t () { 7 P t - > M a k e C o n s i s t e n t ( ) ; 1 8 d i r - > M a k e C o n s i s t e n t ( ) ; , 9 i f ( R e c o m p u t e D e p s ( ) = = E x c e p t io n ) r e t u r n E x c e p t i o n ; 10 v a l i d = T ru e; 11 r e t u r n S u c c e s s ; 12 } ' 13 This procedure recursively asks each of its ancestors to make themselves consistent (line 8-9). At line 10, we know that all the ancestors of the object are consistent. Thereafter, in line 10, all the children of * t h i s are recomputed (including * t h i s ) . Having computed the object correctly, the object is validated and Success is signalled. 6.4.5 RecomputeDeps (SOutcome R o o t : : R eco m p u teD ep s ()) This routine is a part of the R o o t class and gets called automatically by dynamic binding. This calls the R ecom p u te () routine for * t h i s and all its children. 6.4.6 Recom pute (SOutcom e O b j e c t : : R ecom p u te ()) This routine is very similar to M y C o n stru c to r? routine in that it recomputes the object without allocating any of the storage. The difference is that in case a singularity 106 is encountered in M y C o n str u c to r , the arguments to the object are perturbed, whereas here if a singularity is encountered, an exception is signalled and the program returns to some Object*: :M y C o n str u c to r () whose arguments will be perturbed. SOutcom e = { S u c c e s s , E x c e p t i o n } ; SOutcom e P o in t,: : R ecom p u te () { . . . some c o m p u t a t io n i f (1 p a r a l l e l t o m) r e t u r n E x c e p t i o n ; co m p u te t h e i n t e r s e c t i o n ; r e t u r n S u c c e s s ; / K } 6.5 Levels of Granularity When a program is considered at the instruction level (the way we looked at it in Section 4.3), the level of granularity is straightforward and it is clear what is an object and what pointers need to be maintained to store history. However, in an object-oriented design an object is usually considered as composed of several objects. These can arise in the following manner The object either, • is composed of subobjects, A n im a l composed of Head and Trunk. • inherits from another object, Cow inheriting from A n i m a l . Therefore, it is not clear whether the dependency information must be kept at the level of Cow, A n im a l or Trunk (or all of the above). A related issue is the overhead incurred in storing the dependency information. The two factors affecting the design decision are: • The amount of overhead that one has to incur is a constant, irrespective of the size of the object, i.e., if the dependency information is stored in Cow it will cost k units and in Head it will still cost k units. However, the relative cost per unit size for the Head is much higher. Therefore, considering just the cost of the 107 overhead, it is desirable to maintain the pointers at the coarsest level, in our example at the level of Cow. • Let us assume that we maintain the dependency information at the level of A n im a l, and that A n i m a l : : A n im a l () is the relevant constructor. Now, when the object Animal needs to be recomputed, all the previously computed information about the Head and Trunk will be thrown away and all the information will be recomputed from scratch even if these objects did not need to be recomputed This is the problem mentioned in the beginning of this chapter on page 94. Therefore, the finer the level of detail for our dependency maintenance, the m ore information we will be able to reuse. The optimal level of detail at which we maintain the dependency information is best determined experimentally and in our implementation (Boolean operations on planar half-spaces in 2D) it was chosen to be the level of P o i n t s . S o l i d CSG S o l i d BRep ■ £ sg Figure 46: Object Modification vs. Object Construction Another case of relevance is when an object conceptually modifies itself. For example in Figure 46, a solid represented by CSG is being modified to get another Solid object that also has the BRep information. This will be a common scenario in geometric modeling where instead of computing objects from scratch frequently the procedures will just modify the objects. This can once again be considered as a case where our paradigm applies easily if we consider the objects at the lower level of CSG and BRep Sr 108 and maintain the dependency information at that level without maintaining it at the level of Solids. t 109 Chapter 7: Coincidence by Theorem: An Implementation Probabilia... sapientis vita regeretur. (Probabilities direct th i conduct o f the wise man. ) Cicero, De Natura Deorum This chapter contains the description of the ir .plementation of the eoincidence- by-theorem part of the implementation design. This chapter takes the theorem detection part of the approach described in Section 4.3 and transforms it into C++. More specifically this is the piece of code missing from line 7 of the code on page 104. Implementation of class I n t e g e r P is the subject of Section 7.1. The design of exact classes corresponding to the approximate classes is covered in Section 7.2, the new data elements that need to be added are discussed in Section 7.3, and the new routines necessary to test coincidence-by-theorem follow in Section 7.4. 7.1 Exact arithmetic modulo p An exact arithmetic class I n t e g e r P was implemented. The operations +, -, * and / modulo p where p is a large prime number were provided. Note that this is an efficient class because all the objects are of constant size and ■ onstitute a field. Unlike rational representations, only a single number must be stored, 7.2 An Exact Class for an Approximate Class An eye fo r an eye Exodus, Ch. 21, v. 23 For each class existing in the original program (all the classes in Chapter 6) a new class is provided which will correspond to an exact version of the object. This is how the exact decision tree is implemented (Section 4.3.2). In the example given 1 >elo w a new class P o i n t E Is introduced for class P o i n t that duplicates all the operations executed in interval arithmetic in exact arithmetic (using I n t e g e r P ) . Therefore, a parallel hierarchy is constructed that computes everything exactly in modulo p arithmetic. We don’t have the ambiguity problem here because all operations are computed in exact arithmetic. 110 c l a s s P o i n t { p u b l i c : P o i n t s t r a n s l a t e ( d o u b l e s , d o u b l e ) ; p r o t e c t e d : d o u b l e xO, x l ; } c l a s s P o i n t E { p u b l i c : P o in tE S t r a n s l a t e ( I n t e g e r P S, I n t e g e r P s ) ; p r o t e c t e d : I n t e g e r P xO, x l ; } 7.3 Addition of data elements in an Object For each object 0 that has a primitive data element D (a d o u b le ), an additional data element DE (an I n t e g e r P ) is added. This data element is in exact arithmetic and stores the randomized version of D. For example, the class P o i n t is changed as follows: c l a s s P o i n t { 1 p u b l i c : 2 d o u b l e x , y; 3 I n t e g e r P xE, yE; 4 } 5 Line 4 was added to accommodate the storage of randomized, exact values for x and y. Separate storage is important because we will eventually need the original values of x and y . I l l 7.4 Coincidence by Theorem “ Mr. Bond, they have a saying in Chicago: ‘Once is happenstance. Twice is coincidence. The third time it’s enemy action.’ Miami, Sandwich and now Geneva. I propose to wring the truth out o f you.” Ian Fleming, Goldfinger In Chapter 4, we described the handling of coincidence-by-theorem at an abstract level. In the last two sections (Section 7.1 and 7.2) we introduced some of the ingredients necessary to implement the detection of theorems. In this section we show how the different pieces o f the mosaic fit together. The following routines are needed: 7.4.1 Randomize ( v o id O b j e c t : : Randomise ()) This routine is very similar to the P e r t u r b () routine (Section 6.4.3) except this takes place in the exact arithmetic domain. One Randomi z e () routine is implemented for each Ob j e c t in the system. In the example below P o i n t 1 (the point defined by the intersection of two lines) is randomized. Since both the parent pointers (pointing to the half-spaces H [ 0 ] and H [ 1 ] defining the point) are non-primitive objects, the Randomize <) routine is called recursively. v o i d P o i n t 1 : : R andom ize() { 1 H [ 0 ] -> R and om ize(); 2 H [ 1 ] -> R and om ize(); 3 } 4 However, in the case of P o i n t below, the parent pointers don’t exist and the data are primitive. Therefore, random values are assigned to the exact variables in the data part of the object (see Section 7.3 for details on xE vs x). v o i d P o i n t : : R a n d o m iz e () { xE = random () % p; yE = random () % p; } 112 7.4.2 Com puteExact (Ob j e c t E Ob j e c t : : C o m p u teE x a ct ()) Let us assume that we are trying to test whether A is equal to zero by a theorem. Then we have already randomly perturbed all the arguments on which the value of A depends (by calling A -> R a n d o m ize ( ) ). Now we need to recompute A in exact arithmetic modulo p. This is done by calling C om p u teE xact. P o i n t l E & P o i n t l : : C o m p u teE x a ct () { / / Take t h e a l r e a d y r a n d o m i z e d >l i n e s and co m p u te t h e / / i n t e r s e c t i o n p o i n t . • * • co m p u te PE ( t y p e P o in t E ) r e t u r n PE } 7.4.3 MyConstructor ( v o i d O b j e c t : : M y C o n s t r u c t o r ()) This routine was earlier described in Section 6.4.2, but at the time only the history management part of the routine was emphasized. In this section we will elaborate on the theorem-detection part. In Section 6.4.2, we had mentioned that in line 7 of the program (see below, or on page 103), it is quite likely that the coincidence is by theorem and the theorem detection machinery is executed (lines 8-19). v o i d P o i n t 1 : : M y C o n s t r u c t o r () { 1 s t a t i c i n t E i ; / / w i l l b e s e t e q u a l t o 2 i n t R e l i a b i l i t y = 0; / / z e r o b y , d e f a u l t 3 i f (1 p a r a l l e l t o m) / / e x c e p t i o n 4 i f (E i == M axC oinc) 5 / / c h a n c e s h i g h t h a t c o i n c i d e n c e b y t h e o r e m 6 / / . . . c o i n c i d e n c e by t h e o r e m h a n d l i n g 7 w h i l e ( R e l i a b i l i t y < 5) { ; 8 H } 0 ] - > R a n d o m i z e ( ) ; 9 H [ 1 ] - > R a n d o m i z e ( ) ; 10 113 h lE = H [ 0 ] - > C o m p u t e E x a c t { ) ; 11 h£E = H [ 1 ] -> C om puteE K act () ; 12 i f ( h l E . p a r a l l e l E (h2.E) ) 13 R e l i a b i l i t y + + ; 14 e l s e 15 / / w e know t h a t c o i n c i d e n c e was b y 16 / / c o i n c i d e n c e . . . 17 b r e a k ; 18 } 19 e l s e { 20 . . . do t h e u s u a l s t u f f g o t o D; 21 } 22 . . . . o t h e r c o m p u t a t i o n 23 D :; 24 } 25 Let us assume that for the probability in consideration (this is how confident we want to be that the theorem in question is indeed true), we want to check the exact arithmetic computation five times. Initially Reliability is set to 0, and in line 7 we conjecture that the coincidence could be by theorem. So, we randomly perturb H [ 0 ] and H [ 1 ] (lines 9,10; these are the same variables which were P e r tu r b e d to eliminate the ambiguity, see Section 6.4.2); compute the exact versions of the half spaces in h lE and h2E (lines 11, 12); and check whether these two exact half-spaces are parallel to each other (line 13, no ambiguity possible). If the lines are parallel, the value for R e l i a b i l i t y is incremented. Otherwise, we know for certain that the coincidence can be removed. Once the value of Re 1 i a b i 1 i t y reaches five, say, the desired probability is attained and we conclude that the two half-spaces intersect at infinity. The same mechanism can be used to test for theorems in all the if-statements in the program. 114 7.5 Examples The following examples were used to check the coincidence by theorem implementation. The first is a Pappus Theorem Implementation. The theorem is more fully described in section 4.1.1.3 on page 46. According to the theorem, points C, D and E shown in Figure 21 are collinear. To test the theorem, we construct two triangles <aj,C,D> and <a2,C,E> and take their union. The perturbation approach applied blindly will try to separate CD from CE and fail. The program managed to detect that the two lines were coincident by theorem and maintained the truth symbolically. The result of the implementation appears in Figure 47. The next example is an instance of the theorem of Brianchon. This is the point-line dual of Pappus’ Theorem. The program ran successfully on the example. The result appears on Figure 48. 115 •Write [Display Figure 47: Pappus The )rem 116 ; Write Display iLocate Figure 48: Theorem of Brianchon C hapter 8: Conclusions However, the field (of Robust Geometry) is still in the mode o f illustrative examples, with different problems analyzed in very different ways. A general theory o f robustness, applicable to a wide range o f problems, is unavailable. Steve Fortune, 1993 The contributions of this thesis are discussed in Section 8.1. The limitations of the thesis are covered in Section 8.2 and extensions / future work follow in Section 8.3. 8.1 C ontributions The specific contributions of this thesis to the field of geometric modeling in the presence of floating point errors are as follows: • The approach is general, systematic, almost mechanical. The approach can be applied to any program that satisfies the model specified on page 14. • The program run:, fast. It is slower by a factor of about four to five compared to a corresponding floating point program. This is much faster than exact arithmetic (off the shelf exact arithmetic libraries: Integer libraries are slower by a factor of hundred; rational, ten thousand). Note that Yap is aiming for a factor of 10 for his exact arithmetic implementation s. 1 • The program generated by our approach is probabilistically provably correct under the defined model of correctness. Therefore a probability p can be set by the user and the approach generates a solution which is correct with probability p. The probability can be, however, made as high as needed. • Another contribution of the thesis is a new notion of correctness for robust geometric algorithms. This notion (and the accompanying taxonomy) clearly delineates the isst es of validity and correctness, which are mixed up in some of the earlier works. However note that these measures are made based on empirical results made on benchmarks. Therefore, a user can construct examples on which the program will be much slower than the above factor. 118 • A multiprecision interval class was implemented (with Joe Bautista). This class is available and can be used by other programs (not necessarily dealing with the approach propounded in this thesis). All multiprecision implementations known to us (including the one mentioned in the previous paragraph) have a constant overhead for each number stored. This overhead is incurred even if the number to be stored can be stored in a d o u b l e . A zero overhead multiprecision interval class designed (not implemented) as a part of this thesis incurs no overhead in case the number stored is a d o u b l e . A boundary evaluator for 2D Constructive Solid Geometry (CSG) has been implemented. The primitives are half-spaces defined by lines. In addition, vertices generated by the intersection of lines can be used to construct new lines. This demonstrates the applicability of the approach to a complex problem. The program has about 3500 lines of C++ code. All the profile tests (see Figure 36,37, 38 and 39) mn in a few seconds on the program running on a SPARC2. The results of interval growth tests are shown in Figure 41-44. The recomputing overhead in the absence of theorems and the presence of theorems is shown in Table 2. 8.2 Lim itations • The implementation exists only for two dimensional, planar objects. Although conceptually the problem is the same at a 2D -or 3D level, the complexity of the problem increases and the number of variables involved go up significantly. * This implementation supports only planar spaces. This is an important limitation because curved- surfaces handling requires a theoretical leap involving the handling of radicals. This implies an ability to deal with a system of polynomials (i.e. a variety) rather than a single polynomial. ♦ Currently only unidirectional constraints can be handled by the system. Therefore, if we have to declare two lines to be parallel, we have to construct the second line from the first one. There is no way of imposing the constraint on existing entities. The bidirectional constraints involve solutions of systems of 119 equations (polynomials) which reduces to the problem mentioned in the previous entry i.e. the ability to deal with a variety. 8.3 F uture Research Besides overcoming the previous limitations the following problems should be tackled in the future. • Applying the approach to the tolerances present in mechanical designs. • Combining intervals from different sources to reduce interval size. • A probabilistic modeling of intervals. Currently a worst case analysis of intervals is done. This can be replaced with a probabilistic model which considers the distribution of numbers in the interval space. I hate quotations. Tell me what you know. Ralph Waldo Emerson 120 C hapter 9: Bibliography [1] O. Aberth and M. J. Schaefer, “Precise computation using range arithmetic, using C++", ACM Transactions on Mathematical Software, Vol. 18, No. 4,pp. 481-491, December 1992. [2] G. Allen, “Testing the accuracy of solid modelers”, Computer Aided Engineering, Vol. 4, No. 6, pp. 50-54, June 1985. [3] M. Benouamer, D. Michelucei and B. Peroche, “Error-free boundary evaluation using lazy rational arithmetic': A detailed implementation”, Proceedings o f Second Symposium on Solid Modeling and Applications, Montreal Canada, pp. 115-126, May 19-21, 1993. [4] B. Briiderlin, “Detecting ambiguities: An optimistic approach to robustness problems in computational geometry”, Technical Report UUCS-90-003, Computer Science Department, University of Utah, April 1990. [5] B. Briiderlin, “Robust regularized set operations on polyhedra”, Proc. o f Hawaii International Conference on System Sciences, January 1991. [6] S. Cameron, “Efficient Intersection Tests for Objects Defined Constructively” International Journal o f Robotics Research, July 1988. [7] S. Cameron and J. R. Rossignac “Relationship between S-Bounds and Active Zones in Constructive Solid Geometry”, Proceedings o f the International Conference on the Theory and Practice o f Geometric Modeling, October 1988. [8] J. Canny “ The complexity of robot motion planning” , Ph.D. dissertation, Computer Science Department, MIT, 1987. [9] G. B. Dantzig, Linear Programming and Extensions. Princeton, New Jersey: Princeton University Press, 1963. [10] H. Desaulnisrs and N. F. Stewart, “Quasi-rectilinear r-sets”, Proceedings o f the Fourth Canadian Conference on Computational Geometry, St. John’s, Newfoundland, pp. 52-58, August 10-14, 1992. [11] H. Desaulniers and N. F. Stewart “Backward error analysis for floating point operations on rectilinear r-sets”, Technical publication #816, Department d ’informatique et de recherche operationnelle, University de Montreal, J anuary 1993. [12] H. Desaulniers and N. F. Stewart, “Robustness of numerical methods in geometric computation when problem data is uncertain” Computer Aided Design, Vo!, 25, No. 9, September 1993.. 121 [13] T. K. Dey, K. Sugihara and C. Bajaj. Tri angulations in Three Dimensions with Finite Precision Arithmetic, technical report CSD-TR-92-001, CAPO Report CER-92-01, Computer Science Department, Purdue University, January 1992. [14] D. P. Dobkin and D. Silver, “Recipes for geometry and numerical analysis, Part I: An empirical study”, Proceedings o f the Symposium on Computational Geometry, ACM, pp. 93 -105, 1988. [15] D. P. Dobkin and D. Silver, “Applied Computational Geometry: Towards robust solutions of basic problems”, Journal o f Computer and System Sciences, Vol 40, No. 1, pp. 70-87, 1990. [16] H. Edelsbrunner. Algorithms in Combinatorial Geometry. New York: Springer Vedag. 1987. [17] H. Edelsbrunner and E. P. Miicke, “Simulation of Simplicity: A technique to cope with degenerate cases in geometric algorithms”, ACM Annual Symposium, on Computational Geometry, pp. 118-133, June 1988 [18] J. Ely, “The application of variable precision interval arithmetic to the solution of problems in vortex dynamics”, in R. Vichnevetsky and J. J. H. Miller, Eds.; Proceedings o f the 13 th World Congress on Computation and Applied Mathematics, IMACS. Dublin, Ireland 1991, pp. 69-70. [19] I. Emeris and J. Canny, “An Efficient approach to removing geometric degeneracies”, Proceedings o f the Eighth Annual ACM Symposium on Computational Geometry, Berlin, Germany, June 1992. [20] I. Emiris and J. Canny, “A general approach to removing degeneracies”, Proc. 32nd Annual Symposium on Foundations o f Computer Science, pp. 405-413, 1991. [21] S. Fang and B. Briiderlin, “Robustness i;. Geometric Modeling — Tolerance Based Methods”, International Workshop on Computational Geometry CG ’91, Lecture Notes in Computer Science 553, Bern, Switzerland, March 1991. [22] S. Fang, “Robustness in geometric modeling”, Ph.D. Dissertation, Computer Science Department, University of Utah, 1992. [23] S. Fang. Personal Communication, 1993 [24] S. Fortune, “Stable Maintenance of point set triangulations in two dimensions’ , IEEE Annual Symposium on Foundations o f Computer Science, pp. 494-499, 1989. 122 [25] S. Fortune, ‘ ^Progress in Computational Geometry”, chapter 3, in R. Martin, ed., Directions in Geometric Computing. Information Geometers, 1993. [26] S. Fortune and C. J. Van Wyk, “Efficient Exact Arithmetic for Computational Geometry”, Proceedings, o f the 9th Annual Symposium on Computational Geometry, pp. 163-172, May 1993. [27] D. H. Green, and F. F. Yao, “Finite Resolution Computational Geometry”, IEEE Annual Symposium on Foundations o f Computer Science, pp. 143- 152, October 1986. [28] C. M. Hoffmann, J. E. Hopcroft and M. S. Karasick, “Towards implementing robust geometric computations”, ACM Annual Symposium on Computational Geometry, pp. 106-117, June 1988. [29] C. M. Hoffmann, “The Problems of Accuracy and Robustness in Geometric Computation”, IEEE Computer, Vol. 22, Number 3, pp. 31-41, March 1989. [30] J. Hong, “Proving by example and gap theorems”, Proc. 27th IEEE Symposium on Foundations o f Computer Science, pages 107-116, 1986. [31] J. E. Hopcroft and P. J. Kahn, “A paradigm for robust geometric algorithms”, Algorithmica, Vol. 7, pp. 339-380, 1992. [32] W. M. Kahan, “Interval Arithmetic Options in the Proposed IEEE Floating Point Arithmetic Standard”, in K. L. E. Nickel, Ed., Interval Mathematics 1980. New York: Academic Press Co., 1980, pp. 99-128. [33] T. C. Kao and G. D. Knott, “An efficient and numerically correct algorithm for 2D convex hull problem, BIT, 30, pp. 289-300, 1990. [34] M. Karasick, “On the Representation and Manipulations of Rigid Solids” , Ph.D. Dissertation, McGill University, 1989. [35] M. Karasick, D. Lieber and L. R. Nackman, “Efficient Delaunay triangulation using rational arithmetic”, ACM Transactions on Graphics, Vol. 10, No. 1, pp. 71-91, 1990. [36] D. Knuth. The Art o f Computer Programming. Vol II. Semi-Numerical Algorithms. Addison Wesley, 1969. [37] D. H. Laidlaw, W. B. Trombore, and J. P. Hugues, “Constructive Solid Geometry for Polyhedral Objects”, Computer Graphics, (SIGGRAPH ’86 proceedings), Vol. 20, No. 4, pp. 161-180, August 1986. [38] J. C. Latombe, Robot Motion Planning. Boston: Kluwer Academic Publishers, 1991. 123 [39] Z. Li and V. J. Milenkovic, “Constructing strongly convex hulls using exact or rounded arithmetic”, ACM Animal Symposium on Computational Geometry, pp. 234-243, 1990. [40] B. E. Meserve, Fundamental Concepts o f Geometry. New York: Dover Publishers, Inc., 1983, pp. 66. [41] V. J. Milenkovic, “Verifiable Implementations of Geometric Algorithms Using Finite Precision Arithmetic” , Ph.D. Dissertation, Carnegie Mellon University, 1988. [42] V. J. Milenkovic, “Verifiable Implementations of Geometric Algorithms Using Finite Precision Arithmetic”, Artificial Intelligence, Vol. 37, pp. 377-401, 1988. [43] V. J. Milenkovic, “Calculating approximate curve arrangement using rounded arithmetic”, ACM Annual Symposium on Computational Geometry, pp 197-207, 1989. [44] R. E. Moore, “Interval arithmetic and automatic error analysis in digital computing”, Ph.D. Dissertation, Stanford University, October 1962. [45] R. E. Moore, Methods and Applications o f Interval Analysis. Philadelphia: SIAM Studies in Applied Mathematics, 1979. [46] B. F. Naylor, “Binary Space Partition Trees: An Alternative Representation of Polytopes”, Computer Aided Design, Vol. 22, No. 2, March 1990. [47] T. Ottmann, G. Thiemt and C. Ullrich, “Numerical Stability of Geometric Algorithms5 1 , ACM Annual Symposium on Computational Geometry, pp. 119-125, June 1987. [48] F. P. Preparata and M. I. Shamos, Comp utational Geometry - An Introduction, New York : Springer Verlag, 1985. [49] A. A. G. Requicha, “Mathematical Models of Rigid Solid Objects”, TM- 28, Production Automation Project, University of Rochester, November 1977. Also available as Technical Report IRIS-253, Department of Computer S rience e, University of Sou teem California. [50] A. A. G. Requicha, “Toward a theory of geometric tolerancing”, Int. J. Robotics Research, Vol 2, No. 2, pp. 45 -60, Winter 1983. [51] J. R. Rossignac and H. B. Voelcker, “Active Zones in Constructive Solid Geometry for Redundancy and Interference Detection’5 , Technical Memorandum TM-52, Production Automation Project, Rochester University, May 1986. 124 [52] D. Salesin, “Epsilon Geometry: Building Robust Algorithms from Imprecise Computations” , Ph.D. Dissertation, Computer Science Department., Stanford University, 1991. [53] D. Salesin, J. Stolfi and L. Guibas, “Epsilon geometry: Building robust algorithms f rom imprecise computations”, ACM Annual Symposium on Computational Geometry, pp. 208-217,1989. [54] J. T. Schwartz, “Fast Probabilistic Algorithms for Verification of Polynomial Identities”, Journal o f the ACM, Vol 27, No. 4, pp. 701-717, October 1980. [55] M. Segal, “Using tolerances to guarantee valid polyhedral modeling results”, Computer Graphics (SIGGRAPH Proceedings), Vol. 24, No. 4, pp. 105-114, 1990, [56] M. Segal and C. Sequin, “Consistent Calculation for Solid Modeling”, ACM Annual Symposium on Computational Geometry, pp. 29-37, June 1985. [57] M. Segal and C. Sequin, “Partitioning polyhedral objects into nonintersecting parts”, IEEE Computer Graphics and Applications, pp. 53- 67, January 1988. / [58] A. J. Spyridi and A. A. G. Requicha, “Accessibility analysis for the automatic inspection of mechanical parts by coordinate measuring machines”, Proceedings o f the 1990 IEEE International Conference on Robotics and Automation, pp. 1284-1289, May 13-18, 1990. [59] A. J. Spyridi and A. A. G. Requicha, “Accessibility Analysis for Polyhedral Objects”, Proceedings o f European Conference onRobotics and Intelligent Systems (EURISCON), Corfu, Greece, June 23-28, 1991. [60] A. J. Spyridi, Personal Communication, 1992. [61] A. J. K. Stewart, “The Theory and Practice of Robust Geometric Computation or How to Build Robust Solid Modelers” , Ph.D. Dissertation, Technical Report 91-1229, Computer Science Department, Cornell University, 1991. [62] A. J. K. Stewart, “Robust point location in approximate polygons”, 1991 Canadian Conference on Computational Geometry, pp. 179-182, August 1991. [63] K. Sugihara and M. Iri, “A solid modelling system free from topological inconsistency”, Journal o f Information Processing, Vol. 12, No. 4, pp. 380- 393, 1989. 125 [64] Sun OS 4.0, Reference Manual. Last Change 21 October 1989. pages 1090- 1091. [65] R. Tilove, ‘‘ Set Membership Classification: A Unified Approach to Geometric Intersection Problems”, IEEE Transactions on Computers, Vol C-29, No. 10, pp. 874-883, October 1980. [66] R. Tisdale, “Automatic Rounding Error Analysis”, M.S. Dissertation, Department of computer science, University of California at Los Angeles, 1994. [67] C. Yap, “A geometric consistency theorem for a symbolic perturbation scheme”, ACM Annual Symposium on Computational Geometry, pp. 134- 142, June 1988. [68] C. Yap, “Towards Exact Geometric Computation”, July 30 1993, to appear in Proceedings 5th Canadian Conference on Computational Geometry. 126

Abstract (if available)

Linked assets