#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
// Нахождение GCD
long long gcd(long long a, long long b){
if(b == 0)
return a;
return gcd(b, a % b);
}
// Двоичное умножение по модулю для защиты от переполнений
long long mul(long long a, long long b, long long m){
if(b == 1)
return a;
if(b % 2 == 0){
long long t = mul(a, b/2, m);
return (2 * t) % m;
}
return (mul(a, b - 1, m) + a) % m;
}
// Возведение в степень по модулю
long long modpow (long long a, long long b, long long m){
if(b==0)
return 1;
if(b % 2 == 0){
long long t = modpow (a, b/2, m);
return mul(t, t, m) % m;
}
return (mul(modpow (a, b - 1, m), a, m)) % m;
}
int jacobi (int nominator, int denominator)
{
int s = 0;
int remainder, power_of_two, t;
do{
remainder = nominator % denominator;
power_of_two = t = 0;
while(remainder % 2 == 0)
{
power_of_two++;
remainder >>= 1;
}
t = remainder;
s = (s + power_of_two * (denominator * denominator - 1) / 8 + (t - 1) * (denominator - 1) / 4) % 2;
if(t == 1)
return (s) ? -1 : 1;
nominator = denominator;
denominator = t;
}while(t >= 3);
}
int main(int argc, char* argv[])
{
long n; /*n > 2, нечетное натуральное число, которое необходимо проверить на простоту*/
long k; /* параметр, определяющий точность теста */
long a;
long residue;
double t;
int flag = 0;
std::cin >> n;
std::cin >> k;
t = 1 - 1 / pow(2.0, k);
srand(time(NULL));
for (int i=1; i <= k; i++)
{
a = rand() % (n - 2) + 2;
long long p = modpow(a, (n - 1) / 2, n);
int j = jacobi (a, n);
residue = p - n;
if (gcd(a, n) != 1)
{
flag = 1;
break;
}
else if (residue != j)
{
std::cout << "[DEBUG]: p - n = " << (p - n) << std::endl;
std::cout << "[DEBUG]: residue = " << residue << std::endl;
std::cout << "[DEBUG]: modpow = " << p << std::endl;
std::cout << "[DEBUG]: jacobi = " << j << std::endl;
flag = 2;
break;
}
}
if (flag !=0)
{
std::cout << "composite \n";
}
else
std::cout<<"prime with possibility " << t << "\n";
return 0;
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8c3RkbGliLmg+CiNpbmNsdWRlIDx0aW1lLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgovLyDQndCw0YXQvtC20LTQtdC90LjQtSBHQ0QKbG9uZyBsb25nIGdjZChsb25nIGxvbmcgYSwgbG9uZyBsb25nIGIpewoJaWYoYiA9PSAwKQoJCXJldHVybiBhOwoJcmV0dXJuIGdjZChiLCBhICUgYik7Cn0KCi8vINCU0LLQvtC40YfQvdC+0LUg0YPQvNC90L7QttC10L3QuNC1INC/0L4g0LzQvtC00YPQu9GOINC00LvRjyDQt9Cw0YnQuNGC0Ysg0L7RgiDQv9C10YDQtdC/0L7Qu9C90LXQvdC40LkKbG9uZyBsb25nIG11bChsb25nIGxvbmcgYSwgbG9uZyBsb25nIGIsIGxvbmcgbG9uZyBtKXsKCWlmKGIgPT0gMSkKCQlyZXR1cm4gYTsKCWlmKGIgJSAyID09IDApewoJCWxvbmcgbG9uZyB0ID0gbXVsKGEsIGIvMiwgbSk7CgkJcmV0dXJuICgyICogdCkgJSBtOwoJfQoJcmV0dXJuIChtdWwoYSwgYiAtIDEsIG0pICsgYSkgJSBtOwp9CgovLyDQktC+0LfQstC10LTQtdC90LjQtSDQsiDRgdGC0LXQv9C10L3RjCDQv9C+INC80L7QtNGD0LvRjgpsb25nIGxvbmcgbW9kcG93IChsb25nIGxvbmcgYSwgbG9uZyBsb25nIGIsIGxvbmcgbG9uZyBtKXsKCWlmKGI9PTApCgkJcmV0dXJuIDE7CglpZihiICUgMiA9PSAwKXsKCQlsb25nIGxvbmcgdCA9IG1vZHBvdyAoYSwgYi8yLCBtKTsKCQlyZXR1cm4gbXVsKHQsIHQsIG0pICUgbTsKCX0KCXJldHVybiAobXVsKG1vZHBvdyAoYSwgYiAtIDEsIG0pLCBhLCBtKSkgJSBtOwp9CgppbnQgamFjb2JpIChpbnQgbm9taW5hdG9yLCBpbnQgZGVub21pbmF0b3IpCnsKICAgIGludCBzID0gMDsKICAgIGludCByZW1haW5kZXIsIHBvd2VyX29mX3R3bywgdDsKICAgIGRvewogICAgICAgIHJlbWFpbmRlciA9IG5vbWluYXRvciAlIGRlbm9taW5hdG9yOwogICAgICAgIHBvd2VyX29mX3R3byA9IHQgPSAwOwogICAgICAgIHdoaWxlKHJlbWFpbmRlciAlIDIgPT0gMCkKICAgICAgICB7CiAgICAgICAgICAgIHBvd2VyX29mX3R3bysrOyAKICAgICAgICAgICAgcmVtYWluZGVyID4+PSAxOwogICAgICAgIH0KICAgICAgICB0ID0gcmVtYWluZGVyOwogICAgICAgIHMgPSAocyArIHBvd2VyX29mX3R3byAqIChkZW5vbWluYXRvciAqIGRlbm9taW5hdG9yIC0gMSkgLyA4ICsgKHQgLSAxKSAqIChkZW5vbWluYXRvciAtIDEpIC8gNCkgJSAyOwogICAgICAgIGlmKHQgPT0gMSkKICAgICAgICAgICAgcmV0dXJuIChzKSA/IC0xIDogMTsKICAgICAgICBub21pbmF0b3IgPSBkZW5vbWluYXRvcjsKICAgICAgICBkZW5vbWluYXRvciA9IHQ7CiAgICB9d2hpbGUodCA+PSAzKTsKfQoKaW50IG1haW4oaW50IGFyZ2MsIGNoYXIqIGFyZ3ZbXSkKewoJbG9uZyBuOwkJLypuID4gMiwg0L3QtdGH0LXRgtC90L7QtSDQvdCw0YLRg9GA0LDQu9GM0L3QvtC1INGH0LjRgdC70L4sINC60L7RgtC+0YDQvtC1INC90LXQvtCx0YXQvtC00LjQvNC+INC/0YDQvtCy0LXRgNC40YLRjCDQvdCwINC/0YDQvtGB0YLQvtGC0YMqLwoJbG9uZyBrOwkJLyog0L/QsNGA0LDQvNC10YLRgCwg0L7Qv9GA0LXQtNC10LvRj9GO0YnQuNC5INGC0L7Rh9C90L7RgdGC0Ywg0YLQtdGB0YLQsCAqLwoJbG9uZyBhOwoJbG9uZyByZXNpZHVlOwoJZG91YmxlIHQ7CglpbnQgZmxhZyA9IDA7CiAgICAgICAgICAgICAgICAgCglzdGQ6OmNpbiA+PiBuOwoJc3RkOjpjaW4gPj4gazsKICAgIHQgPSAxIC0gMSAvIHBvdygyLjAsIGspOwoKCXNyYW5kKHRpbWUoTlVMTCkpOwogICAgICAgICAgICAgICAgICAgIAoJZm9yIChpbnQgaT0xOyBpIDw9IGs7IGkrKykKCXsKCQlhID0gcmFuZCgpICUgKG4gLSAyKSArIDI7CgkJbG9uZyBsb25nIHAgPSBtb2Rwb3coYSwgKG4gLSAxKSAvIDIsIG4pOwoJCWludCBqID0gamFjb2JpIChhLCBuKTsKCQlyZXNpZHVlID0gcCAtIG47CgoJCWlmIChnY2QoYSwgbikgIT0gMSkKCQl7CgkJCWZsYWcgPSAxOwoJCQlicmVhazsKCQl9CgkJZWxzZSBpZiAocmVzaWR1ZSAhPSBqKQoJCXsKCQkJc3RkOjpjb3V0IDw8ICJbREVCVUddOiBwIC0gbiA9ICIgPDwgKHAgLSBuKSA8PCBzdGQ6OmVuZGw7CgkJCXN0ZDo6Y291dCA8PCAiW0RFQlVHXTogcmVzaWR1ZSA9ICIgPDwgcmVzaWR1ZSA8PCBzdGQ6OmVuZGw7CgkJCXN0ZDo6Y291dCA8PCAiW0RFQlVHXTogbW9kcG93ID0gIiA8PCBwIDw8IHN0ZDo6ZW5kbDsKCQkJc3RkOjpjb3V0IDw8ICJbREVCVUddOiBqYWNvYmkgPSAiIDw8IGogPDwgc3RkOjplbmRsOwoJCQlmbGFnID0gMjsKCQkJYnJlYWs7CgkJfQoJfQogICAgICAgIAogICAgICAgICAgICAgIAoJaWYgKGZsYWcgIT0wKQoJewoJCXN0ZDo6Y291dCA8PCAiY29tcG9zaXRlIFxuIjsKCX0KCWVsc2UKCQlzdGQ6OmNvdXQ8PCJwcmltZSB3aXRoIHBvc3NpYmlsaXR5ICIgPDwgdCA8PCAiXG4iOyAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgIHJldHVybiAwOwp9