libxml2 nanohttp, libcurl, Boost Asio http

by taiji at 2016/06/24 20:29

C/C++ でウェブページの取得をする際、libxml2 ライブラリに含まれている nanohttp は、HTTPS にも圧縮伸長にも対応していないので、結局、実用的ではありません。その上、ドキュメントが充実していないためか、ネット情報にて間違った使用方法が散見されましたので、以下に例示しておくことにしました。ちなみに、nanohttpLocation ヘッダによるリダイレクションに対応していますが、URL に対する相対 Location には対応していないことに注意です。

/*
  nanohtget.cc

  Copyright (C) 2016 Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <libxml/nanohttp.h>

int main(int argc, char *argv[])
{
  int rv = 0;

  if (argc < 2) {
    std::cerr << "usage:\n\t" << argv[0] << " url [filename.out]" << std::endl;
    return 1;
  }
  else if (2 < argc) {
    char *content_type = NULL;

    xmlNanoHTTPInit();
    rv = xmlNanoHTTPFetch(argv[1], argv[2], &content_type);
    std::cout << "Content-Type:\t" << (content_type ? content_type : "") << std::endl;
    std::cout << "Return-Value:\t" << rv << std::endl;
    if (content_type) free(content_type);
    xmlNanoHTTPCleanup();
  }
  else /*if (1 < argc) */{
    void *handler;
    char *content_type = NULL;

    xmlNanoHTTPInit();
    if (!(handler = xmlNanoHTTPOpen(argv[1], &content_type))) {
      std::cerr << "xmlNanoHTTPOpen: failure" << std::endl;
      if (content_type) free(content_type);
      return 1;
    }
    if (xmlNanoHTTPReturnCode(handler) != 200) {
      if (content_type) free(content_type);
      return 1;
    }

    std::cerr <<
      "Location:\t"             << (xmlNanoHTTPRedir(handler) ? xmlNanoHTTPRedir(handler) : argv[1]) << std::endl <<
      "Content-Type:\t"         << (content_type ? content_type : "") <<
      (xmlNanoHTTPEncoding(handler) ? std::string("; charset=") + xmlNanoHTTPEncoding(handler) : "") << std::endl <<
      "Mime-Type:\t"            << (xmlNanoHTTPMimeType(handler) ? xmlNanoHTTPMimeType(handler) : "") << std::endl <<
      "Content-Length:\t"       << xmlNanoHTTPContentLength(handler) << std::endl <<
      "WWW-Authenticate:\t"     << (xmlNanoHTTPAuthHeader(handler) ? xmlNanoHTTPAuthHeader(handler) : "") << std::endl;

    char buf[1024];
    int sz;
    while ((sz = xmlNanoHTTPRead(handler, buf, sizeof(buf) - 1)) > 0) {
      buf[sz] = '\0';
      std::cout << buf;
    }

    if (sz == 0)
      std::cerr << "xmlNanoHTTPRead: Connection closed" << std::endl;
    else if (sz == -1) {
      std::cerr << "xmlNanoHTTPRead: Parameter Error" << std::endl;
      rv = 1;
    }

    if (content_type) free(content_type);
    xmlNanoHTTPClose(handler);
    xmlNanoHTTPCleanup();
  }
  return rv;
}

例えば、以下のようにビルドすればよいでしょう。

$ make CPPFLAGS=-I/usr/include/libxml2 LDLIBS=-lxml2 nanohtget
g++ -I/usr/include/libxml2 nanohtget.cc -lxml2 -o nanohtget

間違えやすいのは xmlNanoHTTPRead にて常にナル終端されるわけではないので、読み込んだバッファを C の文字列として扱うには読み込んだサイズに応じてナル終端する必要があることです。ナル終端された行毎に読み込む xmlNanoHTTPReadLine もあるようですので、そちらを使った方がよいかも知れません。

さて、HTTPS に対応したものとして、libcurl が挙げられます。これを使うと流石に簡単に C/C++ でウェブページの取得が出来ます。やはりこれも、適切なネット情報がないようなので、以下に例示をすることにしました。

/*
  urlfetch.cc

  Copyright (C) 2016 Taiji Yamada <taiji@aihara.co.jp>
*/
#include <cstdio>
#include <iostream>
#include <curl.h>       /* must not be <curl/curl.h> */

int main(int argc, char *argv[])
{
  CURL *curl;
  CURLcode rc;
  FILE *fp = stdout;

  if (argc < 2) {
    std::cerr << "usage:\n\t" << argv[0] << " url [filename.out]" << std::endl;
    return 1;
  }
  else if (2 < argc) {
    if (!(fp = fopen(argv[2], "w"))) {
      std::cerr << "error: fopen() " << argv[2] << std::endl;
      return 1;
    }
  }

  if (!(curl = curl_easy_init())) {
    std::cerr << "error: curl_easy_init()" << std::endl;
    return 1;
  }
  curl_easy_setopt(curl, CURLOPT_URL, argv[1]);
  curl_easy_setopt(curl, CURLOPT_FOLLOWLOCATION, 1);
  curl_easy_setopt(curl, CURLOPT_WRITEDATA, fp);
  curl_easy_setopt(curl, CURLOPT_WRITEHEADER, stderr);

  rc = curl_easy_perform(curl);

  std::cerr << "curl_easy_perform: " << rc << std::endl;

  curl_easy_cleanup(curl);
  if (fp != stdout) fclose(fp);

  return rc == 0 ? 0 : 1;
}

例えば、以下のようにビルドすればよいでしょう。

$ make CPPFLAGS=-I/usr/include/curl LDLIBS=-lcurl urlfetch
g++ -I/usr/include/curl urlfetch.cc -lcurl -o urlfetch

libcurl のインストール状態にも依りますが、このまま HTTPS にも対応しているはずです。

また、リダイレクションは CURLOPT_FOLLOWLOCATION1 にセットしないと有効になりませんので、上記のようになっています。さらに、出力先として上記のようにファイルポインタを指定できます。

次に、Boost Asio を使ったウェブページを取得する方法ですが、まず、URL からホスト名とスキームその他にばらすものがあれば便利なので、まずそれを紹介します。以下の url.hpp です。

/*
  url.hpp

  Copyright (C) 2016 Taiji Yamada <taiji@aihara.co.jp>
*/
#ifndef _URL_HPP_
#define _URL_HPP_
#include <boost/xpressive/xpressive.hpp>

struct url {
  url() {}
  url(const char *arg)
  {
    const boost::xpressive::sregex re =
      !((boost::xpressive::s1 = +~(boost::xpressive::set='/', ':')) >> "://" >>
        !((boost::xpressive::s2 = +~(boost::xpressive::set=':', '@')) >>
          !(':' >> (boost::xpressive::s3 = +~(boost::xpressive::set='@'))) >> '@'
          ) >>
        (boost::xpressive::s4 = *~(boost::xpressive::set='/', ':')) >>
        !(':' >> (boost::xpressive::s5 = *~(boost::xpressive::set='.', '/', ':')))
        ) >>
      (boost::xpressive::s6 = *~boost::xpressive::_n);
    boost::xpressive::smatch what;
    if (boost::xpressive::regex_match(std::string(arg), what, re)) {
      if (what[1]) scheme       = what[1];
      if (what[2]) user         = what[2];
      if (what[3]) passwd       = what[3];
      if (what[4]) host         = what[4];
      if (what[5]) port         = what[5];
      if (what[6]) path         = what[6];
    }
  }
  std::string scheme, host, port, path;
  std::string user, passwd;
  bool internal() const
  {
    return host == "";
  }
  bool externalizable(const url &base_url)
  {
    return internal() && !base_url.internal();
  }
  const url external(const url &base_url)
  {
    url rv(*this);
    if (rv.scheme.empty())      rv.scheme       = base_url.scheme;
    if (rv.user.empty())        rv.user         = base_url.user;
    if (rv.passwd.empty())      rv.passwd       = base_url.passwd;
    if (true)                   rv.host         = base_url.host;
    if (rv.port.empty())        rv.port         = base_url.port;
    std::string base_path = base_url.path;
    size_t i;
    if ((i=base_path.rfind('/')) != std::string::npos)
      base_path = base_path.substr(0, i+1);
    if (rv.path.empty())
      rv.path = base_path;
    else if (rv.path[0] != '/')
      rv.path = base_path + rv.path;
    return rv;
  }
  bool connection_equals(const url &another) const
  {
    return scheme == another.scheme && host == another.host && port == another.port;
  }
  const std::string locator() const
  {
    return (scheme.empty() ? "" : scheme + "://") +
      (user.empty() ? "" : user + (passwd.empty() ? "" : ":" + passwd) + "@") +
      host + (port.empty() ? "" : ":" + port) + path;
  }
};

#endif

Boost Xpressive を使っています。便利です。

さて、ここでは同期通信に限ることにしますが、出来るだけ分かり易く書くと以下のようになりました。

/*
  htget.cc

  Written by Taiji Yamada <taiji@aihara.co.jp>

  Reference(s):
  [1] http://www.boost.org/doc/libs/1_39_0/doc/html/boost_asio/example/http/client/sync_client.cpp
*/
#include <iostream>
#include <boost/asio.hpp>
#include "url.hpp"

int main(int argc, char *argv[])
{
  int rv = 0;

  if (argc < 2) {
    std::cerr << "usage:\n\t" << argv[0] << " url" << std::endl;
    return 1;
  }
  try {
    boost::asio::io_service io_service;
    boost::asio::ip::tcp::socket socket(io_service);

    url location(argv[1]);
    boost::asio::ip::tcp::resolver resolver(io_service);
    boost::asio::ip::tcp::resolver::query query((location.host != "") ? location.host : "localhost",
                                                (location.port != "") ? location.port : location.scheme);
    boost::asio::ip::tcp::endpoint endpoint(*resolver.resolve(query));

    socket.connect(endpoint);

    boost::asio::streambuf request;
    std::ostream request_ostream(&request);
    request_ostream << "GET " << location.path << " HTTP/1.0\r\n"
      "\r\n";
    boost::asio::write(socket, request);

    boost::asio::streambuf response;
    boost::system::error_code error_code;
    while (boost::asio::read(socket, response, boost::asio::transfer_at_least(1), error_code))
      std::cout << &response;

    socket.close();

    if (error_code != boost::asio::error::eof && error_code != boost::asio::error::shut_down)
      throw boost::system::system_error(error_code);
  }
  catch (std::exception &e) {
    std::cout << argv[0] << ": " << e.what() << std::endl;
    rv = 1;
  }
  return rv;
}

boost が /opt/local にインストールされているとして、ビルドは以下のようにします。

$ make CPPFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib LDLIBS='-lboost_system' htget
g++ -I/opt/local/include -L/opt/local/lib htget.cc -lboost_system -o htget

しかし、これでは HTTPS に対応していないので、実用的ではありません。HTTPS に対応したものを Unified diff 形式で表すと以下のようになります。

$ diff -u htget.cc s-htget.cc
--- htget.cc?2016-06-30 16:18:50.000000000 +0900
+++ s-htget.cc?2016-07-01 15:37:05.000000000 +0900
@@ -1,13 +1,15 @@
 /*
-  htget.cc
+  s-htget.cc
 
   Written by Taiji Yamada <taiji@aihara.co.jp>
 
   Reference(s):
   [1] http://www.boost.org/doc/libs/1_39_0/doc/html/boost_asio/example/http/client/sync_client.cpp
+  [2] http://www.boost.org/doc/libs/1_39_0/doc/html/boost_asio/example/ssl/client.cpp
 */
 #include <iostream>
 #include <boost/asio.hpp>
+#include <boost/asio/ssl.hpp>
 #include "url.hpp"
 
 int main(int argc, char *argv[])
@@ -20,7 +22,8 @@
   }
   try {
     boost::asio::io_service io_service;
-    boost::asio::ip::tcp::socket socket(io_service);
+    boost::asio::ssl::context context(io_service, boost::asio::ssl::context::sslv23_client);
+    boost::asio::ssl::stream<boost::asio::ip::tcp::socket> socket(io_service, context);
 
     url location(argv[1]);
     boost::asio::ip::tcp::resolver resolver(io_service);
@@ -28,7 +31,10 @@
                                                (location.port != "") ? location.port : location.scheme);
     boost::asio::ip::tcp::endpoint endpoint(*resolver.resolve(query));
 
-    socket.connect(endpoint);
+    context.set_verify_mode(boost::asio::ssl::context::verify_peer);
+    context.load_verify_file("cacert.pem");
+    socket.lowest_layer().connect(endpoint);
+    socket.handshake(boost::asio::ssl::stream_base::client);
 
     boost::asio::streambuf request;
     std::ostream request_ostream(&request);
@@ -41,7 +47,7 @@
     while (boost::asio::read(socket, response, boost::asio::transfer_at_least(1), error_code))
       std::cout << &response;
 
-    socket.close();
+    socket.lowest_layer().close();
 
     if (error_code != boost::asio::error::eof && error_code != boost::asio::error::shut_down)
       throw boost::system::system_error(error_code);

boost が /opt/local にインストールされているとして、ビルドは以下のようにします。

$ make CPPFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib LDLIBS='-lssl -lcrypto -lboost_system' s-htget
g++ -I/opt/local/include -L/opt/local/lib s-htget.cc -lssl -lcrypto -lboost_system -o s-htget

しかし、これでは HTTP, HTTPS の双方には対応しておらず、リダイレクションにも未対応で、実用には程遠いと言わざるを得ません。教育的価値と言っても最低限、以下のような例示まであると、それなりに使えるものとなると思います。

/*
  shtget.cc

  Written by Taiji Yamada <taiji@aihara.co.jp>

  Reference(s):
  [1] http://www.boost.org/doc/libs/1_39_0/doc/html/boost_asio/example/http/client/sync_client.cpp
  [2] http://www.boost.org/doc/libs/1_39_0/doc/html/boost_asio/example/ssl/client.cpp
*/
#include <iostream>
#include <boost/asio.hpp>
#include <boost/asio/ssl.hpp>
#include <boost/xpressive/xpressive.hpp>
#include "url.hpp"

int main(int argc, char *argv[])
{
  int rv = 0;

  if (argc < 2) {
    std::cerr << "usage:\n\t" << argv[0] << " url" << std::endl;
    return 1;
  }
  try {
    boost::asio::io_service io_service;
    boost::asio::ip::tcp::socket socket(io_service);
    boost::asio::ssl::context context(io_service, boost::asio::ssl::context::sslv23_client);
    boost::asio::ssl::stream<boost::asio::ip::tcp::socket> secure_socket(io_service, context);

    url location(argv[1]);
    const int redirection_limit = 4;
    int redirection_count = 0;
    bool redirection;
    do {
      redirection = false;
      boost::asio::ip::tcp::resolver resolver(io_service);
      boost::asio::ip::tcp::resolver::query query((location.host != "") ? location.host : "localhost",
                                                  (location.port != "") ? location.port : location.scheme);
      boost::asio::ip::tcp::endpoint endpoint(*resolver.resolve(query));

      if (location.scheme != "https")
        socket.connect(endpoint);
      else {
        context.set_verify_mode(boost::asio::ssl::context::verify_peer);
        context.load_verify_file("cacert.pem");
        secure_socket.lowest_layer().connect(endpoint);
        secure_socket.handshake(boost::asio::ssl::stream_base::client);
      }

      boost::asio::streambuf request;
      std::ostream request_ostream(&request);
      request_ostream << "GET " << location.path << " HTTP/1.0\r\n"
        "\r\n";
      if (location.scheme != "https")
        boost::asio::write(socket, request);
      else
        boost::asio::write(secure_socket, request);

      boost::asio::streambuf response;
      boost::system::error_code error_code;
      std::string buffer;
      while ((location.scheme != "https") ?
             boost::asio::read(socket, response, boost::asio::transfer_at_least(1), error_code) :
             boost::asio::read(secure_socket, response, boost::asio::transfer_at_least(1), error_code))
        ;
      buffer = std::string(boost::asio::buffer_cast<const char *>(response.data())).substr(0, response.size());
      std::cout << buffer;

      if (location.scheme != "https")
        socket.close();
      else
        secure_socket.lowest_layer().close();

      if (error_code != boost::asio::error::eof && error_code != boost::asio::error::shut_down)
        throw boost::system::system_error(error_code);

      const boost::xpressive::sregex re =
        boost::xpressive::bos >> "HTTP/1.1" >> +boost::xpressive::_s >> (boost::xpressive::s1 = +boost::xpressive::_d) >> +boost::xpressive::_s >> *~boost::xpressive::_n >> +boost::xpressive::_n >>
        *boost::xpressive::_ >> boost::xpressive::after(boost::xpressive::_n) >> "Location:" >> +boost::xpressive::_s >> (boost::xpressive::s2 = *~boost::xpressive::_n) >> +boost::xpressive::_n >>
        *boost::xpressive::_ >>
        boost::xpressive::eos;
      boost::xpressive::smatch what;
      if (boost::xpressive::regex_match(buffer, what, re) &&
          (std::string(what[1]) == "301" || std::string(what[1]) == "302" || std::string(what[1]) == "303" ||
           std::string(what[1]) == "305" || std::string(what[1]) == "307")) {
        url redirect(std::string(what[2]).c_str());
        location = redirect.externalizable(location) ? redirect.external(location) : redirect;
        redirection = true;
      }
    } while (redirection && redirection_count++ < redirection_limit);
  }
  catch (std::exception &e) {
    std::cout << argv[0] << ": " << e.what() << std::endl;
    rv = 1;
  }
  return rv;
}

boost が /opt/local にインストールされているとして、ビルドは以下のようにします。

$ make CPPFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib LDLIBS='-lssl -lcrypto -lboost_system' shtget
g++ -I/opt/local/include -L/opt/local/lib shtget.cc -lssl -lcrypto -lboost_system -o shtget

気になる点としては、同ホストへのリダイレクションでも再接続する無駄があることです。例示としてはさらに複雑になってしまうので書きませんが、そういう時はソケットを close, connect せずにそのまま使うように改良することもできます。


Boost MPL bitwise not

by taiji at 2016/04/22 15:14

Boost MPL を読んでいて気になったのですが、算術演算の「ビット否定」だけ、どう探しても実装されてないと思います。とは言え、自前で用意するのは非常に簡単。

namespace boost { namespace mpl {

  template <typename T>
  struct compl_ : integral_c<typename T::value_type, ~T::value> {};

} }

こういったものを含む boost/mpl/compl.hpp を用意すればよいでしょう。確認の仕方は、

  BOOST_MPL_ASSERT_RELATION((boost::mpl::compl_<boost::mpl::integral_c<unsigned, 0> >::value), ==, ~0U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::compl_<boost::mpl::integral_c<unsigned, ~0U> >::value), ==, 0);

といった感じです。ちなみに、実装されている他のビット演算の確認の仕方は、

#include <boost/mpl/integral_c.hpp>
#include <boost/mpl/shift_left.hpp>
#include <boost/mpl/shift_right.hpp>
#include <boost/mpl/bitand.hpp>
#include <boost/mpl/bitxor.hpp>
#include <boost/mpl/bitor.hpp>
#include <boost/mpl/and.hpp>
#include <boost/mpl/or.hpp>
#include <boost/mpl/assert.hpp>

必要なヘッダファイルをインクルードした上で、

  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitand_<boost::mpl::integral_c<unsigned, 0>, boost::mpl::integral_c<unsigned, 0> >::value), ==, 0U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitand_<boost::mpl::integral_c<unsigned, 0>, boost::mpl::integral_c<unsigned, 1> >::value), ==, 0U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitand_<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 0> >::value), ==, 0U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitand_<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 1> >::value), ==, 1U);

  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitxor_<boost::mpl::integral_c<unsigned, 0>, boost::mpl::integral_c<unsigned, 0> >::value), ==, 0U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitxor_<boost::mpl::integral_c<unsigned, 0>, boost::mpl::integral_c<unsigned, 1> >::value), ==, 1U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitxor_<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 0> >::value), ==, 1U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitxor_<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 1> >::value), ==, 0U);

  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitor_<boost::mpl::integral_c<unsigned, 0>, boost::mpl::integral_c<unsigned, 0> >::value), ==, 0U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitor_<boost::mpl::integral_c<unsigned, 0>, boost::mpl::integral_c<unsigned, 1> >::value), ==, 1U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitor_<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 0> >::value), ==, 1U);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::bitor_<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 1> >::value), ==, 1U);

  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_left<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 0> >::value), ==, 1U<<0);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_left<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 1> >::value), ==, 1U<<1);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_left<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 2> >::value), ==, 1U<<2);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_left<boost::mpl::integral_c<unsigned, 1>, boost::mpl::integral_c<unsigned, 3> >::value), ==, 1U<<3);

  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<unsigned, ~0U>, boost::mpl::integral_c<unsigned, 0> >::value), ==, ~0U>>0);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<unsigned, ~0U>, boost::mpl::integral_c<unsigned, 1> >::value), ==, ~0U>>1);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<unsigned, ~0U>, boost::mpl::integral_c<unsigned, 2> >::value), ==, ~0U>>2);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<unsigned, ~0U>, boost::mpl::integral_c<unsigned, 3> >::value), ==, ~0U>>3);

  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<signed, ~0>, boost::mpl::integral_c<unsigned, 0> >::value), ==, ~0>>0);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<signed, ~0>, boost::mpl::integral_c<unsigned, 1> >::value), ==, ~0>>1);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<signed, ~0>, boost::mpl::integral_c<unsigned, 2> >::value), ==, ~0>>2);
  BOOST_MPL_ASSERT_RELATION((boost::mpl::shift_right<boost::mpl::integral_c<signed, ~0>, boost::mpl::integral_c<unsigned, 3> >::value), ==, ~0>>3);

というような感じです。


オンライン版各種平均テンプレートライブラリ

by taiji at 2015/12/11 18:53

オンライン版各種平均テンプレートライブラリ(後述するオンライン版各種統計量テンプレートライブラリ statistics_online.hpp などを含みます)の averages_online.hpp を作成しました。これは平均値を一つ前の平均値から求めるもので、リアルタイムで効率的に各種平均を求めるのに適しています。まず、通常どおりに各種平均を求めるための式は以下の通りです。

相加平均(arithmetic average):

\[\mu_A = \dfrac1n\sum_{i=1}^n x_i\]

相乗平均(geometric average):

\[\mu_G = \sqrt[n]{\prod_{i=1}^n x_i} = \exp\left( \dfrac1n\sum_{i=1}^n \log x_i \right)\]

調和平均(harmonic average):

\[\mu_H = \dfrac n{\sum_{i=1}^n \dfrac1{x_i}}\]

一般化平均(generalized average):

\[\mu_{G(m)} = \sqrt[m]{\dfrac1n\sum_{i=1}^n x_i^m}\]

任意の一般化平均(arbitrary generalized average):

\[\mu_{f} = f^{-1}\left(\dfrac1n\sum_{i=1}^n f(x_i)\right)\]

以上より、一つ前の平均値から求めるには、以下の差分方程式に依ります。

相加平均の差分方程式(arithmetic average online):

\[\mu_A{(n)} = \dfrac1n\left\{ (n-1)\mu_A{(n-1)} + x_n \right\}\]

相乗平均の差分方程式(geometric average online):

\[\mu_G{(n)} = \sqrt[n]{\mu_G{(n-1)}^{n-1} x_n} = \exp\left( \dfrac1n\left\{ (n-1)\log(\mu_G{(n-1)}) + \log x_n \right\} \right)\]

調和平均の差分方程式(harmonic average online):

\[\mu_H{(n)} = \dfrac n{\dfrac{n-1}{\mu_H{(n-1)}} + \dfrac1{x_n}}\]

一般化平均の差分方程式(generalized average online):

\[\mu_{G(m)}{(n)} = \sqrt[m]{\dfrac1n\left\{ (n-1)\mu_{G(m)}{(n-1)}^m + x_n^m \right\}}\]

任意の一般化平均の差分方程式(arbitrary generalized average online):

\[\mu_{f}{(n)} = f^{-1}\left({\dfrac1n\left\{ (n-1)f\left(\mu_{f}{(n-1)}\right) + f(x_n) \right\}}\right)\]

また、もともと差分方程式で実装されるべき以下の平均もあります。

指数平滑平均(exponential smoothed average)の差分方程式:

\[\mu_E{(n)} = (1 - \alpha)\mu_E{(n-1)} + \alpha x_n\]

さらに、各種の移動平均は、幅\(T\)の記憶域から除去される最古の\(x_{n-T}\)と追加される最新の\(x_n\)および一つ前の移動平均値から求めることが出来、以下の差分方程式に依ります。

相加移動平均の差分方程式(moving arithmetic online):

\[\mu_{MA}{(n)} = \mu_{MA}{(n-1)} - \dfrac1Tx_{n-T} + \dfrac1Tx_n\]

相乗移動平均の差分方程式(moving geometric online):

\[\mu_{MG}{(n)} = \sqrt[T]{\mu_{MG}{(n-1)}^T\dfrac{x_n}{x_{n-T}}} = \exp\left(\log(\mu_{MG}{(n-1)}) - \dfrac1T\log(x_{n-T}) + \dfrac1T\log(x_n)\right)\]

調和移動平均の差分方程式(moving harmonic online):

\[\mu_{MH}{(n)} = \dfrac1{\dfrac1{\mu_{MH}{(n-1)}} - \dfrac1{Tx_{n-T}} + \dfrac1{Tx_n}}\]

一般化移動平均の差分方程式(moving generalized online):

\[\mu_{MG(m)}{(n)} = \sqrt[m]{\mu_{MG(m)}{(n-1)}^m - \dfrac1Tx_{n-T}^m + \dfrac1Tx_n^m}\]

任意の一般化移動平均(moving arbitrary generalized online)

\[\mu_{Mp}{(n)} = f^{-1}\left({f\left(\mu_{Mp}{(n-1)}\right) - \dfrac1Tf(x_{n-T}) + \dfrac1Tf(x_n)}\right)\]

ちなみに、その tarball には上記の検証の為に、任意のコンテナ(配列、vector, list 等々)で通常どおりに以下の平均を求める average.hpp も用意しました。

  • 相加平均(arithmetic average)
  • 相乗平均(geometric average)
  • 調和平均(harmonic average)
  • 一般化平均(generalized average)
  • 任意の一般化平均(arbitrary generalized average)
  • 指数平滑平均(exponential smoothed average)

これら非オンラインとオンラインとの自乗誤差は 100,000 データにおいて 10E-13 程度に収まっていることは確認できています。

使用例として、その tarball にあるコードを以下に示しておきます。

/*
  averages.cc

  Written by Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <vector>
#include <list>
#include <functional>
#include <cstdlib>      // random(), etc
#include <cassert>
#include "averages.hpp"
#include "averages_online.hpp"

int main()
{
  const double maximum_permissible_error(1e-14);
  const char *result_names[] = {
    "arithmetic mean",
    "geometric mean",
    "harmonic mean",
    "generalized mean",
    "exponential smoothed average",
  };
  std::vector<double> v(1000);
  generic_average_online<double> *online_averages[] = {
    new arithmetic_average_online<double>(),
    new geometric_average_online<double>(),
    new harmonic_average_online<double>(),
    new generalized_average_online<double>(),
    new exponential_smoothed_average_online<double>(),
  };
  srandomdev();
  for (size_t i=0; i<v.size(); ++i) {
    double x((double)random()/((1L<<31)-1));
    v[i] = x;
    for (size_t j=0; j<sizeof(online_averages)/sizeof(online_averages[0]); ++j)
      online_averages[j]->average(x);
  }
  double results[] = {
    arithmetic_average(v.begin(), v.end()),
    geometric_average(v.begin(), v.end()),
    harmonic_average(v.begin(), v.end()),
    generalized_average(v.begin(), v.end()),
    exponential_smoothed_average(v.begin(), v.end()),
  };
  for (size_t j=0; j<sizeof(online_averages)/sizeof(online_averages[0]); ++j) {
    double error(sqrt(pow(online_averages[j]->get_value() - results[j], 2)));
    std::cout << result_names[j] << ":\t" << online_averages[j]->get_value() << "\t" << results[j] << "\t" << error << std::endl;
    assert(error < maximum_permissible_error);
  }

  const char *moving_result_names[] = {
    "moving arithmetic average",
    "moving geometric average",
    "moving harmonic average",
    "moving generalized average",
  };
  std::vector<double> m(v.size());
  generic_average_online<double> *online_moving_averages[] = {
    new moving_arithmetic_average_online<double>(v.size()),
    new moving_geometric_average_online<double>(v.size()),
    new moving_harmonic_average_online<double>(v.size()),
    new moving_generalized_average_online<double>(v.size()),
  };
  for (size_t i=0; i<v.size()*10; ++i) {
    double x((double)random()/((1L<<31)-1));
    m[i % v.size()] = x;
    for (size_t j=0; j<sizeof(online_moving_averages)/sizeof(online_moving_averages[0]); ++j)
      online_moving_averages[j]->average(x);
    if (i < v.size()) continue;
    double moving_results[] = {
      arithmetic_average(m.begin(), m.end()),
      geometric_average(m.begin(), m.end()),
      harmonic_average(m.begin(), m.end()),
      generalized_average(m.begin(), m.end()),
    };
    for (size_t j=0; j<sizeof(online_moving_averages)/sizeof(online_moving_averages[0]); ++j) {
      double error(sqrt(pow(online_moving_averages[j]->get_value() - moving_results[j], 2)));
      if (i+1 == v.size()*10)
        std::cout << moving_result_names[j] << ":\t" << online_moving_averages[j]->get_value() << "\t" << moving_results[j] << "\t" << error << std::endl;
      assert(error < maximum_permissible_error);
    }
  }
  return 0;
}

任意の一般化平均を追加したので、上例にて一致する箇所を入れ替える形で、その使用例を以下に示しておきます。

template <typename T>
struct identity_function : std::unary_function<T, T> {
  T operator()(const T &x) const
  { return x; }
};

template <typename T>
struct reciprocal_function : std::unary_function<T, T> {
  T operator()(const T &x) const
  { return T(1)/x; }
};

template <typename T>
struct square_function : std::unary_function<T, T> {
  T operator()(const T &x) const
  { return x*x; }
};

#define USE_ARBITRARY_GENERALIZED_AVERAGE
#define USE_ARBITRARY_GENERALIZED_AVERAGE_ONLINE

  generic_average_online<double> *online_averages[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE_ONLINE)
    new arithmetic_average_online<double>(),
    new geometric_average_online<double>(),
    new harmonic_average_online<double>(),
    new generalized_average_online<double>(),
#else
    new arbitrary_generalized_average_online<double>(identity_function<double>(), identity_function<double>()),
    new arbitrary_generalized_average_online<double>(log, exp),
    new arbitrary_generalized_average_online<double>(reciprocal_function<double>(), reciprocal_function<double>()),
    new arbitrary_generalized_average_online<double>(square_function<double>(), sqrt),
#endif
    new exponential_smoothed_average_online<double>(),
  };

  double results[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE)
    arithmetic_average(v.begin(), v.end()),
    geometric_average(v.begin(), v.end()),
    harmonic_average(v.begin(), v.end()),
    generalized_average(v.begin(), v.end()),
#else
    generalized_average(v.begin(), v.end(), identity_function<double>(), identity_function<double>()),
    generalized_average(v.begin(), v.end(), log, exp),
    generalized_average(v.begin(), v.end(), reciprocal_function<double>(), reciprocal_function<double>()),
    generalized_average(v.begin(), v.end(), square_function<double>(), sqrt),
#endif
    exponential_smoothed_average(v.begin(), v.end()),
  };

  generic_average_online<double> *online_moving_averages[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE_ONLINE)
    new moving_arithmetic_average_online<double>(v.size()),
    new moving_geometric_average_online<double>(v.size()),
    new moving_harmonic_average_online<double>(v.size()),
    new moving_generalized_average_online<double>(v.size()),
#else
    new moving_arbitrary_generalized_average_online<double>(v.size(), identity_function<double>(), identity_function<double>()),
    new moving_arbitrary_generalized_average_online<double>(v.size(), log, exp),
    new moving_arbitrary_generalized_average_online<double>(v.size(), reciprocal_function<double>(), reciprocal_function<double>()),
    new moving_arbitrary_generalized_average_online<double>(v.size(), square_function<double>(), sqrt),
#endif
  };

    double moving_results[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE)
      arithmetic_average(m.begin(), m.end()),
      geometric_average(m.begin(), m.end()),
      harmonic_average(m.begin(), m.end()),
      generalized_average(m.begin(), m.end()),
#else
      generalized_average(m.begin(), m.end(), identity_function<double>(), identity_function<double>()),
      generalized_average(m.begin(), m.end(), log, exp),
      generalized_average(m.begin(), m.end(), reciprocal_function<double>(), reciprocal_function<double>()),
      generalized_average(m.begin(), m.end(), square_function<double>(), sqrt),
#endif

このように、関数ポインタや関数オブジェクトを指定することで、任意の一般化平均を求めることができるようになっています。但し、averages_online.hpp の方は C++11 または boost::function が必須となります。

さらに、statistics.hpp には、分散、標準偏差、絶対偏差、歪度、尖度などの各種統計量を、任意のコンテナ(配列、vector, list 等々)で以下のように求めるアルゴリズムが実装されています。

      unbiased_variance(v.begin(), v.end()),
      unbiased_standard_deviation(v.begin(), v.end()),
      sample_variance(v.begin(), v.end()),
      sample_standard_deviation(v.begin(), v.end()),
      absolute_deviation(v.begin(), v.end()),
      skewness(v.begin(), v.end()),
      kurtosis(v.begin(), v.end()),
      lag_one_autocorrelation(v.begin(), v.end()),

これらのうちいくつかの統計量についてのオンライン版のアルゴリズムはこの文献によって与えられています。そこで使われている式変形について多少補足しておきます。

平均:

\[\mu=E[X]\]

2次のモーメント:

\[\mu_2=E\left[\left(X-\mu\right)^2\right]=E\left[X^2-2X\mu+\mu^2\right]=E[X^2]-2E[X]\mu+\mu^2=E[X^2]-\mu^2=E[X^2]-E[X]^2=\sigma^2\]

3次のモーメント:

\[\mu_3=E\left[\left(X-\mu\right)^3\right]=E\left[X^3-3X^2\mu+3X\mu^2-\mu^3\right]=E[X^3]-3\mu(E[X^2]-E[X]^2)-\mu^3=E[X^3]-3\mu\mu_2-\mu^3\]

歪度:

\[\gamma=E\left[\left(\dfrac{X-\mu}{\sigma}\right)^3\right]=E\left[\dfrac{(X-\mu)^3}{\sigma^3}\right]=\dfrac{E[(X-\mu)^3]}{\sigma^3}=\dfrac{\mu_3}{\sigma^3}=\dfrac{E[X^3]-3\mu\mu_2-\mu^3}{\mu_2^{3/2}}\]

4次のモーメント:

\[\mu_4=E\left[\left(X-\mu\right)^4\right]=E\left[X^4-4X^3\mu+6X^2\mu^2-4X\mu^3+\mu^4\right]=E[X^4]-4E[X^3]\mu+6E[X^2]\mu^2-3\mu^4=\cdots=E[X^4]-4\mu\mu_3-6\mu^2\mu_2-\mu^4%\]

尖度:

\[\kappa=\dfrac{\mu_4}{\mu_2^2}-3=\dfrac{E\left[\left(X-\mu\right)^4\right]}{{\sigma^2}^2}-3=\dfrac{E[X^4]-4E[X^3]\mu+6E[X^2]\mu^2-3\mu^4}{\sigma^4}-3=\dfrac{E[X^4]-4\mu\mu_3-6\mu^2\mu_2-\mu^4}{\sigma^4}-3\]

さらにこれらの式により、移動平均と同様に移動分散、移動歪度、移動尖度などを求めることが出来ます。

それが実装された statistics_online.hpp の使用例として、その tarball にあるコードを以下に示しておきます。

/*
  statistics.cc

  Written by Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <vector>
#include <list>
#include <cstdlib>      // random(), etc
#include <cassert>
#include "statistics.hpp"
#include "statistics_online.hpp"

int main()
{
  const double maximum_permissible_error(1e-12);
  const char *result_names[] = {
    "mean",
    "variance",
    "standard deviation",
    "biased variance",
    "biased standard deviation",
    "skewness",
    "kurtosis",
    "maximum value",
    "minimum value",
    "geometric mean",
    "harmonic mean",
    "root mean square",
  };
  std::vector<double> v(1000);
  statistics_online<double> statistics_online;
  srandomdev();
  for (size_t i=0; i<v.size(); ++i) {
    double x((double)random()/((1L<<31)-1));
    v[i] = x;
    statistics_online.add(x);
  }
  double results[] = {
    arithmetic_average(v.begin(), v.end()),
    unbiased_variance(v.begin(), v.end()),
    unbiased_standard_deviation(v.begin(), v.end()),
    sample_variance(v.begin(), v.end()),
    sample_standard_deviation(v.begin(), v.end()),
    skewness(v.begin(), v.end()),
    kurtosis(v.begin(), v.end()),
    *max_element(v.begin(), v.end()),
    *min_element(v.begin(), v.end()),
    geometric_average(v.begin(), v.end()),
    harmonic_average(v.begin(), v.end()),
    generalized_average(v.begin(), v.end(), 2),
  };
  double online_results[] = {
    statistics_online.get_mean(),
    statistics_online.get_unbiased_variance(),
    statistics_online.get_unbiased_standard_deviation(),
    statistics_online.get_sample_variance(),
    statistics_online.get_sample_standard_deviation(),
    statistics_online.get_skewness(),
    statistics_online.get_kurtosis(),
    statistics_online.get_maximum_value(),
    statistics_online.get_minimum_value(),
    statistics_online.get_geometric_mean(),
    statistics_online.get_harmonic_mean(),
    statistics_online.get_root_mean_square(),
  };
  for (size_t j=0; j<sizeof(online_results)/sizeof(online_results[0]); ++j) {
    double error(sqrt(pow(online_results[j] - results[j], 2)));
    std::cout << result_names[j] << ":\t" << online_results[j] << "\t" << results[j] << "\t" << error << std::endl;
    assert(error < maximum_permissible_error);
  }

  const char *moving_result_names[] = {
    "moving arithmetic average",
    "moving geometric average",
    "moving harmonic average",
    "moving root mean square",
    "moving sample variance",
    "moving sample standard deviation",
    "moving unbiased variance",
    "moving unbiased standard deviation",
    "moving skewness",
    "moving kurtosis",
    "moving minimum value",
    "moving maximum value",
  };
  std::vector<double> m(v.size());
  moving_statistics_online<double> moving_statistics_online(v.size());
  for (size_t i=0; i<v.size()*10; ++i) {
    double x((double)random()/((1L<<31)-1));
    m[i % v.size()] = x;
    moving_statistics_online.add(x);
    if (i < v.size()) continue;
    double moving_results[] = {
      arithmetic_average(m.begin(), m.end()),
      geometric_average(m.begin(), m.end()),
      harmonic_average(m.begin(), m.end()),
      generalized_average(m.begin(), m.end(), 2),
      sample_variance(m.begin(), m.end()),
      sample_standard_deviation(m.begin(), m.end()),
      unbiased_variance(m.begin(), m.end()),
      unbiased_standard_deviation(m.begin(), m.end()),
      skewness(m.begin(), m.end()),
      kurtosis(m.begin(), m.end()),
      *min_element(m.begin(), m.end()),
      *max_element(m.begin(), m.end()),
    };
    double online_moving_results[] = {
      moving_statistics_online.get_mean(),
      moving_statistics_online.get_geometric_mean(),
      moving_statistics_online.get_harmonic_mean(),
      moving_statistics_online.get_root_mean_square(),
      moving_statistics_online.get_sample_variance(),
      moving_statistics_online.get_sample_standard_deviation(),
      moving_statistics_online.get_unbiased_variance(),
      moving_statistics_online.get_unbiased_standard_deviation(),
      moving_statistics_online.get_skewness(),
      moving_statistics_online.get_kurtosis(),
      moving_statistics_online.get_minimum_value(),
      moving_statistics_online.get_maximum_value(),
    };
    for (size_t j=0; j<sizeof(online_moving_results)/sizeof(online_moving_results[0]); ++j) {
      double error(sqrt(pow(online_moving_results[j] - moving_results[j], 2)));
      if (i+1 == v.size()*10)
        std::cout << moving_result_names[j] << ":\t" << online_moving_results[j] << "\t" << moving_results[j] << "\t" << error << std::endl;
      assert(error < maximum_permissible_error);
    }
  }
  return 0;
}

これら非オンラインとオンラインとの自乗誤差は 1,000 データにおいて 10E-12 程度に収まっていることは確認できています。

加えて、averages_online-20160119.tar.bz2 にて quantiles.hpp, quantiles_online.hpp も実装しました。これはメディアン等を求めるテンプレートライブラリです。しかし、オンライン版は移動平均しか現時点では実装しておりませんので、中途ゆえ解説は省略させて頂きます。使い方はソースコードを見れば一目瞭然だと思います。そして、averages_online-20160129.tar.bz2 にてヒューリスティック手法ですがオンライン版を実装しました。元になるアルゴリズムは、例えばこの文献などから辿れます。

それが実装された quantiles_online.hpp の使用例として、その tarball にあるコードを以下に示しておきます。

/*
  quantiles.cc

  Written by Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <vector>
#include <list>
#include <functional>
#if __cplusplus < 201103L
#include <boost/function.hpp>
#endif
#include <cstdlib>      // random(), etc
#include <cassert>
#if __cplusplus < 201103L
#include <boost/utility.hpp>
#endif
#include "quantiles.hpp"
#include "quantiles_online.hpp"

int main()
{
  const double maximum_permissible_error(1e-12);
  const char *result_names[] = {
    "p=0.25 quantile",
    "p=0.5 quantile",
    "p=0.75 quantile",
    "0-quantile (least)",
    "2-quantile (median)",
    "1-quantile (greatest)",
  };
  std::vector<double> v(1000);
  double probabilities[] = { .25, .5, .75, };
  heuristic_quantiles_online<double> quantiles_online(probabilities);
  srandomdev();
  for (size_t i=0; i<v.size(); ++i) {
    double x((double)random()/((1L<<31)-1));
    v[i] = x;
    quantiles_online.add(x);
  }
  sort(v.begin(), v.end());
  double results[] = {
    quantile_in_sorted_container(v.begin(), v.end(), 0.25),
    quantile_in_sorted_container(v.begin(), v.end(), 0.5),
    quantile_in_sorted_container(v.begin(), v.end(), 0.75),
    least_in_sorted_container(v.begin(), v.end()),
    median_in_sorted_container(v.begin(), v.end()),
    greatest_in_sorted_container(v.begin(), v.end()),
  };
  double online_results[] = {
    quantiles_online.get_quantile(0.25),
    quantiles_online.get_quantile(0.5),
    quantiles_online.get_quantile(0.75),
    quantiles_online.get_minimum_value(),
    quantiles_online.get_median(),
    quantiles_online.get_maximum_value(),
  };
  for (size_t j=0; j<sizeof(online_results)/sizeof(online_results[0]); ++j) {
    double error(sqrt(pow(online_results[j] - results[j], 2)));
    std::cout << result_names[j] << ":\t" << online_results[j] << "\t" << results[j] << "\t" << error << std::endl;
    if (!(j<=2 || j==4)) assert(error < maximum_permissible_error);
  }

  std::vector<double> m(v.size()), s(m.size());
  moving_quantiles_online<double> moving_quantiles_online(v.size());
  for (size_t i=0; i<v.size()*10; ++i) {
    double x((double)random()/((1L<<31)-1));
    m[i % v.size()] = x;
    moving_quantiles_online.add(x);
    if (i < v.size()) continue;
    copy(m.begin(), m.end(), s.begin());
    sort(s.begin(), s.end());
    double moving_results[] = {
      quantile_in_sorted_container(s.begin(), s.end(), 0.25),
      quantile_in_sorted_container(s.begin(), s.end(), 0.5),
      quantile_in_sorted_container(s.begin(), s.end(), 0.75),
      least_in_sorted_container(s.begin(), s.end()),
      median_in_sorted_container(s.begin(), s.end()),
      greatest_in_sorted_container(s.begin(), s.end()),
    };
    double online_moving_results[] = {
      moving_quantiles_online.get_quantile(0.25),
      moving_quantiles_online.get_quantile(0.5),
      moving_quantiles_online.get_quantile(0.75),
      moving_quantiles_online.get_minimum_value(),
      moving_quantiles_online.get_median(),
      moving_quantiles_online.get_maximum_value(),
    };
    for (size_t j=0; j<sizeof(online_moving_results)/sizeof(online_moving_results[0]); ++j) {
      double error(sqrt(pow(online_moving_results[j] - moving_results[j], 2)));
      if (i+1 == v.size()*10)
        std::cout << "moving " << result_names[j] << ":\t" << online_moving_results[j] << "\t" << moving_results[j] << "\t" << error << std::endl;
      assert(error < maximum_permissible_error);
    }
  }
  return 0;
}

さらに、averages_online-20160210.tar.bz2 にて、histogram.hpp, histogram_online.hpp を実装しました。

さらに、averages_online-20160302.tar.bz2 にて、multidimensional_histogram.hpp, multidimensional_histogram_online.hpp を実装しました。

さらに、averages_online-20160309.tar.bz2histogram.hpp, multidimensional_histogram.hpp にて、ヒストグラムからその確率密度に従う乱数を生成する histogram(_probability)_cumulative_distribution(_pseudo) を実装しました。

使い方は、ソースコードをご覧下さい。


初等整数コマンドラインツール群

by taiji at 2015/07/22 09:31

シェルスクリプトやコマンドラインである数が素数か、以下のように調べたい。

$ gmpxx-isprime 1907 && echo prime || echo composite
prime

そうした場合、GNU MP を利用すればそういったコマンドラインプログラムを簡単に作る事ができます。

/*
  isprime.cc

  Copyright (C) 2015 Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <string>
#include <libgen.h>
#include <gmpxx.h>

int usage(std::ostream &os, int argc, char *argv[])
{
  (void)argc;
  os << "Usage:\n\t" << basename(argv[0]) << " [number ...]\n";
  return os == std::cerr ? 1 : 0;
}
int main(int argc, char *argv[])
{
  int a;
  mpz_class n;

  if (argc < 2) return usage(std::cerr, argc, argv);
  for (a=1; a<argc; ++a) {
    n = argv[a];
    if (!mpz_probab_prime_p(n.get_mpz_t(), 25))
      return 1;
  }
  return 0;
}

ここでは GNU MP の C++ API である gmpxx.h を利用しています。このように作成しておけば、数列がすべて素数であるかを判定するプログラムになります。

$ gmpxx-isprime 953 1907 && echo primes || echo composites
primes

他にも以下のようなさまざまなプログラムを用意しました。接頭辞「gmpxx-」がインストールされた時のコマンド名になります。

isprime                 素数判定
iscomposite             合成数判定
primes                  素数を範囲で出力
primes_reverse          素数を範囲で逆順に出力
n_primes                素数を個数で出力
n_primes_reverse        素数を個数で逆順に出力
composites              合成数を範囲で出力
factorize               素因数分解
composite               合成数生成
gcd                     最大公約数
lcm                     最小公倍数
divisors                約数
n_divisors              約数の個数
sum_divisors            約数関数

next_prime              次の素数
previous_prime          前の素数

next_composite          次の合成数
previous_composite      前の合成数

integer_valued_polynomial_progression   整数値多項式の数列生成
prime_polynomial_progression            整数値多項式の素数生成

highlycomposites        高度合成数の数列生成
highlyabundants         高度過剰数の数列生成
superabundants          超過剰数の数列生成

naturals                さまざまな自然数や整数の数列生成

totient_function        トーシェント関数の生成
cototient_function      コトーシェント関数の生成
mebius_function         メビウス関数の生成

primes_util-collate     C/C++ 用素数表作成及び検証

例えば、以下のように高度合成数を出力できます。

$ gmpxx-highlycomposites -u 1000000
1
2
4
6
12
24
36
48
60
120
180
240
360
720
840
1260
1680
2520
5040
7560
10080
15120
20160
25200
27720
45360
50400
55440
83160
110880
166320
221760
277200
332640
498960
554400
665280
720720

そこで、これを素因数分解するには以下のようにします。

$ gmpxx-factorize $(gmpxx-highlycomposites)
1       1^0
2       2^1
4       2^2
6       2^1     3^1
12      2^2     3^1
24      2^3     3^1
36      2^2     3^2
48      2^4     3^1
60      2^2     3^1     5^1
120     2^3     3^1     5^1
180     2^2     3^2     5^1
240     2^4     3^1     5^1
360     2^3     3^2     5^1
720     2^4     3^2     5^1
840     2^3     3^1     5^1     7^1
1260    2^2     3^2     5^1     7^1
1680    2^4     3^1     5^1     7^1
2520    2^3     3^2     5^1     7^1
5040    2^4     3^2     5^1     7^1
7560    2^3     3^3     5^1     7^1
10080   2^5     3^2     5^1     7^1
15120   2^4     3^3     5^1     7^1
20160   2^6     3^2     5^1     7^1
25200   2^4     3^2     5^2     7^1
27720   2^3     3^2     5^1     7^1     11^1
45360   2^4     3^4     5^1     7^1
50400   2^5     3^2     5^2     7^1
55440   2^4     3^2     5^1     7^1     11^1
83160   2^3     3^3     5^1     7^1     11^1
110880  2^5     3^2     5^1     7^1     11^1
166320  2^4     3^3     5^1     7^1     11^1
221760  2^6     3^2     5^1     7^1     11^1
277200  2^4     3^2     5^2     7^1     11^1
332640  2^5     3^3     5^1     7^1     11^1
498960  2^4     3^4     5^1     7^1     11^1
554400  2^5     3^2     5^2     7^1     11^1
665280  2^6     3^3     5^1     7^1     11^1
720720  2^4     3^2     5^1     7^1     11^1    13^1

また、コマンド「composite」は以下の形式の合成数を生成します。

\[c_1\times c_2\cdots\times v_1^{i_1}\times v_2^{i_2}\cdots\times w_1^{j_1}\times w_2^{j_2}\cdots\]

但しここで、\(0<i_i<m, 0\le j_i\le m\)。

$ gmpxx-composite -h
Usage:
        gmpxx-composite [-c number ...] [-v number ...] [-w number ...] [-m times]
$ gmpxx-factorize $(gmpxx-composite -c 2 3 -v 5 7 -w 11 13 -m 2)
210     2^1     3^1     5^1     7^1
1050    2^1     3^1     5^2     7^1
1470    2^1     3^1     5^1     7^2
2310    2^1     3^1     5^1     7^1     11^1
2730    2^1     3^1     5^1     7^1     13^1
7350    2^1     3^1     5^2     7^2
11550   2^1     3^1     5^2     7^1     11^1
13650   2^1     3^1     5^2     7^1     13^1
16170   2^1     3^1     5^1     7^2     11^1
19110   2^1     3^1     5^1     7^2     13^1
30030   2^1     3^1     5^1     7^1     11^1    13^1
80850   2^1     3^1     5^2     7^2     11^1
95550   2^1     3^1     5^2     7^2     13^1
150150  2^1     3^1     5^2     7^1     11^1    13^1
210210  2^1     3^1     5^1     7^2     11^1    13^1
1051050 2^1     3^1     5^2     7^2     11^1    13^1

また、コマンド「totient_function」はオイラーのトーシェント関数\(\phi(n)\):

\[\phi(n) = n\prod_{k=1}^{d}\left(1-\frac{1}{p_k}\right)\]

但しここで、\(p_k\)は\(n\)の互いに素な\(d\)個の素因数。ここで分数があるので有理数の扱いが必要に思われますが、少し工夫すれば素因数分解の過程で割り切れる除算のみで実装できます。

mpz_class totient_function(mpz_class n)
{
  mpz_class f(2), r(n);

  for (; f <= n; mpz_nextprime(f.get_mpz_t(), f.get_mpz_t())) {
    if (n % f == 0) {
      r -= r / f;
      n /= f;
    }
    while (n % f == 0)
      n /= f;
  }
  return r;
}

さて、コマンド「naturals」は相当な説明を加えないとなりません。これは「-type name」オプションで以下のようなさまざまな数列を生成します。

natural                         自然数
odd                             奇数
even                            偶数
singly_even                     単偶数
doublly_even                    複偶数
composite                       合成数
prime                           素数
sophie_germain_prime            ソフィー・ジェルマン素数
interprime                      インター素数
twin_prime                      双子素数
cousin_primes                   いとこ素数
sexy_primes                     セクシー素数
square_free                     非平方数
cube_free                       非立方数
perfect_square                  平方数
perfect_power                   累乗数
powerful                        多冪数

highlycomposite                 高度合成数
highlyabundant                  高度過剰数
superabundant                   超過剰数
abundant                        過剰数
primitive_abundant              原始過剰数
primitive_abundant_have_no_abundant_proper_divisor      自身を除く過剰な約数をもたない原始過剰数*
non_prime_unusual               非素数異常数*

pseudoperfect, semiperfect      擬似完全数
primitive_pseudoperfect, primitive_semiperfect          原始擬似完全数
perfect                         完全数
multiperfect                    倍積完全数
deficient                       不足数
superperfect                    超完全数
hemiperfect                     ヘミ完全数
refactorable                    リファクタブル数
unitary_perfect                 ユニタリ完全数
primary_pseudoperfect           プライマリ擬似完全数
weird                           不思議数
practical                       実際数*
ruth_aaron_pair                 ルース=アーロン・ペア
primorial_ruth_aaron_pair       素数階乗ルース=アーロン・ペア
unusual                         異常数

arithmetic                      算術数*
harmonic_divisor                調和数
round_number                    ラウンド数
lucas_carmichael                ルーカス・カーマイケル数
semiprime                       擬似素数
sphenic                         楔数
amicable                        友愛数
betrothed                       婚約数、準友愛数
augumented_amicable             拡張友愛数
friendly                        フレンドリ数

mersenne                        メルセンヌ数
mersenne_prime                  メルセンヌ素数
wilson_quotient                 ウィルソン商
wilson_prime                    ウィルソン素数

fortunate                       フォーチューン数
hilbert                         ヒルベルト数
pronic                          矩形数
power_of_two                    2の冪数
cullen                          カレン数
cullen_prime                    カレン素数
woodall                         ウッダル数
woodall_prime                   ウッダル素数
regular                         正規数*
factorial_prime                 階乗素数
euclid                          ユークリッド数
euclid_prime                    ユークリッド素数
kummer                          クンマー数
kummer_prime                    クンマー素数
primorial_prime                 素数階乗素数
triangular                      三角数
triangular_pyramidal, tetrahedral       三角錐数、四面体数
pentatope                       五胞体数
pentagonal                      五角数
square_triangular               平方三角数
squared_triangular              平方化三角数*
lazy_caterer                    怠惰な仕出し屋数*
jacobsthal                      ジェイコブザル数
jacobsthal_lucas                ジェイコブザル・ルカ数
kynea                           カイネア数
prime_kynea                     カイネア素数
thabit                          サビット数
prime_thabit                    サビット素数
fermat                          フェルマー数
catalan                         カタラン数
derangement                     完全順列
nontotient                      非トーシェント数
noncototient                    非コトーシェント数
perfect_totient                 完全トーシェント数
highly_totient                  高度トーシェント数
highly_cototient                高度コトーシェント数
odd_semifactrial                奇数の半階乗数
even_semifactrial               偶数の半階乗数

fibonacci                       フィボナッチ数
leonardo                        レオナルド数
lucas                           リュカ数
pell                            ペル数
perrin                          ペラン数
padovan                         パドバン数
motzkin                         モズキン数
schreder_hipparchus             シュレーダー・ヒパルカス数
sylvester                       シルベスター数
telephone                       電話数
alternating_factorial           交互階乗数
hyperfactorial                  超階乗数
exponential_factorial           指数階乗数

ここでアスタリスク「*」付きの日本語名は筆者の命名です。いくつか実装に工夫が必要なものを紹介しましょう。

調和数(harmonic divisor numbers)は自然数\(n\)のうち\(k\)個すべての約数\(d_i\)の調和平均:

\[\frac{k}{\sum_{i=1}^{k}\frac{1}{d_i}}\]

が整数値をもつ\(n\)のことです。ここで分数があるので有理数の扱いが必要に思われますが、最小公倍数を使って工夫すれば最終的に割り切れるか否かで実装できます。

bool is_harmonic_divisor(const mpz_class n)
{
  std::vector<mpz_class> q = divisors(n);
  mpz_class d(0), c(q[0]);
  size_t i;

  for (i=1; i<q.size(); ++i)
    mpz_lcm(c.get_mpz_t(), c.get_mpz_t(), q[i].get_mpz_t());
  for (i=0; i<q.size(); ++i)
    d += c/q[i];
  return (mpz_class(q.size())*c) % d == 0;
}

擬似完全数(pseudoperfect numbers)は自然数\(n\)のうち自身を除くいくつかの約数\(d_i\)の和が元の数に等しい数\(n\)のことです。この「いくつかの」というのが厄介に思われますが、ビット演算で約数にフラグを割り当てすべての組み合わせを検証すれば実装できます。

bool is_pseudoperfect(const mpz_class n)
{
  std::vector<mpz_class> d = proper_divisors(n);
  mpz_class i, c;
  size_t j;

  mpz_ui_pow_ui(c.get_mpz_t(), 2, d.size());
  for (i=1; i<c; ++i) {
    mpz_class s(0);
    for (j=0; j<d.size(); ++j)
      if (mpz_tstbit(i.get_mpz_t(), j))
        s += d[j];
    if (s == n)
      return true;
  }
  return false;
}

さて、ここまで GNU MP を利用した実装を紹介して来ましたが、10桁程度の数値でも速度が遅く、その程度の桁数の整数なら C/C++ の整数で十分ですし、実行速度が全く違います。そこで C/C++ 用の初等整数コマンドラインツール群を用意しました。

但し、素数については実行速度を実用的なものにするため「/var/tmp/isprime_bit_table.bin」を作成します。作成途中で中断すると不完全な素数表となり、正しく動作しません。よって、初回のみ以下のように作成及び検証をじっくりと行うことをお勧めします。

$ gmpxx-prime_util-collate -u 2000000000 && echo valid || echo invalid

これで「/var/tmp/isprime_bit_table.bin」は 1.2GB ほど、Core i7 で、素数表の作成は5分ほどで終わるのですが、GNU MP による検証で1時間ほど掛かります。

さて、GNU MP には mpz_perfect_square_pmpz_perfect_power_p が用意されているので、上記の「平方数(perfect square numbers)」「累乗数(perfect power numbers)」はそれを利用すればよかったのです。しかし、C/C++ のみで実装するとなると、これに変わるものを自前で用意しなければなりません。

こちらの場合、素数表があるので素因数分解が使えますし、以下の最大公約数 boost::math::gcd を利用した実装で事足ります。

template <typename T>
bool is_perfect_square(T n)
{
  T f(2), c, c0(0);
  isprime_bit_table_class isprime_bit_table;

  for (; f <= n; f = isprime_bit_table.nextprime(f)) {
    c = 0;
    while (n % f == 0) {
      ++c;
      n /= f;
    }
    if (c % 2 == 1)
      return false;
    if (c0 != 0) {
      if (boost::math::gcd(c0, c) % 2 == 1)
        return false;
      c0 = boost::math::gcd(c0, c);
    }
    else
      c0 = c;
    if (!(f < isprime_bit_table.get_upper())) throw std::length_error("isprime_bit_table: over the length");
  }
  return true;
}
template <typename T>
bool is_perfect_power(T n)
{
  T f(2), c, c0(0);
  isprime_bit_table_class isprime_bit_table;

  for (; f <= n; f = isprime_bit_table.nextprime(f)) {
    c = 0;
    while (n % f == 0) {
      ++c;
      n /= f;
    }
    if (c == 1)
      return false;
    if (c0 != 0) {
      if (boost::math::gcd(c0, c) == 1)
        return false;
      c0 = boost::math::gcd(c0, c);
    }
    else
      c0 = c;
    if (!(f < isprime_bit_table.get_upper())) throw std::length_error("isprime_bit_table: over the length");
  }
  return true;
}

配布物

参考文献

  1. https://en.wikipedia.org/wiki/Category:Integer_sequences
  2. https://ja.wikipedia.org/wiki/Category:整数の類

ics-tools - Emacs calendar diary 形式と iCalendar 形式のファイルのためのツール群

by taiji at 2013/10/31 23:48

以前、iCalendar 形式のファイルのためのツール群を著したが、加えて、Emacs calendar diary 形式から iCalendar 形式へ変換またはその逆変換する以下のツール群を追加した。

  • ecd2ics
  • ics2ecd

ecd2ics - Emacs calendar diary 形式から iCalendar 形式のファイルを生成する。

例えば、同梱されているデバッグ用の examples/emacs-calendar-diary/diary ファイルから diary.ics ファイルを生成するには、以下のように実行する。

$ iconv -f euc-jp -t utf-8 examples/emacs-calendar-diary/diary | ecd2ics --ical - > diary.ics

ics2ecd - iCalendar 形式のファイルから Emacs calendar diary 形式を生成する。

例えば、同梱されているデバッグ用の examples/emacs-calendar-diary/diary から変換した先の diary.ics ファイルから diary を生成するには、以下のように実行する。

$ ics2ecd diary.ics | iconv -f utf-8 -t euc-jp > diary

成果物は、iCalendar 形式のファイルのためのツール群にて配布。