35 # pragma warning (disable: 4701 4127)
43 : maxit2_(maxit1_ +
Math::digits() + 10)
54 , tolb_(tol0_ * tol2_)
55 , xthresh_(1000 * tol2_)
60 , _ep2(_e2 /
Math::sq(_f1))
77 , _etol2(
real(0.1) * tol2_ /
104 ar = 2 * (cosx - sinx) * (cosx + sinx),
105 y0 =
n & 1 ? *--
c : 0,
y1 = 0;
114 ? 2 * sinx * cosx *
y0
119 unsigned caps)
const {
124 bool arcmode,
real s12_a12,
unsigned outmask,
132 GenPosition(arcmode, s12_a12, outmask,
133 lat2, lon2, azi2, s12, m12, M12, M21, S12);
137 bool arcmode,
real s12_a12,
138 unsigned caps)
const {
145 return GeodesicLine(*
this, lat1, lon1, azi1, salp1, calp1,
146 caps, arcmode, s12_a12);
150 unsigned caps)
const {
155 real a12,
unsigned caps)
const {
160 unsigned outmask,
real& s12,
170 int lonsign = lon12 >= 0 ? 1 : -1;
188 int swapp =
abs(lat1) <
abs(lat2) ? -1 : 1;
194 int latsign = lat1 < 0 ? 1 : -1;
209 real sbet1, cbet1, sbet2, cbet2, s12x, m12x;
228 if (cbet1 < -sbet1) {
230 sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
232 if (
abs(sbet2) == -sbet1)
244 bool meridian = lat1 == -90 || slam12 == 0;
251 calp1 = clam12; salp1 = slam12;
252 calp2 = 1; salp2 = 0;
256 ssig1 = sbet1, csig1 = calp1 * cbet1,
257 ssig2 = sbet2, csig2 = calp2 * cbet2;
260 sig12 =
atan2(
max(
real(0), csig1 * ssig2 - ssig1 * csig2),
261 csig1 * csig2 + ssig1 * ssig2);
264 Lengths(
_n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
266 s12x, m12x,
dummy, M12, M21, Ca);
275 if (sig12 < 1 || m12x >= 0) {
277 if (sig12 < 3 *
tiny_)
278 sig12 = m12x = s12x = 0;
288 real omg12 = 0, somg12 = 2, comg12 = 0;
291 (_f <= 0 || lon12s >=
_f * 180)) {
294 calp1 = calp2 = 0; salp1 = salp2 = 1;
296 sig12 = omg12 = lam12 /
_f1;
297 m12x =
_b *
sin(sig12);
299 M12 = M21 =
cos(sig12);
302 }
else if (!meridian) {
309 sig12 =
InverseStart(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
310 lam12, slam12, clam12,
311 salp1, calp1, salp2, calp2, dnm,
316 s12x = sig12 *
_b * dnm;
319 M12 = M21 =
cos(sig12 / dnm);
321 omg12 = lam12 / (
_f1 * dnm);
337 real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0;
341 for (
bool tripn =
false, tripb =
false;
347 real v =
Lambda12(sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
349 salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
350 eps, domg12, numit <
maxit1_, dv, Ca);
352 if (tripb || !(
abs(
v) >= (tripn ? 8 : 1) *
tol0_))
break;
354 if (
v > 0 && (numit >
maxit1_ || calp1/salp1 > calp1b/salp1b))
355 { salp1b = salp1; calp1b = calp1; }
356 else if (
v < 0 && (numit >
maxit1_ || calp1/salp1 < calp1a/salp1a))
357 { salp1a = salp1; calp1a = calp1; }
358 if (numit < maxit1_ && dv > 0) {
362 sdalp1 =
sin(dalp1), cdalp1 =
cos(dalp1),
363 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
365 calp1 = calp1 * cdalp1 - salp1 * sdalp1;
383 salp1 = (salp1a + salp1b)/2;
384 calp1 = (calp1a + calp1b)/2;
387 tripb = (
abs(salp1a - salp1) + (calp1a - calp1) <
tolb_ ||
388 abs(salp1 - salp1b) + (calp1 - calp1b) <
tolb_);
394 unsigned lengthmask = outmask |
396 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
397 cbet1, cbet2, lengthmask, s12x, m12x,
dummy, M12, M21, Ca);
402 if (outmask &
AREA) {
404 real sdomg12 =
sin(domg12), cdomg12 =
cos(domg12);
405 somg12 = slam12 * cdomg12 - clam12 * sdomg12;
406 comg12 = clam12 * cdomg12 + slam12 * sdomg12;
417 if (outmask &
AREA) {
420 salp0 = salp1 * cbet1,
423 if (calp0 != 0 && salp0 != 0) {
426 ssig1 = sbet1, csig1 = calp1 * cbet1,
427 ssig2 = sbet2, csig2 = calp2 * cbet2,
429 eps = k2 / (2 * (1 +
sqrt(1 + k2)) + k2),
438 S12 =
A4 * (B42 - B41);
443 if (!meridian && somg12 > 1) {
444 somg12 =
sin(omg12); comg12 =
cos(omg12);
449 comg12 > -
real(0.7071) &&
450 sbet2 - sbet1 <
real(1.75)) {
454 real domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
455 alp12 = 2 *
atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
456 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
460 salp12 = salp2 * calp1 - calp2 * salp1,
461 calp12 = calp2 * calp1 + salp2 * salp1;
466 if (salp12 == 0 && calp12 < 0) {
467 salp12 =
tiny_ * calp1;
470 alp12 =
atan2(salp12, calp12);
473 S12 *= swapp * lonsign * latsign;
486 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
487 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
499 real salp1, calp1, salp2, calp2,
501 outmask, s12, salp1, calp1, salp2, calp2,
512 unsigned caps)
const {
513 real t, salp1, calp1, salp2, calp2,
516 0u,
t, salp1, calp1, salp2, calp2,
522 GeodesicLine(*
this, lat1, lon1, azi1, salp1, calp1, caps,
true, a12);
528 real cbet1,
real cbet2,
unsigned outmask,
541 real m0x = 0, J12 = 0,
A1 = 0,
A2 = 0;
558 s12b =
A1 * (sig12 + B1);
562 J12 = m0x * sig12 + (
A1 * B1 -
A2 * B2);
566 for (
int l = 1;
l <=
nC2_; ++
l)
567 Cb[
l] =
A1 * Ca[
l] -
A2 * Cb[
l];
576 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
580 real csig12 = csig1 * csig2 + ssig1 * ssig2;
581 real t =
_ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
582 M12 = csig12 + (
t * ssig2 - csig2 * J12) * ssig1 / dn1;
583 M21 = csig12 - (
t * ssig1 - csig1 * J12) * ssig2 / dn2;
595 if ( !(
q == 0 && r <= 0) ) {
604 disc =
S * (
S + 2 *
r3);
615 u +=
T + (
T != 0 ?
r2 /
T : 0);
621 u += 2 * r *
cos(ang / 3);
626 uv = u < 0 ?
q / (
v - u) : u +
v,
627 w = (uv -
q) / (2 *
v);
655 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
656 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
657 #if defined(__GNUC__) && __GNUC__ == 4 && \
658 (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
672 real sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
674 bool shortline = cbet12 >= 0 && sbet12 <
real(0.5) &&
675 cbet2 * lam12 <
real(0.5);
681 sbetm2 /= sbetm2 +
Math::sq(cbet1 + cbet2);
683 real omg12 = lam12 / (
_f1 * dnm);
684 somg12 =
sin(omg12); comg12 =
cos(omg12);
686 somg12 = slam12; comg12 = clam12;
689 salp1 = cbet2 * somg12;
690 calp1 = comg12 >= 0 ?
691 sbet12 + cbet2 * sbet1 *
Math::sq(somg12) / (1 + comg12) :
692 sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
696 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
698 if (shortline && ssig12 <
_etol2) {
700 salp2 = cbet1 * somg12;
701 calp2 = sbet12 - cbet1 * sbet2 *
702 (comg12 >= 0 ?
Math::sq(somg12) / (1 + comg12) : 1 - comg12);
705 sig12 =
atan2(ssig12, csig12);
713 real y, lamscale, betscale;
724 eps = k2 / (2 * (1 +
sqrt(1 + k2)) + k2);
727 betscale = lamscale * cbet1;
729 x = lam12x / lamscale;
730 y = sbet12a / betscale;
734 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
735 bet12a =
atan2(sbet12a, cbet12a);
740 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
744 betscale =
x < -
real(0.01) ? sbet12a /
x :
746 lamscale = betscale / cbet1;
747 y = lam12x / lamscale;
796 omg12a = lamscale * (
_f >= 0 ? -
x * k/(1 + k) : -
y * (1 + k)/k );
797 somg12 =
sin(omg12a); comg12 = -
cos(omg12a);
799 salp1 = cbet2 * somg12;
800 calp1 = sbet12a - cbet2 * sbet1 *
Math::sq(somg12) / (1 - comg12);
807 salp1 = 1; calp1 = 0;
821 bool diffp,
real& dlam12,
825 if (sbet1 == 0 && calp1 == 0)
832 salp0 = salp1 * cbet1,
835 real somg1, comg1, somg2, comg2, somg12, comg12, lam12;
838 ssig1 = sbet1; somg1 = salp0 * sbet1;
839 csig1 = comg1 = calp1 * cbet1;
847 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
852 calp2 = cbet2 != cbet1 ||
abs(sbet2) != -sbet1 ?
855 (cbet2 - cbet1) * (cbet1 + cbet2) :
856 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
860 ssig2 = sbet2; somg2 = salp0 * sbet2;
861 csig2 = comg2 = calp2 * cbet2;
866 sig12 =
atan2(
max(
real(0), csig1 * ssig2 - ssig1 * csig2),
867 csig1 * csig2 + ssig1 * ssig2);
870 somg12 =
max(
real(0), comg1 * somg2 - somg1 * comg2);
871 comg12 = comg1 * comg2 + somg1 * somg2;
873 real eta =
atan2(somg12 * clam120 - comg12 * slam120,
874 comg12 * clam120 + somg12 * slam120);
877 eps = k2 / (2 * (1 +
sqrt(1 + k2)) + k2);
881 domg12 = -
_f *
A3f(eps) * salp0 * (sig12 + B312);
882 lam12 = eta + domg12;
886 dlam12 = - 2 *
_f1 * dn1 / sbet1;
889 Lengths(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
892 dlam12 *=
_f1 / (calp2 * cbet2);
909 for (
int l = 1;
l <
nC3_; ++
l) {
923 for (
int l = 0;
l <
nC4_; ++
l) {
958 #if GEOGRAPHICLIB_GEODESIC_ORDER/2 == 1
959 static const real coeff[] = {
963 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 2
964 static const real coeff[] = {
968 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 3
969 static const real coeff[] = {
973 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 4
974 static const real coeff[] = {
976 25, 64, 256, 4096, 0, 16384,
979 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
981 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
nA1_/2 + 2,
982 "Coefficient array size mismatch in A1m1f");
985 return (
t + eps) / (1 - eps);
991 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3
992 static const real coeff[] = {
1000 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4
1001 static const real coeff[] = {
1011 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5
1012 static const real coeff[] = {
1024 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6
1025 static const real coeff[] = {
1039 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7
1040 static const real coeff[] = {
1042 19, -64, 384, -1024, 2048,
1056 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8
1057 static const real coeff[] = {
1059 19, -64, 384, -1024, 2048,
1061 7, -18, 128, -256, 4096,
1065 -11, 96, -160, 16384,
1076 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
1078 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
1080 "Coefficient array size mismatch in C1f");
1085 for (
int l = 1;
l <=
nC1_; ++
l) {
1097 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3
1098 static const real coeff[] = {
1106 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4
1107 static const real coeff[] = {
1117 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5
1118 static const real coeff[] = {
1120 205, -432, 768, 1536,
1130 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6
1131 static const real coeff[] = {
1133 205, -432, 768, 1536,
1135 4005, -4736, 3840, 12288,
1145 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7
1146 static const real coeff[] = {
1148 -4879, 9840, -20736, 36864, 73728,
1150 4005, -4736, 3840, 12288,
1152 8703, -7200, 3712, 12288,
1156 -141115, 41604, 92160,
1162 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8
1163 static const real coeff[] = {
1165 -4879, 9840, -20736, 36864, 73728,
1167 -86171, 120150, -142080, 115200, 368640,
1169 8703, -7200, 3712, 12288,
1171 1082857, -688608, 258720, 737280,
1173 -141115, 41604, 92160,
1175 -2200311, 533134, 860160,
1179 109167851, 82575360,
1182 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
1184 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
1186 "Coefficient array size mismatch in C1pf");
1203 #if GEOGRAPHICLIB_GEODESIC_ORDER/2 == 1
1204 static const real coeff[] = {
1208 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 2
1209 static const real coeff[] = {
1213 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 3
1214 static const real coeff[] = {
1216 -11, -28, -192, 0, 256,
1218 #elif GEOGRAPHICLIB_GEODESIC_ORDER/2 == 4
1219 static const real coeff[] = {
1221 -375, -704, -1792, -12288, 0, 16384,
1224 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
1226 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
nA2_/2 + 2,
1227 "Coefficient array size mismatch in A2m1f");
1230 return (
t - eps) / (1 + eps);
1236 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3
1237 static const real coeff[] = {
1245 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4
1246 static const real coeff[] = {
1256 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5
1257 static const real coeff[] = {
1269 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6
1270 static const real coeff[] = {
1284 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7
1285 static const real coeff[] = {
1287 41, 64, 128, 1024, 2048,
1301 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8
1302 static const real coeff[] = {
1304 41, 64, 128, 1024, 2048,
1306 47, 70, 128, 768, 4096,
1310 133, 224, 1120, 16384,
1321 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
1323 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
1325 "Coefficient array size mismatch in C2f");
1330 for (
int l = 1;
l <=
nC2_; ++
l) {
1342 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3
1343 static const real coeff[] = {
1351 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4
1352 static const real coeff[] = {
1362 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5
1363 static const real coeff[] = {
1375 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6
1376 static const real coeff[] = {
1390 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7
1391 static const real coeff[] = {
1407 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8
1408 static const real coeff[] = {
1416 -5, -20, -4, -6, 128,
1427 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
1429 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
1431 "Coefficient array size mismatch in A3f");
1433 for (
int j =
nA3_ - 1;
j >= 0; --
j) {
1444 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3
1445 static const real coeff[] = {
1453 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4
1454 static const real coeff[] = {
1470 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5
1471 static const real coeff[] = {
1493 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6
1494 static const real coeff[] = {
1526 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7
1527 static const real coeff[] = {
1571 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8
1572 static const real coeff[] = {
1606 10, -6, -10, 9, 384,
1616 -7, 20, -28, 14, 1024,
1631 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
1633 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
1635 "Coefficient array size mismatch in C3coeff");
1637 for (
int l = 1;
l <
nC3_; ++
l) {
1638 for (
int j =
nC3_ - 1;
j >=
l; --
j) {
1649 #if GEOGRAPHICLIB_GEODESIC_ORDER == 3
1650 static const real coeff[] = {
1664 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 4
1665 static const real coeff[] = {
1673 4, 24, -84, 210, 315,
1687 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 5
1688 static const real coeff[] = {
1694 1088, -352, -66, 3465,
1696 48, -352, 528, -231, 1155,
1698 16, 44, 264, -924, 2310, 3465,
1704 -896, 704, -198, 10395,
1706 -48, 352, -528, 231, 10395,
1712 320, -352, 132, 17325,
1720 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 6
1721 static const real coeff[] = {
1727 -224, -4784, 1573, 45045,
1729 -10656, 14144, -4576, -858, 45045,
1731 64, 624, -4576, 6864, -3003, 15015,
1733 100, 208, 572, 3432, -12012, 30030, 45045,
1739 5792, 1040, -1287, 135135,
1741 5952, -11648, 9152, -2574, 135135,
1743 -64, -624, 4576, -6864, 3003, 135135,
1749 -8448, 4992, -1144, 225225,
1751 -1440, 4160, -4576, 1716, 225225,
1757 3584, -3328, 1144, 315315,
1765 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 7
1766 static const real coeff[] = {
1772 -4480, 1088, 156, 45045,
1774 10736, -224, -4784, 1573, 45045,
1776 1664, -10656, 14144, -4576, -858, 45045,
1778 16, 64, 624, -4576, 6864, -3003, 15015,
1780 56, 100, 208, 572, 3432, -12012, 30030, 45045,
1786 3840, -2944, 468, 135135,
1788 -10704, 5792, 1040, -1287, 135135,
1790 -768, 5952, -11648, 9152, -2574, 135135,
1792 -16, -64, -624, 4576, -6864, 3003, 135135,
1798 1664, 1856, -936, 225225,
1800 6784, -8448, 4992, -1144, 225225,
1802 128, -1440, 4160, -4576, 1716, 225225,
1808 -2048, 1024, -208, 105105,
1810 -1792, 3584, -3328, 1144, 315315,
1816 3072, -2560, 832, 405405,
1824 #elif GEOGRAPHICLIB_GEODESIC_ORDER == 8
1825 static const real coeff[] = {
1831 20960, -7888, 4947, 765765,
1833 12480, -76160, 18496, 2652, 765765,
1835 -154048, 182512, -3808, -81328, 26741, 765765,
1837 3232, 28288, -181152, 240448, -77792, -14586, 765765,
1839 96, 272, 1088, 10608, -77792, 116688, -51051, 255255,
1841 588, 952, 1700, 3536, 9724, 58344, -204204, 510510, 765765,
1847 -39840, 1904, 255, 2297295,
1849 52608, 65280, -50048, 7956, 2297295,
1851 103744, -181968, 98464, 17680, -21879, 2297295,
1853 -1344, -13056, 101184, -198016, 155584, -43758, 2297295,
1855 -96, -272, -1088, -10608, 77792, -116688, 51051, 2297295,
1859 -928, -612, 3828825,
1861 64256, -28288, 2856, 3828825,
1863 -126528, 28288, 31552, -15912, 3828825,
1865 -41472, 115328, -143616, 84864, -19448, 3828825,
1867 160, 2176, -24480, 70720, -77792, 29172, 3828825,
1871 -16384, 1088, 5360355,
1873 -2560, 30464, -11560, 5360355,
1875 35840, -34816, 17408, -3536, 1786785,
1877 7168, -30464, 60928, -56576, 19448, 5360355,
1881 26624, -8704, 6891885,
1883 -77824, 34816, -6528, 6891885,
1885 -32256, 52224, -43520, 14144, 6891885,
1889 24576, -4352, 8423415,
1891 45056, -34816, 10880, 8423415,
1895 -28672, 8704, 9954945,
1900 #error "Bad value for GEOGRAPHICLIB_GEODESIC_ORDER"
1902 GEOGRAPHICLIB_STATIC_ASSERT(
sizeof(coeff) /
sizeof(
real) ==
1904 "Coefficient array size mismatch in C4coeff");
1906 for (
int l = 0;
l <
nC4_; ++
l) {
1907 for (
int j =
nC4_ - 1;
j >=
l; --
j) {