Исправить ошибки в тесте Миллера и тесте Соловея-Штрассена

// нахождение остатка от деления на n числа number, возведенного в степень power
// доп. информация: http://stackoverflow.com/questions/8496182/calculating-powa-b-mod-n
// unsigned long long можно попробовать заменить на unsigned __int64
unsigned long long mod(unsigned int number, unsigned long long power, unsigned long long n)
{
  unsigned long long res = 1;
  while (power)
  {
    if (power % 2)
      res = (1ull * res * number) % n;
    number = (1ull * number * number) % n;
    power /= 2;
  }
  return res;
}

  unsigned long long i, temp, n = 2971215073ull;
  bool not_prime = 0;
  const count = 58;
  // p - массив простых чисел, их (возможно) должно быть достаточно, поэтому не применяется вычисление значений каждый раз
  // pq - "коэффициенты", их кол-во невелико, можно использовать динамический массив, но для простоты взят "статический"
  unsigned int s = 0, q, m, j, k, p[count] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,\
      131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271}, pq[count];
 
// (1) проверка n на то, является ли оно степенью числа
  for (i = 3; i <= sqrt(n); i+=2)
  {
    temp = n;
    unsigned int power = 0;
    while (temp/i > 0 && temp%i == 0)
    {
      temp /= i;
      power++;
    }
    if (temp == 1 && power > 1)
    {
      ShowMessage("составное (является степенью числа)");
      not_prime = 1;
      break;
    }
  }
 
  if (!not_prime)
  {
// (2) определяем первые m простых чисел p[m]
// функция f(n) = n^0.133 определяется из того, что расширенная гипотеза Римана не доказана
// иначе можно использовать f(n) = 2*log^2(n)
// здесь и далее в коментариях ^ - возведение в степень
    for (m = 0; m < count-1; m++)
      if (p[m]<=pow(n, 0.133) && pow(n, 0.133)<=p[m+1])
      {
        m++;
        break;
      }
// находим s и q, где n-1 = q*2^s, q - нечетное
    temp = n-1;
    while (temp/2 > 0 && temp%2 == 0)
    {
      temp /= 2;
      s++;
    }
    q = (n-1)/pow(2, s);
// i=1, (4)
 
// (3) если i<=m, то i=i+1
    for (i = 1; i <= m; i++)
    {
// (4) не очень понятна запись p[i] | n
// полагаю, что означает она делимость n на p[i], т.е. является ли p[i] делителем n
      if (n % p[i-1] == 0)
      {
        ShowMessage("составное (делится на простое)");
        not_prime = 1;
        break;
      }
// вычисляем p[i]^q mod n, p[i]^(q*2) mod n, ..., p[i]^(q*2^s) mod n, обозначим их pq[i]
      for (j = 0; j <= s; j++)
        pq[j] = mod(p[i-1], q*pow(2, j), n);
// (5) проверяем условие pq[i]^(q*2^s) != 1, если да, то число составное
      if (pq[s] != 1)
      {
        ShowMessage("составное (остаток от деления на pq[s] после возведения в степень не равен 1)");
        not_prime = 1;
        break;
      }
// (6) проверяем условие pq[i]^q == 1, если да, то идем на (3)
      if (pq[0] == 1)
        continue;
// тут тоже мало понятно, предположил, что нужно найти индекс j, для которого pq[j] максимален и не равен 1
      temp = 0;
      for (j = 0; j <= s; j++)
        if (pq[j] != 1 && pq[j] > temp)
        {
          temp = pq[j];
          k = j;  // запоминаем индекс j (в переменной k)
        }
// (7) если pq[j] (j=k, см. чуть выше), то (3)
      if (pq[k] == n-1)
        continue;
 
// (8)
      ShowMessage("составное");
      not_prime = 1;
    }
// (3) случай i>m
    if (!not_prime)
      ShowMessage("простое");
  }

unsigned __int64 nod(unsigned __int64 m, unsigned __int64 n)
{
  while(m!=0 && n!=0)
  {
    if(m>=n)
      m%=n;
    else
      n%=m;
  }
  return m+n;
}

int jacob_char(unsigned __int64 a, unsigned __int64 b)
{
  if (nod(a, b) != 1)
    return 0;
  else
  {
    int r = 1;
    if (a < 0)
    {
      a = -a;
      if (b%4 == 3)
        r = -r;
    }
    do
    {
      int t = 0;
      while (a%2 == 0)
      {
        t += 1;
        a /= 2;
      }
      if (t%2 == 1 && (b%8 == 3 || b%8 == 5))
          r = -r;
      if (a%4 == 3 && b%4 == 3)
        r = -r;
      unsigned __int64 c = a;
      a = b%c;
      b = c;
    }
    while (a != 0);
    return r;
  }
}

  unsigned __int64 n = 45212177, k = 1, i, a;
  bool not_prime = 0;
  randomize();
  for (i = 1; i <= k; i++)
  {
// a должно быть случайным, например, random(n-4)+3;
// думаю, что числа a не должны повторяться и должны быть нечетными
    a = 17;
/*  ShowMessage(a);
    ShowMessage(nod(a, n));
    ShowMessage(mod(a, (n-1)/2, n));
    ShowMessage(jacob_char(a, n));*/
    if (nod(a, n) > 1 || mod(a, 1ull*(n-1)/2, n) != jacob_char(a, n))
    {
      ShowMessage("n составное");
      not_prime = 1;
      break;
    }
  }
  if (!not_prime)
    ShowMessage("n простое с вероятностью " + String(1.-pow(2., -k*1.)));

//---------------------------------------------------------------------------
 
#include "stdafx.h"
#include <iostream>
#include <math.h>
#include <conio.h>
#include <stdlib.h>
#pragma hdrstop
 
using namespace std;
//---------------------------------------------------------------------------
//пример простых чисел 13, 109, 193, 433, 769, 1201, 1453, 2029, 3469, 3889, 4801, 10093, 12289, 13873, 18253, 20173, 21169, 22189, 28813, 37633, 43201, 47629, 60493, 63949, 65713, 69313, 73009, 76801, 84673, 106033, 108301, 112909, 115249
// нахождение остатка от деления на n числа number, возведенного в степень power
 long long mod( int number,  long long power,  long long n)
{
  unsigned long long res = 1;
  while (power)
  {
    if (power % 2)
      res = (res * number) % n;
    number = (number * number) % n;
    power /= 2;
  }
  return res;
}
 
 long NOD(long m, long n)
{
  while(m!=0 && n!=0)
  {
    if(m>=n)
      m%=n;
    else
      n%=m;
  }
  return m+n;
}
 
int jacob_char(int a, int b)
{
  if (NOD(a, b) != 1)
    return 0;
  else
  {
    int r = 1;
    if (a < 0)
    {
      a = -a;
      if (b%4 == 3)
        r = -r;
    }
    do
    {
      int t = 0;
      while (a%2 == 0)
      {
        t += 1;
        a /= 2;
      }
      if (t%2 != 0 && (b%8 == 3 || b%8 == 5))
          r = -r;
      if (a%4 == 3 && b%4 == 3)
        r = -r;
      int c = a;
      a = b%c;
      b = c;
    }
    while (a != 0);
    return r;
  }
}
 
int main(int argc, char* argv[])
{
        long n;/*n>2, нечетное натуральное число, которое необходимо проверить на простоту*/
        long k;/* параметр, определяющий точность теста*/
        long a;
        double t;
        int flag = 0;
                 
 
        cin>>n;
        cin>>k;
        t=1-1/pow(2.0,k);
                    
            for (int i=1; i <= k; i++)
            {
            a=rand()%(n-2)+2;
            if (NOD(a, n) > 1)
            {
                
                flag=1;
                break;
            }
            else if ( mod(a,(n-1)/2,n) <= jacob_char(a, n))
            {
                
                flag=2;
                break;
            }
            }
        
              
  if (flag !=0)
  {
      cout<<"sostavnoe \n";
  }
  else
       cout<<"prostoe s veroyatnost'u = " << t << "\n";
     getch();                                 
   
    
                
    return 0;
}
//---------------------------------------------------------------------------