Skip to content

use stackalloc in RobustPredicates workspace arrays #57

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
uuiitwp1 wants to merge 1 commit into
base: master
Choose a base branch
from
+63 −115
Open

use stackalloc in RobustPredicates workspace arrays #57

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
177 changes: 62 additions & 115 deletions src/Triangle/RobustPredicates.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,38 +28,18 @@ public class RobustPredicates : IPredicates
{
#region Default predicates instance (Singleton)

private static readonly object creationLock = new object();
private static RobustPredicates _default;

/// <summary>
/// Gets the default configuration instance.
/// </summary>
public static RobustPredicates Default
{
get
{
if (_default == null)
{
lock (creationLock)
{
if (_default == null)
{
_default = new RobustPredicates();
}
}
}

return _default;
}
}
public static RobustPredicates Default { get; } = new RobustPredicates();

#endregion

#region Static initialization

private static double epsilon, splitter, resulterrbound;
private static double ccwerrboundA, ccwerrboundB, ccwerrboundC;
private static double iccerrboundA, iccerrboundB, iccerrboundC;
private static readonly double epsilon, splitter, resulterrbound;
private static readonly double ccwerrboundA, ccwerrboundB, ccwerrboundC;
private static readonly double iccerrboundA, iccerrboundB, iccerrboundC;
//private static double o3derrboundA, o3derrboundB, o3derrboundC;

/// <summary>
Expand Down Expand Up @@ -126,7 +106,6 @@ static RobustPredicates()
/// </summary>
public RobustPredicates()
{
AllocateWorkspace();
}

/// <inheritdoc/>
Expand Down Expand Up @@ -422,7 +401,7 @@ public Point FindCircumcenter(Point org, Point dest, Point apex,
/// property. (That is, if e is strongly nonoverlapping, h will be also.) Does NOT
/// maintain the nonoverlapping or nonadjacent properties.
/// </remarks>
private int FastExpansionSumZeroElim(int elen, double[] e, int flen, double[] f, double[] h)
private unsafe int FastExpansionSumZeroElim(int elen, double* e, int flen, double* f, double* h)
{
double Q;
double Qnew;
Expand Down Expand Up @@ -548,7 +527,7 @@ private int FastExpansionSumZeroElim(int elen, double[] e, int flen, double[] f,
/// maintains the strongly nonoverlapping and nonadjacent properties as well. (That is,
/// if e has one of these properties, so will h.)
/// </remarks>
private int ScaleExpansionZeroElim(int elen, double[] e, double b, double[] h)
private unsafe int ScaleExpansionZeroElim(int elen, double* e, double b, double* h)
{
double Q, sum;
double hh;
Expand Down Expand Up @@ -598,7 +577,7 @@ private int ScaleExpansionZeroElim(int elen, double[] e, double b, double[] h)
/// <param name="elen"></param>
/// <param name="e"></param>
/// <returns></returns>
private double Estimate(int elen, double[] e)
private unsafe double Estimate(int elen, double* e)
{
double Q;
int eindex;
Expand Down Expand Up @@ -629,7 +608,7 @@ private double Estimate(int elen, double[] e)
/// the returned value has the correct sign. Hence, this function is usually quite fast,
/// but will run more slowly when the input points are collinear or nearly so.
/// </remarks>
private double CounterClockwiseAdapt(Point pa, Point pb, Point pc, double detsum)
private unsafe double CounterClockwiseAdapt(Point pa, Point pb, Point pc, double detsum)
{
double acx, acy, bcx, bcy;
double acxtail, acytail, bcxtail, bcytail;
Expand All @@ -638,8 +617,8 @@ private double CounterClockwiseAdapt(Point pa, Point pb, Point pc, double detsum
double det, errbound;
// Edited to work around index out of range exceptions (changed array length from 4 to 5).
// See unsafe indexing in FastExpansionSumZeroElim.
double[] B = new double[5], u = new double[5];
double[] C1 = new double[8], C2 = new double[12], D = new double[16];
double* B = stackalloc double[5], u = stackalloc double[5];
double* C1 = stackalloc double[8], C2 = stackalloc double[12], D = stackalloc double[16];
double B3;
int C1length, C2length, Dlength;

Expand Down Expand Up @@ -734,51 +713,82 @@ private double CounterClockwiseAdapt(Point pa, Point pb, Point pc, double detsum
/// the returned value has the correct sign. Hence, this function is usually quite fast,
/// but will run more slowly when the input points are cocircular or nearly so.
/// </remarks>
private double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double permanent)
private unsafe double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double permanent)
{
double* fin1 = stackalloc double[1152];
double* fin2 = stackalloc double[1152];
double* abdet = stackalloc double[64];

double* axbc = stackalloc double[8];
double* axxbc = stackalloc double[16];
double* aybc = stackalloc double[8];
double* ayybc = stackalloc double[16];
double* adet = stackalloc double[32];

double* bxca = stackalloc double[8];
double* bxxca = stackalloc double[16];
double* byca = stackalloc double[8];
double* byyca = stackalloc double[16];
double* bdet = stackalloc double[32];

double* cxab = stackalloc double[8];
double* cxxab = stackalloc double[16];
double* cyab = stackalloc double[8];
double* cyyab = stackalloc double[16];
double* cdet = stackalloc double[32];

double* temp8 = stackalloc double[8];
double* temp16a = stackalloc double[16];
double* temp16b = stackalloc double[16];
double* temp16c = stackalloc double[16];

double* temp32a = stackalloc double[32];
double* temp32b = stackalloc double[32];
double* temp48 = stackalloc double[48];
double* temp64 = stackalloc double[64];
double adx, bdx, cdx, ady, bdy, cdy;
double det, errbound;

double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
double[] bc = new double[4], ca = new double[4], ab = new double[4];
double* bc = stackalloc double[4], ca = stackalloc double[4], ab = stackalloc double[4];
double bc3, ca3, ab3;
int axbclen, axxbclen, aybclen, ayybclen, alen;
int bxcalen, bxxcalen, bycalen, byycalen, blen;
int cxablen, cxxablen, cyablen, cyyablen, clen;
int ablen;
double[] finnow, finother, finswap;
double* finnow, finother, finswap;
int finlength;

double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
double[] aa = new double[4], bb = new double[4], cc = new double[4];
double* aa = stackalloc double[4], bb = stackalloc double[4], cc = stackalloc double[4];
double aa3, bb3, cc3;
double ti1, tj1;
double ti0, tj0;
// Edited to work around index out of range exceptions (changed array length from 4 to 5).
// See unsafe indexing in FastExpansionSumZeroElim.
double[] u = new double[5], v = new double[5];
double* u = stackalloc double[5], v = stackalloc double[5];
double u3, v3;
int temp8len, temp16alen, temp16blen, temp16clen;
int temp32alen, temp32blen, temp48len, temp64len;
double[] axtbb = new double[8], axtcc = new double[8], aytbb = new double[8], aytcc = new double[8];
double* axtbb = stackalloc double[8], axtcc = stackalloc double[8], aytbb = stackalloc double[8], aytcc = stackalloc double[8];
int axtbblen, axtcclen, aytbblen, aytcclen;
double[] bxtaa = new double[8], bxtcc = new double[8], bytaa = new double[8], bytcc = new double[8];
double* bxtaa = stackalloc double[8], bxtcc = stackalloc double[8], bytaa = stackalloc double[8], bytcc = stackalloc double[8];
int bxtaalen, bxtcclen, bytaalen, bytcclen;
double[] cxtaa = new double[8], cxtbb = new double[8], cytaa = new double[8], cytbb = new double[8];
double* cxtaa = stackalloc double[8], cxtbb = stackalloc double[8], cytaa = stackalloc double[8], cytbb = stackalloc double[8];
int cxtaalen, cxtbblen, cytaalen, cytbblen;
double[] axtbc = new double[8], aytbc = new double[8], bxtca = new double[8], bytca = new double[8], cxtab = new double[8], cytab = new double[8];
double* axtbc = stackalloc double[8], aytbc = stackalloc double[8], bxtca = stackalloc double[8], bytca = stackalloc double[8], cxtab = stackalloc double[8], cytab = stackalloc double[8];
int axtbclen = 0, aytbclen = 0, bxtcalen = 0, bytcalen = 0, cxtablen = 0, cytablen = 0;
double[] axtbct = new double[16], aytbct = new double[16], bxtcat = new double[16], bytcat = new double[16], cxtabt = new double[16], cytabt = new double[16];
double* axtbct = stackalloc double[16], aytbct = stackalloc double[16], bxtcat = stackalloc double[16], bytcat = stackalloc double[16], cxtabt = stackalloc double[16], cytabt = stackalloc double[16];
int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
double[] axtbctt = new double[8], aytbctt = new double[8], bxtcatt = new double[8];
double[] bytcatt = new double[8], cxtabtt = new double[8], cytabtt = new double[8];
double* axtbctt = stackalloc double[8], aytbctt = stackalloc double[8], bxtcatt = stackalloc double[8];
double* bytcatt = stackalloc double[8], cxtabtt = stackalloc double[8], cytabtt = stackalloc double[8];
int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
double[] abt = new double[8], bct = new double[8], cat = new double[8];
double* abt = stackalloc double[8], bct = stackalloc double[8], cat = stackalloc double[8];
int abtlen, bctlen, catlen;
double[] abtt = new double[4], bctt = new double[4], catt = new double[4];
double* abtt = stackalloc double[4], bctt = stackalloc double[4], catt = stackalloc double[4];
int abttlen, bcttlen, cattlen;
double abtt3, bctt3, catt3;
double negate;
Expand All @@ -799,13 +809,6 @@ private double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double perm
bdy = (double)(pb.y - pd.y);
cdy = (double)(pc.y - pd.y);

adx = (double)(pa.x - pd.x);
bdx = (double)(pb.x - pd.x);
cdx = (double)(pc.x - pd.x);
ady = (double)(pa.y - pd.y);
bdy = (double)(pb.y - pd.y);
cdy = (double)(pc.y - pd.y);

bdxcdy1 = (double)(bdx * cdy); c = (double)(splitter * bdx); abig = (double)(c - bdx); ahi = c - abig; alo = bdx - ahi; c = (double)(splitter * cdy); abig = (double)(c - cdy); bhi = c - abig; blo = cdy - bhi; err1 = bdxcdy1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); bdxcdy0 = (alo * blo) - err3;
cdxbdy1 = (double)(cdx * bdy); c = (double)(splitter * cdx); abig = (double)(c - cdx); ahi = c - abig; alo = cdx - ahi; c = (double)(splitter * bdy); abig = (double)(c - bdy); bhi = c - abig; blo = bdy - bhi; err1 = cdxbdy1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); cdxbdy0 = (alo * blo) - err3;
_i = (double)(bdxcdy0 - cdxbdy0); bvirt = (double)(bdxcdy0 - _i); avirt = _i + bvirt; bround = bvirt - cdxbdy0; around = bdxcdy0 - avirt; bc[0] = around + bround; _j = (double)(bdxcdy1 + _i); bvirt = (double)(_j - bdxcdy1); avirt = _j - bvirt; bround = _i - bvirt; around = bdxcdy1 - avirt; _0 = around + bround; _i = (double)(_0 - cdxbdy1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - cdxbdy1; around = _0 - avirt; bc[1] = around + bround; bc3 = (double)(_j + _i); bvirt = (double)(bc3 - _j); avirt = bc3 - bvirt; bround = _i - bvirt; around = _j - avirt; bc[2] = around + bround;
Expand Down Expand Up @@ -858,13 +861,13 @@ private double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double perm
return det;
}

errbound = iccerrboundC * permanent + resulterrbound * ((det) >= 0.0 ? (det) : -(det));
det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail))
+ 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
+ ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail))
+ 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
+ ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail))
+ 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
errbound = (iccerrboundC * permanent) + (resulterrbound * (det >= 0.0 ? det : -det));
det += (((adx * adx) + (ady * ady)) * ((bdx * cdytail) + (cdy * bdxtail) - ((bdy * cdxtail) + (cdx * bdytail))))
+ (2.0 * ((adx * adxtail) + (ady * adytail)) * ((bdx * cdy) - (bdy * cdx)))
+ ((((bdx * bdx) + (bdy * bdy)) * ((cdx * adytail) + (ady * cdxtail) - ((cdy * adxtail) + (adx * cdytail))))
+ (2.0 * ((bdx * bdxtail) + (bdy * bdytail)) * ((cdx * ady) - (cdy * adx))))
+ ((((cdx * cdx) + (cdy * cdy)) * ((adx * bdytail) + (bdy * adxtail) - ((ady * bdxtail) + (bdx * adytail))))
+ (2.0 * ((cdx * cdxtail) + (cdy * cdytail)) * ((adx * bdy) - (ady * bdx))));
if ((det >= errbound) || (-det >= errbound))
{
return det;
Expand Down Expand Up @@ -1064,7 +1067,6 @@ private double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double perm
finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother);
finswap = finnow; finnow = finother; finother = finswap;


temp32alen = ScaleExpansionZeroElim(aytbctlen, aytbct, adytail, temp32a);
aytbcttlen = ScaleExpansionZeroElim(bcttlen, bctt, adytail, aytbctt);
temp16alen = ScaleExpansionZeroElim(aytbcttlen, aytbctt, 2.0 * ady, temp16a);
Expand Down Expand Up @@ -1229,7 +1231,6 @@ private double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double perm
finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother);
finswap = finnow; finnow = finother; finother = finswap;


temp32alen = ScaleExpansionZeroElim(cytabtlen, cytabt, cdytail, temp32a);
cytabttlen = ScaleExpansionZeroElim(abttlen, abtt, cdytail, cytabtt);
temp16alen = ScaleExpansionZeroElim(cytabttlen, cytabtt, 2.0 * cdy, temp16a);
Expand All @@ -1243,60 +1244,6 @@ private double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double perm

return finnow[finlength - 1];
}

#region Workspace

// InCircleAdapt workspace:
double[] fin1, fin2, abdet;

double[] axbc, axxbc, aybc, ayybc, adet;
double[] bxca, bxxca, byca, byyca, bdet;
double[] cxab, cxxab, cyab, cyyab, cdet;

double[] temp8, temp16a, temp16b, temp16c;
double[] temp32a, temp32b, temp48, temp64;

private void AllocateWorkspace()
{
fin1 = new double[1152];
fin2 = new double[1152];
abdet = new double[64];

axbc = new double[8];
axxbc = new double[16];
aybc = new double[8];
ayybc = new double[16];
adet = new double[32];

bxca = new double[8];
bxxca = new double[16];
byca = new double[8];
byyca = new double[16];
bdet = new double[32];

cxab = new double[8];
cxxab = new double[16];
cyab = new double[8];
cyyab = new double[16];
cdet = new double[32];

temp8 = new double[8];
temp16a = new double[16];
temp16b = new double[16];
temp16c = new double[16];

temp32a = new double[32];
temp32b = new double[32];
temp48 = new double[48];
temp64 = new double[64];
}

private void ClearWorkspace()
{
}

#endregion

#endregion
}
}
1 change: 1 addition & 0 deletions src/Triangle/Triangle.csproj
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
<AssemblyVersion>1.0.0</AssemblyVersion>
<Version>1.0.0-beta5</Version>
<LangVersion>default</LangVersion>
<AllowUnsafeBlocks>True</AllowUnsafeBlocks>
</PropertyGroup>

<ItemGroup>
Expand Down
Loading