#include /* 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 coeffs; typedef std::vector 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 */