// // Stack Calculator of Big Integers // implemented over GMP // (GNU Multiple Presizion library) // // The program uses C and C++ intarfaces of GMP // for representing big integers // #include #include #include #include #include #include #include #include #include #include #include #include #include #include static void onAdd(); static void onSub(); static void onMul(); static void onDiv(); static void onMod(); static void onPow(); static void onPowMod(); static void onInvMod(); static void onGcd(); static void onEUA(); static void onPush(const char* line); static void onPop(); static void onDup(); static void onExch(); static void onInit(); static void display(); static void printHelp(); static void onShow(); static void onGenerateRandom(); static void onTestPrime(); static void onGeneratePrime(); static void onFactorize(); static bool pollardP1Factorization( const mpz_class& m, mpz_class& factor1, mpz_class& factor2 ); static std::deque stack; // Random state static gmp_randclass randState(gmp_randinit_default); static const int MAX_PRIME_POWER = 500000; static std::vector primes; static std::vector prime_powers; static int num_primes = 0; static int num_prime_powers = 0; static void generate_prime_powers(); int main() { char line[256]; printHelp(); // srand((unsigned int) time(0)); time_t curtime; unsigned long curtim = (unsigned long) time(&curtime); randState.seed(curtim); // gets(line); line[0] = 0; int linelen = scanf("%s", line); while (*line != 0) { linelen = (int) strlen(line); try { if (strcmp(line, "+") == 0) onAdd(); else if (strcmp(line, "-") == 0) onSub(); else if (strcmp(line, "*") == 0) onMul(); else if (strcmp(line, "/") == 0) onDiv(); else if (strcmp(line, "%") == 0 || strcmp(line, "mod") == 0) onMod(); else if (strcmp(line, "^") == 0 || strcmp(line, "pow") == 0) onPow(); else if (strcmp(line, "powmod") == 0) onPowMod(); else if (strcmp(line, "invmod") == 0) onInvMod(); else if (strcmp(line, "gcd") == 0) onGcd(); else if (strcmp(line, "eua") == 0 || strcmp(line, "EUA") == 0) onEUA(); else if (strcmp(line, "rand") == 0) onGenerateRandom(); else if (strcmp(line, "tstprime") == 0) onTestPrime(); else if (strcmp(line, "genprime") == 0) onGeneratePrime(); else if (strcmp(line, "factorize") == 0) onFactorize(); else if (strcmp(line, "=") == 0) display(); else if ( (linelen > 0 && isdigit(line[0])) || (linelen > 1 && (line[0] == '-' || line[0] == '+') && isdigit(line[1])) ) onPush(line); else if (strcmp(line, "pop") == 0) onPop(); else if (strcmp(line, "dup") == 0) onDup(); else if (strcmp(line, "exch") == 0) onExch(); else if (strcmp(line, "init") == 0) onInit(); else if (strcmp(line, "show") == 0) onShow(); else if ( strcmp(line, "") == 0 || (linelen > 0 && (line[0] == 'q' || line[0] == 'Q')) ) break; else printHelp(); } catch (std::exception& e) { printf("Exception: %s\n", e.what()); } // gets(line); if (scanf("%s", line) <= 0) break; } return 0; } static void onAdd() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); stack.push_front(x + y); display(); } static void onSub() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); stack.push_front(x - y); display(); } static void onMul() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); stack.push_front(x * y); display(); } static void onDiv() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); if (y == 0) throw std::runtime_error("Zero division"); stack.push_front(x / y); display(); } static void onMod() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); if (y == 0) throw std::runtime_error("Zero division"); stack.push_front(x % y); display(); } static void onPow() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } long y = mpz_get_si( stack.front().get_mpz_t() ); stack.pop_front(); if (y < 0) throw std::runtime_error("Illegal exponent"); mpz_class x = stack.front(); stack.pop_front(); mpz_class res; mpz_pow_ui(res.get_mpz_t(), x.get_mpz_t(), (unsigned long) y); stack.push_front(res); display(); } static void onPowMod() { if (stack.size() < 3) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class m = stack.front(); stack.pop_front(); mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); mpz_class res; mpz_powm(res.get_mpz_t(), x.get_mpz_t(), y.get_mpz_t(), m.get_mpz_t()); stack.push_front(res); display(); } static void onInvMod() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class m = stack.front(); stack.pop_front(); if (m == 0) throw std::runtime_error("Zero module"); mpz_class x = stack.front(); stack.pop_front(); if (x == 0) throw std::runtime_error("Zero division"); mpz_class res; int invExists = mpz_invert( res.get_mpz_t(), x.get_mpz_t(), m.get_mpz_t() ); if (invExists == 0) { printf("No inverse modulo m.\n"); res = 0; } stack.push_front(res); display(); } static void onGcd() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); mpz_class res; mpz_gcd(res.get_mpz_t(), x.get_mpz_t(), y.get_mpz_t()); stack.push_front(res); display(); } static void onEUA() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class u, v, d; mpz_class y = stack.front(); stack.pop_front(); mpz_class x = stack.front(); stack.pop_front(); mpz_gcdext( d.get_mpz_t(), u.get_mpz_t(), v.get_mpz_t(), x.get_mpz_t(), y.get_mpz_t() ); stack.push_front(u); stack.push_front(v); stack.push_front(d); display(); } static void onPush(const char* line) { mpz_class x(line); stack.push_front(x); } static void onPop() { stack.pop_front(); } static void onDup() { stack.push_front(stack.front()); } static void onExch() { if (stack.size() < 2) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class x = stack.front(); stack.pop_front(); mpz_class y = stack.front(); stack.pop_front(); stack.push_front(x); stack.push_front(y); } static void onInit() { stack.clear(); } static void display() { if (!stack.empty()) { printf("= %s\n", stack.front().get_str().c_str()); } else { printf("stack empty\n"); } } static void onShow() { int d = stack.size(); printf("Depth of stack = %d.", d); if (d > 0) printf(" Stack elements:\n"); else printf("\n"); for (int i = 0; i < d; i++) { printf( "= %s\n", stack[i].get_str().c_str() ); } } static void onGenerateRandom() { if (stack.size() < 1) { stack.clear(); throw std::underflow_error("Stack underflow"); } unsigned long len = mpz_get_ui( stack.front().get_mpz_t() ); stack.pop_front(); if (len > 10000) { printf("Length is too large. Using the value 10000.\n"); len = 10000; } stack.push_front( randState.get_z_bits(len) ); display(); } static void onTestPrime() { if (stack.size() < 1) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class m = stack.front(); stack.pop_front(); int res = mpz_probab_prime_p(m.get_mpz_t(), 25); if (res == 0) { printf("Number is composite.\n"); } else if (res == 1) { printf("Number is probably prime.\n"); } else { assert(res == 2); printf("Number is definitely prime.\n"); } stack.push_front(mpz_class(res)); display(); } static void onGeneratePrime() { if (stack.size() < 1) { stack.clear(); throw std::underflow_error("Stack underflow"); } unsigned long len = mpz_get_ui( stack.front().get_mpz_t() ); stack.pop_front(); if (len <= 1) { len = 2; } else if (len > 10000) { printf("Length is too large. Using the value 10000.\n"); len = 10000; } mpz_class two(2); mpz_class pow2Len; mpz_pow_ui(pow2Len.get_mpz_t(), two.get_mpz_t(), len-1); mpz_class m; bool found = false; int maxAttempts = 100000; int attempt = 0; while (!found && attempt < maxAttempts) { m = randState.get_z_bits(len); if (m < pow2Len) m += pow2Len; if (mpz_divisible_2exp_p(m.get_mpz_t(), 1)) { // even == divisible by 2^1 ++m; } int res = mpz_probab_prime_p(m.get_mpz_t(), 25); if (res != 0) { found = true; } else { ++attempt; } } if (!found) { printf("Could not generate a prime.\n"); m = 0; } stack.push_front(m); display(); } static void onFactorize() { if (stack.size() < 1) { stack.clear(); throw std::underflow_error("Stack underflow"); } mpz_class m = stack.front(); if (m == 0) { printf("Cannot factorize (zero argument).\n"); return; } stack.pop_front(); int res = mpz_probab_prime_p(m.get_mpz_t(), 25); if (res != 0) { // m is probably prime printf(" Probably prime.\n"); stack.push_front(m); stack.push_front(mpz_class(1)); display(); return; } if (num_primes == 0) generate_prime_powers(); int numFactors = 0; int i = 0; while (m >= 1 && i < num_primes && primes[i] <= m) { if (m % primes[i] == 0) { stack.push_front(mpz_class(primes[i])); ++numFactors; printf(" factor %d\n", primes[i]); m /= primes[i]; } else { ++i; // to next prime } } if (m == 1) { stack.push_front(mpz_class(numFactors)); display(); return; } res = mpz_probab_prime_p(m.get_mpz_t(), 25); if (res != 0) { // m is probably prime stack.push_front(m); printf(" factor %s\n", m.get_str().c_str()); ++numFactors; stack.push_front(mpz_class(numFactors)); display(); return; } mpz_class f1, f2; while ( m != 1 && mpz_probab_prime_p(m.get_mpz_t(), 25) == 0 // m is composite ) { if (!pollardP1Factorization(m, f1, f2)) { stack.push_front(m); printf(" factor %s\n", m.get_str().c_str()); ++numFactors; stack.push_front(mpz_class(numFactors)); printf(" Cannot factorize.\n"); display(); return; } /*... printf( " f1 = %s, f2 = %s\n", f1.get_str().c_str(), f2.get_str().c_str() ); ...*/ stack.push_front(f1); printf(" factor %s\n", f1.get_str().c_str()); ++numFactors; m = f2; } if (m != 1) { stack.push_front(m); ++numFactors; printf(" factor %s\n", m.get_str().c_str()); } stack.push_front(mpz_class(numFactors)); display(); } static void printHelp() { printf( "Stack Calculator commands:\n" "\t\tPush to stack\n" "\t+, -, *, /, %%\tAriphmetic operations\n" "\t^\t\tInteger power\n" "\tpowmod\t\tPower modulo m\n" "\tinvmod\t\tInverse modulo m\n" "\tgcd\t\tGreatest Common Divisor\n" "\tEUA\t\tExt.Euclid Alg.: in: m,n; out: u,v,d: um+vn=d=gcd(m,n)\n" "\trand\t\tGenerate random integer (input: length in bits)\n" "\ttstprime\tProbabilistic prime test\n" "\tgenprime\tGenerate random prime (input: length in bits)\n" "\tfactorize\tInteger factorization\n" "\t=\t\tDisplay the stack top\n" "\tpop\t\tRemove the stack top\n" "\tdup\t\tDuplicate the stack top to the stack\n" "\texch\t\tExchange (swap) two upper elements in the stack top\n" "\tshow\t\tShow the stack\n" "\tinit\t\tErase the stack\n" "\tquit\t\tEnd the program\n" ); } static void generate_prime_powers() { // 1. Generate primes <= MAX_PRIME_POWER primes.push_back(2); primes.push_back(3); primes.push_back(5); primes.push_back(7); num_primes = 4; unsigned int p = 11; while (num_primes < MAX_PRIME_POWER && p < MAX_PRIME_POWER) { bool prime = true; for ( int i = 0; prime && i < num_primes && primes[i] * primes[i] <= p; ++i ) { if (p % primes[i] == 0) prime = false; } if (prime) { primes.push_back(p); ++num_primes; //... printf("prime %d\n", p); } p += 2; } num_prime_powers = num_primes; prime_powers.resize(num_primes); unsigned int root = (unsigned int) sqrt((double) MAX_PRIME_POWER); for (int i = 0; i < num_prime_powers; ++i) { unsigned int p = primes[i]; prime_powers[i] = p; if (p <= root) { unsigned int pNext = p*p; while (pNext <= MAX_PRIME_POWER) { prime_powers[i] = pNext; pNext *= p; } } //... printf("prime power %d\n", prime_powers[i]); } //... printf("--------\n", p); } static bool pollardP1Factorization( const mpz_class& m, mpz_class& factor1, mpz_class& factor2 ) { assert(m != 0); if (num_prime_powers == 0) generate_prime_powers(); int res = mpz_probab_prime_p(m.get_mpz_t(), 25); if (res != 0) { // printf("Number is probably prime.\n"); return false; } int maxTries = 4; int maxAttempts = 100; int gcdStep = 16; int step = 0; bool factorized = false; mpz_class b, product, power, d; while (!factorized && step < maxTries) { // Generate a prime base int attempt = 0; while (attempt < maxAttempts) { b = randState.get_z_bits(15); if (mpz_divisible_2exp_p(b.get_mpz_t(), 1)) { // even == divisible by 2^1 ++b; } int res = mpz_probab_prime_p(b.get_mpz_t(), 25); if (res != 0) { break; } else { ++attempt; } } //... printf("base = %s\n", b.get_str().c_str()); product = 1; int primeIdx = 0; int productStep = 0; while (!factorized && primeIdx < num_prime_powers) { mpz_powm_ui( power.get_mpz_t(), b.get_mpz_t(), prime_powers[primeIdx], m.get_mpz_t() ); b = power; ++primeIdx; product *= power - 1; product %= m; ++productStep; if ( productStep >= gcdStep || primeIdx >= num_prime_powers ) { mpz_gcd( d.get_mpz_t(), product.get_mpz_t(), m.get_mpz_t() ); if (d != 1) { factorized = true; factor1 = d; factor2 = m / d; return true; } product = 1; productStep = 0; } } ++step; } return false; }