experimental/Intersection/CubicParameterizationCode.cpp - platform_external_skia - Gitiles

 #include <vector>

 /* Given:
  * Resultant[a*t^3 + b*t^2 + c*t + d - x, e*t^3 + f*t^2 + g*t + h - y, t]
  */

 const char result[] =
 "-d^3 e^3 + c d^2 e^2 f - b d^2 e f^2 + a d^2 f^3 - c^2 d e^2 g + "
 " 2 b d^2 e^2 g + b c d e f g - 3 a d^2 e f g - a c d f^2 g - "
 " b^2 d e g^2 + 2 a c d e g^2 + a b d f g^2 - a^2 d g^3 + c^3 e^2 h - "
 " 3 b c d e^2 h + 3 a d^2 e^2 h - b c^2 e f h + 2 b^2 d e f h + "
 " a c d e f h + a c^2 f^2 h - 2 a b d f^2 h + b^2 c e g h - "
 " 2 a c^2 e g h - a b d e g h - a b c f g h + 3 a^2 d f g h + "
 " a^2 c g^2 h - b^3 e h^2 + 3 a b c e h^2 - 3 a^2 d e h^2 + "
 " a b^2 f h^2 - 2 a^2 c f h^2 - a^2 b g h^2 + a^3 h^3 + 3 d^2 e^3 x - "
 " 2 c d e^2 f x + 2 b d e f^2 x - 2 a d f^3 x + c^2 e^2 g x - "
 " 4 b d e^2 g x - b c e f g x + 6 a d e f g x + a c f^2 g x + "
 " b^2 e g^2 x - 2 a c e g^2 x - a b f g^2 x + a^2 g^3 x + "
 " 3 b c e^2 h x - 6 a d e^2 h x - 2 b^2 e f h x - a c e f h x + "
 " 2 a b f^2 h x + a b e g h x - 3 a^2 f g h x + 3 a^2 e h^2 x - "
 " 3 d e^3 x^2 + c e^2 f x^2 - b e f^2 x^2 + a f^3 x^2 + "
 " 2 b e^2 g x^2 - 3 a e f g x^2 + 3 a e^2 h x^2 + e^3 x^3 - "
 " c^3 e^2 y + 3 b c d e^2 y - 3 a d^2 e^2 y + b c^2 e f y - "
 " 2 b^2 d e f y - a c d e f y - a c^2 f^2 y + 2 a b d f^2 y - "
 " b^2 c e g y + 2 a c^2 e g y + a b d e g y + a b c f g y - "
 " 3 a^2 d f g y - a^2 c g^2 y + 2 b^3 e h y - 6 a b c e h y + "
 " 6 a^2 d e h y - 2 a b^2 f h y + 4 a^2 c f h y + 2 a^2 b g h y - "
 " 3 a^3 h^2 y - 3 b c e^2 x y + 6 a d e^2 x y + 2 b^2 e f x y + "
 " a c e f x y - 2 a b f^2 x y - a b e g x y + 3 a^2 f g x y - "
 " 6 a^2 e h x y - 3 a e^2 x^2 y - b^3 e y^2 + 3 a b c e y^2 - "
 " 3 a^2 d e y^2 + a b^2 f y^2 - 2 a^2 c f y^2 - a^2 b g y^2 + "
 " 3 a^3 h y^2 + 3 a^2 e x y^2 - a^3 y^3";

 const size_t len = sizeof(result) - 1;
 const int factors = 8;

 struct coeff {
     int s; // constant and coefficient sign
     int n[factors]; // 0 or power of a (1, 2, or 3) for a through h
 };

 enum {
     xxx_coeff,
     xxy_coeff,
     xyy_coeff,
     yyy_coeff,
     xx_coeff,
     xy_coeff,
     yy_coeff,
     x_coeff,
     y_coeff,
     c_coeff,
     coeff_count
 };

 typedef std::vector<coeff> coeffs;
 typedef std::vector<coeffs> n_coeffs;

 static char skipSpace(size_t& index) {
     do {
         ++index;
     } while (result[index] == ' ');
     return result[index];
 }

 static char backSkipSpace(size_t& end) {
     while (result[end - 1] == ' ') {
         --end;
     }
     return result[end - 1];
 }

 static void match(coeffs& co, const char pattern[]) {
     size_t patternLen = strlen(pattern);
     size_t index = 0;
     while (index < len) {
         char ch = result[index];
         if (ch != '-' && ch != '+') {
             printf("missing sign\n");
         }
         size_t end = index + 1;
         while (result[end] != '+' && result[end] != '-' && ++end < len) {
             ;
         }
         backSkipSpace(end);
         size_t idx = index;
         index = end;
         skipSpace(index);
         if (!strncmp(&result[end - patternLen], pattern, patternLen) == 0) {
             continue;
         }
         size_t endCoeff = end - patternLen;
         char last = backSkipSpace(endCoeff);
         if (last == '2' || last == '3') {
             last = result[endCoeff - 3]; // skip ^2
         }
         if (last == 'x' || last == 'y') {
             continue;
         }
         coeff c;
         c.s = result[idx] == '-' ? -1 : 1;
         bzero(c.n, sizeof(c.n));
         ch = skipSpace(idx);
         if (ch >= '2' && ch <= '6') {
             c.s *= ch - '0';
             ch = skipSpace(idx);
         }
         while (idx < endCoeff) {
             char x = result[idx];
             if (x < 'a' || x > 'a' + factors) {
                 printf("expected factor\n");
             }
             idx++;
             int pow = 1;
             if (result[idx] == '^') {
                 idx++;
                 char exp = result[idx];
                 if (exp < '2' || exp > '3') {
                     printf("expected exponent\n");
                 }
                 pow = exp - '0';
             }
             skipSpace(idx);
             c.n[x - 'a'] = pow;
         }
         co.push_back(c);
     }
 }

 void cubecode_test();

 void cubecode_test() {
     n_coeffs c(coeff_count);
     match(c[xxx_coeff], "x^3");   // 1 factor
     match(c[xxy_coeff], "x^2 y"); // 1 factor
     match(c[xyy_coeff], "x y^2"); // 1 factor
     match(c[yyy_coeff], "y^3");   // 1 factor
     match(c[xx_coeff], "x^2");    // 7 factors
     match(c[xy_coeff], "x y");    // 8 factors
     match(c[yy_coeff], "y^2");    // 7 factors
     match(c[x_coeff], "x");       // 21 factors
     match(c[y_coeff], "y");       // 21 factors
     match(c[c_coeff], "");        // 34 factors
 #define COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS 0
 #define WRITE_AS_NONOPTIMIZED_C_CODE 0
 #if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
     int count[factors][factors][factors];
     bzero(count, sizeof(count));
 #endif
 #if WRITE_AS_NONOPTIMIZED_C_CODE
     printf("// start of generated code");
 #endif
     for (n_coeffs::iterator it = c.begin(); it < c.end(); ++it) {
         coeffs& co = *it;
 #if WRITE_AS_NONOPTIMIZED_C_CODE
         printf("\nstatic double calc_%c(double a, double b, double c, double d,"
                "\n                     double e, double f, double g, double h) {"
                "\n    return"
                "\n ", 'A' + (it - c.begin()));
         if (co[0].s > 0) {
             printf(" ");
         }
         if (abs(co[0].s) == 1) {
             printf("    ");
         }
 #endif
         for (coeffs::iterator ct = co.begin(); ct < co.end(); ++ct) {
             const coeff& cf = *ct;
 #if WRITE_AS_NONOPTIMIZED_C_CODE
             printf("        ");
             bool firstFactor = false;
             if (ct - co.begin() > 0 || cf.s < 0) {
                 printf("%c", cf.s < 0 ? '-' : '+');
             }
             if (ct - co.begin() > 0) {
                 printf(" ");
             }
             if (abs(cf.s) > 1) {
                 printf("%d * ", abs(cf.s));
             } else {
                 if (ct - co.begin() > 0) {
                     printf("    ");
                 }
             }
 #endif
             for (int x = 0; x < factors; ++x) {
                 if (cf.n[x] == 0) {
                     continue;
                 }
 #if WRITE_AS_NONOPTIMIZED_C_CODE
                 for (int y = 0 ; y < cf.n[x]; ++y) {
                     if (y > 0 || firstFactor) {
                         printf(" * ");
                     }
                     printf("%c", 'a' + x);
                 }
                 firstFactor = true;
 #endif
 #if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
                 for (int y = x; y < factors; ++y) {
                     if (cf.n[y] == 0) {
                         continue;
                     }
                     if (x == y && cf.n[y] == 1) {
                         continue;
                     }
                     for (int z = y; z < factors; ++z) {
                         if (cf.n[z] == 0) {
                             continue;
                         }
                         if ((x == z || y == z) && cf.n[z] == 1) {
                             continue;
                         }
                         if (x == y && y == z && cf.n[z] == 2) {
                             continue;
                         }
                         count[x][y][z]++;
                     }
                 }
 #endif
             }
 #if WRITE_AS_NONOPTIMIZED_C_CODE
             if (ct + 1 < co.end()) {
                 printf("\n");
             }
 #endif
         }
 #if WRITE_AS_NONOPTIMIZED_C_CODE
             printf(";\n}\n");
 #endif
     }
 #if WRITE_AS_NONOPTIMIZED_C_CODE
     printf("// end of generated code\n");
 #endif
 #if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
     const int bestCount = 20;
     int best[bestCount][4];
     bzero(best, sizeof(best));
     for (int x = 0; x < factors; ++x) {
         for (int y = x; y < factors; ++y) {
             for (int z = y; z < factors; ++z) {
                 if (!count[x][y][z]) {
                     continue;
                 }
                 for (int w = 0; w < bestCount; ++w) {
                     if (best[w][0] < count[x][y][z]) {
                         best[w][0] = count[x][y][z];
                         best[w][1] = x;
                         best[w][2] = y;
                         best[w][3] = z;
                         break;
                     }
                 }
             }
         }
     }
     for (int w = 0; w < bestCount; ++w) {
         printf("%c%c%c=%d\n", 'a' + best[w][1], 'a'  + best[w][2],
             'a' + best[w][3], best[w][0]);
     }
 #endif
 #if WRITE_AS_NONOPTIMIZED_C_CODE
     printf("\n");
 #endif
 }

 /* results: variable triplets used 10 or more times:
 aah=14
 ade=14
 aeh=14
 dee=14
 bce=13
 beg=13
 beh=12
 bbe=11
 bef=11
 cee=11
 cef=11
 def=11
 ceh=10
 deg=10
 */
	#include <vector>

	/* Given:
	* Resultant[at^3 + bt^2 + ct + d - x, et^3 + ft^2 + gt + h - y, t]
	*/

	const char result[] =
	"-d^3 e^3 + c d^2 e^2 f - b d^2 e f^2 + a d^2 f^3 - c^2 d e^2 g + "
	" 2 b d^2 e^2 g + b c d e f g - 3 a d^2 e f g - a c d f^2 g - "
	" b^2 d e g^2 + 2 a c d e g^2 + a b d f g^2 - a^2 d g^3 + c^3 e^2 h - "
	" 3 b c d e^2 h + 3 a d^2 e^2 h - b c^2 e f h + 2 b^2 d e f h + "
	" a c d e f h + a c^2 f^2 h - 2 a b d f^2 h + b^2 c e g h - "
	" 2 a c^2 e g h - a b d e g h - a b c f g h + 3 a^2 d f g h + "
	" a^2 c g^2 h - b^3 e h^2 + 3 a b c e h^2 - 3 a^2 d e h^2 + "
	" a b^2 f h^2 - 2 a^2 c f h^2 - a^2 b g h^2 + a^3 h^3 + 3 d^2 e^3 x - "
	" 2 c d e^2 f x + 2 b d e f^2 x - 2 a d f^3 x + c^2 e^2 g x - "
	" 4 b d e^2 g x - b c e f g x + 6 a d e f g x + a c f^2 g x + "
	" b^2 e g^2 x - 2 a c e g^2 x - a b f g^2 x + a^2 g^3 x + "
	" 3 b c e^2 h x - 6 a d e^2 h x - 2 b^2 e f h x - a c e f h x + "
	" 2 a b f^2 h x + a b e g h x - 3 a^2 f g h x + 3 a^2 e h^2 x - "
	" 3 d e^3 x^2 + c e^2 f x^2 - b e f^2 x^2 + a f^3 x^2 + "
	" 2 b e^2 g x^2 - 3 a e f g x^2 + 3 a e^2 h x^2 + e^3 x^3 - "
	" c^3 e^2 y + 3 b c d e^2 y - 3 a d^2 e^2 y + b c^2 e f y - "
	" 2 b^2 d e f y - a c d e f y - a c^2 f^2 y + 2 a b d f^2 y - "
	" b^2 c e g y + 2 a c^2 e g y + a b d e g y + a b c f g y - "
	" 3 a^2 d f g y - a^2 c g^2 y + 2 b^3 e h y - 6 a b c e h y + "
	" 6 a^2 d e h y - 2 a b^2 f h y + 4 a^2 c f h y + 2 a^2 b g h y - "
	" 3 a^3 h^2 y - 3 b c e^2 x y + 6 a d e^2 x y + 2 b^2 e f x y + "
	" a c e f x y - 2 a b f^2 x y - a b e g x y + 3 a^2 f g x y - "
	" 6 a^2 e h x y - 3 a e^2 x^2 y - b^3 e y^2 + 3 a b c e y^2 - "
	" 3 a^2 d e y^2 + a b^2 f y^2 - 2 a^2 c f y^2 - a^2 b g y^2 + "
	" 3 a^3 h y^2 + 3 a^2 e x y^2 - a^3 y^3";

	const size_t len = sizeof(result) - 1;
	const int factors = 8;

	struct coeff {
	int s; // constant and coefficient sign
	int n[factors]; // 0 or power of a (1, 2, or 3) for a through h
	};

	enum {
	xxx_coeff,
	xxy_coeff,
	xyy_coeff,
	yyy_coeff,
	xx_coeff,
	xy_coeff,
	yy_coeff,
	x_coeff,
	y_coeff,
	c_coeff,
	coeff_count
	};

	typedef std::vector<coeff> coeffs;
	typedef std::vector<coeffs> n_coeffs;

	static char skipSpace(size_t& index) {
	do {
	++index;
	} while (result[index] == ' ');
	return result[index];
	}

	static char backSkipSpace(size_t& end) {
	while (result[end - 1] == ' ') {
	--end;
	}
	return result[end - 1];
	}

	static void match(coeffs& co, const char pattern[]) {
	size_t patternLen = strlen(pattern);
	size_t index = 0;
	while (index < len) {
	char ch = result[index];
	if (ch != '-' && ch != '+') {
	printf("missing sign\n");
	}
	size_t end = index + 1;
	while (result[end] != '+' && result[end] != '-' && ++end < len) {
	;
	}
	backSkipSpace(end);
	size_t idx = index;
	index = end;
	skipSpace(index);
	if (!strncmp(&result[end - patternLen], pattern, patternLen) == 0) {
	continue;
	}
	size_t endCoeff = end - patternLen;
	char last = backSkipSpace(endCoeff);
	if (last == '2' \|\| last == '3') {
	last = result[endCoeff - 3]; // skip ^2
	}
	if (last == 'x' \|\| last == 'y') {
	continue;
	}
	coeff c;
	c.s = result[idx] == '-' ? -1 : 1;
	bzero(c.n, sizeof(c.n));
	ch = skipSpace(idx);
	if (ch >= '2' && ch <= '6') {
	c.s *= ch - '0';
	ch = skipSpace(idx);
	}
	while (idx < endCoeff) {
	char x = result[idx];
	if (x < 'a' \|\| x > 'a' + factors) {
	printf("expected factor\n");
	}
	idx++;
	int pow = 1;
	if (result[idx] == '^') {
	idx++;
	char exp = result[idx];
	if (exp < '2' \|\| exp > '3') {
	printf("expected exponent\n");
	}
	pow = exp - '0';
	}
	skipSpace(idx);
	c.n[x - 'a'] = pow;
	}
	co.push_back(c);
	}
	}

	void cubecode_test();

	void cubecode_test() {
	n_coeffs c(coeff_count);
	match(c[xxx_coeff], "x^3"); // 1 factor
	match(c[xxy_coeff], "x^2 y"); // 1 factor
	match(c[xyy_coeff], "x y^2"); // 1 factor
	match(c[yyy_coeff], "y^3"); // 1 factor
	match(c[xx_coeff], "x^2"); // 7 factors
	match(c[xy_coeff], "x y"); // 8 factors
	match(c[yy_coeff], "y^2"); // 7 factors
	match(c[x_coeff], "x"); // 21 factors
	match(c[y_coeff], "y"); // 21 factors
	match(c[c_coeff], ""); // 34 factors
	#define COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS 0
	#define WRITE_AS_NONOPTIMIZED_C_CODE 0
	#if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
	int count[factors][factors][factors];
	bzero(count, sizeof(count));
	#endif
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	printf("// start of generated code");
	#endif
	for (n_coeffs::iterator it = c.begin(); it < c.end(); ++it) {
	coeffs& co = *it;
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	printf("\nstatic double calc_%c(double a, double b, double c, double d,"
	"\n double e, double f, double g, double h) {"
	"\n return"
	"\n ", 'A' + (it - c.begin()));
	if (co[0].s > 0) {
	printf(" ");
	}
	if (abs(co[0].s) == 1) {
	printf(" ");
	}
	#endif
	for (coeffs::iterator ct = co.begin(); ct < co.end(); ++ct) {
	const coeff& cf = *ct;
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	printf(" ");
	bool firstFactor = false;
	if (ct - co.begin() > 0 \|\| cf.s < 0) {
	printf("%c", cf.s < 0 ? '-' : '+');
	}
	if (ct - co.begin() > 0) {
	printf(" ");
	}
	if (abs(cf.s) > 1) {
	printf("%d * ", abs(cf.s));
	} else {
	if (ct - co.begin() > 0) {
	printf(" ");
	}
	}
	#endif
	for (int x = 0; x < factors; ++x) {
	if (cf.n[x] == 0) {
	continue;
	}
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	for (int y = 0 ; y < cf.n[x]; ++y) {
	if (y > 0 \|\| firstFactor) {
	printf(" * ");
	}
	printf("%c", 'a' + x);
	}
	firstFactor = true;
	#endif
	#if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
	for (int y = x; y < factors; ++y) {
	if (cf.n[y] == 0) {
	continue;
	}
	if (x == y && cf.n[y] == 1) {
	continue;
	}
	for (int z = y; z < factors; ++z) {
	if (cf.n[z] == 0) {
	continue;
	}
	if ((x == z \|\| y == z) && cf.n[z] == 1) {
	continue;
	}
	if (x == y && y == z && cf.n[z] == 2) {
	continue;
	}
	count[x][y][z]++;
	}
	}
	#endif
	}
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	if (ct + 1 < co.end()) {
	printf("\n");
	}
	#endif
	}
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	printf(";\n}\n");
	#endif
	}
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	printf("// end of generated code\n");
	#endif
	#if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
	const int bestCount = 20;
	int best[bestCount][4];
	bzero(best, sizeof(best));
	for (int x = 0; x < factors; ++x) {
	for (int y = x; y < factors; ++y) {
	for (int z = y; z < factors; ++z) {
	if (!count[x][y][z]) {
	continue;
	}
	for (int w = 0; w < bestCount; ++w) {
	if (best[w][0] < count[x][y][z]) {
	best[w][0] = count[x][y][z];
	best[w][1] = x;
	best[w][2] = y;
	best[w][3] = z;
	break;
	}
	}
	}
	}
	}
	for (int w = 0; w < bestCount; ++w) {
	printf("%c%c%c=%d\n", 'a' + best[w][1], 'a' + best[w][2],
	'a' + best[w][3], best[w][0]);
	}
	#endif
	#if WRITE_AS_NONOPTIMIZED_C_CODE
	printf("\n");
	#endif
	}

	/* results: variable triplets used 10 or more times:
	aah=14
	ade=14
	aeh=14
	dee=14
	bce=13
	beg=13
	beh=12
	bbe=11
	bef=11
	cee=11
	cef=11
	def=11
	ceh=10
	deg=10
	*/