src/common/peano.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 * Generate an N-dimensional Peano space-filling curve
 * on the interval [0,1).
 */

#include "copyright.h"

extern double   floor();


static double
peanoC(n, t, e)                 /* compute Peano coordinate */
int     n;
double  t, e;
{
        register int    i;
        double  d;

        if (e >= 1./3.)
                return(t);
        d = 0.;
        while (t-d >= 1./3.)
                d += 1./3.;
        i = n;
        while (--i > 0)
                t *= 3.;
        t = 3.*(t - floor(t));
        i = 0;
        while (t >= 1.) {
                t -= 1.;
                i ^= 1;
        }
        t = peanoC(n, t, 3.*e);
        if (i)
                t = 1. - t;
        return(d + t/3.);
}


peano(p, n, t, e)               /* compute Peano point */
register double  p[];
int     n;
double  t, e;
{
        register int    i;
        register int  neg = 0;

        i = n;
        while (i-- > 0) {
                p[i] = peanoC(n, t, e);
                if (neg)
                        p[i] = 1. - p[i];
                t *= 3.;
                while (t >= 1.) {
                        t -= 1.;
                        neg ^= 1;
                }
        }
}
Revision:	2.3
Committed:	Tue Feb 25 02:47:21 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R6P1, rad3R5, rad3R6
Changes since 2.2:	+1 -56 lines
Log Message:	Replaced inline copyright notice with #include "copyright.h"
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.2	static const char RCSid[] = "$Id$";
3	greg	1.1	#endif
4			/*
5			* Generate an N-dimensional Peano space-filling curve
6			* on the interval [0,1).
7	greg	2.2	*/
8
9	greg	2.3	#include "copyright.h"
10	greg	1.1
11			extern double floor();
12
13
14			static double
15			peanoC(n, t, e) /* compute Peano coordinate */
16			int n;
17			double t, e;
18			{
19			register int i;
20			double d;
21
22			if (e >= 1./3.)
23			return(t);
24			d = 0.;
25			while (t-d >= 1./3.)
26			d += 1./3.;
27			i = n;
28			while (--i > 0)
29			t *= 3.;
30			t = 3.*(t - floor(t));
31			i = 0;
32			while (t >= 1.) {
33			t -= 1.;
34			i ^= 1;
35			}
36			t = peanoC(n, t, 3.*e);
37			if (i)
38			t = 1. - t;
39			return(d + t/3.);
40			}
41
42
43			peano(p, n, t, e) /* compute Peano point */
44	greg	1.2	register double p[];
45	greg	1.1	int n;
46			double t, e;
47			{
48			register int i;
49	greg	1.2	register int neg = 0;
50	greg	1.1
51			i = n;
52			while (i-- > 0) {
53			p[i] = peanoC(n, t, e);
54			if (neg)
55			p[i] = 1. - p[i];
56			t *= 3.;
57			while (t >= 1.) {
58			t -= 1.;
59			neg ^= 1;
60			}
61			}
62			}