#------------------------------------------------------------------------------- # INSTRUCTIONS # To patch from mangle1.4.1 to mangle1.4.3: # # cd # patch -p0 < path_to/mangle1.4.1_to_1.4.3.patch # #------------------------------------------------------------------------------- --- CHANGES.orig 2004-03-31 13:17:43.000000000 -0700 +++ CHANGES 2006-09-14 10:01:31.000000000 -0600 @@ -1,6 +1,104 @@ -------------------------------------------------------------------------------- +mangle1.4.3 +14 Sep 2006 + +Thanks to Colin Hill, working with Max Tegmark and Molly Swanson at MIT, +for finding a polygon that defeated mangle. + +1. The problem was numerical round-off, +and the immediate fix was to replace line 71 in split_poly.c + if (area != area_tot) { /* boundary intersects poly1 */ +with + if (area > area_tot) { /* boundary intersects poly1 */ +It should never happen that area < area_tot, since removing a boundary +of a polygon should never decrease its area, but thanks to numerical round-off +the area did decrease. Once that happened, mangle span its wheels. + +2. The killer polygon that was responsible for the numerical failure of item 1 +is discussed below. To improve mangle's defence against such polygons, +I changed the strategy for modifying the tolerance angle to multiple +intersections. Originally, if mangle detected an inconsistency in the +topology of the distribution of vertices around multiple intersections, +then mangle would double the tolerance angle and try again. +The new strategy is to try tightening as well as loosening the tolerance angle. +Mangle now see-saws between and tighter and looser tolerance angles +successive factors of 2 away from the original input tolerance angle. + +3. There was also a bug on line 182 of balkanize.c. + for (i = npoly; i < npolys; i++) { +should have been + for (i = 0; i < npolys; i++) { +This bug became evident as a result of a compiler change. +When tripped, the bug produces a segmentation violation, so if you never +experienced this bug before, it should never have caused an error. + +This is Colin Hill's killer polygon: +polygon 6 ( 5 caps, 1.0 weight, 0.000003265281337 str): + -0.2972995896945569 -0.5791945044073898 -0.7590432662449000 0.5786075961488696 + -0.1092635942528740 0.6422417176858356 -0.7586745303718128 0.5279437897883537 + -0.3927631897959760 -0.5199166454922474 -0.7585669110114012 -0.4794583692304901 + -0.0457867636138749 0.6498588533234551 -0.7586745303718132 -0.5279437897887224 + -0.3467216633109536 -0.5515406106135746 -0.7586745303718130 0.5279437897887302 +It is a long (4.3 degrees) thin (36 arcsec) rectangle, split by a diagonal +which is almost, but not quite, tangent at each end to the long direction. +Each near tangent end is both almost multiply intersecting, and almost kissing. + +If you balkanize the killer polygon with the default tolerance angle +for multiple intersections of 10^-5 arsec, then mangle1.4.1 will discard the +polygon as having zero area: + +% balkanize killer.poly - +---------------- balkanize ---------------- +snap angles: axis 2s latitude 2s edge 2s +multiple intersections closer than 1e-05s will be treated as coincident +... +1 polygons read from j6 +warning from balkanize: following polygons have zero area & are being discarded: + 0 +... + 0 polygons written to output + +Well, the polygon is thin, 36 arcsec wide, so maybe you don't mind losing +the polygon. But with snap tolerance angles of 2 arcsec, you'd think mangle +should keep it. It turns out that the multiple intersections at each end +of the polygon have separations that vary from 10^-8 arcsec in the thin +direction to 80 arcsec in the long direction, and this distribution of +separations conspires so that mangle finds no consistent topology +(satisfying the 64-bit check number) for tolerance angles to multiple +intersections anywhere between 10^-8 arcsec and 40 arcsec. But 40 arcsec +exceeds the 36 arcsec width of the polygon, so mangle concludes that the +area is zero. + +If the tolerance angle for multiple intersections is tightened to +5 x 10^-9 arcsec: + +% balkanize -m5e-9s killer.poly - + +then mangle does find the correct solution. That is, since the polygon +is already a valid polygon, balkanize spits back the killer polygon +unchanged. + +Thus the solution to dealing with this polygon is not to loosen the +tolerance, but to tighten it. The new mangle strategy does this. + +Max Tegmark suggests switching the entire mangle code from double (64-bit) +to quadruple (128-bit) precision. This may be a sensible thing to do +at some point. The need for high precision was always evident. +Inevitably, people are messing with masks that have arcsecond precision. +A 1 x 1 arcsec rectangle has an area of 2 x 10^-11 str, which is 2 x 10^-12 +of the whole sky. This is a factor of 10^4 larger than the effective +precision 2 x 10^-16 of 64-bit floats. The factor of 10^4 is ok, +but with ever more complicated masks appearing, it may no longer be +providing a sufficient safety margin. + +-------------------------------------------------------------------------------- +5 Sep 2005 + +Added '\r' to characters considered blank, so polygon files are dos-compatible. + +-------------------------------------------------------------------------------- mangle1.4.1 -30 Sep 2004 +30 Mar 2004 Many thanks to Michael Blanton, who is working on the SDSS mask, for bringing these two problems to light. --- src/balkanize.c.orig 2004-03-30 12:34:19.000000000 -0700 +++ src/balkanize.c 2006-09-10 09:57:19.000000000 -0600 @@ -8,6 +8,11 @@ #include "manglefn.h" #include "defaults.h" +/* define CARRY_ON_REGARDLESS if you want balkanize() to continue even when the number of polygons hits NPOLYSMAX; + if CARRY_ON_REGARDLESS is defined, then balkanize() will create a possibly incomplete polygon file of polygons */ +#undef CARRY_ON_REGARDLESS +//#define CARRY_ON_REGARDLESS + /* getopt options */ const char *optstr = "dqa:b:t:y:m:s:e:v:p:i:o:"; @@ -174,7 +179,7 @@ msg("balkanizing %d polygons ...\n", np); /* nullify all output polygons */ - for (i = npoly; i < npolys; i++) { + for (i = 0; i < npolys; i++) { polys[i] = 0x0; } @@ -244,7 +249,12 @@ if (n > npolys) { fprintf(stderr, "balkanize: total number of polygons exceeded maximum %d\n", npoly + npolys); fprintf(stderr, "if you need more space, enlarge NPOLYSMAX in defines.h, and recompile\n"); + n = npolys; +#ifdef CARRY_ON_REGARDLESS + break; +#else return(-1); +#endif } } /* copy down non-null polygons */ @@ -262,6 +272,8 @@ n = m + dm; if (dm == 0) break; } + /* too many polygons */ + if (n >= npolys) break; ip++; } msg("added %d polygons to make %d\n", dnp, np); @@ -296,7 +308,12 @@ if (n > npolys) { fprintf(stderr, "balkanize: total number of polygons exceeded maximum %d\n", npoly + npolys); fprintf(stderr, "if you need more space, enlarge NPOLYSMAX in defines.h, and recompile\n"); + n = npolys; +#ifdef CARRY_ON_REGARDLESS + break; +#else return(-1); +#endif } ip++; } --- src/garea.s.f.orig 2003-09-08 14:23:16.000000000 -0600 +++ src/garea.s.f 2006-09-14 09:28:59.000000000 -0600 @@ -26,7 +26,7 @@ C logical warn logical whole real*8 bik,cmi,cmik,cmk,d,darea,dph,ikchk,ikran, - * ph,phm,php,psi,si,xi(3),yi(3) + * ph,phm,php,psi,si,tolin,xi(3),yi(3) c * c * Area of surface of sphere of unit radius bounded by c * 1 - r.rp(i) < cm(i) (if cm(i).ge.0) @@ -57,8 +57,10 @@ c possible multiple intersection when dph < dphmin data dphmin /1.d-8/ c +c input tolerance to multiple intersections + tolin=tol C print *,'--------------------' -c come here with enlarged tolerance to multiple intersections +c come here with modified tolerance to multiple intersections 100 continue C write (*,'(a3,a20,4a24)') C * ' ','x','y','z', @@ -129,6 +131,8 @@ elseif (ni.gt.0) then c find ordering of intersection angles around i circle call findbot(phi,2*np,iord,ni) +C write (*,'("phi")') +C write (*,'(i4,2g24.16)') (j,(phi(k,j),k=1,2),j=1,np) c........contribution to area from each segment of i circle jpl=0 c come here to do another segment @@ -150,7 +154,7 @@ C warn=.true. C print *, C * '*** warning from garea: near multiple intersection at', -C * i,': edge',km,kp,' dph=',dph +C * i,': edge',km,kp,' dph =',dph endif c........segment satisfies conditions nbd=nbd+1 @@ -229,26 +233,18 @@ c write (*,'(i3,4g24.16)') c * (j,(rp(i,j),i=1,3), c * cm(j),j=1,np) -c retry with enlarged tolerance - if (tol.le.0.d0) then - tol=1.d-15 - else - tol=tol*2.d0 - endif +c retry with modified tolerance + call gtol(tol,tolin) goto 100 elseif (tol.gt.0.d0) then C print *,'... from garea: success at tol =',tol endif c--------add/subtract 2*pi's to area retry=garpi(area,iarea,rp,cm,np,whole,nbd0m,nbd0p,nbd,nmult) -c retry with enlarged tolerance +c retry with modified tolerance if (retry.eq.1) then C warn=.true. - if (tol.le.0.d0) then - tol=1.d-15 - else - tol=tol*2.d0 - endif + call gtol(tol,tolin) goto 100 endif c--------done --- src/gspher.s.f.orig 2003-09-08 14:24:48.000000000 -0600 +++ src/gspher.s.f 2006-09-14 09:12:31.000000000 -0600 @@ -28,7 +28,7 @@ logical whole real*8 bik,ci,cmi,cmik,cmk,cti,ctk,ctpsi, * d,darea,dbound(2),dph,dvert(2),ikchk,ikran, - * ph,phi,phii,phm,php,psi,psip,ri,rii,si,sk,sqrt4pi,t, + * ph,phi,phii,phm,php,psi,psip,ri,rii,si,sk,sqrt4pi,t,tolin, * xi(3),yi(3) c * c * Spherical transform of region W of sphere of unit radius bounded by @@ -173,8 +173,10 @@ c possible multiple intersection when dph < dphmin data dphmin /1.d-8/ c +c input tolerance to multiple intersections + tolin=tol C print *,'--------------------' -c come here with enlarged tolerance +c come here with modified tolerance 100 continue c initialise error flag to no error ldegen=.false. @@ -534,11 +536,7 @@ c print *,'*** from gspher: at tol =',tol, c * ', ikchk=',ikchk,' should be 0' c retry with enlarged tolerance - if (tol.le.0.d0) then - tol=1.d-15 - else - tol=tol*2.d0 - endif + call gtol(tol,tolin) goto 100 endif c--------add/subtract 2*pi's to area @@ -550,14 +548,10 @@ i=nint((w(1,1)*sqrt4pi+area)/TWOPI) endif w(1,1)=w(1,1)-i*TWOPI/sqrt4pi -c retry with enlarged tolerance +c retry with modified tolerance if (retry.eq.1) then C warn=.true. - if (tol.le.0.d0) then - tol=1.d-15 - else - tol=tol*2.d0 - endif + call gtol(tol,tolin) goto 100 endif c--------done --- src/gsubs.s.f.orig 2003-09-08 14:25:51.000000000 -0600 +++ src/gsubs.s.f 2006-09-14 09:26:42.000000000 -0600 @@ -1261,3 +1261,29 @@ return end c +c----------------------------------------------------------------------- + subroutine gtol(tol,tolin) + real*8 tol,tolin +c * +c * Modify tolerance tol to multiple intersections. +c * The tolerance tol is changed by successive factors of 2 +c * from the original input tolerance tolin, +c * see-sawing between smaller and larger values. +c * + if (tolin.le.0.d0) then + if (tol.le.0.d0) then + tol=1.d-15 + else + tol=tol*2.d0 + endif +c see-saw tolerance between smaller and larger values + else + if (tol.ge.tolin) then + tol=tolin*tolin/tol/2.d0 + else + tol=tolin*tolin/tol + endif + endif + return + end +c --- src/gvert.s.f.orig 2003-09-08 14:26:01.000000000 -0600 +++ src/gvert.s.f 2006-09-14 09:12:35.000000000 -0600 @@ -26,7 +26,7 @@ integer i,ii,ik,iseg,iv,ive,j,jm,jml,jmu,jp,jpl,jpu,km,kp, * mve,ni,scmi C logical warn - real*8 amve,cmi,dph,ikchk,ph,phm,php,si,xi(3),xv,yi(3),yv,zv + real*8 amve,cmi,dph,ikchk,ph,phm,php,si,tolin,xi(3),xv,yi(3),yv,zv c * c * Vertices, plus points on edges, of the polygon defined by c * 1 - r.rp(i) < cm(i) (if cm(i).ge.0) @@ -120,8 +120,10 @@ c data dphmin /1.d-8/ data big /1.d6/ c +c input tolerance to multiple intersections + tolin=tol C print *,'--------------------' -c come here with enlarged tolerance +c come here with modified tolerance 100 continue c initialise error flag to no error ldegen=.false. @@ -308,12 +310,8 @@ C warn=.true. C print *,'*** from gvert: at tol =',tol, C * ', ikchk=',ikchk,' should be 0' -c retry with enlarged tolerance - if (tol.le.0.d0) then - tol=1.d-15 - else - tol=tol*2.d0 - endif +c retry with modified tolerance + call gtol(tol,tolin) goto 100 elseif (tol.gt.0.d0) then C print *,'... from gvert: success at tol =',tol --- src/gvlim.s.f.orig 2003-09-08 14:26:15.000000000 -0600 +++ src/gvlim.s.f 2006-09-14 09:12:40.000000000 -0600 @@ -28,7 +28,7 @@ integer i,ii,ik,iphi,iseg,iv,j,jm,jml,jmu,jp,jpl,jpu,km,kp,ni,scmi C logical warn real*8 cmi,cmik,dph,ikchk,ph,phin,phif,phm,phmax,phmin,php, - * si,sik,xi(3),xv,yi(3),yv,zv + * si,sik,tolin,xi(3),xv,yi(3),yv,zv c * c * Lifted mostly from gvert and gcmlim. c * @@ -97,7 +97,9 @@ c data dphmin /1.d-8/ data big /1.d6/ c -c come here with enlarged tolerance +c input tolerance to multiple intersections + tolin=tol +c come here with modified tolerance 100 continue c initialise error flag to no error ldegen=.false. @@ -321,12 +323,7 @@ if (ikchk.ne.0.d0) then c print *,'*** from gvlim: ikchk=',ikchk,' should be 0' C warn=.true. -c retry with enlarged tolerance - if (tol.le.0.d0) then - tol=1.d-15 - else - tol=tol*2.d0 - endif + call gtol(tol,tolin) goto 100 endif c--------order vertices right-handedly about polygon --- src/rdangle.c.orig 2003-06-03 10:29:33.000000000 -0600 +++ src/rdangle.c 2005-09-05 09:13:01.000000000 -0600 @@ -20,7 +20,7 @@ */ int rdangle(char *word, char **next, char unit, double *angle) { - const char *blank = " \t\n"; + const char *blank = " \t\n\r"; const char *number = "+-.0123456789eE"; int i, ird, sgn; --- src/rdmask.c.orig 2004-03-31 13:17:43.000000000 -0700 +++ src/rdmask.c 2005-09-05 09:12:51.000000000 -0600 @@ -114,7 +114,7 @@ */ int rdmask(char *name, format *fmt, int npolys, polygon *polys[/*npolys*/]) { - char input[] = "input"; + char *input = "input"; char *line_rest, *word; int ird; int npoly = 0; @@ -205,7 +205,7 @@ */ char *get_keyword(char *str, char **str_rest, format *fmt) { - const char *blank = " \t\n"; + const char *blank = " \t\n\r"; char c; int id, ird, n; @@ -267,7 +267,7 @@ const char *end_fmt = "%d"; const char *unit_fmt = " %c"; char *word; - char blank[] = " \t\n"; + char *blank = " \t\n\r"; int ird, iscan, nholes; size_t word_len; @@ -483,8 +483,8 @@ int get_n(format *fmt, char *lin) { #define WARNMAX 3 - char blank[] = " \t\n"; - char plural[] = "s"; + char *blank = " \t\n\r"; + char *plural = "s"; static int warn = 0; char *ch; int n; @@ -778,7 +778,7 @@ static vertices *vert = 0x0; static int *ev = 0x0; - const char *blank = " \t\n"; + const char *blank = " \t\n\r"; char unit; char *next, *word; int anti, i, iang, iedge, iev, ird, iv, ive, ivert, nedge, nv, nve, nvert, reverse; --- src/split_poly.c.orig 2003-06-15 15:35:31.000000000 -0600 +++ src/split_poly.c 2006-09-12 08:52:50.000000000 -0600 @@ -68,8 +68,7 @@ poly->cm[np1 + ip] = 2.; /* suppress boundary to be tested */ ier = garea(poly, &tol, verb, &area); /* area of intersection sans boundary */ poly->cm[np1 + ip] = cm; /* restore tested boundary */ - if (area != area_tot) { /* boundary intersects poly1 */ - + if (area > area_tot) { /* boundary intersects poly1 */ /* poly2 splits poly1, but do not actually split */ if (!poly3) return(2); @@ -112,6 +111,10 @@ /* poly1 successfully split into poly1 and poly3 */ return(2); + } else if (area < area_tot) { + /* area should be >= area_tot because suppressing a boundary of poly should always increase its area; + but this can happen because of roundoff */ + //fprintf(stderr, "split_poly: area %.16g of polygon with boundary %d suppressed should be >= area %.16g of polygon\n", area, ip, area_tot); } }