test_gmp_c.cpp

test_gmp_c.cpp

/*
  test_gmp_c.cpp

  "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.16.27023\bin\Hostx64\x64\cl.exe"
   -source-charset:utf-8 -execution-charset:utf-8
   -EHsc -Fetest_gmp_c.exe test_gmp_c.cpp
   -I.
   -link
    /LIBPATH:"C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.16.27023\lib\x64"
    /LIBPATH:"C:\Program Files (x86)\Windows Kits\10\Lib\10.0.17763.0\ucrt\x64"
    /LIBPATH:"C:\Program Files (x86)\Windows Kits\10\Lib\10.0.17763.0\um\x64"
    libgmp-10.lib
    libgmpxx-4.lib
    mpir.lib

  del .\test_gmp_c.obj
  test_gmp_c

  gmp limb byte array
  https://lv4.hateblo.jp/entry/2014/01/22/224027
  mpz_t mpq_t mpf_t gmp_randstate_t (3.2 Nomenclature and Types)
  https://gmplib.org/manual/Nomenclature-and-Types
  mpn_ (8 Low-level Functions)
  https://gmplib.org/manual/Low_002dlevel-Functions
  mpz_get_str (5.4 Conversion Functions)
  https://gmplib.org/manual/Converting-Integers
  mp_get_memory_functions
  https://stackoverflow.com/questions/15691477/c-mpir-mpz-t-to-stdstring

  gmp_randinit_default gmp_randseed_ui mpz_urandomb mpz_clear gmp_randclear
  https://qiita.com/ykkhrt/items/ba89fc6e8ef7b76fc401
  how to build (lucasLehmerTest Fast) std::chrono
  https://qiita.com/gis/items/92b7fbf279d3ca6227fe
  functions https://sehermitage.web.fc2.com/etc/gmp.html
  prime factor (cygwin) ulong prec = mpf_get_prec(f); mpf_init2(g, prec);
  https://stdkmd.net/nrr/gmp_ja.htm
  pi FFT AGM https://www.kurims.kyoto-u.ac.jp/~ooura/pi_fft-j.html
  pi Gauss-Legendre algorithm
  1 https://beetreehitsuji.hatenablog.com/entry/2017/01/02/080500 double
  2 https://beetreehitsuji.hatenablog.com/entry/2017/01/02/083408 double repeat
  3 https://beetreehitsuji.hatenablog.com/entry/2017/01/02/121128 gmp fast
  pi https://ja.wikipedia.org/wiki/%E5%86%86%E5%91%A8%E7%8E%87

  (skip)
  gmp.h https://github.com/lattera/glibc/blob/master/stdlib/gmp.h
  gmp.h 621 mingw msys https://packages.msys2.org/package/mingw-w64-x86_64-gmp

  (good)
  mpir-binary https://github.com/ChillMagic/MPIR-Binary (by vc15)
  download MPIR-Binary-master_ChillMagic_38da13f_20180204.zip
   gmp.h gmpxx.h mpir.h mpirxx.h
   mpir.dll (mpir.lib)

  (skip) download mpir-3.0.0.zip
  .../GMP/mpir-3.0.0/build.vc14
  msbuild.bat gc DLL x64 Release -> no exe ?
  edit msbuild.bat (set msbdir="...2017..." from vc15/msbuild.bat)
  msbuild.bat gc DLL x64 Release -> no vc140 ?

  (good)
  download gmp-6.2.1.tar.xz
   (gmp-6.2.1/gmp-h.in gmp-impl.h gmp.pc.in gmpxx.h gmpxx.pc.in)
  gmp.h gmpxx.h (must be build but get binary)
   from MPIR-Binary https://github.com/ChillMagic/MPIR-Binary (by vc15)
  download gimp-2.10.30-setup.exe
  copy libgmp-10.dll libgmpxx-4.dll from ".../GIMP 2/bin"
  to create libgmp-10.lib libgmpxx-4.lib
   export def
   ...\bin\Hostx64\x64\dumpbin.exe /exports libgmp-10.dll > libgmp-10.def
   ...\bin\Hostx64\x64\dumpbin.exe /exports libgmpxx-4.dll > libgmpxx-4.def
   edit libgmp-10.def libgmpxx-4.def
   create lib
   ...\bin\Hostx64\x64\lib.exe /machine:x64 /def:libgmp-10.def
   ...\bin\Hostx64\x64\lib.exe /machine:x64 /def:libgmpxx-4.def
*/

#define UNICODE
#define _UNICODE
#include <wchar.h>

#define _USE_MATH_DEFINES
#include <cmath>

#include <iomanip>
#include <iostream>
#include <sstream>
#include <map>
#include <vector>
#include <string>

#include <stdexcept>
#include <exception>

#include <ctime>
// #include <time.h>
#include <chrono>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// #include <mpirxx.h>
// #include <mpir.h>
// this version gmpxx.h includes mpir.h
#include <gmpxx.h>
#include <gmp.h>

using namespace std;

int test_c()
{
  mpz_t a, b, c, r; // mpz_t mpq_t mpf_t mpn_...(mpz_t, ...)
  mpz_init(a);
  mpz_set_ui(a, 1234567890); // _set_si _set_d
  // mpz_init_set_ui(b, 9876543210); // overflow
  // mpz_set_str(b, "9876543210", 10); // must init before
  mpz_init_set_str(b, "9876543210", 10);
  // mpz_init_set_ui(c, 1234987654321); // overflow
  // mpz_set_str(c, "1234987654321", 10); // must init before
  mpz_init_set_str(c, "1234987654321", 10);
  gmp_printf("a %Zd\n", a);
  gmp_printf("b %Zd\n", b);
  gmp_printf("c %Zd\n", c);

  mpz_set(b, a);
  gmp_printf("b <- a %Zd\n", b);

  mpz_set_str(a, "12394872039846520983745092837458238947528346283745802837468238947520137489572463589", 10);
  mpz_set_str(b, "39846520983745092837458238947528346283745802837468238947520137489572463589", 10);
  mpz_set_str(c, "12394872039846520983745092837458238947528346283745802837468238947520137489572463589", 10);
  gmp_printf("a %Zd\n", a);
  gmp_printf("b %Zd\n", b);
  gmp_printf("c %Zd\n", c);

  // _add _add_ui _sub _sub_ui _mul _mul_ui _mod _mod_ui
  // _tdiv_q _tdiv_q_ui _fdiv_q _fdiv_q_ui _cdiv_q _cdiv_q_ui
  // _neg _abs _fac_ui _mul_2exp _pow_ui _ui_pow_ui _powm _powm_ui
  mpz_add(c, a, b);
  gmp_printf("c %Zd\n", c);

  mpz_init(r); // must init before
  mpz_tdiv_qr(c, r, a, b); // truncate
  gmp_printf("tdiv c %Zd r %Zd\n", c, r);

  mpz_fdiv_r(r, a, b); // floor fdiv r >= 0
  gmp_printf("fdiv r %Zd\n", r);

  mpz_fdiv_q(c, a, b); // floor fdiv c <= tdiv c
  gmp_printf("fdiv c %Zd\n", c);

  mpz_cdiv_r(r, a, b); // ceil cdiv r <= 0
  gmp_printf("cdiv r %Zd\n", r);

  mpz_cdiv_q(c, a, b); // ceil cdiv c >= tdiv c
  gmp_printf("cdiv c %Zd\n", c);

#if 0
  mpz_set_ui(c, 65);
  mpz_out_str(NULL, 10, c); // FILE * == NULL ok, stdout hung up
  mpz_out_raw(NULL, c); // (bin) FILE * == NULL ok, stdout hung up
  FILE *fp = fopen("test_gmp_c_dump.dat", "wb");
  if(fp){
    mpz_out_str(fp, 10, c); // FILE * == fp hung up
    mpz_out_raw(fp, c); // (bin) FILE * == fp hung up
    fclose(fp);
  }
#endif

  mpz_set_ui(c, 65);
  char *s = mpz_get_str(NULL, 10, c);
  if(s){
    fprintf(stdout, "[%s]\n", s);
#if 0
    free(s); // standard free() hungup
#else
    void (*_mp_free)(void *, size_t);
    mp_get_memory_functions(NULL, NULL, &_mp_free); // _alloc, _realloc, _free
    if(_mp_free) _mp_free(s, strlen(s) + 1);
#endif
  }

  mpz_t w, x, y, z;
  mpz_init_set_ui(w, 9);
  mpz_init_set_ui(x, 2);
  mpz_init_set_ui(y, 3);
  mpz_init(z);
  mpz_powm(z, x, y, w);
  gmp_printf("%Zd ^ %Zd = %Zd\n", x, y, z);
  gmp_printf("%Zd is perfect ? : %d\n", w, mpz_perfect_square_p(w));

  mpf_t f, g;
  mpf_init_set_d(f, 5.0);
  mpf_init(g);
  mpf_sqrt(g, f);
  gmp_printf("sqrt(%.*Ff) = %.*Ff\n", 17, f, 17, g);

  mpq_t q;
  mpq_init(q);
  gmp_printf("q hex rational [%#40Qx]\n", q);

#if 0
  mp_limb_t l = mpn ... get ... limb(w);
  gmp_printf("limb %Mu\n", l);
  mp_limb_t *p = &l;
  mp_size_t sz = mpz_size(w);
  gmp_printf("limb array %Nx\n", p, sz);
  mpn ... clear(l);
#endif

  mpq_clear(q);
  mpf_clear(g);
  mpf_clear(f);
  mpz_clear(z);
  mpz_clear(y);
  mpz_clear(x);
  mpz_clear(w);
  mpz_clear(r);
  mpz_clear(c);
  mpz_clear(b);
  mpz_clear(a);
  return 0;
}

int test_cpp()
{
  mpz_class x, y, z;
  x = "1234567890987654321";
  cout << "x = " << x << endl;
  y = 123456789;
  cout << "y = " << y << endl;
  z = x + y;
  cout << "z = x + y = " << z << endl;
  z = x - y;
  cout << "z = x - y = " << z << endl;
  z = x * y;
  cout << "z = x * y = " << z << endl;
  z = x / y;
  cout << "z = x / y = " << z << endl;
  z = x % y;
  cout << "z = x % y = " << z << endl;
  z = x / 1e+9;
  cout << "z = x / 1e+9 = " << z << endl;

  mpf_class f;
  f = x / 1e+9;
  cout << "f = x / 1e+9 = " << f << endl;
  f = 1 / 10.0;
  cout << "f = 1 / 10.0 = " << f << endl;
  f = x / 1e+9;
  cout << "f = x / 1e+9 = " << setprecision(20) << f << endl;
  f = 1 / 10.0;
  cout << "f = 1 / 10.0 = " << setprecision(20) << f << endl;
  return 0;
}

int test_random_mt_c()
{
  gmp_randstate_t state;
  gmp_randinit_default(state);
  time_t t;
  time(&t);
  gmp_randseed_ui(state, t);
  mpz_t r;
  mpz_init(r);
  mpz_urandomb(r, state, 100); // 100bits
  gmp_printf("%Zx\n", r); // hex 25digits
  mpz_clear(r);
  gmp_randclear(state);

  return 0;
}

class genGaussLegendre {
private:
  int precDigits;
  int pd;
  mpf_class a, b, t, p;
public:
  genGaussLegendre(int d) : precDigits(d), pd(d * log2(10)),
    a(1.0, pd), b(1.0 / sqrt(mpf_class{2.0}), pd), t(1.0 / 4.0, pd), p(1.0, pd)
    // initializer must be called with pd before mpf_set_default_prec
  {
    mpf_set_default_prec(pd);
  }
  virtual ~genGaussLegendre(){}
  mpf_class pi(){ return (a + b) * (a + b) / (t * 4.0); }
  void next(){
    mpf_class a_{(a + b) / 2.0};
    mpf_class b_{sqrt(a * b)};
    mpf_class t_{t - p * (a - a_) * (a - a_)};
    mpf_class p_{p * 2.0};
    a = a_, b = b_, t = t_, p = p_;
  }
  int recursion(){ return log2(precDigits); }
  void generate(){ for(int i = recursion(); --i >= 0; ) next(); }
};

class mpPrimeTbl {
private:
  int maxera;
  vector<bool> era;
  vector<mpz_class> tbl;
public:
  mpPrimeTbl(int me) : maxera(me) { init_primes(); }
  virtual ~mpPrimeTbl(){ tbl.clear(); era.clear(); }
  int get_maxera(){ return maxera; }
  bool is_prime(int n){ return (n < 0 || n >= maxera) ? false : !era[n]; }
  size_t nprimes(){ return tbl.size(); }
  mpz_class &nth_prime(int n){ return tbl[n]; }

  size_t init_primes(){
    era.reserve(maxera + 1); // uninitialized
    era[0] = era[1] = true;
    for(int i = 2; i <= maxera; ++i) era[i] = false; // allows [maxera] but NC
    size_t cnt = maxera - 2;
    for(int i = sqrt(maxera); i >= 2; --i)
      for(int j = maxera / i; j >= 2; --j)
        if(!era[i * j]){ era[i * j] = true; --cnt; }
    tbl.reserve(cnt);
    cnt = 0;
    for(int i = 2; i < maxera; ++i)
      if(!era[i]){ mpz_class p; tbl.push_back(p); tbl[cnt++] = i; }
    return tbl.size();
  }
};

int test_prime(mpPrimeTbl &primes)
{
  fprintf(stdout, "%zd primes in %d\n", primes.nprimes(), primes.get_maxera());
#if 0
  for(size_t i = 0; i < primes.nprimes(); ++i)
    gmp_printf(" %Zd", primes.nth_prime(i));
  gmp_printf("\n");
#endif
  return 0;
}

int test_pi(mpPrimeTbl &primes)
{
  mpf_class d = 1.0;
  size_t m = primes.nprimes();
  for(size_t n = 0; n < m; ++n){
    mpz_class &p = primes.nth_prime(n);
    mpf_class q = p * p;
    d *= q / (q - 1);
    mpf_class pi = sqrt(d * 6.0);
    if(!(n % 1000000) || n == m - 1) gmp_printf("pi(%d): %.*Ff\n", n, 17, pi);
  }
  return 0;
}

int main(int ac, char **av)
{
  fprintf(stdout, "sizeof(size_t): %zd\n", sizeof(size_t));

#if 0
  test_c();
#endif
#if 0
  test_cpp();
#endif
#if 1
  test_random_mt_c();
#endif

/* Euler Primes
   maxera  primes pi                  pi[] mpPrimeTbl
       10       4 3.09359216769114542 < 0s <  1s
      100      25 3.13873719507172095 < 0s <  1s
     1000     168 3.14139318478792820 < 0s <  1s
    10000    1229 3.14157723516034506 < 0s <  1s
   100000    9592 3.14159139342024467 < 0s <  1s
  1000000   78498 3.14159254712798260 < 1s <  1s   8MB ?
 10000000  664579 3.14159264438479261 < 2s <  7s  80MB ?
100000000 5761455 3.14159265277861852 < 8s < 82s 800MB ?
*/
  chrono::system_clock::time_point start, end;
  start = chrono::system_clock::now();
  mpPrimeTbl primes(1000000);
  test_prime(primes);
  test_pi(primes);
  end = chrono::system_clock::now();
  cout << chrono::duration_cast<chrono::milliseconds>(end - start).count();
  cout << endl;

/* Gauss Legendre
 digits    prec re. pi
     10      64   3 <    1ms
    100     384   6 <    1ms
    200     704   7 <    1ms
   1000    3328   9 <    1ms
  10000   33280  13 <   10ms
 100000  332224  16 <  130ms
1000000 3321984  19 < 2540ms
*/
  int digits = 1000;
  genGaussLegendre gen(digits);
  gen.generate();
  end = chrono::system_clock::now();
  gmp_printf("%.*Ff\n(recursion: %d)\n", digits, gen.pi(), gen.recursion());
  fprintf(stdout, "(precision: %lld)\n", mpf_get_default_prec());
  cout << chrono::duration_cast<chrono::milliseconds>(end - start).count();
  cout << endl;

  return 0;
}