| Ilya Zakharevich on Sat, 13 Jul 2024 07:25:22 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| Re: Segfault in read() |
On Fri, Jul 12, 2024 at 07:25:14PM +0200, Bill Allombert wrote:
> OK, but what is "pos" supposed to be ?
>
> What dbg_x(pos) gives ?
Oups, I remembered that there was something like dbg_x(), but after
inspecting a couple of dbg_*() functions gave up — too early! I
append the results.
> For what I see, pos is a vector with 5 components, whose first component is a t_LIST
> which is corrupted.
>
> So it is unclear if the bug is in writebin or read.
Now I it seems that I fixed most of MY bugs — only PARI’s shortcomings
and limitations seem to be hit now (but this took about 250 millions
tries of fuzzing!). So I just attach the code (with one “workaround”
commented out — to trigger the bug I was tracing through; it is
irrelevant though to the segfault in question).
Load the code, then do
neg1=neg2=POS1=POS=pos2=[1..2] \\ buckets for debugging
fuzz_convhull3D_avoids0(19,[5],0x1,3,4)
Hope this helps,
Ilya
⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜
[&=0000000003984f20] LIST(lg=3,CLONE):2900000000000003 (subtyp=0,lmax=32):000000000388fb50 000000000388fac0 0000000003a03260 0000000003a033a0 0000000003a03300
1st component = [&=000000000388fb50] VEC(lg=3,CLONE):2300000000000003 000000000388fb88 000000000388fb68
1st component = [&=000000000388fb88] REAL(lg=2):0400000000000002 (0,expo=-129):1fffffffffffff7f
2nd component = [&=000000000388fb68] REAL(lg=4):0400000000000004 (+,expo=-2):5ffffffffffffffe ffffffffffffffff ffffffffffffffff
2nd component = [&=000000000388fac0] VEC(lg=3,CLONE):2300000000000003 000000000388faf8 000000000388fad8
1st component = [&=000000000388faf8] REAL(lg=2):0400000000000002 (0,expo=-128):1fffffffffffff80
2nd component = [&=000000000388fad8] REAL(lg=4):0400000000000004 (-,expo=-1):dfffffffffffffff c000000000000000 0000000000000000
3rd component = [&=0000000003a03260] VEC(lg=3,CLONE):2300000000000003 0000000003a03298 0000000003a03278
1st component = [&=0000000003a03298] REAL(lg=4):0400000000000004 (+,expo=0):6000000000000000 8000000000000000 0000000000000000
2nd component = [&=0000000003a03278] REAL(lg=4):0400000000000004 (-,expo=-2):dffffffffffffffe 8000000000000000 0000000000000000
4th component = [&=0000000003a033a0] VEC(lg=3,CLONE):2300000000000003 0000000003a033d8 0000000003a033b8
1st component = [&=0000000003a033d8] REAL(lg=3):0400000000000003 (-,expo=-128):dfffffffffffff80 a000000000000000
2nd component = [&=0000000003a033b8] REAL(lg=4):0400000000000004 (+,expo=-1):5fffffffffffffff ffffffffffffffff fffffffffffffffb
5th component = [&=0000000003a03300] VEC(lg=3,CLONE):2300000000000003 0000000003a03338 0000000003a03318
1st component = [&=0000000003a03338] REAL(lg=4):0400000000000004 (+,expo=1):6000000000000001 8000000000000000 0000000000000000
2nd component = [&=0000000003a03318] REAL(lg=4):0400000000000004 (-,expo=-1):dfffffffffffffff 8000000000000000 0000000000000000
v3_prod(a,b)=[a[2]*b[3]-a[3]*b[2],a[3]*b[1]-a[1]*b[3],a[1]*b[2]-a[2]*b[1]];
v2_prod(a,b)=/*print(a[1]*b[2]-a[2]*b[1],a,b);*/a[1]*b[2]-a[2]*b[1];
v3_abs2(a)=a[1]^2+a[2]^2+a[3]^2;
v2_abs2(a)=a[1]^2+a[2]^2;
v3_sprod(a,b)=a[1]*b[1]+a[2]*b[2]+a[3]*b[3];
v2_sprod(a,b)=a[1]*b[1]+a[2]*b[2];
convhull_above(Xs,Ys) = { \\ Assume Xs grows. Return x[k], y[k], k, a, b of (y=ax+b) for corners; for the 1st entry a, b are fake
if(#Ys<2, error("Too few points"));
my(k=1 /* , pr_sl = 1+(Ys[2]-Ys[1])/(Xs[2]-Xs[1] */);
Newton_tmp_o=vector(#Ys); \\ XXXX Can't make lexical: need to access from a PARI’s NON-closure
Newton_tmp_o[1]=[Xs[1],Ys[1],1,oo,0];
my(addPt(x,y,K,k)= \\ add one point (index=K) to the known top convex hull of preceding points (with k vertices detected)
for(dk=0,k-1, [neg1[1],neg1[2]]=[x, Newton_tmp_o[k-dk][1]]; my(kk,s=(y-Newton_tmp_o[k-dk][2])/(x-Newton_tmp_o[k-dk][1]),S=Newton_tmp_o[k-dk][4]);
if(S>=s,Newton_tmp_o[kk=(k-dk+(S>s))]=[x,y,K,if(S>s, s, (y-Newton_tmp_o[k-dk-1][2])/(x-Newton_tmp_o[k-dk-1][1])),0];
return(kk)));
error("Panic310")); \\ returns the new value of detected slopes
for(K=2,#Ys, k=addPt(Xs[K],Ys[K],K,k));
for(K=2,k, Newton_tmp_o[K][5] = Newton_tmp_o[K][2]-Newton_tmp_o[K][1]*Newton_tmp_o[K][4]);
Newton_tmp_o[1..k]; \\ format of the entries: [x,y,k,slopeLeft]; the leftmost slopeLeft is Min unless #==1 and PreferMax>=0 (when it is Av(Min,Max) or Max)
}
addhelp(convhull_above, "Returns the corners of the top part of the convex hull; Xs are assumed to grow (strictly). The format of a corner is [xK,yK,K,s,b] (with yK=s xK+b), the reported slope s is of the edge to the left of the corner.")
\\ tstA(n,L0=0,L1=1)=my(v=vector(n,N,random(20)));r=convhull_above([1..n],v);[v,r]
\\ install("gcmp","iGG") \\ as cmp(), but does numeric comparison on numbers
\\ DANGER: 0e-16 < 1e-30 is false! So actually one should check !x || x<e0 !!!
\\ 1) Skip starting-earlier entries; 2)
convhulls_maxsum(P,e0) = { \\ P[1] and P[2] each list (distinct) vertices of graphs of two piecewise-linear functions (in ↗-order of X-coordinates).
if(#(P[1]) && #(P[2]),, return([-oo])); \\ Returns [max(f1+f2),P_index(1|2),offset in indexed, (larger of two, if applicable) offset in the other] or [-oo].
\\ invariants: work indices the “graph” with the smaller max-of-processed Xs; mx takes into account points up to this smaller max
my(Ind=[1,1],diffs=abs(P[1][1][1]-P[2][1][1])>e0, work = if(P[1][1][1]<P[2][1][1],1,2), v, vP, vv, yy, ind, mx=[-oo]); \\ mx=if(cmp,-oo,);
if(diffs, while(Ind[work] < #(P[work]) && P[work][Ind[work]++][1]<P[3-work][1][1],), \\ Skip where f1+f2 is undefined (may break the invariant)
if(#(P[work])==1, return([P[1][1][2]+P[2][1][2],1,1,1]), Ind[work]++)); \\ ensure the invariant
if(P[work][Ind[work]][1] < P[3-work][1][1]-e0, return([mx[1],work,Ind[work],1])); \\ outside the range of intersection
work = 3 - work; \\ restore the invariant;
while(1, \\ the smaller P[work][ind] is not yet counted into mx, and is between two last Xs in the other P
vP = (v = P[3-work])[(ind = Ind[3-work])-1]; vv = P[work][Ind[work]];
\\ print("P="P"; work="work"; Ind="Ind"; ind="=ind);
\\ print("dist="P[work][Ind[work]][1] - v[ind][1],"; between: "P[work][Ind[work]][1],", "v,"; slope=",(P[3-work][ind-1][2]-v[ind][2])/(P[3-work][ind-1][1] - v[ind][1])", val="=v[ind][2]);
\\ print(" dist*slope=",(P[work][Ind[work]][1] - v[ind][1])*(P[3-work][ind-1][2]-v[ind][2])/(P[3-work][ind-1][1] - v[ind][1]));
\\ DANGER: linear interpolation between (x₁,y₁) and (x₂,y₂) is very tricky with PARI arithmetic. It is easy to get answers
\\ “outside” of the range [y₁,y₂] (or “partially outside” — when precision is taken into account).
\\ First of all, values ≈ 0 appearing due to rounding errors would have only 1-word precision; likewise for (or values near
\\ ∞ obtained by projective transforms. So interpolating with y₂ ≈ ∞ when choosing x₂ AS REFERENCE, when x≳x₁ would
\\ essentially do y₂-Δy, which is ∼0 with the ABSOLUTE error of this “y₂ ≈ ∞”. Since the relative error of y₂ is
\\ 1 word, and |1/y₂| is ∼ precision, this absolute error is expected to be 1<<((n-1)*WORDWID) with n words of precision).
\\ (This makes the range of about Δx>>WORDWID near x₁ “dangerous”.
\\ Second: even when x₁ is the reference point and x≈x₁, it may happen that the difference is approx0, as in 0.e-38). Then
\\ bigneg*approx0 becomes bigapprox0, so looses the info about the sign (why it may have “big range”). Hence
\\ neg+bigneg*approx0 may be >=0. This seems to only affect the case when x-x₁≡0 in PARI.
\\ Third: the erroneous almost-vertical slopes are ∼1/lsb, and the error in calculation of x⸣s is ∼lsb. This gives the
\\ error in y as ∼1 — unacceptable. So we must ignore small x-x₁ and small x-x₂ as well… (Assuming that the caller
\\ ensures that small Δx mean small Δy.)
yy = vv[2] + if(!(yy=vv[1] - vP[1]) || yy<e0, /* print(1); v = v[ind]; */ vP[2],
!(yy=vv[1] - (v = v[ind])[1]) || yy>-e0, /* print(2); */ v[2],
abs(vv[1]-vP[1]) < abs(vv[1]-v[1]),
/* print(3); */ vP[2] + (vP[2]-v[2])/(vP[1] - v[1])*(vv[1] - vP[1]),
/* print(4); */ v[2] + (vP[2]-v[2])/(vP[1] - v[1])*(vv[1] - v[1]));
\\ print(" yy="yy", mx="mx", vtx=",P[work][Ind[work]][2]", vP[2]="vP[2]" distP=",(vv[1] - vP[1])" dy=",(vP[2]-v[2])/(vP[1] - v[1])*(vv[1] - vP[1]));
if(mx[1]<yy, mx=[yy,work,Ind[work],Ind[3-work]]);
Ind[work]++;
if(Ind[work] > #(P[work]), return(mx)); \\ all in “work” accounted; the last indexed in 3-work is out-of-range
if(P[work][Ind[work]][1] > P[3-work][ind][1], work=3-work); \\ restore the invariant;
);
return(mx);}
addhelp(convhulls_maxsum, "Assume P[1] and P[2] each list (distinct) vertices (in increasing order of X-coordinates) of graphs of two piecewise-linear functions f and g. Returns the array [max(f+g),W,indA,indB] with the last 3 entries indicating where this max is achieved: on the indA'th vertex of P[W], and either on indB'th vertex of the other P[n], or on the edge before tindB'th vertex. max is -oo if the domains do not intersect; then indA and indB are for vertices with the closest x-coordinate. e0 gives the leeway in checking that the domains do not intersect.")
/* do_CK() = {
ck_convhulls_maxsum([[], []],[-oo]);
ck_convhulls_maxsum([[], [[1,2]]],[-oo]);
ck_convhulls_maxsum([[[1,2]], []],[-oo]);
ck_convhulls_maxsum([[[1,2]], [[2,3]]],[-oo,1,1,1]);
ck_convhulls_maxsum([[[3,2]], [[2,3]]],[-oo,2,1,1]);
ck_convhulls_maxsum([[[1,2]], [[1,3]]],[5,1,1,1]);
ck_convhulls_maxsum([[[1,2]], [[1,3],[5,7]]],[5,1,1,2], 0b1010);
ck_convhulls_maxsum([[[1,2],[3,5],[4,-3]], [[1,3],[5,7]]],[10,1,2,2]);
ck_convhulls_maxsum([[[1,2],[3,5],[5,4]], [[1,3],[4,7],[6,6]]],[11+1/2,2,2,3]);
ck_convhulls_maxsum([[[-10,-10],[1,2]], [[1,3]]],[5,2,1,2],0b100);
ck_convhulls_maxsum([[[-10,-10],[1,2]], [[1,3],[5,7]]],[5,2,1,2],0b1100);
ck_convhulls_maxsum([[[-10,-10],[1,2],[3,5],[4,-3]], [[1,3],[5,7]]],[10,1,3,2]);
ck_convhulls_maxsum([[[-10,-10],[1,2],[3,5],[5,4]], [[1,3],[4,7],[6,6]]],[11+1/2,2,2,4]);
ck_convhulls_maxsum([[[-10,10],[1,2]], [[1,3]]],[5,2,1,2],0b100);
ck_convhulls_maxsum([[[-10,10],[1,2]], [[1,3],[5,7]]],[5,2,1,2],0b1100);
ck_convhulls_maxsum([[[-10,10],[1,2],[3,5],[4,-3]], [[1,3],[5,7]]],[10,1,3,2]);
ck_convhulls_maxsum([[[-10,10],[1,2],[3,5],[5,4]], [[1,3],[4,7],[6,6]]],[11+1/2,2,2,4]);
ck_convhulls_maxsum([[[-10,10],[0,3],[1,2]], [[1,3]]],[5,2,1,3],0b100);
ck_convhulls_maxsum([[[-10,10],[0,3],[1,2]], [[1,3],[5,7]]],[5,2,1,3],0b1100);
ck_convhulls_maxsum([[[-10,10],[0,3],[1,2],[3,5],[4,-3]], [[1,3],[5,7]]],[10,1,4,2]);
ck_convhulls_maxsum([[[-10,10],[0,3],[1,2],[3,5],[5,4]], [[1,3],[4,7],[6,6]]],[11+1/2,2,2,5]);
\\ ckNear_convhulls_maxsum(1e-30,1e-33); \\ ??? NOT WORKING YET. (Handle in the caller!)
}
*/
ck1t_convhulls_maxsum(P,o,f,e0) = my(OO = convhulls_maxsum(P, e0)); if(((o[1]!=OO[1]) && (abs(o[1]-OO[1])>1e-30)) || (!f && o!=OO && o[2..4] != OO[2..4] && !(o[3..4]==OO[3..4] && OO[3..4]==[1,1])), print("convhulls_maxsum:" P " -> " OO ". Expecting: " o". diff(exp,rc)="o[1]-OO[1]));
ck1_convhulls_maxsum(P,o,f,e0) = ck1t_convhulls_maxsum(P,o,f,e0); ck1t_convhulls_maxsum(1.*P,concat([if(o[1]>-oo,1.*o[1],o[1])],o[2..#o]),f,e0); ck1t_convhulls_maxsum([1.*P[1],P[2]],concat([if(o[1]>-oo,1.*o[1],o[1])],o[2..#o]),f,e0); ck1t_convhulls_maxsum([P[1],1.*P[2]],concat([if(o[1]>-oo,1.*o[1],o[1])],o[2..#o]),f,e0);
ck2_convhulls_maxsum(P,o,f,e0) = ck1_convhulls_maxsum(P,o,bitand(f,0b1),e0); ck1_convhulls_maxsum([P[2],P[1]],if(#o<2, o, [o[1],3-o[2],o[3],o[4]]),bitand(f,0b10),e0);
ck_convhulls_maxsum(P,o,f,e0) = my(V); ck2_convhulls_maxsum(P,o,f,e0); \
ck2_convhulls_maxsum([vector(#(P[1]),k,[6-(V=P[1][#(P[1])+1-k])[1],V[2]]),vector(#(P[2]),k,[6-(V=P[2][#(P[2])+1-k])[1],V[2]])],if(#o<2, o, [o[1],o[2],#(P[o[2]])+1-o[3],#(P[3-o[2]])+if(o[4]==1,1,2)-o[4]]),f>>2,e0);
ckNear_convhulls_maxsum(e0,st) = foreach([[1,2],[2,1]], v1, foreach([[3,4],[4,3]], v2, forvec(X1=[[1,4],[1,4]], forvec( X2=[[1,4],[1,4]], ck1t_convhulls_maxsum([[[1+X1[1]*st,v1[1]],[1+X1[2]*st,v1[2]]],[[1+X2[1]*st,v2[1]],[1+X2[2]*st,v2[2]]]],[6],0b1111,e0),2),2)));
\\ May modify elements of Vs. (“the external normal” unimplemented yet ??? The normal exists if 0 is on boundary too, but we ignore this)
convhull3D_avoids0(Vs,dbg,eps=2^-floor(0.66*bitprecision(1.)), e0=eps) = { \\ Given a collection Vs of non-0 3D vectors, returns the external normal if 0 is not in the convex hull, otherwise 0
if(!#Vs, return(1));
my( pos=List(), neg=List(), haveZ, Zpos=List(), Zneg=List(), ZZpos, ZZneg, Zslope=-oo, Znegmin=[oo], Zposmin=Znegmin, Znegmax=[-oo], Zposmax=Znegmax,
perm, v, vv, s1, s2, xRefl=1, yRefl=1, asX=1, asY=2, s, numcmp(x,y)=if(x==oo||y==oo,if(x<y,-1,x>y),x-y)); \\ temporaries, and parameters of linear transformations
e0 *= vecmax(apply(vecmax,abs(Vs))); \\ numcmp(x,y)=if(x<y,-1,x>y): ??? lex() and gcmp() are wrong for avoiding x-y==0
until(haveZ,
for(k=1,#Vs, if((v=Vs[k][3])>e0, listput(~pos,Vs[k]), v<-e0, listput(~neg,Vs[k]), \\ ① Need an “arbitrary” hyperplane! Assume ≈ to v₃=0:
(v=Vs[k][2])>e0, listput(~Zpos,Vs[k]), v<-e0, listput(~Zneg,Vs[k]), \\ but rotated a bit along the 2nd coordinate:
(v=Vs[k][1])>e0, ZZpos++, v<-e0, ZZneg++, return(0))); \\ and rotated a bit along the 1st:
if(dbg, print("pos="pos", Zpos="Zpos", neg="neg", Zneg="Zneg", ZZ(+,-)=",[ZZpos,ZZneg]));
if(#Zpos || #Zneg || ZZpos || ZZneg, haveZ++, \\ ② Ensure that one of the projectivizations is at ∞
#pos && #neg, for(k=1,#Vs, v=Vs[k][2]/Vs[k][3]; if(v > Zslope, vv=k; Zslope = v));\\ Find the minimal rotation to put on ∞-line
if(Zslope < .5, s1 = 1/2-Zslope; Zslope = 1/2; for(k=1,#Vs, Vs[k][2]+=Vs[k][3]*s1)); \\ Make V₃=const shear to 2nd if needed
s2 = Zslope = 1/Zslope;
for(k=1,#Vs, Vs[k][3]-=Vs[k][2]*s2); \\ Make V₂=const shear to 3rd
pos=List(); neg=List(); Zpos=List(); Zneg=List(); ZZpos=ZZneg=0; Vs[vv][3] = 0, \\ Due to rounding, it may be not 0
return(if(#pos, [0,0,-1], [0,0,1]))); \\ Triggers before the linear transforms only
);
if(ZZpos && ZZneg, return(0)); \\ , ZZpos=(!!ZZpos)-!!ZZneg, listput(~Zpos, [ZZpos,0,0])); \\ XXXX Would not survive translation to v below???
\\ ③ Find the convex hull of Zpos and -Zneg (BUG: XXXX should be Zneg !!! ???)
\\ (instead of shifting neg in the direction of Zneg·∞ in (?), can shift pos in the direction of -Zneg·∞)
for(k=1,#Zpos, vv = Zpos[k]; v = vv[1]/vv[2]; if(v<Zposmin[1], Zposmin=[v,vv]); if(v>Zposmax[1], Zposmax=[v,vv]));
for(k=1,#Zneg, vv = Zneg[k]; v = vv[1]/vv[2]; if(v<Znegmin[1], Znegmin=[v,vv]); if(v>Znegmax[1], Znegmax=[v,vv]));
if(dbg, print("Zposmin="Zposmin", Zposmax="Zposmax", Znegmin="Znegmin", Znegmax="Znegmax));
if(max(Znegmin[1],Zposmin[1]) <= if((v = min(Znegmax[1],Zposmax[1])) > -oo, v+eps, v) \\ ranges (approximately) intersect (max may be ∞, min: -∞)
|| (Znegmin[1] < Zposmin[1] && #Zpos && ZZneg) || (Znegmax[1] > Zposmax[1] && #Zpos && ZZpos), return(0)); \\ Z's ∪ ZZ's span everything (we take into account possibility of ∞)
if(#pos && #neg, , return(1)); \\ Since now we know that Z's do not span 0
\\ ≥1 of Zneg/Zpos may be empty (then ZZneg/ZZpos matter); if not, then “effectively” either ZZneg or ZZpos are “as if present”; 8 cases???
\\ This is too many cases; we reduce the complexity by doing suitable reflections:
\\ Ⓐ if ZZneg is “in the convex hull of Z’s”, multiply x’s by -1; then there are two cases where Zpos is empty: ⓐ with Zneg≠∅
\\ ⓑ with ZZpos present. in the first case multiply y’s by -1; in the second exchange x’s and y’s.
\\ (Only Z*max/min are updated now — the rest is postponed via xRefl, yRefl, asX/asY.)
if(#Zneg && !#Zpos, yRefl = -1; Zposmin = -Znegmax; Zposmax = -Znegmin; Znegmin = oo; Znegmax = -oo; Zneg=[]; Zpos=[oo]); \\ Reflect ↕ if only neg is present
if(ZZneg || (#Zneg && #Zpos && Znegmax[1] > Zposmax[1]), \\ Reflect ↔ if “ZZneg is in the convex hull
xRefl = -1; v = Zposmin; Zposmin = -Zposmax; Zposmax = -v; v = Znegmin; Znegmin = -Znegmax; Znegmax = -v; ZZpos = ZZneg; ZZneg = 0);
if(ZZpos && !(#Zneg || #Zpos), \\ Exchange x and y if only ZZpos is present (we exchange 1 and 2 in vecsort() below, etc.)
asX = 2; asY = 1; [xRefl, yRefl] = [yRefl, xRefl]; \\ Zposmin = [0, [0,1]]; Zposmax = [0, [0,1]]), \\ Now ⟪0, Z*⟫ intersects y>0 and does not contain (-1,0).
Zposmin=Zposmax=[0,[0,1]];
for(k=1, #pos, pos[k] /= pos[k][3]);
for(k=1, #neg, neg[k] /= neg[k][3]),
\\ else (have Zneg or Zpos)
if(#Zneg, Zposmax = [v = Znegmax[1] - Zposmin[1], [-v,-1]], ZZpos, Zposmax = [oo,[1,0]] ); \\ Now Zposmin is “the left boundary of ⟪0,Z's⟫”; find the right boundary (after “shearing” the left boundary to 0)
s = Zposmin[1]*xRefl*yRefl; \\ For doing shear before reflections
for(k=1, #pos, pos[k] = [pos[k][1]-s*(v=pos[k][2]), v]/pos[k][3]);
for(k=1, #neg, neg[k] = [neg[k][1]-s*(v=neg[k][2]), v]/neg[k][3]);
if(#Zneg, , if(Zposmax[1]<oo, Zposmax[1] -= Zposmin[1])); \\ ??? SHOULD NOT WE MERGE IT ABOVE???
);
\\ After-③: not in the hull ⇔ can separate projectivizations of neg and pos + ℝ₊·⟪Zposmin,Zposmax⟫ by a hyperplane
\\ ④ The orientation of this hyperplane is clear since ⟪Zposmin,Zposmax⟫ is not empty!
\\ my([xxRefl, yyRefl] = [xRefl, yRefl]);
\\ if(asX==2, [xxRefl, yyRefl] = [yyRefl, xxRefl]); \\ for arguments to convhull_above
if(dbg, writebin("pos.wrbin",pos); dbg_x(pos); POS[1]=pos;
print("asX="asX", pos="pos", neg="neg", xRefl="xRefl", yRefl="yRefl,if(#pos>1,[pos2[1],pos2[2]]=pos[1..2]); if(#neg>1,[neg2[1],neg2[2]]=neg[#neg-1..#neg];Str(", lst neg x diff="neg[#neg][1]-neg[#neg-1][1]," da="neg[#neg-1][1]-0.35," db="neg[#neg][1]-0.35), "")));
\\ ??? Need to repeat twice to avoid the bug with using subtraction for a sorting function
neg = vecsort(vecsort(neg,(x,y)->numcmp(x[asY],y[asY]),4*(yRefl==1)),(x,y)->numcmp(x[asX],y[asX]),8+4*(xRefl==-1));
\\ neg = vecsort(vecsort(neg,(x,y)->numcmp(x[asY],y[asY]),4*(yRefl==1)),(x,y)->numcmp(x[asX],y[asX]),8+4*(xRefl==-1));
pos = vecsort(vecsort(pos,(x,y)->numcmp(x[asY],y[asY]),4*(yRefl==-1)),(x,y)->numcmp(x[asX],y[asX]),8+4*(xRefl==-1)); \\ sort lexicographically in order [“asX”,-“asY”] and [“asX”,“asY”] and keep the first entry
\\ pos = vecsort(vecsort(pos,(x,y)->numcmp(x[asY],y[asY]),4*(yRefl==-1)),(x,y)->numcmp(x[asX],y[asX]),8+4*(xRefl==-1)); \\ sort lexicographically in order [“asX”,-“asY”] and [“asX”,“asY”] and keep the first entry
if(dbg, print("srt pos="pos", neg="neg", xRefl="xRefl", yRefl="yRefl,if(#neg>1,[neg2[1],neg2[2]]=neg[#neg-1..#neg];Str(", lst neg x diff="neg[#neg][1]-neg[#neg-1][1]," da="neg[#neg-1][1]-0.35," db="neg[#neg][1]-0.35), "")));
if(#neg > 1, neg = convhull_above(xRefl*vector(#neg,k,neg[k][asX]),yRefl*vector(#neg,k, neg[k][asY])), neg = [[xRefl*neg[1][asX], yRefl*neg[1][asY]]]);
if(#pos > 1, pos = convhull_above(xRefl*vector(#pos,k,pos[k][asX]),yRefl*vector(#pos,k,-pos[k][asY])), pos = [[xRefl*pos[1][asX], -yRefl*pos[1][asY]]]); \\ Negated!
if(Zposmax[1] == -oo, print(" !!! Error")); \\ ??? REMOVE THIS
if(dbg, print("Zposmax="Zposmax);
print("hull pos="pos", neg="neg));
\\ ??? With fuzz_convhull3D_avoids0(83,[3],1,2,5) on the last step pos[2][2] is twice as large as before???
\\ ??? Why eps below saves the problematic cases, and hurts others ??? fuzz_convhull3D_avoids0(892,[3],1,2,2) ⇒ small >0
if( Zposmax[1], \\ > eps || Zposmax[1] < -eps, \\ We have not yet taken into account the slope vv extending pos “to the right” to infinity (cannot trust slope of ultra-short in-clusters-gaps):
vv = if(Zposmax[1] < oo, -1/Zposmax[1], 0); \\ (Zposmax is 1/slope) the part-of-pos-at-∞ stretches from slope=∞ to slope = vv
v = 1; my(mx = pos[1][2]); \\ The part to the right of v is consumed by the hull with slope vv
for(k=2,#pos, if(pos[k][4]<=vv, break); \\ Take into account short interval-of-x only if y is increasing on them
if(pos[k][1]-pos[k-1][1]>e0, mx=pos[v=k][2], if(pos[k][4]>0, v=k); mx=max(mx,pos[k][2]))); \\ ??? WAS: && (pos[k][1]-pos[k-1][1])>e0
if(pos[v][1] < neg[#neg][1], pos = concat(pos[1..v], [[neg[#neg][1], mx + (neg[#neg][1] - pos[v][1])*vv]])));
if(dbg, print("addR pos="pos", neg="neg", v="v));
v=1; vv=neg[1][2]; for(k=2,#neg, if(neg[k][1]-neg[k-1][1]>=e0, break); v=k; vv=max(vv,neg[k][2]));
for(k=1, v, neg[k][2] = vv); \\ “Inside”, taking convex hulls handles “clumps of nearby x-coordinates”. Do this at start too.
v=0; vv=neg[#neg][2]; for(k=1,#neg-1, if(neg[#neg-k+1][1]-neg[#neg-k][1]>=e0, break, v=k; vv=max(vv,neg[#neg-k][2])));
for(k=0, v, neg[#neg-k][2] = vv); \\ And at end.
v=1; vv=pos[1][2]; for(k=2,#pos, if(pos[k][1]-pos[k-1][1]>=e0, break); v=k; vv=max(vv,pos[k][2]));
for(k=1, v, pos[k][2] = vv); \\ “Inside”, taking convex hulls handles “clumps of nearby x-coordinates”. Do this at start too.
if(dbg, print("cmprs pos="pos", neg="neg));
v=0; vv=pos[#pos][2]; for(k=1,#pos-1, if(pos[#pos-k+1][1]-pos[#pos-k][1]>=e0, break, v=k; vv=max(vv,pos[#pos-k][2])));
for(k=0, v, pos[#pos-k][2] = vv); \\ And at end.
if(dbg, print("cmprs pos="pos", neg="neg", delta1x="pos[1][1]-neg[1][1]", v="v", vv="vv));
if((v = convhulls_maxsum(vv = [pos,neg], e0))[1] >= -e0, if(dbg,print("v="v)); 0, \\ The hulls intersect
#v == 1, error("Bug: the x-ranges of pos or neg empty"), \\ Should not happen
1, if(dbg,print("v="v)); 1, \\ ??? The rest not implemented yet
#v, error("Non-intersecting cases not yet implemented ???") \\ NOT YET: x-coordinates do not overlap
\\ ???
); \\ ??? No normal yet
}
ckNO_convhull3D_avoids0(n,m,dbg,N=2^31,f=1.) = { \\ 0 on success
my(o, randV(D,lim=2,pos=0,N)=vector(D,k,lim*(if(pos,0,-1/2/(1+(N==2)))+random(N)/N)));
my(seed = f*vector(n,nn,randV(3,,,N)), mix = f*randV(n,5,1,N), cap=-sum(nn=1,n,seed[nn]*mix[nn]), rest = f*vector(m,mm,randV(3,10,,N)), perm=numtoperm(n+m+1,random((n+m+1)!)));
if(n==0, cap=[0,0,0]); \\ empty sum would give 0
convhull3D_avoids0(vecextract(concat(seed,concat([cap],rest)),perm),dbg);
}
ckYES_convhull3D_avoids0(n,dbg,N=2^31,f=1.) = { \\ non-0 on success
my(NNN=if(N==-2,N,max(N,10000)), randV(D,lim=2,pos=0,N)=vector(D,k,lim*(if(pos,0,-1/2/(1+(N==2)))+random(N)/N)), v = f*randV(3,,,NNN), vv);
while(v3_sprod(v,v)<0.1-1e-31, v = f*randV(3,,,NNN)); \\ N=2 cannot generate vv-samples with 10000
convhull3D_avoids0(vector(n,nn, vv=f*[0,0,0]; while(v3_sprod(v,vv)<=1e-28, vv = f*randV(3,,,N)); vv),dbg);
}
fuzz_convhull3D_avoids0(cnt=100,Ns=concat([2..18],[2^31]),dbg0,pre_=4,post_=8,NN=12,e_cnt) = {
if(type(cnt)==type([]), foreach(cnt,c, fuzz_convhull3D_avoids0(1,Ns,dbg0,pre_,post_,NN,c-1)); return());
my(dbg=bitand(dbg0, 0x0F)); \\ use 0x16 to debug aborts
for(k=if(dbg0,cnt,1)+e_cnt,cnt+e_cnt,
if(k%1000, , print1(if(k%10000,".","%")));
foreach(Ns, N,
if(dbg!=1,
for(nn=if(dbg,NN,0),NN, foreach([1,1.], f, setrand(k); if(!ckYES_convhull3D_avoids0(nn,dbg, N, f), print("failYES rand="k" NN="nn",\tunits=1/"N", f="f), bitand(dbg0,0x10), print("OK rand="k" NN="nn",\tunits=1/"N", f="f)))));
if(dbg!=2,
for(pre=if(dbg,pre_,0),pre_, for(post=if(dbg,post_,0), post_, foreach([1,1.], f,
setrand(k); if(ckNO_convhull3D_avoids0(pre,post,dbg,N,f), print("failNO rand="k" pre="pre" post="post",\tunits=1/"N", f="f), bitand(dbg0,0x10), print("OK rand="k" pre="pre" post="post",\tunits=1/"N", f="f))))))));
}
ck__(fn,args,expc) = my(rc=call(eval(fn),args)); if(expc!=rc, print("!!! "fn"( "args" ) -> " rc ". Expecting: " expc));
/* ck__("convhull3D_avoids0", [[[1,2,3],[-2,-4,-6]]], 0)
for(pre=0,1, for(post=0, 1, for(k=1,10, setrand(k); if(ckNO_convhull3D_avoids0(pre,post), print("fail rand="k" pre="pre" post="post)))))
*/
/* ; 0≤eps≪1 applied heuristically only???
\\ n = listpop(~neg); p = listpop(~pos); \\ Move the hyperplane=line ⟪ℝ˟n,ℝ˟p⟫ to become infinity; p to go to (0,+oo), n to (-oo,0)
n = neg[1]; p = pos[1]; \\ Move the hyperplane=line ⟪ℝ˟n,ℝ˟p⟫ to become infinity; p to go to (0,+oo), n to (-oo,0)
n /= sqrt(v3_sprod(n,n)); p /= sqrt(v3_sprod(p,p));
if(perp = v3_prod(n, p), , return(1)); \\ proportional, in opposite ½-spaces
my(perp_nrm, pos_mintan, pos_maxtan, t);
perp /= (perp_nrm = sqrt(v3_sprod(perp,perp))); \\ perp × “•” rotates “•” 90° in the plane ⟪p,n⟫
n_perp = v3_prod(perp,n);
for(k=1,#pos, if(abs(v3_sprod(perp,pos[k]))<eps,
if( abs(t = v3_sprod(pos[k],n_perp)) < eps, return(1));
t = v3_sprod(pos[k],n)/t; \\ (co???)tangent of angle with n
listput(~on_pln, pos[k]))); \\ Actually, we only need the convex hull of these vectors???
*/
v3_projdir(v,p,Cmp=v3_prod([0,0,1],p))={\\ expresses how close is the direction of the projection of v along p to the projection of rel
my(rotPr=v3_prod(v,p)); v3_sprod(rotPr,Cmp)/sqrt(v3_abs2(rotPr)); \\ the 3rd argument should be v3_prod(rel,p) with the default rel=[0,0,1]
} \\ larger is closer; assumes all v's are on the same side of the plane <p,rel>
highest_ang_indices(v1,v2,vs,s=1,e=#vs,rel=[0,0,1])={\\ returns a set of indices of vs which lead to the highest angle of slope of the plane w.r.t. rel
my(pr=v1-v2, ang, Cmp=v3_prod(rel,pr), ang0=-oo, inds);
for(k=s,e, ang = v3_projdir(vs[k]-v1,pr,Cmp);
if(ang>=ang0,
if(ang>ang0,
inds=List([k]); ang0=ang,
listput(inds,k))
)
); inds;
}
v2_in_triang(p,p1,p2,p3,sn=v2_prod(p3-p1,p2-p1),chkN=1, chkOp=1)={ \\ returns 0 if p is in △ p1 p2 p3; assumes p,p3 are on the same side of (p1 p2)
if(chkN && sn*v2_prod(p3-p1,p-p1)<0,return(-1)); \\ Did not reach v yet
if(chkOp && sn*v2_prod(p3-p2,p2-p)<0,return(1)); \\ beyond the edge (p2 p3)
0; \\ found
}
v3_in_hat_of_convhull_ind(vs,cont=0,minI=1,maxI=#vs, omit=1, v=sum(k=minI,maxI-omit,vs[k]/(maxI-minI)),do1=1)={ \\ assumes 2D projections of vs are vertices of a convex poly P (in the cyclic order), v is in P.
\\ Returns true if v is in a top face of the convex hull which (unless cont) contains vN as a vertex. Assumes v is on the required side of v1→vs[ind[1]]. (and of ???)
while(do1, do1 = 0; \\ if cont, finds the actual face over v
my(rel = vector(maxI-minI,k, vs[minI-1+k]-vs[maxI]), VV=vs[minI]-vs[maxI-1],
Rel = vector(#rel,k, rel[k]/v2_prod(rel[k],VV)), \\ secs = Rel × rel[1] grows from 0 to 1
secs = vector(#rel,k, v2_prod(Rel[k],rel[1]) ), hull=convhull_above(secs,vector(#Rel,k,Rel[k][3])), \\ conic hull above
sn = v2_prod(vs[minI+1]-vs[maxI],vs[minI]-vs[maxI]));
for(k=2, #hull, my(pl=v2_in_triang(v,vs[maxI], vs[minI-1+hull[k-1][3]], vs[minI-1+hull[k][3]],
sn, k<#hull, hull[k][3]>1+hull[k-1][3]));
if(pl>0, if(!cont, return(0),
do1=1; minI = minI-1+hull[k-1][3]; maxI = minI-1+hull[k][3]; break),
!pl, return([minI-1+hull[k-1][3], minI-1+hull[k][3],maxI])));
if(!do1, error("Panic311"), /* print("Redo: ",[minI,maxI, hull]) */ )
);
}
v3_extrapolate_ind(ps,i,i1,i2,i3)={ \\ find the height of the plane through the points ps[iK] over the point ps[i]
my(d2=ps[i2]-ps[i1], d3=ps[i3]-ps[i1], d=ps[i]-ps[i1], pr=v2_prod(d2,d3));
v2_prod(d,d3)/pr*d2[3] + v2_prod(d2,d)/pr*d3[3] + ps[i1][3];
}
v3_coeff_ind(ps,i1,i2,i3)={ \\ find the slopes of the plane through the points ps[iK] (returns a “normalized” normal vector)
my(nrml=v3_prod(ps[i2]-ps[i1], ps[i3]-ps[i1]));
-nrml/nrml[3];
}
linlog_est_from_below(xs,ys,x,ytarget,cont=0)={ \\ xs positive, len≥3, without limit of leeway
my(len=1+#xs, vs=vector(len,k,if(k==len,[x,log(x),-ytarget],[xs[k],log(xs[k]),-ys[k]])),
res = v3_in_hat_of_convhull_ind(vs,cont)); \\ if !cont, true if the expected value at x is ≥ ytarget
if(cont && res[3]<=#xs,
if(cont>0, return([-v3_extrapolate_ind(vs,#vs,res[1],res[2],res[3]), res]),
return([-v3_coeff_ind(vs,res[1],res[2],res[3]), res])));
["OK",res];
} \\ on data x^C on [1..4] vs 10 underestimates (OK: 67 of 100) for C=2, and overestimates by 7% (√11.2 instead of √10) for C=½
powlog_pow_est_from_above(xs,ys)={ \\ |xs| with many > 1, ys ≠ 0, len≥3. Find est
my(inds = [k|k<-[1..#xs], abs(xs[k])>1 && ys[k]!=0]);
if(#inds < 3, return([]));
xs = [log(abs(xs[inds[k]])) | k<-[1..#inds]];
ys = [log(abs(ys[inds[k]])) | k<-[1..#inds]];
my(len=#inds, vs=vector(len,k,[xs[k],log(xs[k]),ys[k]]),
res = v3_in_hat_of_convhull_ind(vs,1,,,0)); \\ find the exact face over the c.o.m. of all
v3_coeff_ind(vs,res[1],res[2],res[3]);
} \\ 0.47 for √(x+1), 1.88 for (x+1)² on [1..4]
powlog_est_from_below(xs,ys,x,ytarget,cont=0,withPow=1)={ \\ xs positive, len≥3, without limit of leeway
my(Pow0 = if(withPow,powlog_pow_est_from_above(xs,ys), []));
Pow0 = if(#Pow0, /* print("Pow="Pow0[1]); */ Pow0[1], 1); \\ may give up on estimating the power
my(len=1+#xs, vs=vector(len,k,if(k==len,[x^Pow0,log(x),-ytarget],[xs[k]^Pow0,log(xs[k]),-ys[k]])),
res = v3_in_hat_of_convhull_ind(vs,cont)); \\ if !cont, true if the expected value at x is ≥ ytarget
if(cont && res[3]<=#xs,
if(cont>0, return([-v3_extrapolate_ind(vs,#vs,res[1],res[2],res[3]), res]),
return([-v3_coeff_ind(vs,res[1],res[2],res[3]), res])));
["OK",res];
} \\ for (x+1)² on [1..4] vs 10 overestimates 0.3% (121.4 vs. 121), and underestimates (OK) by 1/2000 for √(x+1)
\\ Badness: 3: the most probable estimate is <ytarget at x (“not good enough”)
\\ 2: ¬3, but can make bad with quality remaining ≤leeway
\\ 1: ¬2, ¬3, but can make bad with quality remaining ≤best_possible + leeway_extra
\\ 0: every fit which is ≤leeway or ≤best_possible + leeway_extra is good enough
\\ In “1”, “2”, there is a good “chance” that the actual data is actually “good at x”; ⇒ no incentive to skip 2x and go to 4x etc.
\\ With “3”, may want to calculate the trend to find how many subdivisions we may skip (Is not needed with
lw_linlog_est_from_below(xs,ys,x,ytarget,leeway,leeway_extra)={ \\ xs positive, len≥3, without limit of leeway
my(len=1+#xs, vs=vector(len,k,if(k==len,[x,log(x),-ytarget],[xs[k],log(xs[k]),-ys[k]])),
v=sum(k=1,#xs,vs[k])/#xs, VV, res = v3_in_hat_of_convhull_ind(vs,1,,,,v), \\ ¬abort early; find face over v
bad = if(res[3]<=#xs, 3, \\ do not even need to use leeway \\ interpolate at v:
my(cf = -v3_coeff_ind(vs,res[1],res[2],res[3]), V = res[1][3]+cf[1]*(v[1]-res[1][1])+cf[2]*(v[2]-res[1][2]));
if(V>v[3]-leeway, 2, \\ Decreasing “quality” by leeway, can make bad
res = v3_in_hat_of_convhull_ind(vs,1,,#xs,,v)); \\ ignore vs[#vs]; ??? should restrict to the hat-of-vs[#vs] only???
cf = -v3_coeff_ind(vs,res[1],res[2],res[3]); VV = res[1][3]+cf[1]*(v[1]-res[1][1])+cf[2]*(v[2]-res[1][2]);
if(V>VV-leeway_extra, 1, 0) \\ VV = the best fit to data at v, decreasing quality by leeway_extra can make bad
)); \\ Somewhat similar to “conditional probability”: an unprobable event “VV” already happened; what extra is needed to get “V”?
\\ TBC……………
[bad, res]; \\ returning cf makes no sense: it is needed with bad==3, but is calculated only with bad<3
}
/*
It seems that the criterion below interprets “badly approximated by the given law” as an evidence of “low probability” of WHAT???
So if the “shifted center of mass” below is inside the convex hull, one should essentially shift it more to the boundary of
the convex hull. If inside face, then the plane of this face gives “the best” main part matching the data, so the estimate
at the extrapolation point should be given by the intersection with this plane.
==============
Assume that we know that a function φ has a form ε(t)+ο(t) (with ε an “expected” form, and ο≥0 a “random overshoot”), and that ο(t)
is distributed by the exponential law with the distance of decay D(t).
Suppose that we know that ε∈⟨f₁,…,fₙ⟩ and that we know φ(tₖ)≕φₖ; in addition to φ, apply the same convention to ε, ο and D. Given
this, how is φ(T) distributed? Let vₖ≔(f₁(tₖ),…,fₙ(tₖ)) and write ε≕c₁f₁+…+cₙfₙ and c≔(c₁,…,cₙ). Then we now that ⟨c,vₖ⟩ ≤ φₖ, and
the probability of a particular value of ⟨c,vₖ⟩ is proportional to exp(D(tₖ)(⟨c,vₖ⟩-φₖ). Hence the total probability is proportional
to exp ∑ₖD(tₖ)(⟨c,vₖ⟩-φₖ) ≕ exp(⟨c,V⟩-Φ), with V≔∑ₖD(tₖ)vₖ (and Φ≔∑ₖD(tₖ)φₖ — but we do not need this).
CONCLUSION: we got a polyhedral region U in ℝⁿ with c satisfying ⟨c,vₖ⟩ ≤ φₖ, and the exponential measure with density proportional
to exp⟨c,V⟩. If we want to estimate ε(T) ≝ ⟨c,w⟩ with w≔(f₁(T),…,fₙ(T)), we need to intersect U with hyperplanes ⟨c,w⟩=const and
inspect how the function ⟨c,V⟩ looks on the intersections.
In the simplest cases U is unbounded as well as its intersections with hyperplanes ⟨c,w⟩=const. Consider the point P∈U where ⟨c,V⟩
takes its maximal value M on U. Call the hyperplane ⟨c,w⟩=const passing through P as Π; it is clear (???) that the corresponding
value W of ⟨c,w⟩ is the most probable one.
Consider the half-space where ⟨c,w⟩ is smaller than W.
In the simplest cases there is no vertex of U inside this half-space. Then all the sections of U by hyperplanes ⟨c,w⟩=const are
just translations of each other in the direction d of a certain edge of P; hence their relative probability is described by α⟨d,V⟩
when the translation is by αd. Hence to find the translation where the probability decreases exp -N times, it is enough to find
the minimum of ⟨c,w⟩ on {⟨c,V⟩=M-N}∩U.
In fact, this reduces to two problems in dimension n-1: first, find m≔max⟨c,V⟩ under the condition that ⟨c,w⟩ is the target value at
T. Then check whether intersection of U with ⟨c,V⟩ = m - N is empty. Dualizing the first problem, we consider lines passing below
the given collection of points Qₖ, and want the intersection with the given vertical line ℓ to be as high as possible (so need to
intersect ℓ with the hull of Qₖ). For the second problem, to see whether the given point R is inside the convex hull of other
points Qₖ, draw projectivizations of P-Qₖ (with the red/blue color indicating the direction w.r.t. a given plane hyperplane Π);
check whether one can separate red points from blue ones by a hyperplane. (One may assume that one blue and one red point are at
infinity.)
Although this is a cheating, use this approach as a proxy in the general case as well. ???
Actually, if we know the intersections, it is probably easy to calculate the honest integrals of exp.
============= Older
Suppose that we know that ε∈⟨f,g,h⟩ and φ(tₖ)≕φₖ, and we
want to understand whether it is very probable that ε(T)≥ε° (with probability ≥1-δ and a given δ≪1).
We are going to use a proxy for probability: the largest probability density of a possible combination of ο(tₖ).
By a linear transformation of f,g,h we may assume that f(T)=1, g(T)=0, and h(T) has a zero of higher order than for g(T). Below,
s(t) ≔ h(t)/g(t); so s(T)=0. Define τ as s⁻¹; in other words, τ(s(t)) ≡ t. Let F(s) ≔ f∘τ, G(s) ≔ g∘τ. Then s·G(s) = h∘τ. Let
sₖ ≔ s(tₖ).
Seems that we assume that functions f, g, h do not change sign at the range of arguments we are interested in.
============= Estimation for ε(T)=ε₀
Now E≔ε∘τ∈⟨F,G,s·G⟩ (hence E(0)=ε₀), and E=ε₀F+ℓG for a certain linear affine function ℓ. So let Φ(s) ≔ (φ∘τ - ε₀F)/G. We see that
Φ-ℓ≥0, and at s its values are distributed exponentially with the distance of decay D(t(s))/G(s).
Hence the density of probability at the given point is exp -∑ₖ (φₖ - ε₀f(tₖ) - ℓ(sₖ)g(tₖ))/D(tₖ) if all the summands are
non-negative, and 0 otherwise. Writing ℓ(s) = αs + β we see that the dependence on α and β is linear, so the expression above
may be written as exp -(X+Y ℓ(S₀)+Zε₀) for a certain S₀. Hence if we want high enough density, we want X+Y ℓ(S₀)+Zε₀ ≥ Δ₀ (with
Δ₀ depending on the threshold δ), as well as ℓ(sₖ) ≤ (φₖ - ε₀f(tₖ))/g(tₖ).
This (and the condition ε(T)≥ε₀) defines a polytope on the space (ε₀,α,β). ???
ℓ(sₖ) + ε₀σₖ ≤ φₖ/g(tₖ) with σₖ ≔ f(tₖ)/g(tₖ)
ℓ(S₀) + ε₀Σ₀ ≥ (Δ₀-X)/Y with Σ₀ = Z/Y (assuming Y>0)
So: want a plane separating ??? points (sₖ,σₖ) and c(S₀,Σ₀) (with a suitable c ≔ φₖY/(g(tₖ)(Δ₀-X)) ???) with a bound on ε₀. In
other words, given a collection S of points in ℝ³ and a point P with an angle with a vertex at P, find whether there is a plane
passing through P (“inside” this angle) such that S is on one side of this plane.
Given a subset S₀ of a sphere and an arc on the square, find a big circle intersecting the arc and not separating points of S.
In other words, the convex hull S₁ of S₀ and its opposite -S₁ should not contain the whole arc.
May it be that (as for a fixed ε₀) this is a problem of finding an intersection of a ray (from a center of mass???) to the
boundary of the convex hull???
================ Joining the estimate in the extrapolation point
If we want to find whether the value in the extrapolation point can be less than something, this just adds another point where
the extrapolation should pass below it. (So the only difference is that its weight in the center of mass is 0.)
This reduces the problem to the linear programming problem: find the intersection of the convex hull of points with the given line.
However, if the given shift of the center of mass turns out to be in the non-extended hull, we need to move it to the boundary
of the smaller hull. Then we need to check whether it is inside the larger hull.
================ Convex hull in a convex-projection case
Assume that the points xₖ are on edges of (infinite) vertical convex prism (the preimage of vertices pₚ of a convex polygon P w.r.t.
the projection along a vector e₃ — directed up). Then the segments [xₖxₗ] which are preimages of the sides of P are edges of the
convex hull; such “lateral” edges split the boundary of the hull into the “top” and the “bottom” halves. The face on the “top side”
adjacent to a lateral edge [xₖxₗ] can be found by looking for the vertex (or vertices) xₘ such that the direction of their
projection along this edge is closest to the projection of e₂.
Or with the maximal ((xₘ-xₗ)×(xₖ-xₗ),e₃×(xₖ-xₗ))/|(xₘ-xₗ)×(xₖ-xₗ)|.
Or use (xₘ,-e₃×(xₖ-xₗ)×(xₖ-xₗ))/???. ???
Knowing the face adjacent to an edge, it splits the “top side” into several parts, and the argument above is applicable to all the
parts. Hence one can inductively finish finding the top side.
In applications, P is known in advance. Can it help with calculations???
================ Convex hull in an almost-convex-projection case
Suppose that in addition to xₖ as above, we add another point X about which it is known that it is X₀+te₃ with X₀ inside the
convex hull of {xₖ}, and t>0. Then in the calculation above one can additionally check whether the plane through xₖ¸ₗ¸ₘ passes
above X. If it does, then the found face “survives” the addition of X. The “non-surviving” faces form a part of the convex hull
of {xₖ} which is replaced by “a conic part” adjacent to X..
This conic part may have non-△ faces only if X is on the (continuation of) certain faces of the hull of {xₖ}.
Intersections of a vertical segment Σ on the boundary of the prism (or an “extended” convex prism) with planes between X and the
hull of {xₖ} are determined by the non-surviving faces of the hull of {xₖ} and the “new” faces including X. ???
So it seems we are ony interested in the intersection of the hull with the plane passing through X and Σ???
================ Convexity
In applications, the points Pₖ are proportional to vectors (1, x, log x) with x>0. Due to convexity of log, these are on a boundary
of a convex cone.
At infinity, the hull of this cone has a straight part connecting (0:0:-1) and (0:1:0). Near (0:0:-1) the cone is exponentially
close to the direction (ε:0:1) (so has a dihedral angle). Near (0:1:0) the cone is t/log t-close to the opposite direction to
the connecting line. (So it is not smooth, but does not have a dihedral angle there!)
Applying a projective transformation sending 0 to -∞·e₃ (±???) reduces the problem to one above.
OLD_v2_in_convhull_ind(v,vs,i1,iN,ind,sn)={ \\ assumes vs are vertices of a convex poly (in the cyclic order), v is inside. i1,ind,iN ↗
my(); \\ Returns 0 if v is inside vs[i1,ind,iN]. Assumes v is on the required side of vs[i1,ind[1]].
if(sn*v2_prod(vs[iN]-vs[ind[#ind]],v-vs[i1])<0,return([-2])); \\ Assumes sign(sn) == sign(v2_prod(vs[3]-vs[1],vs[2]-vs[1]))
\\ \\ if(ind[1]-i1>1 && sn*v2_prod(v-vs[i1],vs[ind[1]]-vs[i1])<0,return(0));
for(k=1,#ind-1,
if(ind[k+1]-ind[k] > 1 && sn*v2_prod(v-vs[ind[k]],vs[ind[k+1]]-vs[ind[k]])<0,return([k]));
);
0; \\ success
}
v2_in_convhull_ind(v,vs,v1,vN,i1,ind,sn)={ \\ assumes vs are vertices of a convex poly (in the cyclic order), v is inside. index(v1),ind,index(vN) ↗
\\ Returns 0 if v is inside v1, vs[ind], vN. Assumes v is on the required side of v1→vs[ind[1]].
\\ print(v"---"ind); print(v1" ",[vs[k]|k<-ind]""vN); \\ Assumes sign(sn) == sign(v2_prod(vs[ind[1]]-v1,vs[ind[2]]-v1))
if(sn*v2_prod(vs[ind[#ind]]-vN,v-vN)>0,return(-1)); \\ Did not reach v yet
if(ind[1]-i1 > 1 && sn*v2_prod(v-v1,vs[ind[1]]-v1)>0,return(1)); \\ beyond the first edge
for(k=1,#ind-1,
if(ind[k+1]-ind[k] > 1 && sn*v2_prod(v-vs[ind[k]],vs[ind[k+1]]-vs[ind[k]])>0,return(k+1)); \\ beyond this edge
);
0; \\ found
}
v3_in_hat_of_convhull_ind(vs,v=sum(k=1,#vs-1,vs[k])/(#vs-1))={ \\ assumes 2D projections of vs are vertices of a convex poly P (in the cyclic order), v is in P.
\\ Returns true if v is in a top face of the convex hull which contains vN as a vertex. Assumes v is on the required side of v1→vs[ind[1]].
my(edge=1,e=#vs-1,vN=vs[#vs],inds,IN, sn=v2_prod(vs[2]-vs[1],vs[3]-vs[1]));
while(edge<e, inds = highest_ang_indices(vs[edge],vN,vs,edge+1,e);
if(0==(IN=v2_in_convhull_ind(v,vs,vs[edge],vN,edge,inds,sn)), return([edge,inds])); \\ found
if(IN > 0, return(0)); \\ “beyond” the face.
edge = inds[#inds]);
return(0);
}
if(res[3]<=#xs, \\ not good enough
if(!"estim", return([-v3_extrapolate_ind(vs,#vs,res[1],res[2],res[3]), res]), \\ DONOT: return [estimate, face]
my(cf = -v3_coeff_ind(vs,res[1],res[2],res[3]));
return([cf[1],cf[2],res[1][3]-cf[1]*res[1][1]-cf[2]*res[1][2], res])), \\ BUT: [plane-coeffs, face]
my(V=-v3_extrapolate_ind(vs,#vs,res[1],res[2],res[3])); \\ a candidate; check leeways
);
["OK",res];
*/