primes.c++

#include <cassert>
#include <cstdio>
#include <cstring>

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>

using namespace std;


class Number {
  private:
    int n;
    unsigned *data;
    Number ( int n );
  public:
    Number ( int n, const string &decimal );
    Number ( const Number &number );
    Number& operator = ( const Number &number );
    ~Number ();

    explicit operator bool () const;
    bool equals ( unsigned value ) const;
    bool operator == ( const Number &b ) const;
    bool operator != ( const Number &b ) const;
    bool operator <= ( const Number &b ) const;
    bool operator < ( const Number &b ) const;
    bool operator >= ( const Number &b ) const;
    bool operator > ( const Number &b ) const;

    int operator % ( int b ) const;
    Number& operator /= ( int b );

    Number operator + ( const Number &b ) const;
    Number operator - ( const Number &b ) const;
    Number operator * ( const Number &b ) const;
    Number operator % ( const Number &b ) const;

    Number reduce () const;
    Number shift () const;
    Number negate () const;

    static Number random ( const Number &limit );
    Number limit () const;

    unsigned long long debug_value () const;
    string debug_str () const;
    void debug ( const char *name, bool full = false ) const;
};


class Modulo {
  private:
    Number n, r, n1;
    Number mm_shift ( const Number &t ) const;
  public:
    Modulo ( const Number &n );

    Number one_ () const;
    Number mul_ ( const Number &a_, const Number &b_ ) const;
    Number pow_ ( Number a_, Number b ) const;
};


static int cnt[10] = {0};

const Number one (10, "1");
const Number two (10, "2");

Modulo::Modulo ( const Number &_n) : n (_n), r (_n), n1 (_n) {
  r = n.limit () % n;
  n1 = one;
  Number t = n;
  for (int i = 1; i < 319; i++) {
    n1 = (n1 * t).reduce ();
    cnt[1]++;
    t = (t * t).reduce ();
    cnt[1]++;
  }
  n1 = n1.negate ();
}

Number Modulo::mm_shift ( const Number &t ) const {
  Number t2 = t.reduce ();
  t2 = (t2 * n1).reduce ();
  cnt[3]++;
  Number tt = t + t2 * n;
  cnt[3]++;
  Number result = tt.shift ();
  if (result >= n) {
    result = result - n;
  }
  return result;
}

Number Modulo::one_ () const {
  return r;
}

Number Modulo::mul_ ( const Number &a_, const Number &b_ ) const {
  return mm_shift (a_ * b_);
}

Number Modulo::pow_ ( Number a_, Number b ) const {
  Number r_ = r;
  for (; b; b /= 2) {
    if (b % 2) {
      cnt[2] += 3;
      r_ = mm_shift (r_ * a_);
    }
    cnt[2] += 3;
    a_ = mm_shift (a_ * a_);
  }
  return r_;
}

vector <int> primes_small;

bool check_prime ( const Number &p ) {
  if (p < two)
    return false;
  for (auto x : primes_small) {
    if (p % x != 0)
      continue;
    return p.equals (x);
  }

  Modulo m (p);
  Number one_ = m.one_ ();

  Number u = p - one;
  int k = 0;
  while (u % 2 == 0)
    k++, u /= 2;

  for (int step = 0; step < 100; step++) {
    Number x_ = Number::random (p - one - one) + one, y_ (p);
    x_ = m.pow_ (x_, u);
    if (x_ == one_)
      continue;
    for (int i = 0; i < k; i++) {
      y_ = m.mul_ (x_, x_);
      if (y_ == one_ && !(x_ == p - one_))
        return false;
      x_ = y_;
      if (x_ == one_)
        break;
    }
    if (!(x_ == one_))
      return false;
  }
  return true;
}

int main () {
  const int PRIMES_LIMIT = 300;
  static int prime_dividor[PRIMES_LIMIT];
  memset (prime_dividor, 0, sizeof (prime_dividor));
  for (int i = 2; i < PRIMES_LIMIT; i++) {
    if (prime_dividor[i] == 0) {
      primes_small.push_back (i);
      prime_dividor[i] = i;
    }
    for (int j = 0; j < (int)primes_small.size () && primes_small[j] <= prime_dividor[i] && i * primes_small[j] < PRIMES_LIMIT; j++)
      prime_dividor[i * primes_small[j]] = primes_small[j];
  }

  srand (301703);
  string sa, sb;
  cin >> sa >> sb;
  Number a (10, sa);
  Number b (10, sb);

  int ans = 0;
  for (; a <= b; a = a + one) {
    if (check_prime (a)) {
      fprintf (stderr, "[debug] prime: %s\n", a.debug_str ().c_str ());
      ans++;
    }
  }
  cout << ans << endl;

  for (int i = 0; i < 10; i++)
    fprintf (stderr, "cnt[i=%d] = %d\n", i, cnt[i]);

  return 0;
}


Number::Number ( int _n ) : n (_n) {
  data = new unsigned[n];
  // memset (data, 0, sizeof (data[0]) * n);
}

Number::Number ( int _n, const string &decimal ) : n (_n) {
  data = new unsigned[n];
  memset (data, 0, sizeof (data[0]) * n);
  for (auto digit : decimal) {
    unsigned long long carry = digit - '0';
    for (int i = 0; i < n; i++) {
      carry = carry + (unsigned long long)data[i] * 10;
      data[i] = carry & 0xffffffff;
      carry >>= 32;
    }
    assert (carry == 0);
  }
}

Number::Number ( const Number &number ) : n (number.n) {
  data = new unsigned[n];
  memcpy (data, number.data, sizeof (data[0]) * n);
}

Number& Number::operator = ( const Number &number ) {
  assert (n == number.n);
  memcpy (data, number.data, sizeof (data[0]) * n);
  return *this;
}

Number::~Number () {
  delete[] data;
}

Number::operator bool () const {
  for (int i = 0; i < n; i++)
    if (data[i])
      return true;
  return false;
}

bool Number::equals ( unsigned value ) const {
  for (int i = 1; i < n; i++)
    if (data[i])
      return false;
  return data[0] == value;
}

bool Number::operator == ( const Number &b ) const {
  assert (n == b.n);
  for (int i = n - 1; i >= 0; i--)
    if (data[i] != b.data[i])
      return false;
  return true;
}

bool Number::operator <= ( const Number &b ) const {
  assert (n == b.n);
  for (int i = n - 1; i >= 0; i--)
    if (data[i] != b.data[i])
      return data[i] < b.data[i];
  return true;
}

bool Number::operator < ( const Number &b ) const {
  assert (n == b.n);
  for (int i = n - 1; i >= 0; i--)
    if (data[i] != b.data[i])
      return data[i] < b.data[i];
  return false;
}

bool Number::operator >= ( const Number &b ) const {
  assert (n == b.n);
  for (int i = n - 1; i >= 0; i--)
    if (data[i] != b.data[i])
      return data[i] > b.data[i];
  return true;
}

int Number::operator % ( int b ) const {
  unsigned long long result = 0;
  for (int i = n - 1; i >= 0; i--) {
    result = ((result << 32) + data[i]) % b;
  }
  return result;
}

Number& Number::operator /= ( int b ) {
  unsigned long long carry = 0;
  for (int i = n - 1; i >= 0; i--) {
    carry = (carry << 32) + data[i];
    data[i] = carry / b;
    carry %= b;
  }
  return *this;
}

Number Number::operator + ( const Number &b ) const {
  assert (n == b.n);
  Number result (n);
  unsigned long long carry = 0;
  for (int i = 0; i < n; i++) {
    carry = carry + data[i] + b.data[i];
    result.data[i] = carry & 0xffffffff;
    carry >>= 32;
  }
  assert (carry == 0);
  return result;
}

Number Number::operator - ( const Number &b ) const {
  assert (n == b.n);
  Number result (n);
  unsigned long long carry = 0;
  for (int i = 0; i < n; i++) {
    carry = (1ll << 32) + data[i] - b.data[i] - carry;
    result.data[i] = carry & 0xffffffff;
    carry >>= 32;
    carry = !carry;
  }
  assert (carry == 0);
  return result;
}

Number Number::operator * ( const Number &b ) const {
  cnt[0]++;
  assert (n == b.n);
  static long long tmp[200];
  assert (2 * n <= 200);
  memset (tmp, 0, sizeof (tmp[0]) * 2 * n);
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++) {
      unsigned long long t = (unsigned long long)data[i] * b.data[j];
      tmp[i + j] += t & 0xffffffff;
      tmp[i + j + 1] += t >> 32;
    }

  Number result (2 * n);
  unsigned long long carry = 0;
  for (int i = 0; i < 2 * n; i++) {
    carry += tmp[i];
    result.data[i] = carry & 0xffffffff;
    carry >>= 32;
  }
  assert (carry == 0);
  return result;
}

Number Number::operator % ( const Number &b ) const {
  assert (n >= b.n);
  Number tmp = (*this);

  int ib = b.n, it = tmp.n;
  while (ib > 0 && b.data[ib - 1] == 0)
    ib--;
  assert (ib);

  while (true) {
    while (it > 0 && tmp.data[it - 1] == 0)
      it--;
    if (it < ib)
      break;
    bool less = false;
    for (int i = 0; i < ib; i++) {
      if (tmp.data[it - i - 1] != b.data[ib - i - 1]) {
        less = tmp.data[it - i - 1] < b.data[ib - i - 1];
        break;
      }
    }
    if (it == ib && less)
      break;

    unsigned long long k = tmp.data[it - 1] / b.data[ib - 1];
    int shift = 0;
    if ((k == 1 && less) || k == 0) {
      k = ((unsigned long long)tmp.data[it - 1] << 32) + tmp.data[it - 2];
      k /= b.data[ib - 1];
      shift = 1;
    }
    if (k > 1)
      k /= 2;

    assert (k != 0);

    unsigned long long carry = 0;
    for (int i = 0; i < ib; i++) {
      int j = i + it - shift - ib;
      unsigned long long value = tmp.data[j] - k * b.data[i] - carry;
      tmp.data[j] = value;
      carry = (tmp.data[j] - value) >> 32;
    }
    if (shift) {
      assert (tmp.data[it - 1] >= carry);
      tmp.data[it - 1] -= carry;
    }
  }

  Number result (b.n);
  memcpy (result.data, tmp.data, sizeof (result.data[0]) * b.n);
  return result;
}

Number Number::reduce () const {
  Number result (n / 2);
  for (int i = 0; i < n / 2; i++)
    result.data[i] = data[i];
  return result;
}

Number Number::shift () const {
  Number result (n / 2);
  for (int i = 0; i < n / 2; i++)
    result.data[i] = data[i + n / 2];
  return result;
}

Number Number::negate () const {
  if (!*this)
    return *this;
  Number result (n);
  for (int i = 0; i < n; i++) {
    result.data[i] = 0xffffffff - data[i];
  }
  return result + one;
}

Number Number::random ( const Number &limit ) {
  Number result (limit.n);
  bool less = false;
  for (int i = limit.n - 1; i >= 0; i--) {
    result.data[i] = rand () ^ (rand () << 16);
    if (less)
      continue;
    if (limit.data[i] < 0xffffffff)
    result.data[i] %= (limit.data[i] + 1);
    if (result.data[i] < limit.data[i])
      less = true;
  }
  return result;
}

Number Number::limit () const {
  Number result (n * 2);
  memset (result.data, 0, sizeof (result.data[0]) * 2 * n);
  result.data[n] = 1;
  return result;
}

unsigned long long Number::debug_value () const {
  unsigned long long value = 0;
  for (int i = n - 1; i >= 0; i--)
    value = (value << 32) + data[i];
  return value;
}

string Number::debug_str () const {
  string number = "";
  Number tmp = *this;
  while (tmp) {
    number += '0' + (tmp % 10);
    tmp /= 10;
  }
  reverse (number.begin (), number.end ());
  return number;
}

void Number::debug ( const char *name, bool full ) const {
  unsigned long long value = 0;
  for (int i = n - 1; i >= 0; i--)
    value = (value << 32) + data[i];
  string number = "";
  Number tmp = *this;
  while (tmp) {
    number += '0' + (tmp % 10);
    tmp /= 10;
  }
  reverse (number.begin (), number.end ());
  fprintf (stderr, "[debug] %s = %s\n", name, number.c_str ());
  if (full) {
    for (int i = n - 1; i >= 0; i--)
      fprintf (stderr, "    [data %d] = 0x%08x\n", i, data[i]);
  }
}