ref: e2426e3277266e230b352bb843c4f58a91a0eb2e
parent: 3582b9cc27bf2a83800c00b4a84873037e99642d
author: Sigrid Solveig Haflínudóttir <sigrid@ftrv.se>
date: Wed Jul 31 20:59:41 EDT 2024
qslv: sort and remove duplicates on the go
--- a/rast.c
+++ b/rast.c
@@ -75,47 +75,54 @@
return v + 1;
}
-static int
+static Sval *
qslv(Sval *qs, Sval a, Sval b, Sval c)
{- Sval p, q, d;
- int n;
+ Sval p, q, d, t;
+ int i;
- if(unlikely(is₀(a))){- if(is₀(b))
- return 0;
- qs[0] = -c/b;
- return qs[0] > 0 && qs[0] < 1;
+#define ins(q) \
+ if(q > 0 && q < 1){ \+ i = -1; \
+ if((t = qs[i]) < q){ \+ *qs++ = q; \
+ }else if(q < t){ \+ qs++; \
+ do{ \+ qs[i--] = t; \
+ qs[i] = q; \
+ t = qs[i-1]; \
+ }while(q < t); \
+ } \
}
- p = b/(2.0*a);
- d = p*p - c/a;
- if(d < ε)
- return 0;
-
- d = sqrt(d);
- n = 0;
- if((q = -d-p) > 0 && q < 1)
- qs[n++] = q;
- if((q = d-p) > 0 && q < 1)
- qs[n++] = q;
- return n;
+ if(unlikely(is₀(a))){+ if(!is₀(b)){+ q = -c/b;
+ ins(q);
+ }
+ }else{+ p = b/(2.0*a);
+ d = p*p - c/a;
+ if(d > ε){+ d = sqrt(d);
+ q = -d-p;
+ ins(q);
+ q = d-p;
+ ins(q);
+ }
+ }
+#undef ins
+ return qs;
}
static int
-Svalcmp(void *a_, void *b_)
-{- Sval a = *(Sval*)a_, b = *(Sval*)b_;
- return (a > b) - (a < b);
-}
-
-static int
qKL(SegQ *s₀, int jj, Sval px, Sval py, Sval *K, Sval *L)
{- Sval qs[1+4*2+1];
+ Sval *q, *e, qs[1+1+4*2+1];
Sval a, b, c, w, v;
+ int i, n, r;
SegQ s, z;
- int i, j, n, r;
/* transform */
s.p0.x = s₀->p0.x*jj - px;
@@ -129,43 +136,34 @@
#define within(v) ((w = e(v, x)) >= -ε && w <= 1+ε && (w = e(v, y)) >= -ε && w <= 1+ε)
/* clip to quadtree cell */
- n = 0;
+ q = qs;
+ *q++ = -1;
if(s.p0.x >= 0 && s.p0.x <= 1 && s.p0.y >= 0 && s.p0.y <= 1)
- qs[n++] = 0;
-
+ *q++ = 0;
c = s.p0.x;
a = c - 2*s.p1.x + s.p2.x;
b = 2*(s.p1.x - c);
- n += qslv(qs+n, a, b, c);
- n += qslv(qs+n, a, b, c-1);
+ q = qslv(q, a, b, c);
+ q = qslv(q, a, b, c-1);
c = s.p0.y;
a = c - 2*s.p1.y + s.p2.y;
b = 2*(s.p1.y - c);
- n += qslv(qs+n, a, b, c);
- n += qslv(qs+n, a, b, c-1);
- qsort(qs, n, sizeof(Sval), Svalcmp);
-
+ q = qslv(q, a, b, c);
+ q = qslv(q, a, b, c-1);
if(s.p2.x >= 0 && s.p2.x <= 1 && s.p2.y >= 0 && s.p2.y <= 1)
- qs[n++] = 1;
- j = 0;
- for(i = 0; i < n; i++){- v = qs[i];
- if(within(v))
- qs[j++] = v;
- }
-#undef within
+ *q++ = 1;
- /* remove duplicates */
n = 1;
- for(i = 1; i < j; i++){- v = qs[i];
- if(v != qs[i-1])
+ for(e = qs+1; e < q; e++){+ v = *e;
+ if(within(v))
qs[n++] = v;
}
+#undef within
/* calculate K & L on subsections */
r = 0;
- for(i = 0; i < n-1; i++){+ for(i = 1; i < n-1; i++){a = qs[i];
b = qs[i+1];
@@ -257,7 +255,7 @@
Cxy(SegQ *s, int ns, int jj, int px, int py, Sval *c, u64int *m, u64int tm)
{int (*f)(SegQ*, int, Sval, Sval, Sval*, Sval*);
- Sval K[4][2], L[4][2], q[6], j;
+ Sval K[4][2], L[4][2], x₀, x₁, x₂, y₀, y₁, y₂, psz;
u64int tx₀, tx₁, z, all;
u8int w;
int i;
@@ -272,20 +270,20 @@
tx₀ = 1ULL<<(px*jj + py);
tx₁ = 1ULL<<((px+1)*jj + py);
}
- j = 1.0/(Sval)jj;
- q[0] = j*px;
- q[1] = j*py;
- q[2] = q[0]+j;
- q[3] = q[1]+j;
- q[4] = q[2]+j;
- q[5] = q[3]+j;
+ psz = 1.0/(Sval)jj;
+ x₀ = psz*px;
+ y₀ = psz*py;
+ x₁ = psz+x₀;
+ y₁ = psz+y₀;
+ x₂ = psz+x₁;
+ y₂ = psz+y₁;
all = 0;
for(i = 0; i < ns; i++, s++){if((m[i] & tm) == 0 ||
- s->p0.x <= q[0] && s->p1.x <= q[0] && s->p2.x <= q[0] ||
- s->p0.x >= q[4] && s->p1.x >= q[4] && s->p2.x >= q[4] ||
- s->p0.y <= q[1] && s->p1.y <= q[1] && s->p2.y <= q[1] ||
- s->p0.y >= q[5] && s->p1.y >= q[5] && s->p2.y >= q[5])
+ s->p0.x >= x₂ && s->p1.x >= x₂ && s->p2.x >= x₂ ||
+ s->p0.x <= x₀ && s->p1.x <= x₀ && s->p2.x <= x₀ ||
+ s->p0.y >= y₂ && s->p1.y >= y₂ && s->p2.y >= y₂ ||
+ s->p0.y <= y₀ && s->p1.y <= y₀ && s->p2.y <= y₀)
continue;
K[0][0] = K[0][1] = 0;
@@ -298,12 +296,11 @@
L[3][0] = L[3][1] = 0;
z = 0;
f = s->p1.x == s->p2.x && s->p1.y == s->p2.y ? lKL : qKL;
-
w =
- (s->p0.x <= q[2] || s->p1.x <= q[2] || s->p2.x <= q[2])<<0 |
- (s->p0.x >= q[2] || s->p1.x >= q[2] || s->p2.x >= q[2])<<1 |
- (s->p0.y <= q[3] || s->p1.y <= q[3] || s->p2.y <= q[3])<<2 |
- (s->p0.y >= q[3] || s->p1.y >= q[3] || s->p2.y >= q[3])<<3;
+ (s->p0.x <= x₁ || s->p1.x <= x₁ || s->p2.x <= x₁)<<0 |
+ (s->p0.x >= x₁ || s->p1.x >= x₁ || s->p2.x >= x₁)<<1 |
+ (s->p0.y <= y₁ || s->p1.y <= y₁ || s->p2.y <= y₁)<<2 |
+ (s->p0.y >= y₁ || s->p1.y >= y₁ || s->p2.y >= y₁)<<3;
if(tx₀ == 0){if((w & 5) == 5) z |= f(s, jj, px+0, py+0, K[0], L[0]);