『Binary Hacks Rebooted―低レイヤの世界を探検するテクニック89選』に込めた思い

2024-08-28 19:01:51 JST · 元記事

はじめに

2024年8月28日に『Binary Hacks Rebooted―低レイヤの世界を探検するテクニック89選』という本の販売が開始されました。 著者の一人として、この本の紹介をします。 お持ちでない方は、Amazon や、楽天、お近くの書店などで購入が可能です。 購入していただけると大変嬉しいです。 どうぞよろしくおねがいします。

どのような本か

『Binary Hacks Rebooted』は2006年に発売された『Binary Hacks』の再構築版です。 再構築版と聞くと『Binary Hacks』とどれくらい異なるのが気になると思いますが、ほぼすべてが異なります。 詳細な目次はこちら で確認できますが、『Binary Hacks』からそのまま引き継いだ Hack は一つもなく、全てが新しい Hack です。

『Binary Hacks Rebooted』は低レイヤプログラミングの技術を幅広く、かつ深く解説した書籍です。 「低レイヤプログラミング」という単語は広い意味を持ちますが、この書籍ではELF フォーマットの解析からファームウェアの作成、浮動小数点例外まで幅広い内容を扱っています。 そして、各トピックについて日本語はもちろん、英語でも類似する解説が少ないと思われる内容を紹介しています。 カーネル／VM探検隊 に参加したことがある方は、「読むカーネル／VM探検隊」のような感覚で読んでいただけるかと思います。

この本に込めた思い

私達はこの本を2つのことを目標に書いていました。

1つ目は「低レイヤプログラミングを面白いと感じてもらう」ことです。 細かいことはわからなくても構いません。 「あーこんなことできるんだ、おもしろ」と思ってもらえれば大変嬉しいです。 あなたが好きなように「探検」してみてください。

低レイヤプログラミングは一見とっつきにくい分野ですが、一度どのように動いているのかを理解すると、様々な面白いことができるようになります。 例えば、自分だけのファームウェアを作ってChromebook を動かすこともできますし、DWARF を使って四則演算を行うこともできます。 いつかどこかで私と会った時に、Hack の内容で盛り上がれることを楽しみにしています。

2つ目は「低レイヤプログラミングで困っている方の助けになる」ことです。 低レイヤプログラミングを実際に行ってみると、不可解なバグに多数直面します。 そのうちの多くは把握していない仕様、挙動によるものが多く、このようなバグを解決するためには幅広い分野の知識が必要です。 ELF ハックの多くは私が過去に実際に知りたかった知識をまとめたものが多く、同様のバグで困っている方の助けになればと思っています。

Thread Local Storage (TLS) の解説はその一例です。 TLS は非常に多くのプログラムで人知らず使われている機能ですが、低レイヤプログラミングを実際に行っているとその仕組みを理解するように迫られることがあります。 ところが、TLS についての詳細な解説は日本語でも英語でもあまり見かけません。 「#10 TLSのしくみを理解する 」は誰よりも3年前の過去の自分に当てて書いた Hack です。 この内容を3年前に知っていたら何ヶ月もデバッグに時間を費やすことはなかったでしょう。

私と同様に意味不明な挙動に悩まされている方が、世の中にはきっとたくさんいるはずです。 segmentation fault (core dumped) を見飽きたときに、この本を開いてみてください。

どんなふうに読んでほしいか

この本は最初のページから順番に読んでいく必要はありません。 いくつかの続きものの Hack を除いて各 Hack は独立して読めるようになっています。　 目次を見て面白そうなところを読んでみたり、サイコロを振って出た目の Hack を読んでみたりしてください。 どんな方でもパラパラとめくっているうちに意外な気づきがあるはずです。

サンプルコードを実際に動かして見てから文書をしっかり読んでみるのも良いかもしれません。 サンプルコードはhttps://github.com/oreilly-japan/binary-hacks-rebooted/にあります。

おわりに

この本を書くために膨大な時間を費やし、各種の仕様書、Linux や glibc のソースコードを何度も読み返しました。 自信を持って面白い本だと言えます。 この本を読んで、一人でも「低レイヤプログラミングを面白い」と感じてもらえる方がいれば幸いです。

最後に、もし、あなたが2042年ごろにこの文書を読んでおり、『Binary Hacks』の再再構築版が出したいと思っているなら、ぜひ私に連絡してください。 喜んで協力したいです。

2024年8月28日

河田 旺

gcc ビルトイン関数の呼び出し

2022-01-12 21:50:53 JST · 元記事

ダイナミックリンカのソースコード(elf/rtld.c)を見ていると興味深いコメントがあった。

  /* Partly clean the `bootstrap_map' structure up.  Don't use
     `memset' since it might not be built in or inlined and we cannot
     make function calls at this point.  Use '__builtin_memset' if we
     know it is available.  We do not have to clear the memory if we
     do not have to use the temporary bootstrap_map.  Global variables
     are initialized to zero by default.  */

「memsetはビルドインでないもしくはインライン化されていない可能性があるため、ここでは呼び出してはいけない」。このコメントはダイナミックリンカの起動直後にあるため、再配置が終わっていない関数を呼び出すとアドレスが埋まっていないGOTを参照してしまい、Segmentation Faultが発生することを危惧したものだと考えられる。

ということはgccのビルトイン関数は必ずインライン化されなければならない。実際gccのビルトイン関数の説明には以下のような記述がある。

With the exception of built-ins that have library equivalents such as the standard C library functions discussed below, or that expand to library calls, GCC built-in functions are always expanded inline and thus do not have corresponding entry points and their address cannot be obtained. Attempting to use them in an expression other than a function call results in a compile-time error.

ビルトイン関数を使う小さなプログラムでも実験的に確かめることができる。

> cat main.c
#include <stdio.h>

int main() {
    int c = __builtin_clz(0x10101);
    printf("c = %d\n", c);
    return 0;
}
> gcc main.c -static -o main
> objdump -d main | grep "<main>" -A 16
0000000000401ce5 <main>:
  401ce5:       f3 0f 1e fa             endbr64 
  401ce9:       55                      push   %rbp
  401cea:       48 89 e5                mov    %rsp,%rbp
  401ced:       48 83 ec 10             sub    $0x10,%rsp
  401cf1:       c7 45 fc 0f 00 00 00    movl   $0xf,-0x4(%rbp)
  401cf8:       8b 45 fc                mov    -0x4(%rbp),%eax
  401cfb:       89 c6                   mov    %eax,%esi
  401cfd:       48 8d 3d 00 33 09 00    lea    0x93300(%rip),%rdi        # 495004 <_IO_stdin_used+0x4>
  401d04:       b8 00 00 00 00          mov    $0x0,%eax
  401d09:       e8 92 ee 00 00          callq  410ba0 <_IO_printf>
  401d0e:       b8 00 00 00 00          mov    $0x0,%eax
  401d13:       c9                      leaveq 
  401d14:       c3                      retq   
  401d15:       66 2e 0f 1f 84 00 00    nopw   %cs:0x0(%rax,%rax,1)
  401d1c:       00 00 00 
  401d1f:       90                      nop

__builtin_clzはインライン化されているようだ。

LD_AUDITとGlobal Offset Table

2022-01-08 22:05:26 JST · 元記事

GNU製のLinux向け動的リンカld-linux.so には環境変数 LD_AUDIT 経由で使える監視APIがあります。このAPIを使えば、ld-linux.soの様々な挙動にフックして、その挙動を監視・干渉することができます。また、興味深いことにLD_AUDITを有効にするとGOT(Global Offset Table)周りの挙動も変化します。

PLT(Procedure Linkage Table)やGOTについて知らない場合はこちらの記事が分かりやすいです。この記事のglibcはすべてコミット 137ed5ac44を参照しています。

動作例

LD_AUDITを使ってltrace(1)まがいのものを作ってみます。以下の3つのファイルを作ります。なお、このサンプルはLinux glibcで、LD_AUDIT機能による関数トレースを参考にしました。

audit.c

#define _GNU_SOURCE
#include <link.h>
#include <stdio.h>

unsigned int la_version(unsigned int version) { return LAV_CURRENT; }

unsigned int la_objopen(struct link_map* map, Lmid_t lmid, uintptr_t* cookie) {
    printf("[LD_AUDIT] %s loaded\n", map->l_name);
    return LA_FLG_BINDTO | LA_FLG_BINDFROM;
}

ElfW(Addr) la_x86_64_gnu_pltenter(ElfW(Sym) * sym, unsigned int ndx,
                                  uintptr_t* refcook, uintptr_t* defcook,
                                  La_x86_64_regs* regs, unsigned int* flags,
                                  const char* symname, long* framesizep) {
    printf("[LD_AUDIT] %s entering\n", symname);
    return sym->st_value;
}

hello.c

#include <stdio.h>

int main() {
    puts("Hello World!");
    puts("Hello World!");
    puts("Hello World!");
    return 0;
}

run.sh

#! /bin/bash -eux

gcc -o libaudit.so -shared -fpic -fPIC -Wl,-soname,libaudit.so audit.c
gcc -o hello hello.c

run.sh を実行すると以下のような結果が得られて、ltraceっぽいことができているのがわかります。

> ./run.sh
+ gcc -o libaudit.so -shared -fpic -fPIC -Wl,-soname,libaudit.so audit.c
+ gcc -o hello hello.c
+ LD_AUDIT=libaudit.so
+ ./hello
[LD_AUDIT]  loaded
[LD_AUDIT] /lib64/ld-linux-x86-64.so.2 loaded
[LD_AUDIT] /lib/x86_64-linux-gnu/libc.so.6 loaded
[LD_AUDIT] _dl_find_dso_for_object entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] __tunable_get_val entering
[LD_AUDIT] puts entering
Hello World!
[LD_AUDIT] puts entering
Hello World!
[LD_AUDIT] puts entering
Hello World!
+ LD_BIND_NOW=1
+ LD_AUDIT=libaudit.so

解説

la_versionはLD_AUDITを使うために必ず必要な関数、la_objopenは新しいshared objectがロードされたときに呼び出される関数であり、今回は深入りしません。

la_x86_64_gnu_pltenterは PLT(Procedure Linkage Table)がGOT(Global Offset Table)を埋めるときに呼び出される関数であり、具体的にはここで呼び出されます。

この実行結果で興味深いのは [LD_AUDIT] puts entering が複数回出力されていることです。一度PLTのエントリがGOTに正しいアドレスを埋めた後は、putsの呼び出しはlibcのアドレスを埋める部分を通らず、libc中のputs本体に直接jmpしla_x86_64_gnu_pltenterを呼び出すとことはできないはずです。

ここから考えるとLD_AUDITが設定されているときはGOTにアドレスを埋める処理が行われない、と推測できます。glibcの実装を確認してみましょう。

• LD_AUDIT によって設定されたla_x86_64_gnu_pltenterはdl-runtime.cのL435-L439で呼び出されます。
• これは dl-runtime.cの_dl_profile_fixupの中から呼ばれています。
• また、_dl_profile_fixupにはGOTの値を更新する処理がありません。
• 一方、dl-runtime.cには _dl_fixupもあり、LD_AUDITの有無で _dl_profile_fixupと _dl_fixupを切り替えているのだと推察できます。
• _dl_fixupは dl-trampoline.hの_dl_runtime_resolveから、_dl_profile_fixupは dl-trampoline.hの_dl_runtime_profile からそれぞれ呼び出されています。
• これら２つの関数のどちらかが dl-machine.hでgot[2]に設定されます。
• この分岐を決定するのはelf_machine_runtime_setupのprofile引数で dynamic-link.hのELF_DYNAMIC_RELOCATEが呼び出し元です。
• このマクロはdl-reloc.cで呼び出されています。
• dl-reloc.cの_dl_relocate_objectの最初の方でconsider_profilingがLD_AUDITの有無に従って設定されています。よって、LD_AUDITが設定されているときは _dL_profile_fixupが呼び出され、GOTが更新されないことがわかります。

また、dl-machine.hにはprofiling extensionが使用されている場合、GOTを書き換えない旨のコメントがあります。

/ The got[2] entry contains the address of a function which gets called to get the address of a so far unresolved function and jump to it. The profiling extension of the dynamic linker allows to intercept the calls to collect information. In this case we don't store the address in the GOT so that all future calls also end in this function. /

参考

• PLTとGOTってなんだっけ
• Linux glibcで、LD_AUDIT機能による関数トレース
• glibc 137ed5ac44

GNU_IFUNCとは何か

2021-11-07 10:05:40 JST · 元記事

GNU_IFUNC(GNU indirect function)とは同一の関数の複数の実装から、ロード時に最適な実装を選択する仕組みです¹。通常、同一の関数の複数の実装が存在した場合、動的リンカ(ld-linux.so)が最初に見つけた実装が選択されますが、その選択するための基準を開発者が決定できるところに、この機能の意義があります。

GNU_IFUNCを使った簡単な例を以下に示します。resolve_foo で foo の実装が選択されていることがわかります。

> cat foo.c 
#include <stdio.h>

extern void foo();
void foo_default() { printf("foo_default\n"); }
void foo_1() { printf("foo_1\n"); }
void foo_2() { printf("foo_2\n"); }

void foo() __attribute__((ifunc("resolve_foo")));

static void *resolve_foo(void) {
    if (0)
        return foo_1;
    else if (42 == 41 + 1)
        return foo_2;
    else
        return foo_default;
}
> cat main.c 
void foo();

int main() {
    foo();
    return 0;
}
> gcc -shared -fpic -fPIC foo.c -o libfoo.so
> gcc -o main main.c libfoo.so
> ./main
foo_2

glibcでは memmove, memset, memcpy, strcmp, strstr などの関数でCPUの機能に応じて高度に最適化された実装をロード時に選択するために使われています²。例えば x86-64のmemcpy の場合はAVX512等のSIMD拡張が使えるかどうかを判定し、最適な実装を選択するコードが入っています。³

/*  glibc/sysdeps/x86_64/multiarch/ifunc-impl-list.c  */
 IFUNC_IMPL (i, name, __memcpy_chk,
          IFUNC_IMPL_ADD (array, i, __memcpy_chk,
                  CPU_FEATURE_USABLE (AVX512F),
                  __memcpy_chk_avx512_no_vzeroupper)
          IFUNC_IMPL_ADD (array, i, __memcpy_chk,
                  CPU_FEATURE_USABLE (AVX512VL),
                  __memcpy_chk_avx512_unaligned)
          IFUNC_IMPL_ADD (array, i, __memcpy_chk,
                  CPU_FEATURE_USABLE (AVX512VL),
                  __memcpy_chk_avx512_unaligned_erms)
          IFUNC_IMPL_ADD (array, i, __memcpy_chk,
                  CPU_FEATURE_USABLE (AVX),
                  __memcpy_chk_avx_unaligned)

最後にGNU_IFUNCの実装選択関数を悪用してみます。glibcのGNU_IFUNCに関する再配置情報を処理する部分を見ると、ロード時に実装を選択するための関数が呼ばれていることがわかります⁴。

/* glibc/sysdeps/x86_64/dl-machine.h */
else if (__glibc_unlikely (r_type == R_X86_64_IRELATIVE))
    {
      ElfW(Addr) value = map->l_addr + reloc->r_addend;
      if (__glibc_likely (!skip_ifunc))
    value = ((ElfW(Addr) (*) (void)) value) ();
      *reloc_addr = value;
    }

つまり、__attribute__(constructor)⁵ と同様に main() の前に関数を呼び出す手段として使えそうです。やってみましょう。

> cat foo.c 
#include <stdio.h>

extern void foo();
void foo_default() { printf("foo_default\n"); }
void foo_1() { printf("foo_1\n"); }
void foo_2() { printf("foo_2\n"); }

void foo() __attribute__((ifunc("resolve_foo")));

static void *resolve_foo(void) {
    printf("resolve_foo\n");
    if (0)
        return foo_1;
    else if (42 == 41 + 1)
        return foo_2;
    else
        return foo_default;
}
> cat main.c 
#include <stdio.h>

void foo();

int main() {
    printf("beginning of main()\n");
    foo();
    return 0;
}
> gcc -shared -fpic -fPIC foo.c -o libfoo.so
> gcc -o main main.c libfoo.so
> LD_BIND_NOW=1
> ./main
resolve_foo
beginning of main()
foo_2

R_X86_64_IRELATIVE が main() 起動後に遅延して解決されるのを防ぐために LD_BIND_NOW=1 が必要でしたが、無事 main() の前に resolve_foo を呼び出せたことがわかります。

setjmp, longjmpで使われるjmp_bufには何がどのように保存されているのか

2020-12-31 18:43:39 JST · 元記事

setjmp, longjmpはC言語で実行コンテキストを保存し、保存したコンテキストに復帰するために使われる関数(orマクロ)である。g++の例外実装の一つはこのsetjmp, longjmpを用いて実現されている。

この記事では 実行コンテキスト とは何を指すのか? また、その保存方法の実装について調べる。

jmp_buf

setjmpの実装はglibcの setjmp/setjmp.h にある。

typedef struct __jmp_buf_tag jmp_buf[1];

/* Store the calling environment in ENV, also saving the signal mask.
   Return 0.  */
extern int setjmp (jmp_buf __env) __THROWNL;

このjmp_bufはsetjmp/bits/types/struct___jmp_buf_tag.h で定義されている。

/* Calling environment, plus possibly a saved signal mask.  */
struct __jmp_buf_tag
  {
    /* NOTE: The machine-dependent definitions of `__sigsetjmp'
       assume that a `jmp_buf' begins with a `__jmp_buf' and that
       `__mask_was_saved' follows it.  Do not move these members
       or add others before it.  */
    __jmp_buf __jmpbuf;     /* Calling environment.  */
    int __mask_was_saved;  /* Saved the signal mask?  */
    __sigset_t __saved_mask;    /* Saved signal mask.  */
  };

x86-64の場合、__jmp_bufはsysdeps/x86/bits/setjmp.hで定義されている。__jmp_bufは8個のint型整数からなり、更に64bitかどうかで場合分けがなされていることがわかる。

# if __WORDSIZE == 64
typedef long int __jmp_buf[8];
# elif defined  __x86_64__
__extension__ typedef long long int __jmp_buf[8];
# else
typedef int __jmp_buf[6];
# endif

setjmp

setjmp自体の実装はsysdeps/x86_64/setjmp.Sにある。Linux x86-64では第一引数を%rdiで渡すので、第一引数として渡された構造体にレジスタの値を保存していることがわかる。保存されているのは、%rbx, %r12, %r13, %r14, %r15, %rsp, PCである。しかし、%rsp, %rbp, PCはPTR_MANGLEで処理されてから保存されている。PTR_MANGLEとは何だろうか? また、なぜ必要なのだろう?

ENTRY (__sigsetjmp)
    /* Save registers.  */
    movq %rbx, (JB_RBX*8)(%rdi)
#ifdef PTR_MANGLE
# ifdef __ILP32__
    /* Save the high bits of %rbp first, since PTR_MANGLE will
       only handle the low bits but we cannot presume %rbp is
       being used as a pointer and truncate it.  Here we write all
       of %rbp, but the low bits will be overwritten below.  */
    movq %rbp, (JB_RBP*8)(%rdi)
# endif
    mov %RBP_LP, %RAX_LP
    PTR_MANGLE (%RAX_LP)
    mov %RAX_LP, (JB_RBP*8)(%rdi)
#else
    movq %rbp, (JB_RBP*8)(%rdi)
#endif
    movq %r12, (JB_R12*8)(%rdi)
    movq %r13, (JB_R13*8)(%rdi)
    movq %r14, (JB_R14*8)(%rdi)
    movq %r15, (JB_R15*8)(%rdi)
    lea 8(%rsp), %RDX_LP    /* Save SP as it will be after we return.  */
#ifdef PTR_MANGLE
    PTR_MANGLE (%RDX_LP)
#endif
    movq %rdx, (JB_RSP*8)(%rdi)
    mov (%rsp), %RAX_LP /* Save PC we are returning to now.  */
    LIBC_PROBE (setjmp, 3, LP_SIZE@%RDI_LP, -4@%esi, LP_SIZE@%RAX_LP)
#ifdef PTR_MANGLE
    PTR_MANGLE (%RAX_LP)
#endif
    movq %rax, (JB_PC*8)(%rdi)

#ifdef SHADOW_STACK_POINTER_OFFSET
# if IS_IN (libc) && defined SHARED && defined FEATURE_1_OFFSET
    /* Check if Shadow Stack is enabled.  */
    testl $X86_FEATURE_1_SHSTK, %fs:FEATURE_1_OFFSET
    jz L(skip_ssp)
# else
    xorl %eax, %eax
# endif
    /* Get the current Shadow-Stack-Pointer and save it.  */
    rdsspq %rax
    movq %rax, SHADOW_STACK_POINTER_OFFSET(%rdi)
# if IS_IN (libc) && defined SHARED && defined FEATURE_1_OFFSET
L(skip_ssp):
# endif
#endif
#if IS_IN (rtld)
    /* In ld.so we never save the signal mask.  */
    xorl %eax, %eax
    retq
#else
    /* Make a tail call to __sigjmp_save; it takes the same args.  */
    jmp __sigjmp_save
#endif
END (__sigsetjmp)

PTR_MANGLEは sysdeps/unix/sysv/linux/x86_64/sysdep.h に定義がある。

#  define PTR_MANGLE(reg)    xor %fs:POINTER_GUARD, reg;           \
                rol $2*LP_SIZE+1, reg
#  define PTR_DEMANGLE(reg) ror $2*LP_SIZE+1, reg;                \
                xor %fs:POINTER_GUARD, reg

Linux x86-64では%fsレジスタがThread Local Storageのベースポインタとして使われることを思い出すと、このコードは

• スレッド固有の値(%fs:POINTER_GUARD)とのxorを取り
• $2*LP_SIZE+1だけレジスタのビットを回転させる

ことがわかる。PTR_MAGNLEの詳細は Pointer Encryptionにある。

なぜsetjmp, longjmpに際してポインタを暗号化するのだろうか? 最近のglibcではatexit関数やjmp_bufを狙った攻撃は効かない (PTR_MANGLE)やPointer Subterfugeによれば、書き込み可能な領域に生のEIPを保存すると攻撃対象になるからだそうだ。

結論

• setjmpは%rbx, %r12, %r13, %r14, %r15, %rsp, PCをjmp_bufに保存する
• このとき、生のEIPを保存するとそのデータ自体が攻撃目標になるため、PTR_MANGLEで暗号化する

検索エンジンを自作する夢

2020-03-06 19:25:27 JST · 元記事

そろそろ学生も終わるので、1年前の失敗と反省を書き留めておこうと思います。私達(@a_kawashiroと@gky360)は2018年度の未踏プロジェクトで分野限定型検索エンジンを複数組み合わせた分散型検索エンジンとして採用され、9ヶ月間を検索エンジンの作成に費やしました。この記事ではプロジェクトの紹介と、作成した検索エンジンが満足な性能を出せなかった原因について述べます。

膨大なインターネットの情報の海の中で目的の情報に辿り着くためには、検索エンジンが必要不可欠です。しかし世の中のメジャーな検索エンジンは大企業が提供してくれているものばかりで、ユーザが自分で管理・設定できるものはほとんどありません。それなら一つ作ってみるか、ということで始めたのがこのプロジェクトでした。

検索エンジンは大まかに言って3つの機能で構成されます。

• クロール: インターネット上からウェブページをダウンロードして情報を蓄積する
• 検索: 蓄積した情報のなかからユーザのクエリにマッチするものを探し出す
• 表示: マッチするものを適切に並び替えて表示する

どれも最低限の実装をするならそれほど難しくないので、検索エンジンを作るのはそれほど難しくないように思われます。

しかし、現在のインターネットはあまりに巨大です¹。千里眼やMondouが作られた時代はせいぜい300万サイトしかありませんでしたが、今は16億サイトを超えていて素直な方法で検索エンジンを作ってもまともな結果が帰ってこないのは明らかです。

Googleのような巨大なデータセンターを使えればよいのですが、お金に余裕がなかったので、インターネット上の個人のサーバをうまく協力させて一つの検索エンジンを構成することにしました。このような検索エンジンは分散検索エンジンと呼ばれ、yacyなどがよく知られています。yacyは完全なpeer to peer型の分散検索エンジンですが、私達は

1. 数多くのwebページを分野ごとに分類し各分野ごとの検索エンジンを作り、
2. 分野ごとの検索エンジンを複数つなぎ合わせることでより大きな検索エンジンを作る

という設計を採用することにしました。この設計のメリットは、

• 各分野ごとの検索エンジンはその分野に特化することで、大きなリソースを使わずに精度の高い検索サービスを提供でき、
• 更にそれらを複数つなぎ合わせることでより広範囲の検索サービスを段階的に構築できる

の2点です。完全なpeer to peer型の分散検索エンジンに比べて、ある程度の階層構造があることで検索精度の改善が期待できます。

この設計を元に実装したのが、kearchです。kearchには2つの検索エンジンの実装が含まれます。一つは専門検索エンジン(各分野ごとの検索エンジン)、もう一つはメタ検索エンジン(専門検索エンジンをつなぎ合わせるための検索エンジン)です。

ユーザからのクエリはメタ検索エンジンによって適切な専門検索エンジンに割り振られ、ユーザはその専門検索エンジンの結果をメタ検索エンジンを通して受け取ります。例えばこの検索エンジンの構成で、"pip3 install flask"と検索するとPythonの専門検索エンジンサーバにクエリが飛んでPythonに関する検索結果が返ってきます。

2019年1月から3月にかけて、私達はさくらクラウド上で上の図の検索エンジンを実際に運用してみました。確かに専門検索エンジンもメタ検索エンジンも動くのですが実用には程遠い検索性能でした。その理由は以下の2つです。

• 専門検索エンジンの数が足りない
• 専門検索エンジンのクロール性能が低い

まず、1つ目の「専門検索エンジンの数が足りない」ですが、これはつまり、人間が日常的に検索する分野の範囲が広すぎるということです。ブラウザの履歴を確認してもらうとわかるのですが、1人の人間が1ヶ月で検索する分野は非常に広大です。明らかに4台では足りませんでした。この問題に対する解決策は未だに思いついていません。

2つ目の「専門検索エンジンのクロール性能が低い」は具体的にはwebページの分類性能が低いということです。現在のkearchでは該当分野の単語リストを使ってwebページを分類するのですが、単語リストだけでは十分な性能が出ません。該当分野と全く関係ないwebページがクロールされていることもよくありました。

2つ目の問題にはいくつかの解決策が考えられます。例えば、該当分野の論文を読み込ませて単語リストを拡充したり、信頼できるwebページからのホップ数を考慮することで分類精度を改善できるでしょう。またFessなどのエンタープライズサーチ向けの製品を専門検索エンジンとして使えるように改造するのも解決策の一つだと思います。

最後に、専門検索エンジンについては様々な問題を感じる一方、メタ検索エンジンについてはさほど問題を感じませんでした。「小さな検索エンジンをつなぎ合わせてより大きな検索エンジンを作る」という基本設計は間違っていないと今でも思っています。

formalized-egison -- Egisonの型安全性の証明に向けて

2019-12-23 18:02:13 JST · 元記事

概要

Egisonのパターンマッチの動作を定理証明支援系Coq上で定義し、「Egisonの処理系が項Mを値vに評価する」という命題をCoq上で表現できるようにしました。

動機

Egisonの中核的な機能はユーザが拡張可能なパターンマッチ機構です。例えば以下に示すように、ペアに対してその順序を無視してパターンマッチすることが可能です。

(define $unordered-pair
  (matcher
    {[<pair $ $> [something something]
      {[[$x $y] {[x y] [y x]}]}]
     [$ something
      {[$tgt {tgt}]}]}))

(match-all [1 2] unordered-pair {[<pair $a $b> [a b]]})
; ===> {[1 2] [2 1]}

しかし、その拡張性のためにこのパターンマッチの動作は大変複雑です。Egi, Nishiwaki (2018)はEgisonのパターンマッチの動作をbig-step styleの操作的意味論として次のように定義しています。

型健全性などのEgisonの性質を証明する際はこの操作的意味論を使うことになりますが、規則の種類、数が多いので証明を手で書くと間違えてしまいそうです。そこで今回の記事ではEgi, Nishiwaki (2018)の操作的意味論をCoqに書き直しました。

formalized-egison

formalized-egisonはEgi, Nishiwaki (2018)によるEgisonの操作的意味論をCoqに書き直したものです。意味論を定義している部分はEgison.vの以下の部分です。

Inductive eval : (env * tm * tm)-> Prop :=
  | evarin : forall i Gamma t, (Gamma i) = Some t -> eval (Gamma, (tvar i), t)
  | evarout : forall i Gamma, (Gamma i) = None -> eval (Gamma, (tvar i), (tvar i))
  | eint : forall i e, eval (e, (tint i), (tint i))
  | etpl : forall Gamma t1 t2 v1 v2, eval (Gamma, t1, v1) -> eval (Gamma, t2, v2) -> eval (Gamma, ttpl t1 t2, ttpl v1 v2)
  | ecll : forall e ts vs, same_length_list ts vs ->
                      Forall eval (map (fun tpl => let '(t,v) := tpl in (e,t,v)) (zip ts vs)) -> eval (e, (tcll ts), (tcll vs))
  | epair : forall e t1 t2 v1 v2, eval (e, t1, v1) -> eval (e, t2, v2) -> eval (e, (tpair t1 t2), (tpair v1 v2))
  | esm : forall Gamma, eval (Gamma, tsm, tsm)
  (* | emtc : forall Gamma (ts: (list (pptn * tm * (list (dptn * tm))))), eval Gamma ((tmtc ts), (tmtc vs)) *)
  | etplmtc : forall Gamma t1 t2 v1 v2, eval (Gamma, t1, v1) -> eval (Gamma, t2, v2) -> eval (Gamma, ttplmtc t1 t2, ttplmtc v1 v2)
  | etmal : forall Gamma M N p L v_v v m_m m_e Delta_v, same_length_list Delta_v v_v -> eval (Gamma, M,v) -> evalmtc Gamma N [(m_m, m_e)] -> evalms3 [[([(p,m_m,m_e,v)], Gamma, empty)]] Delta_v ->
                                                  Forall eval (map (fun t => let '(d,v) := t in (Gamma @@ d, L, v)) (zip Delta_v v_v)) ->
                                                  eval (Gamma, (tmal M N (p, L)), tcll v_v)

  with evalmtc : env -> tm -> list (tm * env) -> Prop :=
  | emtcsm : forall Gamma, evalmtc Gamma tsm [(tsm, Gamma)]
  | emtcmtc : forall Gamma l, evalmtc Gamma (tmtc l) [((tmtc l), Gamma)]
  | emtctpl : forall Gamma m1 m2 n1 n2, eval (Gamma, (ttplmtc m1 m2), (ttplmtc n1 n2)) -> evalmtc Gamma (ttplmtc m1 m2) [(n1, Gamma); (n2, Gamma)]

  with evaldp : dptn -> tm -> option env -> Prop :=
  | edpvar : forall z v, value v -> evaldp (dpvar z) v (Some (z |-> v))
  | edppair : forall p1 p2 v1 v2 g1 g2,
      value v1 -> value v2 -> evaldp p1 v1 (Some g1) -> evaldp p2 v2 (Some g2) ->
      evaldp (dppair p1 p2) (tpair v1 v2) (Some (g1 @@ g2))
  | edpfail : forall t p1 p2, not (is_tpair t) -> evaldp (dppair p1 p2) t None

  with evalpp : pptn -> env -> ptn -> option ((list ptn) * env) -> Prop :=
  | eppdol : forall g p, evalpp ppdol g p (Some ([p], empty))
  | eppvar : forall i g m v, eval (g, m, v) -> evalpp (ppvar i) g (pval m) (Some ([], (i |-> v)))
  | epppair : forall pp1 pp2 p1 p2 g pv1 pv2 g1 g2,
                evalpp pp1 g p1 (Some (pv1,g1)) -> evalpp pp2 g p2 (Some (pv2,g2)) ->
                evalpp (pppair pp1 pp2) g (ppair p1 p2) (Some ((pv1 ++ pv2), (g1 @@ g2)))
  | eppvarfail : forall y g p, not (is_pval p) -> evalpp (ppvar y) g p None
  | epppairfail : forall pp1 pp2 p g, not (is_ppair p) -> evalpp (pppair pp1 pp2) g p None

  with evalms1 : ((list ms) * option env * option (list ms)) -> Prop :=
  | ems1nil : evalms1 ([], None, None)
  | ems1anil : forall sv g d, evalms1 ((([],g,d)::sv), (Some d), (Some sv))
  | ems1 : forall p m mg v av g d sv avv d1,
        evalma (g @@ d) (p,m,mg,v) avv d1 ->
        evalms1 ((((p,m,mg,v)::av, g, d)::sv), None, (Some ((map (fun ai => (ai ++ av, g, d @@ d1)) avv) ++ sv)))

  with evalms2 : (list (list ms)) -> (list env) -> (list (list ms)) -> Prop :=
  | ems2 : forall svv gvv svv1 gvv1 svv2,
      same_length_list3 svv gvv svv1 ->
      Forall evalms1 (zip3 svv gvv svv1) ->
      (filtersome gvv) = gvv1 ->
      (filtersome svv1) = svv2 ->
      evalms2 svv gvv1 svv2

  with evalms3 : (list (list ms)) -> (list env) -> Prop :=
  | ems3nil : evalms3 [[]] []
  | ems3 : forall svv gv svv1 dv gdv, evalms2 svv gv svv1 -> evalms3 svv1 dv -> gdv = gv ++ dv ->
                             evalms3 svv gdv

  with evalma : env -> ma -> list (list ma) -> env -> Prop :=
  | emasome : forall x g v d, evalma g (pvar x, tsm, d, v) [[]] (x |-> v)
  | emappfail : forall p g pp m sv pv d v avv g1,
      evalpp pp g p None -> evalma g (p,(tmtc pv),d,v) avv g1 ->
      evalma g (p,tmtc ((pp,m,sv)::pv),d,v) avv g1
  | emadpfail : forall p g pp m dp n sv pv d v pv1 d1 avv g1,
      evalpp pp g p (Some (pv1, d1)) -> evaldp dp v None ->
      evalma g (p, tmtc ((pp,m,sv)::pv),d,v) avv g1 ->
      evalma g (p, tmtc ((pp,m,(dp,n)::sv)::pv),d,v) avv g1
  | ema : forall p Gamma pp M dp N sigma_v Delta v phi1_v p1_v Delta1 Delta2 v1_vv m1_v,
      evalpp pp Gamma p (Some (p1_v, Delta1)) ->
      evaldp dp v (Some Delta2) ->
      eval (Delta @@ Delta1 @@ Delta2, N, tcll v1_vv) ->
      evalmtc Gamma M m1_v ->
      evalma Gamma (p, tmtc ((pp,M,(dp,N)::sigma_v)::phi1_v), Delta, v)
             ((map (fun tpl => match tpl with
                              | (ttpl v11 v12) => map (fun t => let '(v1, (m1, Gamma1), p1) := t in (p1,m1,Gamma1,v1)) (zip3 [v11;v12] m1_v p1_v)
                              | v11 => map (fun t => let '(v1, (m1, Gamma1), p1) := t in (p1,m1,Gamma1,v1)) (zip3 [v11] m1_v p1_v)
                              end
                   ) v1_vv)) empty.

最も重要なのはeval (Gamma, M, v)で「環境Gammaのもとで項Mが値vに評価される」という命題です。例えば、冒頭のプログラム中の(match-all [1 2] unordered-pair {[<pair $a $b> [a b]]})が{[1 2] [2 1]}に評価される、というのは次のようなCoqの命題として表すことができます。

Definition unordered_pair: tm :=
    (tmtc [(pppair ppdol ppdol, ttplmtc tsm tsm,
            [(dppair (dpvar "x") (dpvar "y"), (tcll [(ttpl (tvar "x") (tvar "y")); ttpl (tvar "y") (tvar "x")]))]);
           (ppdol, tsm,
            [(dpvar "tgt", tcll [tvar "tgt"])])]).

Definition match_all_example: tm :=
    (tmal (tpair (tint 1) (tint 2)) unordered_pair (ppair (pvar "a") (pvar "b"),ttpl (tvar "a") (tvar "b"))).

Theorem unordered_pair_example : eval (empty, match_all_example, tcll [ttpl (tint 1) (tint 2);ttpl (tint 2) (tint 1)]).

また、Egisonの操作的意味論上でこの評価結果が正しいということも証明できます。 つまりunordered_pair_exampleを証明することができます。証明は長くなるので省略しますが、Egison.vの230行目以降にあります。

今後の課題

動機の項で述べたようにformalized-egisonの最終的な目標はEgisonの型安全性の証明です。Coqで表現すると

Theorem type_soundness:
    forall Gamma M T, is_typed Gamma M T => exists v, eval Gamma M v /\ is_typed Gamma v T.

となります。ただしis_typed Gamma M Tは環境Gammaの下で項Mに型Tがつくと読みます。 今後はtyped-egisonの型付け規則をCoqで書き直してis_typedを定義した上でtype_soundnessを証明する予定です。

MNISTを可能な限り高速に分類する

2019-07-21 18:54:53 JST · 元記事

概要

• MNISTの分類をする学習済みモデルを軽量化し、
• 更にSIMD命令を使った高速化を行うことで、
• onnxruntimeより高速なMNISTの分類が可能になりました。(シングルスレッドで)

前回までのあらすじ

この記事は前回の続きです。 前回はMNISTを分類する学習済みニューラルネットワークから不要な枝を削除し、軽量化した学習済みモデルを走らせる専用のランタイムsonnxを作ってMNIST分類の高速化を試みました。 今回はSIMD命令とマルチスレッド化による最適化でMNIST分類速度の限界に挑みます。

今回の目標タイムは既存のONNXランタイムonnxruntimeです。 この記事の各実行時間の計測は10回行い、その平均と分散を求めました。 onnxruntimeと最適化無しのsonnxの実行時間は、OS: Ubuntu19.04、CPU: Core i7 8th Gen、 メモリ: 16GiB、GPU無しの環境下で以下のようになりました。¹

機械学習の推論過程、特にMNISTを分類する場合、においてボトルネックになるのは重み行列を入力ベクトルに乗算する処理です。 前回は重み行列の中で絶対値の小さい要素を無視することで計算の効率化を図り、行列の要素の80%を無視して約5倍の高速化に成功しました。 しかし、80%を削減してもなお乗算はボトルネックでした。

具体的にはsonnx.cppのこの部分がボトルネックになります。

for(int i=0;i<n;i++){
        ret[B_row[i]] += B_scale[i] * x[B_column[i]];
}

SIMDによる高速化

まず、sonnxをSIMDで高速化してみます。 Core i7 8th GenではAVX2命令セットが使えるので256bitの演算を一度に行うことができ、 今回は32bit浮動小数点数で計算しているので最大8倍の高速化が見込めます。

しかし、元のコードはメモリへの間接参照を２つ(ret[B_row[i]]とx[B_column[i]])含んでおりそのままではSIMD化するのが困難です。 まずret[B_row[i]]への書き込みは同じB_row[i]の値を持つものをまとめて計算し、メモリへの書き込みを一度に行います。 x[B_column[i]]からの読み出しはAVX2のgather命令を使ってSIMD化します。 完成したコードがこれです。

float r = 0;
__m256 vr = _mm256_setzero_ps();
for(cur;cur<n1;cur+=8){
    __m256i vc = _mm256_loadu_si256((__m256i*)(&B_column_p[cur]));
    __m256 vx = _mm256_i32gather_ps(x_p, vc, 4);
    __m256 vs = _mm256_load_ps(&B_scale_p[cur]);
    vr = _mm256_fmadd_ps(vs, vx, vr);
}
for(cur;cur<n;cur++){
    r += B_scale_p[cur] * x_p[B_column_p[cur]];
}
__attribute__((aligned(32))) float t[8] = {0};
_mm256_store_ps(t, vr);
ret[B_row_p[cur-1]] += r + t[0] + t[1] + t[2] + t[3] + t[4] + t[5] + t[6] + t[7];

実行時間はこんな感じになりました。

ウェルチのt検定を用いて検定すると²有意水準5%で帰無仮説が棄却され、確かにonnxruntimeより高速に推論できています。 onnxruntime(シングルスレッド)と大きな差がつかないのはgather命令が遅いからでしょうか。

マルチスレッドによる高速化

(SIMDなしの)マルチスレッド化はどうでしょうか。 オリジナルのコードはメモリへの書き込みが一箇所(ret[B_row[i]]への書き込み)しかないので、この部分だけ複数スレッドで同時に書き込まないようにすればデータレースは起こりません。 スレッド数は手元で最適なものを探索した結果4にしています。 ソースコードはここにあります。

void CompressedGemm::calc_partially(const int index, const vector<float> &x){
    int m = B_nrows_threads[index].size();
    int cur = 0;
    for(int i=0;i<m;i++){
        int n = B_nrows_threads[index][i];
        float r = 0;
        for(int j=cur;j<cur+n;j++){
            r += B_scale_threads[index][j] * x[B_column_threads[index][j]];
        }
        ret[B_row_threads[index][cur]] += r;
        cur += n;
    }
}

vector<float> CompressedGemm::calc(const vector<float> &x){
    ret = C;
    vector<thread> ths;
    
    for(int i=0;i<n_thread;i++){
        ths.push_back(thread(&CompressedGemm::calc_partially, this, i, x));
    }
    for(int i=0;i<n_thread;i++){
        ths[i].join();
    }
    return ret;
}

結果はこんな感じです。

完敗です。

SIMDとマルチスレッドの併用

では最後にSIMDとマルチスレッディングを併用してみます。 ソースコードはここにあります。

あまり効果がありませんね...複数コアからの結果をまとめるのに時間がかかっているのかもしれません。

まとめ

影響力の小さい要素を削除して学習済みのモデルを軽量化した上でSIMDとマルチスレッド化を用いて推論の高速化を試みました。 マルチスレッド環境ではonnxruntimeに勝てませんでしたが、シングルスレッドではonnxruntimeより高速に推論できました。

参考

• 実行時間の測定データ

WebAssemblyで自作言語用のGCを書く

2019-07-11 18:28:15 JST · 元記事

前回の続きです。

概要

前回の記事ではWebAssemblyを出力するMLのサブセットコンパイラを作りました。しかし、WebAssemblyにはGabage Collection (GC) が未だに実装されていないため(2019/7/8時点)、メモリ管理は全て自分で行う必要があります。前回は超適当mallocで間に合わせていたのですが、今回はmalloc, freeを含むGCをWebAssembly(手書き)で作りました。

GCを含むMLコンパイラのソースコードはここで、GC本体はここです。

WebAssemblyの基礎知識

WebAssemblyはスタックマシンで実行される言語であり、各命令は引数をスタックのトップから取り出し結果をスタックのトップに積みます。トップ以外のスタックの中身を自由に見ることはできません¹。また、スタックの他にもリニアメモリも持ち、i32.store, i32.loadといった命令でアクセスすることができます。リニアメモリは以下のような特徴を持ちます。

• アドレスは0から始まる
• アドレスは32bit
• 整数又は浮動小数点数をstore/loadできる

GCの実装

GCを書く前にまずmallocとfreeが必要です。今回はfree listでブロックを管理する簡単なものを書きました。各ブロックは下図のような構造になっていて、次のブロックのアドレス、ブロックサイズ、フラグ、確保したメモリを保持します。フラグにはfree listの終端かどうかの情報とブロックが使用中かどうかの情報が入っています。それぞれ4byteずつ使うのでmallocのためのヘッダだけで12byteも使います。贅沢ですね。

+-------------------------------------------------------+
| next block address | block size | flags  | contents   |
| 4bytes             | 4bytes     | 4bytes | 4*n bytes  |
+-------------------------------------------------------+

mallocはfree listを先頭から走査して要求されたサイズ以上の空きブロックがあれば確保、最後まで走査して該当するブロックがなければ末尾に新たにブロックを作ります。freeは指定されたブロックの状態を空きブロックに変えるだけです。空きブロック同士を結合する処理はめんどくさいのでやっていません。

このmallocとfreeを使ってGCを作ります。『ガベージコレクション』によれば、GCには次の3種類があります。

• マーク・アンド・スイープ
• コピーGC
• 参照カウント

このうち、マーク・アンド・スイープとコピーGCは生きているオブジェクトを探索するのにルートセットが必要です。ルートセットは一般にレジスタ、スタック、グローバル変数から構成されます。ところが、WebAssemblyではスタックの中身を自由に見ることができないためルートセットが必要なGCは採用できません²。このため今回は参照カウント方式でGCを実装しました。

gc_malloc を呼び出すとmallocが指定されたサイズ+GCのヘッダサイズ分のメモリをリニアメモリ上に確保し、GCの管理情報を書き込んだ上でそのアドレスを返します。GCの管理領域は下図のようになっていて、確保したメモリのサイズ、参照カウントの値、フラグ、オブジェクトのためのメモリ領域から構成されています。フラグには参照カウント時に行う深さ優先探索用のビットと格納しているオブジェクトが値(整数値又は浮動小数点数)か否かの情報が含まれます。

+----------------------------------------------------+
| memory size | reference count | flags  | contents  |
| 4bytes      | 4bytes          | 4bytes | 4*n bytes |
+----------------------------------------------------+

これがmallocのcontentsの中に入っているので実際にはこのようになります。

+---------------+------------+--------+--------------------------------------------------------+
| next malloc   | malloc     | malloc | malloc contents                                        |
| block address | block size | flags  | +----------------------------------------------------+ |
| 4bytes        | 4bytes     | 4bytes | | memory size | reference count | flags  | contents  | |
|               |            |        | | 4bytes      | 4bytes          | 4bytes | 4*n bytes | |
|               |            |        | +----------------------------------------------------+ |
+---------------+------------+--------+--------------------------------------------------------+  

参照カウントはコンパイラが変数のスコープに応じてオブジェクトの参照カウントを増減させるコード(gc_increase_rc, gc_decrease_rc)を挿入することで操作します。基本的にオブジェクトの作成時に参照カウントが1増加し、そのオブジェクトを指す変数のスコープが終了したところで参照カウントが1減らします。参照カウントが0になるとそのオブジェクトは解放され、またそのオブジェクトが保持していたアドレスが指すオブジェクトの参照カウントを1減らします。この操作は再帰的に行われます。ただし、解放されたオブジェクトが値だった場合はこの再帰的な操作は行われません。

↓gc_increase_rc, gc_decrease_rcはこんな感じで挿入される

(i32.const 1)
(i32.store)))
(; let fun_fib_0 end ;)
(call $gc_increase_rc)
(get_local $val_b_6)
(call $gc_decrease_rc)
(drop)
(get_local $val_a_5)
(call $gc_decrease_rc)
(drop)

基本的に参照カウントの操作は変数のスコープと連動させればよいのですが、関数終了時には特殊な処理が必要になります。関数内のローカル変数が指すオブジェクトは関数の終端で解放されますが、関数の戻り値が開放されるオブジェクトを含む場合、すでに解放されたオブジェクトが関数の戻り値に含まれることになります。このようなバグを防ぐためには、関数内のローカル変数を開放する処理を行う前に戻り値オブジェクトの参照カウントをインクリメントしておく必要があります。

その他にも配列の要素を書き換えるときは前のオブジェクトの参照カウントを減らしてから、新しく代入するオブジェクトの参照カウントをインクリメントする必要があるなど、細かいところでいろいろハマりました。

結果

フィボナッチ数列を第10項まで計算するプログラムでプログラム終了時のメモリ使用量をGC有りの場合と無しの場合で結果を比較してみると、以下のようになりました。使用メモリの4分の3ぐらいを解放できているのでそれなりに(僕が)満足できる結果です。

let rec fib n = 
  if n < 3 
    then 1
    else 
        let a = fib (n - 1) in
        let b = fib (n - 2) in
        a + b in
let rec print x =
  if 10 < x
    then 1
    else
      let a = fib x in
      print_i32 a;
      print (x+1) in
print 0

参考リンク

• WebAssemblyでGCを実装する
  GCの管理領域の設計はほとんどこの記事のものをそのまま使っています。
• WebAssemblyでGC
  bitmap式のMark and SweepのGCを書く話。

余談

• コンパイル時にWebAssemblyで書いたGCをリンクする必要があるのですが、スマートな方法が見つからずテキスト表現のまま無理やりリンクしています。lldを使えばいいのではないか?という情報を頂いたのですが、WebAssemblyのテキスト表現をELF形式に変換する方法がわかりませんでした。
• WebAssemblyの直書きは辛い...バグると何も言わずに止まっちゃうし。
• 公式の人は早くGCを組み込んで欲しい。やっぱりスタックを自由に走査できない環境でGCを書くのは無理がある。

ONNXファイルから不要な枝を削ってMNISTの推論を高速化してみる

2019-03-07 20:13:04 JST · 元記事

この記事の中のソースコードは全てhttps://github.com/akawashiro/sonnxにあります。

概要

• ニューラルネットワークから要らなそうな枝を80%削除しても精度が変わらなかった
• ONNXの中身をいじるのが大変だった
• onnxruntimeには勝てなかった

背景

機械学習の学習済みモデルを小さなデバイスで動かす、というのが最近流行っているそうです。機械学習では、学習には大きな計算コストがかかりますが、推論はそれほど大きな計算コストがかかりません。このため、学習だけを別のコンピュータで行っておいて、実際の推論は小さなデバイスで行うということが可能です。

ただし、推論だけでもそれなりに計算資源が必要です。そこで、学習済みのモデルの高速化が重要になります。Raspberry Piに搭載されているGPUを使うIdeinとか有名です。

僕も学習済みモデルの推論を高速化できそうな方法を思いついたので実験してみました。

アイデア

今回はMNISTを分類する学習済みモデルを高速化します。今回使用するモデルは次の図のようなものです。画像は28*28(=784)pxなので入力は784個、出力は各数字の確率なので10個あり、中間層が2つ挟まっています。各層間は全結合しており、活性化関数としてReluを使います。

このモデルを教師データを使って学習すると、枝の重みが変わってこんな感じになります。

僕のアイデアは学習後のネットワークから重みの小さい枝を取り去ってもちゃんと動くんじゃないか、というものです。重みの小さい枝を取り去るとこんな感じになります。

枝の本数が少なくなれば、推論を高速化できそうな気がします。

アイデアの裏付け

実際に学習後のモデルで赤丸で囲った部分の重みの分布を確認します。

分布はこのようになっています。

重みが0の部分が非常に多いです。まず、学習済みモデルから重みが0の枝を削除しても推論結果に影響しないはずです。また、グラフが左右対称になっているので、絶対値の小さい順に削除していけば、各パーセプトロンへの入力はそれほど変化しない気がします。

手法

やることは非常に単純です。ニューラルネットワーク中の各層間枝の重みの統計を取り、重み上位何%かを残して残りを削除します。つまり、何らかの方法で学習済みのモデルから枝の重みのデータを取り出し、枝をカットし、さらに加工後のモデルのデータを使って推論できるようにします。

幸い今はONNXという良いものがあります。ONNXとは学習済みモデルのデータを出力する形式の一つで、多くのフレームワークが対応しています。

今回はChainerで書いたモデルからONNXデータを出力し、そのデータを加工することにしました。また加工後のデータはC++で書いた俺俺ONNXランタイムに実行してもらうことにしました。

纏めると、
1. Chainerでニューラルネットワークを書いて学習する
2. 学習済みのニューラルネットワークからONNXデータを出力する
3. ONNXデータを俺俺ONNXランタイムに読み込んで加工、実行する

となります。一つづつ何をやるのかを説明します。

1. Chainerでニューラルネットワークを書いて学習する

2. 学習済みのニューラルネットワークからONNXデータを出力する

1と2は簡単です。onnx-chainerを使えばすぐにできます。

python3 learn_mnist.py

でmnist.onnxというファイルができます。

3. ONNXデータを俺俺ONNXランタイムに読み込んで加工、実行する

ここが大変でした。ディープラーニングフレームワークからのONNXモデルの出力は多くの人が試しているのですが、出力したONNXモデルをチューニングしようとする人はほとんどいないようです。

3.1 ONNXデータを解析する

とりあえずnetronというONNXの可視化ツールでmnist.onnxを可視化してみました。

GemmはGeneral matrix multiplyの略です。各Gemmノードは行列BとベクトルCを持ち、ベクトルxを入力としてBx+Cを出力します。Reluは活性化関数です。

Reluはmax(0,x)で定義されている関数のでONNXから抽出する必要は無く、各Gemmノードの行列BとベクトルCの情報だけをを抽出できれば良いです。

今回は各GemmノードのBとCをテキストファイルとして抽出します。

> python3 analyze_mnist_onnx.py

は次のファイルを出力します。

• *************_matrix.txt
mnist.onnxの全てのGemmノードのB行列とC行列
• *************_matrix.png
  各行列の中の重みの分布
• mnist.onnx.json
  ONNXファイルをJSONに変換したもの
• mnist_train.txt
• mnist_test.txt
  MNISTをC++から読み込みやすくするためにテキストファイルに変換したもの

各*************_matrix.txtがどのGemmノードに対応するのかはmnist.onnx.jsonを睨むとわかります(←ここ超不親切)。

3.2 抽出したデータを加工、実行する

g++ -O3 -mtune=native -march=native -mfpmath=both sonnx.cpp && ./a.out

で出力した重みのデータを読み込み、推論を実行します。sonnx.cppは簡単なONNXランタイムになっており、mnist.onnxから抽出した行列のデータとMNISTの画像データのテキストファイルから推論を行います。デフォルトではmnist_test.txtの10000枚について推論を行います。

出力はこのようになります。

accuracy, time
0.9817000031, 15.93887043
compress_ratio, accuracy, time
0, 0.9817000031, 36.66616821
0.05, 0.9817000031, 34.61940384
0.1, 0.9818000197, 32.64707184
0.15, 0.9818000197, 30.78090286
0.2, 0.9817000031, 29.01914406
..........................

2行目は読み込んだデータを加工せずで隣接行列表現で保持したときの推論の精度と計算時間、4行目以降は読み込んだデータを加工したときの推論の圧縮率(compress_ratio)、精度と計算時間です。圧縮率が0のときは全く加工しないのと同じ、圧縮率が0.3のときは枝の重みの絶対値が小さい30%を削除したときの結果です。圧縮率が1になると全ての枝が削除されます。

データを加工して枝を削除するとき、sonnx.cppでは行列データを非零成分のインデックスとその値の組の配列として保持しています。これは枝を削除したときに保持するデータ量が減らし、高速に計算するためです。

結果

圧縮率を変化させたときの精度と計算時間のグラフです。圧縮率を0.8まで上げても推論の精度が変わっていません。これは学習済みモデル中の8割の枝を削除しても、推論精度が保てるということです。驚きですね!

一方、計算時間はほぼ枝の数に反比例しています。これは予想通りでした。

表にしてみます。推論精度を1%犠牲にするだけで2倍も高速化できています。やったね!

※ 圧縮率0のときの値は重みのデータを隣接行列表現で保持したときのものです。

おまけ

最後にonnxruntimeで元の学習済みモデルを実行したときの時間を計ってみます。

> python3 analyze_mnist_onnx.py
Accuracy rate =  0.9817 , time =  3.285153865814209

え、3秒?? 普通に考えると16秒ぐらいになるはずです。なぜこんなに速いんだ...

理由は色々あると思いますが、僕が思いつくのは以下の2つです。

• AVXなどのSIMD拡張命令を使っている
  (僕もAVXを使おうとしましたが、scatter命令が僕のPCで使えないので諦めました。)
• CPUのキャッシュに乗るようなプログラムになっている

関連記事・研究

• 深層学習のモデル圧縮・高速化に関する論文80本ノック
  学習済みのモデルを圧縮する論文がまとめられています。このブログ記事の内容はリンク先の「スパースなモデル」の章に該当します。
• モデルアーキテクチャ観点からのDeep Neural Network高速化
  DNNの高速化に関するスライドです。Pruningの章がこのブログ記事に該当します。
• Reducing the Model Order of Deep Neural Networks Using Information Theory
  上のリンク中の論文で特にこの記事に近いものです。このブログ記事でやっているのは、この論文のmagnitude-based pruningです。
• Deep Compression: Compressing Deep Neural Networks with Pruning, Trained Quantization and Huffman Coding
  CPUでの実行速度もオリジナルのモデルより速くなってるっぽくてすごい。何をどう実装したのか知りたい。
• keras-surgeon
  kerasのモデルをpruningするためのソフトウェア。

まとめ

• ONNXランタイムを自作したらナイーブな実装の二倍の速度で推論ができた
• onnxruntimeには勝てなかった

高階関数を基本的な関数の合成で作った関数でQuickCheckする

2018-12-25 22:25:58 JST · 元記事

高階関数をQuickCheckでテストしてみる

QuickCheckを知っていますか?
QuickCheckと言うのはHaskellのデータ駆動型のテスト用ライブラリで
テストしたい関数を指定するとその引数に合わせて適当なテストデータを生成してくれます。

では高階関数(関数を引数に取る関数)をQuickCheckでテストすると
どんなテストデータ(関数)を生成してくれるのでしょうか。
ちょっと確認してみましょう。

-- QCTest.hs
import Test.QuickCheck
import Test.QuickCheck.Function

prop :: Fun Integer Integer -> Bool
prop f = apply f 10 == apply f 20

main = quickCheck prop

% stack runghc QCTest.hs
*** Failed! Falsifiable (after 2 tests and 12 shrinks):
{20->0, _->1}

...なんだこの関数は?
入力と出力の組が列挙されてるだけやんけ。
僕の思ってる関数と違うんですが。
map (+2)とかtailとかそういうやつを生成してほしいなぁ。

まともなテスト用関数を生成する

QuickCheckのテストデータが気に入らないので、まともな(主観)テストデータを生成してみました。
今回生成するテストデータは[Int] -> [Int]の関数に限定します。

そもそも関数を生成するってどうやればいいんでしょうか?
QuickCheckのように入力と出力の対を列挙するのは筋が悪い気がします。
mapや(+)、tailといった基本的な関数を合成して関数を生成したいのです。
我々がプログラミングするときってそうやりますよね。
入力と出力の組を列挙する人はあんまりいないと思います。

[Int] -> [Int]型の関数を列挙してみます。

tail
reverse
map (* 10)
map (+ 10)

tailやreverseはそもそも[Int] -> [Int] 型です。
mapは(Int -> Int) -> [Int] -> [Int] 型なので
(+ 10) :: Int -> Intを渡すと[Int] -> [Int]型になります。

こんな風にどの関数(例 map)にどの関数(例 (+ 10))を適応すれば
どんな型(例 [Int] -> [Int])になるかを規則として書き出します。
型付け規則っぽく書くとこんな感じです。

ちょっとわかりにくいかも知れません。

[Int] -> [Int]型の関数を2つ作ってみました。

左はmap (+ 10)のような関数、右はmap (+ 10 (* 5))のような関数を表しています。
一番下の型にある関数の規則を持ってきて、上に必要な型を書いて
さらにその必要な型をもつ関数の規則を持ってきて...という感じです。

とりあえずこの規則を使って関数を生成できるようになりました。
[Int] -> [Int]型の関数を生成した例がこちらです。

% stack exec qcfun-exe
> Test datum are following.
> ["tail","reverse","(map (* 75))","tail","reverse"]
...

今見せた例はこのようなデータ型で表現されています。

data QProg = QMap1 QProg | QMap2 QProg QProg 
             | QTail1 QProg | QTail | QRev | QMult QProg | QAdd QProg | QRand Int deriving (Eq)

instance Show QProg where
  show (QMap1 p) = "(map " ++ show p ++ ")"
  show (QMap2 p1 p2) = "(map " ++ show p1 ++ " " ++ show p2 ++ ")"
  show (QTail1 p) = "(tail " ++ show p ++ ")"
  show (QTail) = "tail"
  show QRev = "reverse"
  show (QMult p) = "(* " ++ show p ++ ")"
  show (QAdd p) = "(* " ++ show p ++ ")"
  show (QRand i) = show i

QProg型がプログラム(関数)を表す型です。
このQProg型に対して適当なShowインスタンスを定義してあげると関数っぽく見えるようになります。

作ったテストデータ用の関数で高階関数をテストする

次に作った関数を使って高階関数をテストしてみましょう。
今回テストする対象のプログラムはこれです。

func f = reverse . f
prop f = (func f [1,2,3]) == [1,2,3]

関数fを引数に取る簡単な高階関数funcの性質をテストします。

しかし作ったテストデータ用の関数はすべてQProg型です。
つまりQProg型をなんとかして[Int] -> [Int]型に変換したうえで
テストしたい高階関数に適用する必要があります。

今回はTemplate Haskellを使ってQProg型を[Int]->[Int]型のHaskellの関数に変換しました。
Template Haskellというのはプログラムの中でプログラムを作るためのツールです。
Lensなんかで使われています。 動作例がこちらです。

% stack exec qcfun-exe
...
Test results are following.
[(False,"tail"),(True,"reverse"),(False,"(map (* 75))"),(False,"tail"),(True,"reverse")]

テストデータは先頭から順に["tail","reverse","(map (* 75))","tail","reverse"]だったので正常に動いているようです。

まとめ

高階関数をまともな(主観)関数でテストできるようにしました。
しかし技術的な制約のためテストデータをコンパイル時に生成するので
テストデータを完全にランダムに生成することはできませんでした。

参考

http://haskell.g.hatena.ne.jp/mr_konn/20111218/1324220725

付録

githubへのリンクです
https://github.com/akawashiro/qcfun
以下のコマンドで実行できます。

git clone https://github.com/akawashiro/qcfun.git
stack build
stack exec qcfun-exe

ソースコードはTemplate Haskell的な事情で2つにわかれています。

WebAssemblyを出力するMinCamlコンパイラを実装しました

2018-10-31 21:14:24 JST · 元記事

概要

WebAssemblyを出力するMinCamlコンパイラml2wasmをフルスクラッチで実装しました。
github.com

マンデルブロ集合を計算するこんな↓感じのMinCamlのソースコードが

こんな↓感じのWebAssemblyに変換されて

実行して、適切にプロットするとこんな↓感じになります。

導入

WebAssemblyとは

WebAssemblyとはブラウザで実行可能な低級プログラミング言語であり、仮想的なスタックマシン上で動作します。 また、JavaScriptに比べて構文解析が容易であり高速に動作します。 WebAssemblyは最近のほとんどのブラウザで動かすことができます。*1

;; WebAssembly のテキスト表現の例
(module
  (func (export "add") (result i32)
    (i32.add
      (i32.const 1)
      (i32.const 3))))

MinCamlとは

MinCamlとは東北大学の住井先生*2が設計したプログラミング言語MLのサブセットです。 言語仕様が小さく、処理系を制作するのが比較的容易なのが特徴です。 住井先生は未踏事業としてMinCamlのコンパイラも制作されています。 このコンパイラは読みやすくドキュメントもしっかりしているので、コンパイラを作ってみたい人には非常にオススメです。 ml2wasmの開発でもこのコンパイラを大いに参考にしました。

(* MinCaml の例 *)
let rec add = 1 + 3 in add

コンパイラの作成

今回はMinCamlのソースコードをWebAssemblyのテキスト表現に変換するコンパイラml2wasmを作成しました。 基本的にml2wasmはオリジナルのMinCamlコンパイラと同じ構造ですが、WebAssemblyがスタックマシンであることを利用していくつかの変換の工程が省略されています。

ml2wasmの実装の詳細

ml2wasmでは次の5段階の工程でMinCamlをWebAssemblyに変換します。

1. パース
2. アルファ変換
3. 型推論
4. クロージャ変換
5. コード生成

特に説明のない部分はオリジナルのMinCamlコンパイラと同じです。 今回はHaskellで実装しましたが、もともとの実装言語がOCamlでHaskellとそれほど変わらないので、その点は苦労しませんでした。 以下、WebAssembly特有のコンパイラ実装事情をいくつか説明します。

オリジナルのMinCamlコンパイラには3と4の間にK正規化、4と5の間にレジスタ割当という処理がありましたが、ml2wasmではこれら2つの処理を省略しています。 K正規化はWebAssemblyでは命令のオペランドとしてスタックの先頭の値を使い、オペランドをレジスタに格納する必要がないため、 レジスタ割当はWebAssemblyにはレジスタが存在しないため、それぞれ省略できました。

逆に、WebAssembly特有の難しい工程としてはコード生成があります。 通常、関数型言語のコンパイラではプログラムをクロージャ変換し、関数から自由変数を削除します。 その後コード生成の工程で関数呼び出しは間接ジャンプ(MIPSならjr命令)に変換されます。

この間接ジャンプをWebAssemblyではcall_indirect命令で実現します。

;; call_indirect命令はこんな感じで呼び出される。
...
(i32.load)
(call_indirect (param f32) (param f32) (param i32) (result i32)))))
...
;; (param f32) (param f32) (param i32) (result i32) ← 型が必要

ここで問題になるのがcall_indirect命令が呼び出す関数の型をオペランドとして要求する点です。 これはクロージャ変換後の中間言語に型情報を残しておく必要があることを意味します。 コンパイラを書き始めたときにはこのことに気づかなかったので、あとから型情報を追加するのに苦労しました。

他には公式のリファレンスがどこにあるのかイマイチわからない問題もあります。 結局これを使ったのですが、 URLに公式感がなく最後までこのリファレンスが公式のものなのかはわからずじまいでした。

結果

マンデルブロ集合を出力するプログラムをコンパイルできるようになりました! やったね!

マンデルブロ集合を出力するMinCamlのソースコード
mandelbrot.ml · GitHub
コンパイル結果
mandelbrot.wast · GitHub

ちなみにマンデルブロ集合っていうのはこんな↓やつです。

今回作ったコンパイラでコンパイルし実行した結果(を適切にプロットしたもの)は↓

まあ大体あってますね。点の数が以上に少ないのは大量の点をプロットしようとするとすぐにメモリ不足に陥るからです。

今後の課題

ml2wasmで出力したWebAssemblyコードはすぐにメモリ不足に陥ります。 理由は単純でGCがないからです。 WebAssemblyでGCを実装する試みとしては κeenさんのWebAssemblyでGC | κeenのHappy Hacκing Blog がありますが、今回はめんどくさいのでメモリは確保しっぱなしです。

そのうちWebAssemblyにデフォルトでGCが入るとかいう噂も聞くので、そのときにGCを入れればいいやと思っています。

参考リンク

• WebAssembly
  WebAssemblyの公式ページ
• 速攻MinCamlコンパイラ概説
  MinCaml
• WebAssemblyでGC | κeenのHappy Hacκing Blog
  κeenさんによるWebAssemblyでのGCの実装
• wasm-reference-manual/WebAssembly.md at master · sunfishcode/wasm-reference-manual · GitHub
  WebAssemblyの仕様書(?)。多分公式ではないが一番役にたった

ネットワークデバイスドライバを一からビルドしてインストールした

2018-09-28 09:44:50 JST · 元記事

概要

新し目のコンピュータにDebianを入れたら、ネットワークデバイスドライバが入ってなくて結構大変だった。

経緯

最近、自分の使っていたコンピュータ(6年前くらいのデスクトップ)に限界を感じ始めたので、日本橋で新しいものを買ってくることにした。購入に当たっては、換装が難しいCPUをCorei7 8700kに決定し、フィーリングで選ぶことにした。日本橋を端から端まで歩き回って、結局LM-iH700XD1-EX4を買った。

Windows10をDebian streachで上書きして起動すると、ネットワークドライバが入ってないことがわかった。マザーボードが新しすぎるのでドライバがDebianに取り込まれていないのだろう。こういう場合は自前でドライバをビルドする必要がある。

NIC(Network Interface Card)の型番を調べる

まず、NICの型番を調べる。商品のスペックを確認したが、「LAN:10/100/1000BASE-T LAN(オンボード)」としか書いていない。マザーボードの方には「マザーボード:Intel B360 Micro ATX LGA1151」と書いてあった。オンボードのNICなのだから、たぶん「Intel B360」のほうが関係しているのだろうと考えられる。ここから「Intel B360 network driver」などで検索して一時間ほど頑張ったが、なんの成果も得られなかった。

次はLinux側からNICを調べてみる。

% lspci
(省略)
00:1f.6 Ethernet controller: Intel Corporation Device 15bc (rev 10)

なるほど、「Intel Corporation Device 15bc linux driver」で検索すれば良さそうである。検索すると、Linux Kernel Driver DataBaseがトップに出てきた。「15bc 」でページ検索すると

vendor: 8086 ("Intel Corporation"), device: 15bc ("Ethernet Connection (7) I219-V")

が出てくる。「I219-V linux driver」で検索してみると「e1000e」の最新のドライバを入れればいいらしい?

そこで、もうどこで見たのか忘れたけど Intel Ethernet Drivers and Utilities の最新版をインストールすればいいらしい。(検索しすぎて記憶がないorz)

makeを入れる

というわけでe1000e-3.4.2.1.tar.gzをビルドする。 ネットワークに繋がらないので別のパソコンでファイルをダウンロードしてUSBメモリでコピーした。

展開して、make installすると...makeコマンドが入ってないと言われた。 Debianのインストール元のlive USBをマウント後、apt install build-essentialで入れた。 最初から入れといてほしい...

デバイスドライバをビルドする

e1000e-3.4.2.1.tar.gzを展開して、make installするとカーネルのヘッダファイルが足りないと言われた。

% uname -a
Linux gley 4.9.0-8-amd64 #1 SMP Debian 4.9.110-3+deb9u4 (2018-08-21) x86_64 GNU/Linux

4.9.0-8-amd64のヘッダファイルが必要なようだ。別のパソコンで該当のdebファイルを落としてきてdpkg -iで入れる。ちなみにcommonとamd64の両方が必要だった。

その後、再度make install。 再起動するとネットワークに繋がっていた。

ネットワークデバイスドライバをインストールできた〜〜〜〜嬉しい〜〜〜〜

— a_kawashiro (@a_kawashiro) September 27, 2018

まとめ

新し目のパソコンでLinuxを動かすのは大変。 ただドライバをビルドしてインストールする作業は初めてで勉強になった。 もうやりたくない。

ちなみにオーディオデバイスドライバが入っていないので音は出ない。 そのうち頑張る。

このマザーボード新しすぎて音も出せないwww

— a_kawashiro (@a_kawashiro) September 27, 2018

プログラミング言語Egisonの型システムを設計するインターンをした

2018-09-22 07:25:07 JST · 元記事

概要

2018年3月11日から4月14日まで楽天技術研究所でEgisonの型システムの設計 & 型検査器を書くインターンをしていた。 もうインターンから半年も経ってしまったが、何も書かないよりはましというわけでブログ記事を書くことにした。

実装した型検査器はここにある。 まだまだ不完全ではあるが、一応動くのでよかったら試してみてほしい。

Egisonとは

Egisonとは江木さんという方が開発しているLisp likeなプログラミング言語で、非常に強力なパターンマッチ機構を持つ言語である。 江木さんは2007年からほぼ一人でEgisonを開発し続けており、その歴史はここで見ることができる。 最近ではパターンマッチ機構を活かして、数式処理システムの制作に使われている(江木さんがつかっている)。

インターンが始まるまで

楽天技術研究所では，現代数学を簡潔に記述するためのプログラミング言語EgisonをHaskellで開発する方(アルバイトやインターンなど)を募集しています．
Haskellで言語実装することに興味ある方は，ぜひDMください！

— Egison (@Egison_Lang) 2017年11月13日

採用が決まった後は、Egisonが微分幾何の記述に使われるということなので局面と曲線の微分幾何を読んでいた。 研究室でPPL2018に行ったら隣の席に江木さんがいて、そのとき初めて江木さんに会った。

インターン初週

インターン初日、東京で疲れてしまったのか、高熱と吐き気がして会社を休んだ。 社会人失格かと思ったが、次の日出社したら「大丈夫ですか」と気遣ってくれて嬉しかった。 後でTwitterを見た感じ、東京で風邪が流行っていたらしい。

インターン初週はEgisonの実装を読むことが主な業務だった。 Egisonは巨大な構文木インタープリタとして実装されており、相次ぐ機能の拡張によってエライことになっていることが分かった。

インターン2週目

江木さんがEgisonに型検査を入れたいというので、Egisonのサブセットを抽出、定式化をして、そのサブセットに対して型システムを設計しようという話になった。 Egisonにはmatcherやprimitive patternなどの通常の言語には無い構文要素があり、ちょっと大変だった。

この週の金曜日、江木さんが東大の出身研究室(萩谷研)を訪問するというのでくっついて行った。 萩谷研では、Egisonの意味論をどう定義するかという話をしていた(僕は議論する様子を眺める係をしていた)。 この議論は後々論文になり、APLAS 2018に採択されている。

インターン3週目

2週目に設計した型システムを型検査器として実装した。 基本的にはHM型の制約解消アルゴリズムで実装したが、 型付けできない式を受理するためにワイルドカード型を入れたりした(Egisonはもともと型なしの言語なので、型付けできないが動く式が基本ライブラリ中に存在する)。 当時は知らなかったのだが、あとで確認してみると漸進的型付けの型検査と同じような実装になっていた。

インターン4週目

引き続き型検査器を実装した。 matcherを正しく型付けできたときがとても嬉しかったのを覚えている。

↑インターン中の食事の一例

インターン5週目

成果発表のためのスライド制作&発表練習をした。 スライドが全直しになったりしていい経験になった。 この経験が後々、未踏の二次審査で役に立ったような気がする。 最終発表の後で江木さんと二人でもんじゃ焼きを食べにいった。

感想

インターンを通して基本的に江木さんと二人だけで仕事をしていたが、江木さんがとても良い方だったのであまり苦ではなかった。 また週に一度kotapikuさんや梅崎さんがアルバイトに来て、数学の話をしたりするのが楽しかった。 食堂も大変美味しいのでついつい食べすぎてしまい、かなり太った。

ISUCON8予選に参加した

2018-09-17 01:32:10 JST · 元記事

ISUCON8の予選にokeigoさん、mayokoさんとチーム「しょラーさんのおかげ」として参加していた。 最終的な得点は8000点くらいで予選落ちだった。悔しい。

チームメイトの参加記

brookbach.com

予選まで

ISUCON出てみたいけどその分野について全く知らない

— a_kawashiro (@a_kawashiro) July 15, 2018

ご検討ください https://t.co/fOPEOwDecN

— しょラー (@shora_kujira16) July 15, 2018

ISUCONの参加が決定してからは手元でISUCON7の予選、ISUCON6の予選を立てて練習していた。 N+1クエリ問題について初めて知った。 あとsystemd、MySQL、nginx、pythonのプロファイラの設定も初めてでいろいろ勉強になった。

予選の一週間前にokeigoさんがISHOCON2をセットアップしてくれて、そのサーバを使って3人で練習した。 3人とも離れたところに住んでいるので、Discordで音声通話、Slackでチャット、GitHubでコードを共有することになった。 言語は3人とも読めるPythonを使うことにした。

予選当日

朝

6時ぐらいに目が覚めた。いつもより6時間くらい早起き。

10時〜

予選スタート。 設定ファイルやpythonのアプリの実装をデプロイするスクリプトを書いた。 h2oをnginxに切り替えるのはokeigoさんがやってくれた。 この間にmayokoさんがアプリを読んでいた。

11時〜

とりあえずベンチマークを取った。 /favicon.icoが異常に重かったのだが、okeigoさんが設定のtypoを見つけて直したら解消した。

13時〜

僕はapp.pyの中でsheetをランダムに並び替えているところを小手先の小技で高速化しようとしていた。 が、あまり早くならなかった。 mayokoさんとokeigoさんがapp.pyのget_eventの改善をやっていた。

14時〜

mayokoさんの実装が終わった。 マージしてベンチを取ったら5000点ぐらい出て、学生3位になったので正直勝ったと思った。

15時〜

okeigoさんがredisを使ってキャッシュを実装していた。 venv環境にどうやってredisをインストールするのかがわからずに悪戦苦闘した。 結局バイナリをコピーして対応したらしい。SUGOI

16時〜

キャッシュの実装が動いたので、そろそろ3台構成にしましょうか、という話になった。 が、ポートの開放とかでいろいろ詰まった。 僕はあんまり何もできなかった。

17時〜

17:30ぐらいに3台構成が動いた。 あとは何度がベンチマークをとって終了。 最終的には7900点ぐらいでたはず。

18時〜

お祈りした。

真のお祈りタイム #isucon

— a_kawashiro (@a_kawashiro) September 16, 2018

19時ごろ

結果が出た。予選を通過することはできなかった。

燃え尽きたので起き上がれない

— a_kawashiro (@a_kawashiro) September 16, 2018

反省

• okeigoさんの負担が重すぎた。もうちょっとサーバの設定とか練習しとけばよかった。
• 役割分担がきちんと決まっていなかった。アプリ2人、インフラ1人というのがバランスが良さそうだった。
• 複数台構成におけるサーバ間の通信についてわかっていなかった。ISUCONのような特殊な環境ではファイヤーウォールを切って、全サービスが0.0.0.0:portをリッスンしてしまえば良いと思った。
• サービスとして登録するpythonスクリプトに新しくライブラリ(redis)を追加する処理についてわかっていなかった。
• プロファイルをちゃんと見ていなかった。ログの可視化スクリプトを作成しておくべきだった。

まとめ

チカラが足りなかった。悔しい。

Vimで初めてコマンドを定義してみた

2018-02-02 15:25:23 JST · 元記事

書かないと忘れちゃうから…

これだけ

:command! MT :MerlinTypeOf
:MT

Vimをコンパイルしてインストールする

2017-11-23 20:55:08 JST · 元記事

使っているプラグインの関係で+python,+lua,+perlなVimが必要になりました。
しかし、Ubuntuのパッケージで条件を満たすものが見当たりません。
そこで自分で一からコンパイルしてインストールしました。

僕は自分でコンパイルしたソフトウェアを$HOME/.local/以下にインストールしているので、
--prefix=$HOME/.localとしています。
必要に応じて適当に変更してください。

git clone https://github.com/vim/vim
cd vim
./configure --prefix=$HOME/.local/ --with-features=huge --enable-multibyte --enable-gpm --enable-cscope --enable-fail-if-missing --enable-perlinterp --enable-python3interp --enable-rubyinterp
make install

あとは$HOME/.local/bin/vimにパスを通すだけです。

nvimに移行したほうがいいのかなぁ。

githubのバックアップを取る

2017-11-23 02:41:39 JST · 元記事

年末も近くなってきたのでgithubが落ちた時に備えてリポジトリのバックアップを取ることにしました。
以下のスクリプトを実行するとgithub以下にすべてのリポジトリがバックアップされます。

YOUR_ACCESS_TOKENはSign in to GitHub · GitHubで作れます。
Select scopesのところでreposの欄をすべてチェックしてトークンを作成してください。

#! /bin/bash

token=YOUR_ACCESS_TOKEN
curl "https://api.github.com/user/repos?access_token=${token}&per_page=100&affiliation=owner" 2>/dev/null -o tmp
repos=`cat tmp | grep ssh_url | cut -d" " -f 6 | sed -e 's/"//g' | sed -e 's/,//g'`
mkdir github
cd github
for r in ${repos}; do
    name=`echo ${r} | sed -e "s:git.*/::" -e "s:[.]git::"`
    git clone ${r}
    cd ${name}
    for branch in `git branch -a | grep remotes | grep -v HEAD | grep -v master `; do
        git branch --track ${branch#remotes/origin/} $branch
    done
    cd ..
done
rm tmp

これでgithubが落ちても安心ですね。

Freeモナドについて今更勉強してみた

2017-10-20 00:01:23 JST · 元記事

5年ほど前、Freeモナドという言葉がHaskell界隈でよく見られました。
僕のフォローしている人たちが次々に解説記事をアップしていたように思います。
5年たった今になって、Freeモナドについて興味が湧いてきたので記事を書いてみることにしました。

Freeモナドの定義に行く前にFree型について説明します。

-- The kind of f is * -> * like Maybe, [].
-- The kind of a is * like Int, Char, [Char].
data Free f a = Pure a | Free (f (Free f a)) 

ここでfはカインド * -> * を持つ型コンストラクタであり、あとで述べるようにFunctorのインスタンスです。
aはカインド * のただの型です。
fが[]、aがIntの場合について調べてみましょう。

free1 :: Free [] Int
free1 = Pure 10

free2 :: Free [] Int
free2 = Free [Free [Free [Free []]]]

free3 :: Free [] Int
free3 = Free [Pure 10]

これら3つの定数関数はすべてFree [] aの型を持ちます。
直感的に言うとリストもリストのリストのリストもただの整数もすべて同じ型になっています。
これだけではFree型をどう使っていいのかわかりません。

これをモナドにしてみましょう。

-- functor is typeclass which "fmap" is defined
-- (<$>) :: Functor f => (a->b) -> f a -> f b
-- (<$>) = fmap
instance Functor f => Functor (Free f) where
  fmap f (Pure a) = Pure (f a)
  fmap f (Free fa) = Free (fmap (fmap f) fa)

instance Functor f => Monad (Free f) where
    return = Pure
    Pure a >>= k = k a
    Free fm >>= k = Free (fmap (>>= k) fm)

モナドになると何が良いのでしょうか?
簡単なゲームを作る例を見ながら、Freeモナドの利点を確認しましょう。

このゲームは

1. 入力を受け取る
2. 文字列を出力する
3. 一秒待つ

の3つの要素の組み合わせで構成されています。
ゲームを表すデータ型がこれです。

data Game cont = GetInput (String -> cont) | Render String cont | Wait cont

instance Functor Game where
  fmap f (GetInput g) = GetInput (f . g)
  fmap f (Render str cont) = Render str (f cont)
  fmap f (Wait cont)     = Wait (f cont)

各コンストラクタは動作を行ったあとの続きの計算を持ちます。
例えば、GetInput fだったら入力を受け取ってfに渡すことで
続きの計算contを計算し、contを実行します。

ここでGameデータ型をFreeモナドと組み合わせて使う利点を説明します。

game1 :: Free Game ()
game1 = Free (Wait (Free (Render "Hello" (Pure ()))))

game1' :: Free Game ()
game1' = (Free (Wait (Pure ()))) >>  (Free (Render "Hello" (Pure ())))

一秒待って"Hello"と出力するゲームを2つの方法で記述してみました。
game1は愚直な記述、game1'は>>演算子を用いた記述です。

実はこの2つは全く同じものです。
game1'を定義にしたがって展開してみます。

(Free (Wait (Pure ()))) >>  (Free (Render "Hello" (Pure ())))
==> (Free (Wait (Pure ()))) >>= (\x -> (Free (Render "Hello" (Pure ()))))
==> Free (fmap (>>= (\x -> (Free (Render "Hello" (Pure ()))))) (Wait (Pure ())))
==> Free (Wait ((>>= (\x -> (Free (Render "Hello" (Pure ()))))) (Pure ())))
==> Free (Wait ((Pure ()) >>= (\x -> (Free (Render "Hello" (Pure ()))))))
==> Free (Wait (Free (Render "Hello" (Pure ()))))

同じになりましたね。
モナドの枠組みを使えばGameデータ型同士を簡単に合成できることがわかります。

input :: (String -> Free Game ()) -> Free Game ()
input f = Free (GetInput f)

render :: String -> Free Game ()
render str = Free (Render str (Pure ()))

wait :: Free Game ()
wait = Free (Wait (Pure ()))

miniGame = do
  render "Please input your name."
  input (\x -> render $ "Your name is " ++ x)
  wait
  render "End"

miniGame関数ではdo記法で3種類のGameコンストラクタを合成しています。
あとはFree Game ()をIO ()に変換する、つまり実際の入出力を行う部分を書けば完成です。

runGame :: Free Game () -> IO ()
runGame (Pure ()) = return ()
runGame (Free (Wait cont)) = runGame cont
runGame (Free (GetInput f)) = getLine >>= \s -> runGame (f s)
runGame (Free (Render str cont)) = putStrLn str >> runGame cont

main :: IO ()
main = runGame miniGame

コンパイルしてmain関数を実行すると名前を聞いてくれます。(ソースコードはこの記事の末尾にあります。)

結局Freeモナドを使うとどんな利点があるのでしょうか。
この記事の例ではゲームの構造と入出力が分離できています。
つまり、テストがしやすいということが言えます。
またゲームの構造だけ作っておいてあとで入出力の部分を差し替えたりするのも簡単です。
Freeモナドを使った2Dゲーム用のライブラリというのもあります。

参考

Freeモナドって何なのさっ！？ - capriccioso String Creating(Object something){ return My.Expression(something); }

そろそろFreeモナドに関して一言いっとくか - モナドとわたしとコモナド

以下ソースコード

module Main where

-- The kind of f is * -> * like Maybe, [].
-- The kind of a is * like Int, Char, [Char].
data Free f a = Pure a | Free (f (Free f a)) 

-- functor is typeclass which "fmap" is defined
-- (<$>) :: Functor f => (a->b) -> f a -> f b
-- (<$>) = fmap
instance Functor f => Functor (Free f) where
  fmap f (Pure a) = Pure (f a)
  fmap f (Free fa) = Free (fmap (fmap f) fa)

-- applicative is typeclass which derived from functor
-- and offer "pure" and <*>
-- (<*>) :: f (a -> b) -> f a -> f b
instance Functor f => Applicative (Free f) where
  pure = Pure
  Pure a <*> Pure b = Pure $ a b
  Pure a <*> Free mb = Free $ fmap a <$> mb
  Free ma <*> b = Free $ (<*> b) <$> ma

instance Functor f => Monad (Free f) where
    return = Pure
    Pure a >>= k = k a
    Free fm >>= k = Free (fmap (>>= k) fm)

-- free1 :: Num a => Free f a
free1 :: Free [] Int
free1 = Pure 10

free2 :: Free [] Int
free2 = Free [Free [Free [Free []]]]

free3 :: Free [] Int
free3 = Free [Pure 10]

data Game cont = GetInput (String -> cont) | Render String cont | Wait cont

instance Functor Game where
  fmap f (GetInput g) = GetInput (f . g)
  fmap f (Render str cont) = Render str (f cont)
  fmap f (Wait cont)     = Wait (f cont)

game1 :: Free Game ()
game1 = Free (Wait (Free (Render "Hello" (Pure ()))))

game1' :: Free Game ()
game1' = (Free (Wait (Pure ()))) >> (Free (Render "Hello" (Pure ())))

input :: (String -> Free Game ()) -> Free Game ()
input f = Free (GetInput f)

render :: String -> Free Game ()
render str = Free (Render str (Pure ()))

wait :: Free Game ()
wait = Free (Wait (Pure ()))

miniGame = do
  render "Please input your name."
  input (\x -> render $ "Your name is " ++ x)
  wait
  render "End"

runGame :: Free Game () -> IO ()
runGame (Pure ()) = return ()
runGame (Free (Wait cont)) = runGame cont
runGame (Free (GetInput f)) = getLine >>= \s -> runGame (f s)
runGame (Free (Render str cont)) = putStrLn str >> runGame cont

main :: IO ()
main = runGame miniGame

型レベル自然数を使って型付き行列を作ってみた

2016-12-06 21:24:31 JST · 元記事

#動機 研究で行列計算ライブラリを使っているのですが、numpyはなぜか異なる型の行列同士を足せてしまいます。 このせいでデバッグがめちゃくちゃ大変でした。

こんな単純な間違いはコンパイル時に弾いてくれればいいのにと思ったので、違った型の行列は足せないようなちょっとしたライブラリを作ってみました。

※行列を安全に計算したい人のためには

Numeric.LinearAlgebra.Static

↑のライブラリが既に存在しました。不勉強ですみません。この記事は「型レベル自然数ってこんなに簡単に(僕でも)使えるんだ」、程度に読んでください。

型レベル変数

最近のGHCでは、自然数をパラメーターとするデータ型を定義できる。固定長リストの長さを型に持たせるとか、行列の大きさを型に持たせるとか、そういうことができる

引用元:Haskell (GHC) の型レベル自然数とリフレクションを試してみる

行列の大きさの情報を型に持たせることができる。 →やってみました

仕様

• 行列の足し算・掛け算をサポートする・
• 足し算/掛け算できないものを計算しようとしたらコンパイル時にエラーにする

実装

{-# LANGUAGE DataKinds, KindSignatures #-}
{-# LANGUAGE TemplateHaskell #-}

module SizeFixedMatrix where
import GHC.TypeLits
import Data.List.Split
import Data.List

data SizeFixedMatrix (h :: Nat) (w :: Nat) = MakeMartix [[Integer]] deriving (Eq,Show)

mat1 = MakeMartix [[1,2],[2,3]] :: (SizeFixedMatrix 2 2)
mat2 = MakeMartix [[2,9],[1,3],[4,3]] :: (SizeFixedMatrix 3 2)
mat3 = MakeMartix [[5,9],[3,3],[0,3]] :: (SizeFixedMatrix 3 2)


add :: SizeFixedMatrix p q -> SizeFixedMatrix p q -> SizeFixedMatrix p q
add (MakeMartix m1) (MakeMartix m2) = MakeMartix (map (\(a,b) -> map (\(c,d) -> c+d) (zip a b)) (zip m1 m2))

mul :: SizeFixedMatrix p q -> SizeFixedMatrix q r -> SizeFixedMatrix p r
mul (MakeMartix m1) (MakeMartix m2) = MakeMartix $ splitEvery (length (head m2)) [sum $ zipWith (*) v w | v <- m1, w <- transpose m2]

…これだけです。拍子抜けするくらい簡単に作れました。もう少し凝ったものを作りたかったのですが、締切に間に合いそうにありませんでした。

一応ちょこっとだけ解説しておくと

SizeFixedMatrix p q

が縦p行横q列の行列の型になります。

不適切な型同士でaddしたりmulしたりするとコンパイルエラーになるので安心してプログラミングできます。

結果

$stack ghci
(省略)
*Main FiniteField Lib SizeFixedMatrix> mat1
MakeMartix [[1,2],[2,3]]
*Main FiniteField Lib SizeFixedMatrix> mat2
MakeMartix [[2,9],[1,3],[4,3]]
*Main FiniteField Lib SizeFixedMatrix> mat3
MakeMartix [[5,9],[3,3],[0,3]]

足し算

*Main FiniteField Lib SizeFixedMatrix> add mat1 mat2

<interactive>:4:10: error:
    • Couldn't match type ‘3’ with ‘2’
      Expected type: SizeFixedMatrix 2 2
        Actual type: SizeFixedMatrix 3 2
    • In the second argument of ‘add’, namely ‘mat2’
      In the expression: add mat1 mat2
      In an equation for ‘it’: it = add mat1 mat2

mat1は2行2列の行列型、mat2は3行2列の行列型なので足せないと言われてしまいます

*Main FiniteField Lib SizeFixedMatrix> add mat2 mat3
MakeMartix [[7,18],[4,6],[4,6]]

掛け算

*Main FiniteField Lib SizeFixedMatrix> mul mat1 mat2

<interactive>:8:10: error:
    • Couldn't match type ‘3’ with ‘2’
      Expected type: SizeFixedMatrix 2 2
        Actual type: SizeFixedMatrix 3 2
    • In the second argument of ‘mul’, namely ‘mat2’
      In the expression: mul mat1 mat2
      In an equation for ‘it’: it = mul mat1 mat2
*Main FiniteField Lib SizeFixedMatrix> 

mat1*mat2(行列積)は存在しません

*Main FiniteField Lib SizeFixedMatrix> mul mat2 mat1
MakeMartix [[20,31],[7,11],[10,17]]

感想

型レベル変数などという難しい名前の機能をはじめて使いました。案外簡単に使えたので、他にも応用できると思います。

参考文献

• Haskell (GHC) の型レベル自然数とリフレクションを試してみる

• hstackを教えてもらったツイート

@igrep @a_kawashiro hmatrixのNumeric.LinearAlgebra.Staticがそんな感じです。 https://t.co/e66NSsTNDr

— Masahiro Sakai (@masahiro_sakai) December 15, 2016

CTF入門記

2015-12-12 00:31:26 JST · 元記事

11月の終わりからCTFをかじり始めました。この熱意がいつまでつづくのか分からないのですがとりあえずこの二週間くらいでやったことを書いてみようと思います。

絶望

別にPCを初めて触るわけでもないし、ksnctfのポイントの低い問題ならなんとかなるやろ…

わからなすぎてペンペン草も生えない

— a_kawashiro (@a_kawashiro) 2015, 11月 30

舐めてました。難しすぎる。というわけで梅田の丸善ジュンク堂で

をお買い上げてしまいました。その後も

POSTリクエストの中身をいじって送信するだけのことができない

— a_kawashiro (@a_kawashiro) 2015, 12月 1

SQLインジェクションって存在するんだ(無知)

— a_kawashiro (@a_kawashiro) 2015, 12月 4

などと(初心者なので)当然すぎる知識不足を露呈し続けました。

n番煎じのyesod入門ーとりあえずプロジェクトを作ってみる

2013-12-16 00:24:12 JST · 元記事

cabal install yesod-platform yesod-bin

cabal install cabal-dev
mkdir -p ~/haskell/yesod
cd /haskell/yesod
cabal-dev install yesod-platform yesod-bin

./cabal-dev/bin/yesod init

cd myFirstYesod/
../cabal-dev/bin/yesod --dev devel

cabal: Cannot find the program 'ghc' at 'yesod-ghc-wrapper' or on the path

cabal-dev/bin

cabal-dev/lib

Resolving dependencies...
Configuring myFirstYesod-0.0.0...
cabal: At least the following dependencies are missing:
aeson -any,
conduit >=1.0,
fast-logger >=2.0,
hamlet ==1.1.*,
以下略

Haskellでインベーダーゲームを作ってみた

2013-09-12 15:27:03 JST · 元記事

夏休みがあまりに暇なのでインベーダーゲームまがいの物を作ってみました。
ゲーム全体の状態を持つ変数をIORefで包み
一定時間毎に更新することでゲームの状態を変更・保持しています。
(IORefを使うとIOモナドの中で変数が更新できます。（実はよく分かってないのですが…））
　　
基本的な四角形の描画、テクスチャの貼り付けは↓のコピペです。
HaskellとOpenGLを使って四角形を描画しテクスチャを貼り付ける - 春まで冬眠します
　　
スクリーンショット↓

　　　

import Graphics.UI.GLUT hiding (Bitmap)
import qualified Graphics.Rendering.OpenGL.GL as GL
import System.Exit ( exitFailure, exitWith, ExitCode(ExitSuccess) )
import Data.IORef
import Codec.Picture.Repa
import Data.Array.Repa as R hiding (reshape , map )
import qualified Data.Array.Repa.Repr.ForeignPtr as RF
import Foreign.ForeignPtr
import Data.Word()
import Control.Applicative
import Control.Monad
import Unsafe.Coerce

data Fighter = Fighter { xfighter::GLfloat , yfighter::GLfloat }
data Bullet = Bullet { xbullet::GLfloat , ybullet::GLfloat ,lbullet::Int  }
data Enemy = Enemy { xenemy::GLfloat , yenemy::GLfloat , lenemy::Int }
data Flame = Flame { xflame::GLfloat , yflame::GLfloat , lflame::Int }
data State = Start | Main deriving (Eq)
data Game = Game { fighter::Fighter,bullets::[Bullet],enemies::[Enemy],flames::[Flame],state::State,ftex::TextureObject,btex::TextureObject,etex::TextureObject,fltex::TextureObject,sttex::TextureObject }

cz :: GLfloat
cz = 50 -- character size

timerInterval :: Timeout
timerInterval = 40

inite h = map (\i->(Enemy i h 1)) [100,170..500]

main :: IO ()
main = do 
    
    --GLUTの初期化
    initialDisplayMode $= [RGBAMode, DoubleBuffered]
    initialWindowSize $= Size 640 480
    initialize "" []

    --ウィンドウを作る
    createWindow "invader"

    --テクスチャを読み込む
    texbullet <- loadTextureFromFile "./bullet.jpg"
    texfighter <- loadTextureFromFile "./fighter.jpg"
    texenemy <- loadTextureFromFile "./enemy.jpg"
    texflame <- loadTextureFromFile "./flame.jpg"
    texstart <- loadTextureFromFile "./start.jpg"
    texture Texture2D $= Enabled
    GL.blend $= GL.Enabled
    GL.blendFunc $= (GL.SrcAlpha, GL.OneMinusSrcAlpha)
    
    game <- newIORef (Game (Fighter 50 50) [(Bullet (1000) 50 0),(Bullet (1000) 50 0),(Bullet (1000) 50 0)] (concat [inite 450,inite 380,inite 310]) [] Start texfighter texbullet texenemy texflame texstart)
    

    displayCallback $= display game
    reshapeCallback $= Just reshape
    keyboardMouseCallback $= Just (keyboardProc game)
    addTimerCallback timerInterval $ timerProc (display game)
    mainLoop

processBullet :: [Enemy] -> [Bullet] -> [Bullet]
--processBullet es bs = ((map f1) . (map f2) . (map f3)) bs
processBullet es bs = map (f1 . f2 . f3) bs
    where f1 b = b{ybullet=ybullet b+4}
          f2 b = if (ybullet b)>505 then b{lbullet=0} else b
          f3 b = if any (\e->g1 (lenemy e)>0&&g1 (xbullet b-xenemy e)<cz/2&&g1 (ybullet b-yenemy e)<cz/2) es
                    then b{lbullet=0}
                    else b
          g1 f = if f>0 then f else -f

processEnemy :: [Bullet] -> [Enemy] -> [Enemy]
--processEnemy bs es = ((map f4).(map f1).(map f2)) es
processEnemy bs es = map (f4 . f1 . f2) es
    where f1 e = if lenemy e`div`20`mod`2==0 then e{xenemy=xenemy e-2} else e{xenemy=xenemy e+2}
          f2 e = if any (\b->(lbullet b)>0&&g1 (xbullet b-xenemy e)<cz/2&&g1 (ybullet b-yenemy e)<cz/2) bs
                    then e{lenemy=0}
                    else e     
          g1 f = if f>0 then f else -f
          f4 e = if lenemy e>0 then e{yenemy=yenemy e-0.2,lenemy=lenemy e+1} else e

processFlame :: [Bullet] -> [Enemy] -> [Flame] -> [Flame]
processFlame bs es fs = ((Prelude.++ f4) . (filter f2) . (map f1)) fs
    where f1 f = f{lflame=lflame f-1}
          f2 f = (lflame f) > 0
          f3 (b,e) = if (lbullet b)>0&&(lenemy e)>0&&
                        g1 (xbullet b-xenemy e)<cz/2&&g1 (ybullet b-yenemy e)<cz/2
                        then [Flame {xflame=xenemy e,yflame=yenemy e,lflame=10}]
                        else []
          f4 = concat (map f3 g2)
          g1 f = if f>0 then f else -f
          g2 = concat (map (\b-> (map (\e->(b,e)) es)) bs)
          
processGame :: Game -> Game
processGame g = if state g == Main 
                   then g{bullets=processBullet (enemies g) (bullets g),
                              enemies=processEnemy (bullets g) (enemies g),
                              flames=processFlame (bullets g) (enemies g) (flames g)}
                   else g     

display :: IORef Game -> IO ()
display game = do
    g <-readIORef game
    modifyIORef game processGame
    clearColor $= Color4 0.0 0.0 0.0 0.0
    clear [ColorBuffer]
    loadIdentity
    
    renderGame g
    swapBuffers

renderGame :: Game -> IO()
renderGame g = if state g == Main 
    then do
         mapM_ (renderBullet (btex g)) (bullets g)
         mapM_ (renderEnemy (etex g)) (enemies g)
         mapM_ (renderFlame (fltex g)) (flames g)
         render (xfighter (fighter g)) (yfighter (fighter g)) (ftex g)
    else do
         currentColor $= Color4 1 1 1 1 
         textureBinding Texture2D $= Just (sttex g)
         preservingMatrix $ do
           -- translate (Vector3 0 0 (0::Double))
             renderPrimitive Quads $ do
                 texCoord2f (TexCoord2 0 0)
                 vertex3f (Vertex3 0 0 0.0)
                 texCoord2f (TexCoord2 0 1)
                 vertex3f (Vertex3 0 480 0.0)
                 texCoord2f (TexCoord2 1 1)
                 vertex3f (Vertex3 640 480 0.0)
                 texCoord2f (TexCoord2 1 0)
                 vertex3f (Vertex3 640 0 0.0)        

renderBullet :: TextureObject -> Bullet -> IO ()
renderBullet texbullet b = if lbullet b > 0 then render (xbullet b) (ybullet b) texbullet
                                            else return ()     

renderEnemy :: TextureObject -> Enemy -> IO ()          
renderEnemy texenemy e = if lenemy e > 0 then render (xenemy e) (yenemy e) texenemy
                                         else return ()     

renderFlame :: TextureObject -> Flame -> IO()
renderFlame texflame f = if lflame f > 0 then render (xflame f) (yflame f) texflame
                                         else return ()     

render :: GLfloat -> GLfloat -> TextureObject -> IO ()
render x y tex = do
    currentColor $= Color4 1 1 1 1 
    textureBinding Texture2D $= Just tex
    preservingMatrix $ do
        -- translate (Vector3 0 0 (0::Double))
        renderPrimitive Quads $ do
            texCoord2f (TexCoord2 0 0)
            vertex3f (Vertex3 (x-cz/2) (y-cz/2) 0.0)
            texCoord2f (TexCoord2 0 1)
            vertex3f (Vertex3 (x-cz/2) (y+cz/2) 0.0)
            texCoord2f (TexCoord2 1 1)
            vertex3f (Vertex3 (x+cz/2) (y+cz/2) 0.0)
            texCoord2f (TexCoord2 1 0)
            vertex3f (Vertex3 (x+cz/2) (y-cz/2) 0.0)

texCoord2f = texCoord :: TexCoord2 GLfloat -> IO ()
vertex3f = vertex :: Vertex3 GLfloat -> IO ()
          

--タイマが呼ばれるたびにactを繰り返す
timerProc act = do
    act
    addTimerCallback timerInterval $ timerProc act

--ウィンドウのサイズが変更された時の処理
reshape :: Size -> IO()
reshape (Size w h)=do
    viewport $= (Position 0 0, (Size w h)) --ウィンドウ全体を使う
    matrixMode $= Projection
    loadIdentity
    ortho 0.0 640.0 0.0 480.0 (-1000.0) 1000.0
    matrixMode $= Modelview 0 

keyboardProc :: IORef Game -> Key -> t -> t1 -> t2 -> IO ()
keyboardProc game ch st _ _
    | ch == SpecialKey KeyLeft  = modifyIORef game (\g->g{fighter=changeX (-2) (fighter g)})
    | ch == SpecialKey KeyRight = modifyIORef game (\g->g{fighter=changeX 2 (fighter g)})
    | ch == Char ' ' = space (readIORef game)
    | otherwise            = return ()  
    where changeX dx f = if 0+cz<=xfighter f+dx&&xfighter f+dx<=640-cz 
                            then  f { xfighter = xfighter f + dx }
                            else f
          shot = modifyIORef game (\g->g{bullets=shotProcess (fighter g) (bullets g)})
          space g' = do
              g <- g'
              if state g == Main 
                  then shot 
                  else modifyIORef game (\g->g{state=Main})
    
          
shotProcess :: Fighter -> [Bullet] -> [Bullet]
shotProcess f bs = if (all (\b->ybullet b-50>yfighter f) bs) && (any (\b->lbullet b==0) bs)
                       then concat [[(head (filter (\b->lbullet b==0) bs)){xbullet=xfighter f,ybullet=yfighter f,lbullet=1}],
                                    (tail (filter (\b->lbullet b==0) bs)),
                                    (filter (\b->lbullet b/=0) bs)]
                       else bs


-- テクスチャをファイルから読み込む
-- 何をやっているのか未だによくわからない
-- FreeGameのサンプルをちょっといじったもの
loadTextureFromFile :: FilePath -> IO GL.TextureObject
loadTextureFromFile path = do
    content <- delay <$> (flipVertically.imgData) <$> either error id <$> (readImageRGBA path)
    let (Z :. width :. height :. _) = R.extent content
    [tex] <- GL.genObjectNames 1
    GL.textureBinding GL.Texture2D GL.$= Just tex
    GL.textureFilter Texture2D $= ((Nearest, Nothing), Nearest)
    fptr <- liftM RF.toForeignPtr $ R.computeP $ content
    withForeignPtr fptr
        $ GL.texImage2D Nothing GL.NoProxy 0 GL.RGBA8 (GL.TextureSize2D (gsizei width) (gsizei height)) 0
        . GL.PixelData GL.RGBA GL.UnsignedInt8888
    return tex

gsizei :: Int -> GL.GLsizei
{-# INLINE gsizei #-}
gsizei x = unsafeCoerce x

Codeforces Round #176 div 2

2013-03-24 08:35:41 JST · 元記事

久しぶりに（一年ぐらい）Codeforcesに出ました。
正解したのはＡ問題のみというひどい結果でした。
Ｂ問題はコーナーケースで落ち、Ｃ問題以降は
手も足も出ませんでした。
もう一度初心に帰って勉強する必要性を感じました。

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

using namespace std;

bool check(vector<string> &vs){
     for(int i=0;i<3;i++)
	  for(int j=0;j<3;j++)
	       if(vs[i][j]==vs[i+1][j]&&vs[i][j]==vs[i][j+1]&&vs[i][j]==vs[i+1][j+1])
		    return true;
     return false;
}

char opposite(char c){
     if(c=='#')return '.';
     else return '#';
}

int main(){
     vector<string> vs;
     for(int i=0;i<4;i++){
	  string s;
	  cin>>s;
	  vs.push_back(s);
     }
     string ans="NO";
     if(check(vs))
	  ans="YES";
     for(int i=0;i<4;i++)
	  for(int j=0;j<4;j++){
	       vs[i][j]=opposite(vs[i][j]);
	       if(check(vs))
		    ans="YES";
	       vs[i][j]=opposite(vs[i][j]);
	  }
     cout<<ans<<endl;
     return 0;
}

#include <iostream>
#include <algorithm>

using namespace std;
typedef long long LL;

int main(){
     LL n,k;
     cin>>n>>k;
     if(n==1){
	  cout<<0<<endl;
	  return 0;
     }
     if(n<=k){
	  cout<<1<<endl;
	  return 0;
     }
     LL u=k-1,d=0;
     for(int i=0;i<100;i++){
	  LL m=(u+d)/2;
	  LL a=(k-m)*(k+m-1)/2;
	  if(n-1<=a)
	       d=m;
	  else
	       u=m;
     }
     // cout<<"u="<<u<<" d="<<d<<endl;
     if(n-1<=(k-u)*(k+u-1)/2)
	  cout<<k-u<<endl;
     else if(n-1<=(k-d)*(k+d-1)/2)
	  cout<<k-d<<endl;
     else
	  cout<<-1<<endl;
     return 0;
}

Marathon Match 78 の記録

2013-03-20 20:48:25 JST · 元記事

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <string>
#include <vector>
#include <queue>
#include <map>
#include <set>
#include <algorithm>
#include <functional>
#include <cstring>
#include <sstream>

#define FOR(i,k,n)  for (int i=(k); i<(int)(n); ++i)
#define REP(i,n)    FOR(i,0,n)
#define FORIT(i,c)	for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)
#define SZ(i) ((int)i.size())
#define pb          push_back
#define mp          make_pair
#define ALL(X)      (X).begin(),(X).end()
typedef long long LL;

// #define DFS_DEPTH 11
#define MAX_ANS_LENGTH 20200
#define INF 100000000

using namespace std;

int dfsDepth=10;
int mv[4][2]={{-1,0},{1,0},{0,-1},{0,1}};
int c2m(char c){if(c=='U')return 0;if(c=='D')return 1;if(c=='L')return 2;if(c=='R')return 3;return INF;}
int iabs(int n){if(0<=n)return n;else return -n;}
vector<string> dia;
int H,W;
int isVisited[110][110];
int nVisit[110][110];
int nTouch[110][110];
int bfs[110][110];
string udlr="UDLR";
string snake;
int point;
pair<pair<int,int>,string> ans;
int ystart,xstart,yhead,xhead,ygoal,xgoal;
class BIT{
public:
     int bit[MAX_ANS_LENGTH+10];
     int sum(int i){
	  i++;
	  int s=0;
	  while(0<i){
	       s+=bit[i];
	       i-=i&-i;
	  }
	  return s;
     }
     void add(int i,int x){
	  i++;
	  while(i<=MAX_ANS_LENGTH+5){
	       bit[i]+=x;
	       i+=i&-i;
	  }
     }
};
BIT turnCount;    

int rorlturn(char s,char t){//left turn = 1 right turn = -1
     if(s=='U'){if(t=='L')return 1;if(t=='R')return -1;}
     if(s=='L'){if(t=='D')return 1;if(t=='U')return -1;}
     if(s=='D'){if(t=='R')return 1;if(t=='L')return -1;}
     if(s=='R'){if(t=='U')return 1;if(t=='D')return -1;}
     return 0;
}

bool isInLoop(const int w,const int v,const char n,const int y,const int x,const char m,const int nv){
     if(y==0||y==H||x==0||x==W)return false;
     //右回りのループを検出
     int yy=y,xx=x;
     if(m=='U')xx--;
     if(m=='D')xx++;
     if(m=='L')yy++;
     if(m=='R')yy--;
     if(0<isVisited[yy][xx]&&turnCount.sum(nv)-turnCount.sum(nVisit[yy][xx]+1)<=-3)
	  return true;
     //左回りのループを検出
     yy=y,xx=x;
     if(m=='U')xx++;
     if(m=='D')xx--;
     if(m=='L')yy--;
     if(m=='R')yy++;
     if(0<isVisited[yy][xx]&&turnCount.sum(nv)-turnCount.sum(nVisit[yy][xx]+1)>=3)
	  return true;
     int yf=w+mv[c2m(n)][0],xf=v+mv[c2m(n)][1];
     if(yf==0||yf==H||xf==0||xf==W)return false;
     //右回りのループを検出
     yy=yf,xx=xf;
     if(n=='U')xx--;
     if(n=='D')xx++;
     if(n=='L')yy++;
     if(n=='R')yy--;
     if(0<isVisited[yy][xx]&&turnCount.sum(nv-1)-turnCount.sum(nVisit[yy][xx]+1)<=-3&&rorlturn(n,m)==-1)
     	  return true;
     //左回りのループを検出
     yy=yf,xx=xf;
     if(n=='U')xx++;
     if(n=='D')xx--;
     if(n=='L')yy--;
     if(n=='R')yy++;
     if(0<isVisited[yy][xx]&&turnCount.sum(nv-1)-turnCount.sum(nVisit[yy][xx]+1)>=3&&rorlturn(n,m)==1)
     	  return true;
     return false;
}

pair<pair<int,int>,string> DFS(const int y,const int x,const string s,const int p,const char l){
     if(dfsDepth==SZ(s)){
	  return make_pair(make_pair(p,iabs(ygoal-y)+iabs(xgoal-x)),s);
     }else{
	  vector<pair<pair<int,int>,string> > res;
	  for(int i=0;i<4;i++){
	       int yy=y+mv[i][0];
	       int xx=x+mv[i][1];
	       if(yy<0||H<yy||xx<0||W<xx)
		    continue;
	       if(ygoal==H&&xgoal==W&&yy<xx)
	       	    continue;
	       if(ygoal==H&&xgoal==W&&yy==xx&&i==0)
	       	    continue;
	        if(ygoal==H&&xgoal==W&&y==x+1&&i==2)
	       	    continue;
	       vector<int> tpos;
	       string ss=s+udlr[i];
	       char ll=udlr[i];
	       if(i==0){
		    if(0<x) tpos.push_back((y-1)*1000+x-1);
		    if(x<W) tpos.push_back((y-1)*1000+x);
	       }else if(i==1){
		    if(0<x) tpos.push_back(y*1000+x-1);
		    if(x<W) tpos.push_back(y*1000+x);
	       }else if(i==2){
		    if(0<y) tpos.push_back((y-1)*1000+x-1);
		    if(y<H) tpos.push_back(y*1000+x-1);
	       }else{
		    if(0<y) tpos.push_back((y-1)*1000+x);
		    if(y<H) tpos.push_back(y*1000+x);
	       }
	       int pp=p;
	       for(int j=0;j<SZ(tpos);j++){
		    int yyy=tpos[j]/1000,xxx=tpos[j]%1000;
		    if(nTouch[yyy][xxx]==dia[yyy][xxx]-'0')   pp--;
		    if(nTouch[yyy][xxx]+1==dia[yyy][xxx]-'0') pp++;
		    nTouch[yyy][xxx]++;
	       }
	       isVisited[yy][xx]++;
	       int oldnVisit=nVisit[yy][xx];
	       nVisit[yy][xx]=SZ(snake)+SZ(ss);
	       turnCount.add(SZ(snake)+SZ(ss),rorlturn(l,ll));
	       if((x==0&&i==0)||(y==H&&i==2)||(x==W&&i==1)||(y==0&&i==3)){
		    ;
	       }else if(ygoal==yy&&xgoal==xx){
		    if(2<snake.size()+ss.size()&&ans.first.first<point+pp)
			 ans=make_pair(make_pair(point+pp,iabs(ygoal-yy)+iabs(xgoal-xx)),snake+ss);
	       }else if(1<isVisited[yy][xx]){
		    ;
	       }else{
		    if(!isInLoop(y,x,l,yy,xx,udlr[i],SZ(snake)+SZ(ss))){
		     	 res.push_back(DFS(yy,xx,ss,pp,ll));
		    }
	       }
	       for(int j=0;j<SZ(tpos);j++){
		    int yyy=tpos[j]/1000,xxx=tpos[j]%1000;
		    nTouch[yyy][xxx]--;
	       }
	       isVisited[yy][xx]--;
	       nVisit[yy][xx]=oldnVisit;
	       turnCount.add(SZ(snake)+SZ(ss),-1*rorlturn(l,ll));
	  }
	  if(SZ(res)==0)
	       return make_pair(make_pair(-INF,iabs(ygoal-y)+iabs(xgoal-x)),"");
	  else{
	       sort(res.begin(),res.end());
	       return res.back();
	  }
     }
}

string solve(void){
     memset(turnCount.bit,0,sizeof(turnCount));
     turnCount.add(100,1);
     if(H<=20){
	  dfsDepth=12;
     }else if(H<=30){
	  dfsDepth=11;
     }else if(H<=40){
	  dfsDepth=10;
     }else if(H<=50){
	  dfsDepth=10;
     }else if(H<=60){
	  dfsDepth=9;
     }else if(H<=70){
	  dfsDepth=9;
     }else if(H<=80){
	  dfsDepth=9;
     }else if(H<=90){
	  dfsDepth=8;
     }else if(H<=100){
	  dfsDepth=8;
     }

     point=0;
     FOR(i,0,H)FOR(j,0,W)if(dia[i][j]=='0')point++;
     memset(isVisited,0,sizeof(isVisited));
     memset(nVisit,-1,sizeof(nVisit));
     memset(nTouch,0,sizeof(nTouch));
     ystart=0;xstart=0;isVisited[0][0]=1;nVisit[0][0]=0;
     yhead=1;xhead=0;isVisited[1][0]=1;nVisit[1][0]=1;
     ygoal=H;xgoal=W;
     snake="D";
     ans=make_pair(make_pair(-INF,iabs(ygoal-yhead)+iabs(xgoal-xhead)),"");
     nTouch[0][0]=1;
     if(dia[0][0]=='0')point--;
     if(dia[0][0]=='1')point++;
     for(int time=0;time<MAX_ANS_LENGTH/dfsDepth/2;time++){
	  pair<pair<int,int>,string> res=DFS(yhead,xhead,"",0,snake[SZ(snake)-1]);
	  if(res.first.first==-INF)break;
	  point+=res.first.first;
	  snake+=res.second;
	  FOR(i,0,SZ(res.second)){
	       char ch=res.second[i];
	       if(ch=='U')yhead--;
	       if(ch=='D')yhead++;
	       if(ch=='L')xhead--;
	       if(ch=='R')xhead++;
	       isVisited[yhead][xhead]++;
	       nVisit[yhead][xhead]=SZ(snake)-SZ(res.second)+i+1;
	       turnCount.add(SZ(snake)-SZ(res.second)+i,rorlturn(snake[SZ(snake)-SZ(res.second)+i-1],snake[SZ(snake)-SZ(res.second)+i]));
	  } 
     }
     fprintf(stderr,"ans.first.first=%d ans.second=\"%s\"\n",ans.first.first,ans.second.c_str());
     if(ans.first.first==-INF){
	  stringstream tss;
	  tss<<0<<" "<<0<<" "<<snake<<endl;
	  return tss.str();
     }

     pair<pair<int,int>,string> fans=ans;

     memset(isVisited,0,sizeof(isVisited));
     memset(nTouch,0,sizeof(nTouch));
     memset(nVisit,-1,sizeof(nVisit));
     memset(turnCount.bit,0,sizeof(turnCount.bit));
     yhead=0,xhead=0;
     isVisited[yhead][xhead]++;
     nVisit[yhead][xhead]=0;
     FOR(i,0,SZ(fans.second)){
	  char ch=fans.second[i];
	  if(ch=='U')yhead--;
	  if(ch=='D')yhead++;
	  if(ch=='L')xhead--;
	  if(ch=='R')xhead++;
	  isVisited[yhead][xhead]++;
	  nVisit[yhead][xhead]=i+1;
	  if(0<i)
	       turnCount.add(i,rorlturn(fans.second[i-1],fans.second[i]));
     }
     point=0;
     for(int y=0;y<H;y++)
	  for(int x=0;x<W;x++){
	       nTouch[y][x]+=(isVisited[y][x]&&isVisited[y+1][x]);
	       nTouch[y][x]+=(isVisited[y+1][x]&&isVisited[y+1][x+1]);
	       nTouch[y][x]+=(isVisited[y+1][x+1]&&isVisited[y][x+1]);
	       nTouch[y][x]+=(isVisited[y][x+1]&&isVisited[y+1][x]);
	       point+=(nTouch[y][x]==dia[y][x]-'0');
	  }

     snake=fans.second+"U";
     ystart=H;xstart=W;isVisited[H][W]=1;
     turnCount.add(SZ(snake)-1,rorlturn('R','U'));
     yhead=H-1;xhead=W;isVisited[H-1][W]=1;nVisit[H-1][W]=snake.size();
     ygoal=0;xgoal=0;
     ans=make_pair(make_pair(-INF,iabs(ygoal-yhead)+iabs(xgoal-xhead)),"");
     for(int time=0;time<MAX_ANS_LENGTH/dfsDepth/2;time++){
	  pair<pair<int,int>,string> res=DFS(yhead,xhead,"",0,snake[SZ(snake)-1]);
	  if(res.first.first==-INF)break;
	  point+=res.first.first;
	  snake+=res.second;
	  FOR(i,0,SZ(res.second)){
	       char ch=res.second[i];
	       if(ch=='U')yhead--;
	       if(ch=='D')yhead++;
	       if(ch=='L')xhead--;
	       if(ch=='R')xhead++;
	       isVisited[yhead][xhead]++;
	       nVisit[yhead][xhead]=SZ(snake)-SZ(res.second)+i+1;
	       turnCount.add(SZ(snake)-SZ(res.second)+i,rorlturn(snake[SZ(snake)-SZ(res.second)+i-1],snake[SZ(snake)-SZ(res.second)+i]));
	  }
     }

     fprintf(stderr,"ans.first.first=%d ans.second=\"%s\"\n",ans.first.first,ans.second.c_str());
     fprintf(stderr,"snake=\"%s\"\n",snake.c_str());
    
     if(ans.first.first==-INF){
	  stringstream tss;
	  tss<<0<<" "<<0<<" "<<snake<<endl;
	  return tss.str();
     }

     stringstream ss;
     ss<<0<<" "<<0<<" "<<ans.second;
     return ss.str();
}

class FixTheFence{
public:
     string findLoop(const vector<string>& __diagram){
	  dia=__diagram;
	  H=SZ(dia);
	  W=SZ(dia[0]);
	  return solve();
     }
};

HaskellとOpenGLを使って四角形を描画しテクスチャを貼り付ける

2013-03-04 17:29:58 JST · 元記事

こんな感じになります↓

結構苦労しました。誰かの役に立ってくれると嬉しいです。

参考サイト
　　いろいろなサンプルがある：Index of /GLUT/examples/RedBook

　　　　今回役に立ったのは Tex〜.hs

　　視野の設定：Lichu's_Base
　　テクスチャの読み込み・設定：https://github.com/fumieval/free-game/blob/master/Graphics/FreeGame/Backends/GLFW.hs
　　

import Graphics.UI.GLUT hiding (Bitmap)
import qualified Graphics.Rendering.OpenGL.GL as GL
import System.Exit ( exitFailure, exitWith, ExitCode(ExitSuccess) )
import Data.IORef
import Codec.Picture.Repa
import Data.Array.Repa as R hiding (reshape)
import qualified Data.Array.Repa.Repr.ForeignPtr as RF
import Foreign.ForeignPtr
import Data.Word
import Control.Applicative
import Control.Monad
import Unsafe.Coerce

-- テクスチャをファイルから読み込む
-- 何をやっているのか未だによくわからない
-- FreeGameのサンプルをちょっといじったもの
loadTextureFromFile :: FilePath -> IO GL.TextureObject
loadTextureFromFile path = do
    content <- delay <$> (flipVertically.imgData) <$> either error id <$> (readImageRGBA path)
    let (Z :. width :. height :. _) = R.extent content
    [tex] <- GL.genObjectNames 1
    GL.textureBinding GL.Texture2D GL.$= Just tex
    GL.textureFilter Texture2D $= ((Nearest, Nothing), Nearest)
    fptr <- liftM RF.toForeignPtr $ R.computeP $ content
    withForeignPtr fptr
        $ GL.texImage2D Nothing GL.NoProxy 0 GL.RGBA8 (GL.TextureSize2D (gsizei width) (gsizei height)) 0
        . GL.PixelData GL.RGBA GL.UnsignedInt8888
    return tex

gsizei :: Int -> GL.GLsizei
{-# INLINE gsizei #-}
gsizei x = unsafeCoerce x

--タイマの間隔
timerInterval = 40

main = do 
    --回転の角度を初期化
    rot <- newIORef 0.0
    arg <- newIORef 10.4
    
    --GLUTの初期化
    initialDisplayMode $= [RGBAMode, DoubleBuffered]
    initialWindowSize $= Size 640 480
    initialize "" []

    --ウィンドウを作る
    createWindow "DrawPictureWithOpneGL"

    --テクスチャを読み込む
    tex <- loadTextureFromFile "読み込みたい画像ファイルへのパス"
    texture Texture2D $= Enabled
    GL.blend $= GL.Enabled
    GL.blendFunc $= (GL.SrcAlpha, GL.OneMinusSrcAlpha)


    --表示に使うコールバック関数の指定
    displayCallback $= display rot arg tex
    
    --ウィンドウのサイズが変更された時に呼ぶコールバック関数の指定
    reshapeCallback $= Just reshape
    
    --キーボードやマウスのコールバック
    keyboardMouseCallback $= Just (keyboardProc arg)
    
    --タイマを作る
    addTimerCallback timerInterval $ timerProc (display rot arg tex)

    --GLUTのメインループに入る
    mainLoop

display rot arg texname = do
    --回転させる
    w<-readIORef arg 
    --w <- readIORef hoge
    modifyIORef rot (+w)
    r <- readIORef rot
 
    --背景を黒にする
    clearColor $= Color4 0.0 0.0 0.0 0.0
    clear [ColorBuffer]
    
    --単位行列を読み込む
    loadIdentity
    
    -- texCoord,vetexは多くの型について定義されているので
    -- オーバーロードに関する問題を避けるために以下の関数が必要
    let texCoord2f = texCoord :: TexCoord2 GLfloat -> IO ()
        vertex3f = vertex :: Vertex3 GLfloat -> IO ()

    -- 表示
    -- 描写する直前に色を指定するようだ
    -- Color4 赤 緑 青 透明度(０で完全に透明、１で完全に不透明）
    currentColor $= Color4 1 1 1 1 

    -- テクスチャの設定
    textureBinding Texture2D $= Just texname

    preservingMatrix $ do
        translate (Vector3 100 100 (0::Double))
        
        -- rotate 度数法での角度
        rotate r (Vector3 0.0 0.0 1.0 :: Vector3 GLdouble)
        
        -- 描画
        renderPrimitive Quads $ do
　　　　　　-- ここから
            texCoord2f (TexCoord2 0 0)
            vertex3f (Vertex3 (-50) (-50) 0.0)
            -- ここまでで一つの頂点の設定（テクスチャに関するもの、場所に関するもの）
            texCoord2f (TexCoord2 0 1)
            vertex3f (Vertex3 (-50) 50 0.0)
            texCoord2f (TexCoord2 1 1)
            vertex3f (Vertex3 50 50 0.0)
            texCoord2f (TexCoord2 1 0)
            vertex3f (Vertex3 50 (-50) 0.0)
    
    --バッファの入れ替え
    swapBuffers

--タイマが呼ばれるたびにactを繰り返す
timerProc act = do
    act
    addTimerCallback timerInterval $ timerProc act

--ウィンドウのサイズが変更された時の処理
reshape (Size w h)=do
    viewport $= (Position 0 0, (Size w h)) --ウィンドウ全体を使う
    
    --ビューボリュームの設定-
    matrixMode $= Projection
    loadIdentity

    -- 視野のセッティング
    ortho 0.0 640.0 0.0 480.0 (-1000.0) 1000.0
    
    -- これ大事（理由はよくわからない）
    matrixMode $= Modelview 0 

--キー入力の処理
keyboardProc arg ch state _ _
    -- Qが押されたら終了 
    | ch     == Char 'q'    = exitWith ExitSuccess
    -- Aが押されたら回転速度を二倍にする
    | ch    == Char 'a'     = modifyIORef arg (*(2))
    -- Sが押されたら回転速度を二分の1にする
    | ch    == Char 's'     = modifyIORef arg (*(0.5))
    -- それ以外なら回転の方向を変える
    | state    == Down        = modifyIORef arg (*(-1))    
    | otherwise            = return ()

HaskellとOpenGLをつかって四角形を描画する

2013-03-04 16:05:20 JST · 元記事

うまく動いたら雛形にでも使ってやってください。



参考サイト
　　　GLUTによる「手抜き」OpenGL入門
　　　プログラミング/Haskell/HSDL - Flightless wing

import Graphics.UI.GLUT
import Graphics.Rendering.OpenGL.GLU
import System.Exit
import Data.IORef

--タイマの間隔
timerInterval = 40

main = do 
    --回転の角度を初期化
    rot <- newIORef 0.0
    arg <- newIORef 10.4
    
    --GLUTの初期化
    initialDisplayMode $= [RGBAMode, DoubleBuffered]
    initialWindowSize $= Size 640 480
    initialize "" []

    --ウィンドウを作る
    createWindow "GLUTpractice"

    --表示に使うコールバック関数の指定
    displayCallback $= display rot arg
    
    --ウィンドウのサイズが変更された時に呼ぶコールバック関数の指定
    reshapeCallback $= Just reshape
    
    --キーボードやマウスのコールバック
    keyboardMouseCallback $= Just (keyboardProc arg)
    
    --タイマを作る
    addTimerCallback timerInterval $ timerProc (display rot arg)

    --GLUTのメインループに入る
    mainLoop

display rot arg= do
    --回転させる
    w<-readIORef arg 
    --w <- readIORef hoge
    modifyIORef rot (+w)
    r <- readIORef rot
 
    --背景を黒にする
    clearColor $= Color4 0.0 0.0 0.0 0.0
    clear [ColorBuffer]
    
    --単位行列を読み込む
    loadIdentity
    
    --表示
    currentColor $= Color4 0 0 1 0  -- 描写する直前に色を指定するようだ
    preservingMatrix $ do
        translate (Vector3 100 100 (0::Double))
        rotate r (Vector3 0.0 0.0 1.0 :: Vector3 GLdouble) -- rotate 度数法での角度
        {-- renderPrimitive Quads $ mapM_ vertex [
                    Vertex3 0.10 0.10 0.0,
                    Vertex3 (-0.10) 0.10 0.0,
            Vertex3 (-0.10) (-0.10) 0.0,
            Vertex3 0.10 (-0.10) 0.0 :: Vertex3 GLfloat] --}       
        renderPrimitive Quads $ mapM_ vertex [
                    Vertex3 50 50 0.0,
                    Vertex3 (-50) 50 0.0,
            Vertex3 (-50) (-50) 0.0,
            Vertex3 50 (-50) 0.0 :: Vertex3 GLfloat]
    
    --バッファの入れ替え
    swapBuffers

--タイマが呼ばれるたびにactを繰り返す
timerProc act = do
    act
    addTimerCallback timerInterval $ timerProc act

--ウィンドウのサイズが変更された時の処理
reshape (Size w h)=do
    viewport $= (Position 0 0, (Size w h)) --ウィンドウ全体を使う
    
    --ビューボリュームの設定-
    matrixMode $= Projection
    loadIdentity
    ortho 0.0 640.0 0.0 480.0 (-1000.0) 1000.0
    -- ortho2D 0.0 640.0 0.0 480.0
    
    matrixMode $= Modelview 0 -- これ大事
    
    {- --少し後ろから撮影
    lookAt (Vertex3 0.0 0.0 (600.0)) (Vertex3 0.0 0.0 0.0) (Vector3 0.0 1.0 0.0)
    matrixMode $= Modelview 0 -}
    

--キー入力の処理
keyboardProc arg ch state _ _
    | ch     == Char 'q'    = exitWith ExitSuccess        --qが押されたら終了
    | ch    == Char 'a'     = modifyIORef arg (*(2))
    | ch    == Char 's'     = modifyIORef arg (*(0.5))
    | state    == Down        = modifyIORef arg (*(-1))    --それ以外なら回転の方向を変える
    | otherwise            = return ()

cabal install を使う前に

2012-10-28 07:12:44 JST · 元記事

apt-cache search パッケージ名
とするとよいかも

AOJ 0090 Overlaps of Seals

2012-10-13 02:35:14 JST · 元記事

cabal install monadius したら結構ハマった

2012-04-03 23:28:20 JST · 元記事

JOI2012春合宿

2012-03-28 10:24:17 JST · 元記事

一日目：１２５
二日目：１１
三日目：１０
四日目：７０

Haskellで文字列圧縮(エンコーダーのみ）

2012-03-03 18:26:00 JST · 元記事

import Data.Map as M
import Data.List
import Data.Char

allChar = Prelude.map chr [0..127]

genDict str = fromList $ Prelude.map func $ zip [1..] $ Prelude.map snd $ reverse $ sort $ zip (Prelude.map (countChar str) allChar) allChar
	where 
		func (a,b) = (b,makeCodeLengthN a)
		makeCodeLengthN 1 = "1"
		makeCodeLengthN n = '0' : makeCodeLengthN (n-1)
		countChar [] _ = 0
		countChar (x:xs) c = if x==c
			then 1 + countChar xs c
			else countChar xs c

translate [] d = return []
translate (x:xs) d = do
	y <- M.lookup x d
	ys <- translate xs d
	return (y++ys)

bitToChar str = func (reverse str) 0
	where
		func [] a = chr (33+a)
		func (x:xs) a = if x=='1'
			then func xs (a*2+1)
			else func xs (a*2)

split6 str = if length str < 6
	then [str]
	else fst (splitAt 6 str) : split6 (snd (splitAt 6 str) )

compress str = do
	a <- translate str $ genDict str
	return ( Prelude.map bitToChar (split6 a) )

Haskellでheadコマンド

2012-03-02 17:33:59 JST · 元記事

Haskellで作るものが思いつかないのでheadコマンドを実装してみた。

main = do
	cs <- getContents
	putStr $ unlines $ take 10 $ lines cs

std::setでデータ同士の間に順位付けができないと無視されてしまう件

2012-02-25 00:54:18 JST · 元記事

#include <stdio.h>
#include <set>
using namespace std;
struct Data{
	int a,b;
};
Data makeD(int aa,int bb){
	Data d;
	d.a=aa,d.b=bb;
	return d;
}
bool operator < (const Data &d,const Data &e){
	return d.a<e.a;
}
int main(){
	set<Data> se;
	se.insert(makeD(0,1));
	se.insert(makeD(0,2));
	printf("se.size()==%d\n",se.size());
	for(set<Data>::iterator it=se.begin();it!=se.end();it++)
		printf("a=%d b=%d\n",(*it).a,(*it).b);
	return 0;
}

を実行すると

se.size()==1
a=0 b=1

となってしまう。
これはsetが(0,1)と(0,2)の違いを分からないからである。

se.size()==2
a=0 b=1
a=0 b=2

このようにするには

bool operator < (const Data &d,const Data &e){
        if(d.a==e.a)           //<=
               return d.b<e.b; //<=
	return d.a<e.a;
}

こんな風に書き換えると良い。

HaskellでBrainfuckの処理系を実装する

2012-02-19 16:34:22 JST · 元記事

-- Brainfuck.hs
import Data.Char

main = do 
    cs <- getContents
    putStrLn $ runBrainfuck cs
runBrainfuck p = display $ runOneStep $ Machine p 0 [0,0..] 0 ""
runOneStep (Machine a b c d e) = if (length a) == b then Machine a b c d e
                                                    else runOneStep $ incrementProgramPointer func
    where func = case (a!!b) of
                    '>' -> incrementPointer (Machine a b c d e)
                    '<' -> decrementPointer (Machine a b c d e)
                    '+' -> incrementMemory (Machine a b c d e)
                    '-' -> decrementMemory (Machine a b c d e)
                    '.' -> putCharToDisplay (Machine a b c d e)
                    ']' -> jumpEnd (Machine a b c d e)
                    '[' -> jumpStart (Machine a b c d e)
                    otherwise -> (Machine a b c d e)
incrementProgramPointer (Machine a b c d e) = Machine a (b+1) c d e
data Machine = Machine { program::String,programPointer::Int,memory::[Int],memoryPointer::Int,display::String   }
    deriving Show 
incrementPointer (Machine a b c d e) = Machine a b c (d+1) e
decrementPointer (Machine a b c d e) = Machine a b c (d-1) e
putCharToDisplay (Machine a b c d e) = Machine a b c d (e++[chr (c!!d)])
incrementMemory (Machine a b c d e) = Machine a b (changeMemory c 0 d (+ 1))d e
decrementMemory (Machine a b c d e) = Machine a b (changeMemory c 0 d (subtract 1))d e
jumpStart (Machine a b c d e) = if (c!!d)==0 then Machine a (func 0 b) c d e
                                             else Machine a b c d e 
    where func depth pos = case (a!!pos) of
                            '[' ->  func (depth+1) (pos+1)
                            ']' ->  if depth==1 then pos
                                                else func (depth-1) (pos+1)
                            otherwise   -> func depth (pos+1)
jumpEnd (Machine a b c d e) = if(c!!d)==0 then Machine a b c d e
                                          else Machine a (func 0 b) c d e
    where func depth pos = case (a!!pos) of
                            '[' ->  if depth==1 then pos
                                                else func (depth-1) (pos-1)
                            ']' -> func (depth+1) (pos-1)
                            otherwise   -> func depth (pos-1)
changeMemory [] i k f = []
changeMemory (x:mem) i k f = if i==k    then f x:mem
                                        else x:changeMemory mem (i+1) k f

情報オリンピックの過去問の入出力をチェックするのに便利だったコマンド

2012-02-18 21:55:17 JST · 元記事

nkf -Lu -overwrite TEST.txt

とやるとファイル中のCRLFが全てLFになって便利だった。

JOI2007 春合宿 1日目 2番

2012-02-18 20:11:29 JST · 元記事

Factorial: 階乗

2012-02-18 06:17:16 JST · 元記事

#include <stdio.h>
#include <vector>
#include <map>
using namespace std;
typedef long long LL;
vector<pair<LL,int> > primeFactor(LL num){ //first->prime,second->pow{
	if(num==1){
		vector<pair<LL,int> > r;
		r.push_back(make_pair(1,1));
		return r;
	}
	else{
		vector<pair<LL,int> > r;
		for(LL i=2;i*i<=num;i++){
			if(num%i==0){
				int p=0;
				while(num%i==0)
					num/=i,p++;
				r.push_back(make_pair(i,p));
			}
		}
		if(num!=1)
			r.push_back(make_pair(num,1));
		return r;
	}
}
int main(){
	LL n;
	scanf("%lld",&n);
	vector<pair<LL,int> > ps=primeFactor(n);
	LL res=0;
	for(int i=0;i<(int)ps.size();i++){
		LL a=0,b=ps[i].second;
		while(0<b){
			a+=ps[i].first;
			for(LL t=ps[i].first;a%t==0;t*=ps[i].first)
				b--;
		}
		res=max(res,a);
	}
	printf("%lld\n",res);
	return 0;
}

AOJ 1203

2012-02-17 23:21:23 JST · 元記事

• 問題文要約

記号を含む文字列が与えられるので、記号を除いて小文字を大文字に変換した後その文字列が含む回文を列挙せよ。ただしcXcとなる回文が存在した場合Xを出力してはならない。

• 解法

dp[i][j]=(i文字目からj文字目までが回文となっているかどうか) でDPする。

#include <string>
#include <iostream>
#include <set>
#include <string.h>
using namespace std;
int dp[1030][1030];
int main(){
	string s;
	while(getline(cin,s)){
		string in;
		for(int i=0;i<(int)s.size();i++)
			if('A'<=s[i]&&s[i]<='Z')	in+=s[i];
			else if('a'<=s[i]&&s[i]<='z')	in+=s[i]-'a'+'A';
		int n=in.size();
		memset(dp,0,sizeof(dp));
		set<string> ans,ng;
		for(int i=0;i<n;i++)
			dp[i][i]=1,dp[i][i+1]=1;
		for(int l=2;l<=n;l++)
			for(int i=0;i+l<=n;i++)
				if(dp[i+1][i+1+l-2]==1&&in[i]==in[i+l-1]){
					dp[i+1][i+l-2]=0;
					ng.insert(in.substr(i+1,l-2));
					dp[i][i+l]=1;
					ans.insert(in.substr(i,l));
				}
		for(set<string>::iterator it=ans.begin();it!=ans.end();it++)
			if(3<=(*it).size()&&ng.find(*it)==ng.end())
				cout<<*it<<" ";
		cout<<endl;
	}
	return 0;
}

AOJ 1204

2012-02-17 03:04:01 JST · 元記事

• 解法

最後のnクロックの状態と次にpipelineに入れられる仕事の番号を使って拡張ダイクストラをする。int dist[21][1<<20]とかやるとMLEするのでsetを使う。

//AOJ1204
#include <stdio.h>
#include <vector>
#include <string>
#include <string.h>
#include <queue>
#include <map>
#include <set>
#define INF (1<<20)
#define W 10
using namespace std;
struct State{
	int time,next,line;
};
State makeS(int t,int nn,int l){
	State s;
	s.time=t,s.next=nn,s.line=l;
	return s;
}
bool operator < (const State &s,const State &t){
	return s.time>t.time;
}
int n;
int used[5];
char as[5][30];
int pi[30];
int check(int line){
	memset(used,0,sizeof(used));
	for(int i=n-1;0<=i;i--)
		if(line&(1<<i))
			used[pi[n-1-i]]++;
	for(int i=0;i<5;i++)
		if(1<used[i])
			return 0;
	return 1;
}
int main(){
	while(1){
		scanf("%d",&n);
		if(n==0)	break;
		for(int i=0;i<5;i++)scanf("%s",as[i]);
		for(int i=0;i<n;i++)for(int j=0;j<5;j++)if(as[j][i]=='X')pi[i]=j;
		int ans;
		priority_queue<State> que;
		set<pair<int,int> > se;
		que.push(makeS(0,0,0));
		se.insert(make_pair(0,0));
		while(!que.empty()){
			State s=que.top();que.pop();
			if(s.next==10&&s.line==1){ans=s.time;break;}
			int l=s.line>>1,ok=check(l);
			if(!ok)	continue;
			if(s.next==W){
				if(se.find(make_pair(W,l))==se.end()){
					se.insert(make_pair(W,l));
					que.push(makeS(s.time+1,W,l));
				}
			}else{
				int ll=l|(1<<(n-1));
				if(check(ll)){
					if(se.find(make_pair(s.next+1,ll))==se.end()){
						se.insert(make_pair(s.next+1,ll));
						que.push(makeS(s.time+1,s.next+1,ll));
					}
				}
				if(se.find(make_pair(s.next,l))==se.end()){
					se.insert(make_pair(s.next,l));
					que.push(makeS(s.time+1,s.next,l));
				}
			}
		}
		printf("%d\n",ans);
	}
}

JOI2012本選5番

2012-02-13 21:51:10 JST · 元記事

ダイクストラして最大全域木求めてLCAする全部入り問題。

#include <vector>
#include <queue>
#include <stdio.h>
#include <algorithm>
#include <string.h>
#include <queue>
#define MAX_V 100000+10
#define MAX_Q 100000+10
#define INF (1<<28)
#define ROOT 1
using namespace std;
struct Edge{
	int to,cost;
};
Edge makeE(int t,int c){
	Edge e;e.to=t,e.cost=c;
	return e;
}
bool operator < (const Edge &e,const Edge &f){
	return e.cost>f.cost;
}
struct Biedge{
	int a,b,cost;
};
Biedge makeB(int aa,int bb,int c){
	Biedge e;e.a=aa,e.b=bb,e.cost=c;
	return e;
}
bool operator < (const Biedge &e,const Biedge &f){
	return e.cost>f.cost;
}
class Unionfind{
	public:
		vector<int> data;
		Unionfind(int n){
			for(int i=0;i<n;i++)
				data.push_back(i);
		}
		int root(int a){
			//printf("root(%d)\n",a);
			if(data[a]==a)	return a;
			else		return data[a]=root(data[a]);
		}
		void set(int a,int b){
			if(root(a)==root(b))
				return;
			else
				data[root(a)]=root(b);
		}
};
int fest[MAX_V],dist[MAX_V],qa[MAX_Q],qb[MAX_Q],qr[MAX_Q],depth[MAX_V],score[32][MAX_V],parent[32][MAX_V];
vector<Biedge> biedge;
vector<Edge> G[MAX_V];
vector<Edge> MG[MAX_V];
void dfs(int v,int p,int s,int d){
	parent[0][v]=p;
	score[0][v]=s;
	depth[v]=d;
	for(int i=0;i<(int)MG[v].size();i++)
		if(MG[v][i].to!=p)
			dfs(MG[v][i].to,v,MG[v][i].cost,d+1);
}
void init(int n){
	dfs(ROOT,-1,dist[ROOT],0);
	for(int k=0;k+1<32;k++)
		for(int i=0;i<n;i++){
			if(parent[k][i]<0)	parent[k+1][i]=-1;
			else{
				score[k+1][i]=min(score[k][parent[k][i]],score[k][i]);
				parent[k+1][i]=parent[k][parent[k][i]];
			}
		}
}
int lca(int u,int v){
	int res=INF;
	if(depth[u]>depth[v])	swap(u,v);
	for(int k=0;k<32;k++)
		if((depth[v]-depth[u])>>k&1){
			res=min(res,score[k][v]);
			v=parent[k][v];
		}
	//printf("lca-m\n");fflush(stdout);
	if(u==v)	return res;
	for(int k=31;0<=k;k--)
		if(parent[k][u]!=parent[k][v]){
			res=min(res,score[k][u]);
			u=parent[k][u];
			res=min(res,score[k][v]);
			v=parent[k][v];
		}
	res=min(res,score[0][u]);
	res=min(res,score[0][v]);
	return res;
}
int main(){
	int n,m,k,q;
	scanf("%d%d%d%d",&n,&m,&k,&q);
	for(int i=0;i<m;i++){
		int a,b,c;
		scanf("%d%d%d",&a,&b,&c);
		G[a].push_back(makeE(b,c));
		G[b].push_back(makeE(a,c));
		biedge.push_back(makeB(a,b,0));
	}
	for(int i=0;i<k;i++)
		scanf("%d",fest+i);
	for(int i=0;i<q;i++)
		scanf("%d%d",qa+i,qb+i);
	for(int i=0;i<MAX_V;i++)
		dist[i]=INF;
	priority_queue<Edge> que;
	for(int i=0;i<k;i++){
		que.push(makeE(fest[i],0));
		dist[fest[i]]=0;
	}
	while(!que.empty()){
		int t=que.top().to,c=que.top().cost;
		que.pop();
		if(dist[t]<c)
			continue;
		for(int i=0;i<(int)G[t].size();i++){
			int tt=G[t][i].to,cc=c+G[t][i].cost;
			if(cc<dist[tt]){
				que.push(makeE(tt,cc));
				dist[tt]=cc;
			}
		}
	}
	for(int i=0;i<m;i++)
		biedge[i].cost=min(dist[biedge[i].a],dist[biedge[i].b]);
	sort(biedge.begin(),biedge.end());
	memset(qr,-1,sizeof(qr));
	Unionfind un(n+1);
	for(int i=0;i<m;i++){
		if(un.root(biedge[i].a)!=un.root(biedge[i].b)){
			un.set(biedge[i].a,biedge[i].b);
			MG[biedge[i].a].push_back(makeE(biedge[i].b,biedge[i].cost));
			MG[biedge[i].b].push_back(makeE(biedge[i].a,biedge[i].cost));
		}
	}
	init(n+1);
	for(int i=0;i<q;i++)
		printf("%d\n",lca(qa[i],qb[i]));
	return 0;
}

JOI2012本選４番

2012-02-13 19:17:54 JST · 元記事

講義資料 - いもす研 (imos laboratory)
↑参照。
累積和を使ったコード

#include <stdio.h>
int nails[5010][5010];
int main(){
	int n,m;
	scanf("%d%d",&n,&m);
	for(int i=0;i<m;i++){
		int a,b,x;
		scanf("%d%d%d",&a,&b,&x);
		nails[a][b]++,nails[a][b+1]--;
		nails[a+x+1][b]--,nails[a+x+2][b+1]++;
		nails[a+x+1][b+x+2]++,nails[a+x+2][b+x+2]--;
	}
	for(int i=0;i<5010;i++){
		int a=0;
		for(int j=0;j<5010;j++){
			a+=nails[i][j];
			nails[i][j]=a;
		}
	}
	for(int i=0;i<5010;i++){
		int a=0;
		for(int j=0;j<5010;j++){
			a+=nails[j][i];
			nails[j][i]=a;
		}
	}
	for(int i=0;i<5010;i++){
		int a=0;
		for(int j=0;i+j<5010;j++){
			a+=nails[i+j][j];
			nails[i+j][j]=a;
		}
	}
	int ans=0;
	for(int i=0;i<5010;i++)
		for(int j=0;j<5010;j++)
			if(0<nails[i][j])
				ans++;
	printf("%d\n",ans);
	return 0;
}

最大値の伝搬を使ったコード（ただし後ろから）。

#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
int nails[5010][5010];
int main(){
	int n,m;
	scanf("%d%d",&n,&m);
	memset(nails,0,sizeof(nails));
	for(int i=0;i<m;i++){
		int a,b,x;
		scanf("%d%d%d",&a,&b,&x);
		nails[a+x][b+x]=max(nails[a+x][b+x],x+1);
	}
	int ans=0;
	for(int y=n;0<=y;y--)
		for(int x=n;0<=x;x--){
			if(0<nails[y][x])
				ans++;
			if(0<x)
				nails[y][x-1]=max(nails[y][x-1],nails[y][x]-1);
			if(0<x&&0<y)
				nails[y-1][x-1]=max(nails[y-1][x-1],nails[y][x]-1);
		}
	printf("%d\n",ans);
	return 0;
}

JOI2012本選3番

2012-02-13 19:15:16 JST · 元記事

時間がSにかぶるときだけ気をつけてdpすればOK。

#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
int n,S,T;
int A[3010],B[3010];
int dp[3010][3010];
int rec(int store,int time){
	if(store==n)
		return 0;
	else if(dp[store][time]!=-1)
		return dp[store][time];
	else{
		int res=rec(store+1,time);
		if(time<S&&S<time+B[store]){
			if(S+B[store]<=T){
				res=max(res,rec(store+1,S+B[store])+A[store]);
			}
		}else{
			if(time+B[store]<=T){
				res=max(res,rec(store+1,time+B[store])+A[store]);
			}
		}
		return dp[store][time]=res;
	}
}
int main(){
	scanf("%d%d%d",&n,&T,&S);
	for(int i=0;i<n;i++)
		scanf("%d%d",A+i,B+i);
	memset(dp,-1,sizeof(dp));
	int ans=rec(0,0);
	printf("%d\n",ans);
	return 0;
}

JOI2012本選2番

2012-02-13 19:10:09 JST · 元記事

dp[ap=Aの何枚目か][bp=Bの何枚目か]でDPした。

#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
int dp[5010][5010];
int ac[5010],bc[5010];
int rec(int ap,int bp){
	if(ap<0||bp<0)
		return 0;
	else if(dp[ap][bp]!=-1)
		return dp[ap][bp];
	else{
		int r=rec(ap-1,bp);
		if(ac[ap]==bc[bp])
			r=max(r,rec(ap-1,bp-1)+1);
		return dp[ap][bp]=r;
	}
}
int main(){
	int an,bn;
	scanf("%d%d",&an,&bn);
	for(int i=0;i<an;i++)
		scanf("%d",ac+i);
	for(int i=0;i<bn;i++)
		scanf("%d",bc+i);
	memset(dp,-1,sizeof(dp));
	int ans=0;
	for(int i=0;i<an;i++)
		for(int j=0;j<bn;j++)
			ans=max(ans,rec(i,j));
	printf("%d\n",ans);
	return 0;
}

JOI2012本選1番

2012-02-13 19:09:52 JST · 元記事

解説を聞いてから自分のやっていたことがRunLength圧縮と同じであると気づいた。

#include <stdio.h>
#include <vector>
using namespace std;
char in[10000000+10];
int main(){
	scanf("%s",in);
	vector<int> ch,sq;
	int c=in[0],s=1;
	for(int i=1;in[i]!=0&&in[i]!='\n';i++){
		if(c!=in[i]){
			ch.push_back(c);
			sq.push_back(s);
			c=in[i];
			s=1;
		}else{
			s++;
		}
	}
	ch.push_back(c);
	sq.push_back(s);
	int ans=0;
	for(int i=0;i+2<(int)ch.size();i++){
		if(ch[i]=='J'&&ch[i+1]=='O'&&ch[i+2]=='I'&&sq[i+1]<=min(sq[i],sq[i+2]))
			ans=max(ans,min(sq[i],sq[i+1]));
	}
	printf("%d\n",ans);
	return 0;
}

JOI2012本選まとめ

2012-02-13 19:03:31 JST · 元記事

↑の大会に参加してきました。
実力を限界まで発揮できたと思います。
合宿は問題を楽しむだけにしようと思います。

JOI2011本選5番

2011-08-30 00:37:20 JST · 元記事

kについて二分探索した後、その判定を頑張ってO(n)で出来るようにする。

#include <stdio.h>
#include <utility>
#include <set>
#include <algorithm>
#include <vector>
#define MAX_V 1000100

using namespace std;

int n;
int ai[MAX_V],can[MAX_V],bi2ai[MAX_V];
pair<int,int> bi[MAX_V];

int ok(int k)
{
	if(k==0)
		return 1;
	double sumA=0;
	for(int i=0;i<n;i++)
		can[i]=1;
	for(int i=0;i<k;i++)
		sumA+=ai[i];
	int lastA=k-1;
	for(int b=0;b+k<=n;b++)
	{
		if(sumA/(double)k<=(double)bi[b].first)
			return 1;
		int tP=bi2ai[b];
		can[tP]=0;
		if(tP<=lastA)
		{
			sumA-=ai[tP];
			lastA++;
			if(lastA>=n)
				return 0;
			while(!can[lastA])
				lastA++;
			sumA+=ai[lastA];
		}
	}
	return 0;
}

int main()
{
	scanf("%d",&n);
	for(int i=0;i<n;i++)
	{
		int a,b;
		scanf("%d%d",&a,&b);
		ai[i]=a;
		can[i]=1;
		bi2ai[i]=-1;
		bi[i]=make_pair(b,-a);
	}
	sort(ai,ai+n);
	sort(bi,bi+n);
	for(int i=0;i<n;i++)
	{
		int a=-bi[i].second;
		int tB=lower_bound(ai,ai+n,a)-ai;
		int tE=upper_bound(ai,ai+n,a)-ai;
		int tP;
		for(tP=tB;tP<tE;tP++)
			if(can[tP])
				break;
		can[tP]=0;
		bi2ai[i]=tP;
	}
	int b=0,c=n;
	while(1<c-b)
	{
		int m=(b+c)/2;
		if(ok(m))
			b=m;
		else
			c=m;
	}
	if(ok(c))
		printf("%d\n",c);
	else
		printf("%d\n",b);
	return 0;
}

SuperCon2011に参加しました

2011-08-28 17:06:35 JST · 元記事

２０１１年８月２２日(月) 〜 ２６日(金)にかけて大阪大学豊中キャンパスで
SuperCon2011にチームMKとして参加してきました。
　
今回の課題は「https://www.cp.cmc.osaka-u.ac.jp/supercon/naklon.html」。
リンクをクリックして表示されるのはSIZE＝８のものですが
実際の本選はSIZE＝２５０という巨大なくろんなので全探索はできません。
　
それではどうするか。
僕達のチーム（主に僕）は２つの方針を考えました。
１：１０００手ほど乱択しそれを何度か繰り替えして最もよいものを選ぶことを何度か繰り返す
２：５手ほど全探索して評価関数で評価して最も良い手を選ぶ。
　
１は相方に実装してもらったのですがあまり結果が芳しくなかったので
僕は２を頑張っていました。
　
３日目の昼頃に評価関数を得点÷落ちた玉に変えた所
非常に結果が良かったので最後までその方針で行きました
　
結果は５位でした。
２位から５位まではほとんど同じ方針だったということなので悔しいです

Haskellで数式をパースする

2011-08-15 12:00:27 JST · 元記事

使い方
　:load で読み込んだ後
　myEquation "1+2*3"

import Text.ParserCombinators.Parsec

myDigit :: Parser Char
myDigit = oneOf ['0'..'9']

myNumberString :: Parser String
myNumberString = do
	c<-myDigit
	do
			cs <- myNumberString
			return (c:cs)
		<|> return [c]

myNumber :: Parser Integer
myNumber = do
	s<-myNumberString
	return (read s)

myPlusMinus :: Parser Char
myPlusMinus = oneOf "+-"

myMultiDiv :: Parser Char
myMultiDiv = oneOf "*/"

myCalc '+' left right = (left+right)
myCalc '-' left right = (left-right)
myCalc '*' left right = (left*right)
myCalc '/' _ 0 = 0
myCalc '/' left right = (left `div` right)

myEquation :: Parser Integer
myEquation = do{ left<-myFactor; rest left }
		where
			rest l = do{ op<-myPlusMinus
					;r<-myFactor
					;rest (myCalc op l r)
					}
				<|> return l
{-<|> do{ op<-myPlusMinus
		;right<-myEquation
		;return (myCalc op x right
		}-}

myFactor :: Parser Integer
--myFactor = do{ x<-myTerm; rest x}
myFactor = do {left<-myTerm; rest left }
		where
			rest l = do{ op<-myMultiDiv
					;r<-myTerm
					;rest (myCalc op l r)
					}
				<|> return l
				{-
		left<-myFactor
		op<-myMultiDiv
		right<-myTerm
		return (myCalc op left right)
	<|>	myTerm-}

myTerm :: Parser Integer
myTerm = do
		char '('
		x<-myEquation
		char ')'
		return x
	<|> myNumber

{-

myExpr2List String->[String]
myExpr2List "" = ""
myExpr2List x:xs = if myPlusMinus 

myFactor :: Parser Integer
myFactor =
	do
		char '('
		x<-myExpr
		char ')'
		return x
	<|>	myExpr

myExpr :: Parser Integer
myExpr = do
		left<-myNumber
		do
				op<-myPlusMinus
				right<-myFactor
				return (myCalc op left right)
			<|> return left
	<|> myFactor
-}

HaskellでHTMLをパースする

2011-08-15 11:25:18 JST · 元記事

使い方
　hoge.hsに下のコードを保存する。
　runghc "hoge.hs" "target.html" ※ダブルクオーテーションは必要

　このままでは空行が大量に含まれているので
　runghc "hoge.hs" "target.html" | sed '/^ *$/d'
　すると良い

--sed '/^ *$/d'
import Text.ParserCombinators.Parsec
import System.Environment

data MyTag = MakeMyTag {tagName::String,tagContent::[String]} deriving (Eq,Show)

testMyTag = MakeMyTag "aa" ["bb","cc"]

data MyHTMLContent = MakeContentChild { myLeftTag::MyTag,myContentMain::[MyHTMLContent],myRightTag::MyTag }
	| MakeTagSelf MyTag
	| MakePlainText { myText::String }
	| MakeComment
	deriving (Eq,Show)

testHTMLContent1 = (MakePlainText "dd")
testHTMLContent2 = (MakeTagSelf testMyTag)
testHTMLContent3 = (MakeContentChild testMyTag [testHTMLContent1,testHTMLContent2] testMyTag)

outputHTML (MakeComment) n = ""
outputHTML (MakePlainText "") n = ""
outputHTML (MakePlainText s) n = (sTab n++s++"\n")
outputHTML (MakeTagSelf x) n = ""
outputHTML (MakeContentChild l con r) n = sCon con (n+1)
		
sTab 0 = ""
sTab n = "  "++sTab (n-1)

sCon [] _ = "" 
sCon (x:xs) n = (outputHTML x n++sCon xs n)

main = do
	inStr<-getArgs
	tmp <- parseFromFile myHTML (head inStr)
	case tmp of
		Left err -> do {print "error"; }
		Right xs -> do {putStrLn (func xs); }

func :: [MyHTMLContent]->String
func [] = ""
func (x:xs) = (outputHTML x 0)++(func xs)

isNotEqual c = if c=='='
	then False
	else True

myDigit :: Parser Char
myDigit = oneOf ['0'..'9']

myNumberString :: Parser String
myNumberString = do
	c<-myDigit
	do
			cs <- myNumberString
			return (c:cs)
		<|> return [c]

myTagContentLeft :: Parser String
myTagContentLeft = many1 (noneOf "=>/")
-- myTagContentLeft = many1 (satisfy (\c -> c /= '='))
--myTagContentLeft = many1 (alphaNum <|> char ':')
{-where
		rest = do{ x<-satisfy (\c -> c /= '=')
				; xs<-rest
				; return (x:xs) 
				}
			<|> return ""
-}

myTagContentRight :: Parser String
myTagContentRight = do
		lD<-char '"'
		mD<-inC
		rD<-char '"'
		return ([lD]++mD++[rD])
	<|> myNumberString
		where

			inC = do{ x<-satisfy (\c -> c /= '"')
					; xs<-inC
					; return (x:xs)
					}
				<|> return ""

myTagContent :: Parser String
myTagContent = do
	left<-myTagContentLeft
	char '='
	right<-myTagContentRight
	return (left++"="++right)

myTagName :: Parser String
myTagName = many1 alphaNum

myEndTag :: Parser MyTag
myEndTag = do
	string "</"
	many space
	name<-myTagName
	many space
	char '>'
	return (MakeMyTag name [])

myBeginTag :: Parser MyTag
myBeginTag = do
	char '<'
	name<-myTagName
	many space
	tagContent<-many p
	char '>'
	return (MakeMyTag name tagContent)
		where
			p = do
				r<-myTagContent
				many space
				return r

mySelfTag :: Parser MyTag
mySelfTag = do
	char '<'
	name<-myTagName
	many space
	tagContent<-many p
	string "/>"
	return (MakeMyTag name tagContent)
		where
			p = do
				r<-myTagContent
				many space
				return r

myComment :: Parser MyHTMLContent
myComment = do
		string "<!--"
		rest "aa"
		char '>'
		return MakeComment
	-- <?> "myComment"
	
	where 
		rest before =
			if before == "--"
			then return ()
			else do
				x<-anyChar
				rest (tail before++[x])
				
myContent :: Parser MyHTMLContent
myContent = do {p<-many1 (satisfy (\c-> c/='<')); return ( MakePlainText p)}
	<|> myBeginAndEndTag
		<|> do  {t<-try mySelfTag; return (MakeTagSelf t)}
		<|> try myComment
		<?> "myContent"

{--myContent = try( do {b<-myBeginTag; l<-myContentList; e<-myEndTag; return  (MakeContentChild b l e)})
		<|> try( do  {t<-mySelfTag; return (MakeTagSelf t)})
		<|>  do {p<-many1 (satisfy (\c-> c/='<')); return ( MakePlainText p)}--}

myBeginAndEndTag :: Parser MyHTMLContent
myBeginAndEndTag = do
	b<-try myBeginTag
	l<-myContentList
	e<-myEndTag
	if tagName b == tagName e 
		then return (MakeContentChild b l e)
		else fail $ "tag don't match "++tagName b++" and "++tagName e
	
myContentList :: Parser [MyHTMLContent]
myContentList = (do	{ r<-myContent; l<-myContentList; return ([r]++l)})
			<|> do return []
--myContentList = try (do	{r<-myContent; many myComment; l<-myContentList; return ([r]++l)})

myHTML :: Parser [MyHTMLContent]
myHTML = do
	x <- myContentList
	eof
	return x

ubuntu11.04でmp600を使う

2011-07-31 09:34:48 JST · 元記事

環境：xubuntu11.04+canon mp600

キヤノン：ダウンロード｜ScanGear MP Ver.1.00 for Linux
ここから共通ファイルとmp600用のファイルをrpm形式でダウンロードして
キヤノンよりダウンロードした２つのファイルを~/canonに置いた。
　　
sudo apt-get install alien
　　　
sudo apt-get install libglib2.0
　　　
sudo apt-get install libgtk2.0
　　　
sudo apt-get install libpng3
　　　
cd canon
　　　
sudo alien -c *.rpm
　　　
sudo dpkg -i scangearmp-common_1.00-3_i386.deb
　　　
sudo dpkg -i scangearmp-mp600_1.00-2_i386.deb
　　　
sudo ln -s /usr/lib/i386-linux-gnu/libtiff.so.4 /usr/lib/libtiff.so.3
　　　
sudo /etc/init.d/cups restart

プリンターの電源を入れパソコンに接続した
　　
動いた
　　　
参考にしたサイト
Ubuntu日本語フォーラム / プリンター動作報告
の#116、Unpaさんの投稿

やり方が悪いのかもしれませんが、ここで
/etc/cups/printers.confを強引に書き換えないと、プリンタが動きませんでした。
　　
sudo vi /etc/cups/printers.conf
「DeviceURI usb://Canon/MP600」の行の内容を
「DeviceURI cnij_usb:/dev/usb/lp0」に書き換え

ここは必要なかった

JavaScriptでインベーダーゲームを作るーその１

2011-07-12 21:58:35 JST · 元記事

JavaScriptの実行されるタイミングがよく分からない。

SuperCon2011予選

2011-06-22 22:10:59 JST · 元記事

PriorityQueueを使って全ての地点までのk最短路を求めた後にメモ化再帰。
ソースコードの後半は指定されていたテンプレート＋入力データのコピーだけ。

/* SuperCon 2011 予選問題Ｃ用テンプレート
   ・解答プログラムはこのテンプレートに従って作成すること．
　 ・解答プログラムは１つのファイルで，チーム名.c という名前にすること．
　 ・入力の方式は，あらかじめ入力ファイル（例：input_sample.txt）を作っ
　　 ておき，実行時にファイル名を指定する方式です．
*/

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

/* ↓以下の範囲は変更可能 */

#include <limits.h>

#define HEAPTABLESIZE 4*210*210*210

int heapData[HEAPTABLESIZE];
int heapTable[HEAPTABLESIZE][3];
int heapSize=0;
int heapTableSize=HEAPTABLESIZE;
int heapTableCur=0;

void heapPop(int*,int*,int*);
void heapPush(int,int,int);
int heapEmpty();

void heapDown(int);
void heapUp(int);

int rec(int x,int y,int o);
void solve(int _h,int _w,int _k,int *ans_len,int *ans_num);

int k,w,h;
int mv[4][2]={{0,1},{1,0},{0,-1},{-1,0}};
int dist[210][210][210];
int filled[210][210];
int MAP[200+1][200+1][4];
int cache[210][210][210];

int heapEmpty()
{
	return (heapSize==0);
}

void heapPop(int *c,int *x,int *y)
{
	if(heapSize==0)
	{
		*c=-1,*x=-1,*y=-1;
	}
	else if(heapSize==1)
	{
		int r=heapData[1];
		heapData[heapSize]=-1;
		heapSize=0;
		*c=heapTable[r][0],*x=heapTable[r][1],*y=heapTable[r][2];
	}
	else
	{
		int r=heapData[1];
		heapData[1]=heapData[heapSize];
		heapSize--;
		heapDown(1);
		*c=heapTable[r][0],*x=heapTable[r][1],*y=heapTable[r][2];
	}
}

void heapPush(int c,int x,int y)
{
	heapTableCur=(heapTableCur+1)%heapTableSize;
	heapTable[heapTableCur][0]=c;
	heapTable[heapTableCur][1]=x;
	heapTable[heapTableCur][2]=y;

	heapSize++;
	heapData[heapSize]=heapTableCur;

	heapUp(heapSize);
}

void heapUp(int k)
{
	int p=k/2;
	if(k==1)
		;
	else if(heapTable[heapData[k]][0]<heapTable[heapData[p]][0])
	{
		int t=heapData[k];
		heapData[k]=heapData[p];
		heapData[p]=t;
		heapUp(p);
	}
	else
		;
}

void heapDown(int k)
{
	int n1=2*k,n2=2*k+1;
	if(heapSize<n1)
	{
		;
	}
	else if(heapSize<n2)
	{
		if(heapTable[heapData[n1]][0]<heapTable[heapData[k]][0])
		{
			int t=heapData[k];
			heapData[k]=heapData[n1];
			heapData[n1]=t;
			heapDown(n1);
		}
		else
		{
			;
		}
	}
	else
	{
		if(heapTable[heapData[n1]][0]<heapTable[heapData[n2]][0] && 
				heapTable[heapData[n1]][0]<heapTable[heapData[k]][0])
		{
			int t=heapData[k];
			heapData[k]=heapData[n1];
			heapData[n1]=t;
			heapDown(n1);
		}
		else if(heapTable[heapData[n2]][0]<=heapTable[heapData[n1]][0] && 
				heapTable[heapData[n2]][0]<heapTable[heapData[k]][0])
		{
			int t=heapData[k];
			heapData[k]=heapData[n2];
			heapData[n2]=t;
			heapDown(n2);
		}
		else
		{
			;
		}
	}	
}

int rec(int x,int y,int o)
{
	int pp;

	int *r=&cache[x][y][o];
	if(*r!=-1)
		;
	else if(x==0 && y==0 && o==0)
		*r=1;
	else
	{
		*r=0;
		for(pp=0;pp<4;pp++)
			if(MAP[x][y][pp]!=0)
			{
				int xx=x+mv[pp][0],yy=y+mv[pp][1];
				int cc=dist[x][y][o]-MAP[x][y][pp];
				int bo=0,ce=k-1;
				while(1<ce-bo)
				{
					int me=(bo+ce)/2;
					if(dist[xx][yy][me]<cc)
						bo=me;
					else
						ce=me;
				}
				if(dist[xx][yy][bo]==cc)
					*r+=rec(xx,yy,bo);
				else if(dist[xx][yy][ce]==cc)
					*r+=rec(xx,yy,ce);
				else
					;
			}
	}
	return *r;
}

void solve(int _h,int _w,int _k,int *ans_len,int *ans_num)
{
	int pp,qq,rr;
	int c,x,y;
	int cc,xx,yy;

	k=_k,w=_w,h=_h;
	h++,w++;

	for(pp=0;pp<210;pp++)
		for(qq=0;qq<210;qq++)
		{
			for(rr=0;rr<210;rr++)
				dist[pp][qq][rr]=INT_MAX;
			filled[pp][qq]=0;
		}

	heapPush(0,0,0);

	while(!heapEmpty())
	{
		heapPop(&c,&x,&y);

		if(filled[x][y]==k ||
			 	(0<filled[x][y] && dist[x][y][filled[x][y]-1]==c))
		{
			;
		}
		else
		{
			dist[x][y][filled[x][y]]=c;
			filled[x][y]++;

			//if(x!=w-1 || y!=h-1)
				for(pp=0;pp<4;pp++)
					if(MAP[x][y][pp]!=0)
					{
						cc=c+MAP[x][y][pp];
						xx=x+mv[pp][0];
						yy=y+mv[pp][1];

						//if(!(filled[xx][yy]==k || (0<filled[xx][yy] && dist[xx][yy][filled[xx][yy]-1]==cc)))
						//{
							heapPush(cc,xx,yy);
						//}
					}
		}
	}

	memset(cache,-1,sizeof(cache));
	if(dist[w-1][h-1][k-1]!=INT_MAX)
	{
		*ans_len=dist[w-1][h-1][k-1];
		*ans_num=rec(w-1,h-1,k-1);
	}
	else
	{
		*ans_len=0;
		*ans_num=0;
	}
}

/* ↑上記の範囲は変更可能 */

int main(int argc, char** argv) {
  int answer_length = -1; /* この変数に k 番目に長い経路の長さを代入してください */
  int answer_number = -1; /* この変数に k 番目に長い経路の総数を代入してください */
  int m, n, k;
  int D[200+1][200+1][4];
  char* problem_file;
  clock_t start, end;
  FILE* fp;

  int i, x, y;
  char buf[0xffff];
  char* p;

  if (argc <= 1) {
    fprintf(stderr, "Enter the input file.\n");
    exit(EXIT_FAILURE);
  }

  problem_file = argv[1];
  fp = fopen(problem_file, "r");
  if (fp == NULL) {
    fprintf(stderr, "Cannot open %s.\n", problem_file);
    exit(EXIT_FAILURE);
  }

  p = fgets(buf, 0xffff, fp);
  assert(p != 0);

  m = atoi(strtok(buf, ", \n"));
  n = atoi(strtok(NULL, ", \n"));
  k = atoi(strtok(NULL, ", \n"));
  assert(m > 0 && m <= 200);
  assert(n > 0 && n <= 200);
  assert(k > 0 && k <= 200);
  for (y = 0; y <= n; y++) {
    p = fgets(buf, 0xffff, fp);
    assert(p != 0);
    p = strtok(buf, ", \n");
    for (i = 0; i < m; i++) {
      int length = atoi(p);
      assert(length >= 0 && length <= 20);
      D[i][y][1] = length;
      D[i+1][y][3] = length;
      p = strtok(NULL, ", \n");
    }
    D[0][y][3] = 0;
    D[m][y][1] = 0;
  }
  for (x = 0; x <= m; x++) {
    p = fgets(buf, 0xffff, fp);
    assert(p != 0);
    p = strtok(buf, ", \n");
    for (i = 0; i < n; i++) {
      int length = atoi(p);
      assert(length >= 0 && length <= 20);
      D[x][i][0] = length;
      D[x][i+1][2] = length;
      p = strtok(NULL, ", \n");
    }
    D[x][0][2] = 0;
    D[x][n][0] = 0;
  }

  fclose(fp);

  /* 入力した情報を画面に出力する（デバッグ等に使って下さい）
　　 提出時は削除しないで，このようにコメントアウトすること
  printf("The input graph\n");
  for (i = 2*n; i >= 0; i--) {
    y = i / 2;
    if (i % 2 == 0) {
      printf("+");
      for (x = 0; x < m; x++) {
        assert(D[x][y][1] == D[x+1][y][3]);
        printf("%d+", D[x][y][1]);
      }
      assert(D[0][y][3] == 0);
      assert(D[m][y][1] == 0);
    } else if (i > 0) {
      for (x = 0; x <= m; x++) {
        assert(D[x][y][0] == D[x][y+1][2]);
        assert(y != n-1 || D[x][y+1][0] == 0);
        printf("%d ", D[x][y][0]);
      }
    } else {
      for (x = 0; x <= m; x++) {
        assert(D[x][y][2] == 0);
      }
    }
    printf("\n");
  }
  printf("\n");
  */

  /* 時間計測用（デバッグ等に使って下さい）
　　 提出時は削除しないで，このようにコメントアウトすること
  start = clock();
  */

/* ↓以下の範囲は変更可能 */
	int pp,qq,rr;
	for(pp=0;pp<201;pp++)
		for(qq=0;qq<201;qq++)
			for(rr=0;rr<4;rr++)
				MAP[pp][qq][rr]=D[pp][qq][rr];
	solve(n,m,k,&answer_length,&answer_number);

/* ↑上記の範囲は変更可能 */


  /* 時間計測用（デバッグ等に使って下さい）
　　 提出時は削除しないで，このようにコメントアウトすること
  end = clock();
  printf("%s, %f, %d, %d\n", problem_file, (double)(end - start) / CLOCKS_PER_SEC, answer_length, answer_number);
  */

  printf("%s, %d, %d\n", problem_file, answer_length, answer_number);
  return 0;
}

SRM370 div1medium

2011-06-21 21:04:00 JST · 元記事

通信範囲を１づつ大きくしながらメモ化再帰。

#include <algorithm>
#include <vector>

using namespace std;

class ConnectTheCities 
{
	public:

		int n,dist,range;
		vector<int> pos;
		int cache[105][55];

		int rec(int bP,int k)
		{
			//printf("bP=%d k=%d\n",bP,k);

			int &r=cache[bP][k];
			if(r!=-1)
				;
			else if(k==n)
			{
				if(dist-range<=bP)
					r=0;
				else
					r=1000000;
			}
			else
			{
				r=1000000;
				for(int nP=bP;nP<=min(bP+range,dist);nP++)
					r=min(r,abs(nP-pos[k])+rec(nP,k+1));
			}
			return r;
		}		

		int minimalRange(int distance, int funds, vector <int> position) 
		{
			dist=distance;
			pos=position;
			sort(pos.begin(),pos.end());
			n=pos.size();
			pos.push_back(distance);

			//printf("n=%d\n",n);

			for(range=1;range<=distance;range++)
			{
				memset(cache,-1,sizeof(cache));
				if(rec(0,0)<=funds)
					return range;
			}
		}
};

無題

2011-06-20 00:00:00 JST · 元記事

ややこしすぎて、なんで通ったのか分からない。

#include <algorithm>
#include <string>
#include <vector>
#include <iostream>
#include <sstream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <set>
#include <map>
#include <climits>
#include <cfloat>

using namespace std;
int K,gN;
int E[40];
int G[16][16];
double C[40][40];
double cache[16][1<<16];

void init()
{
	for(int i=0;i<40;i++)
		for(int j=0;j<40;j++)
			if(i<j)
				C[i][j]=0;
			else if(j==0)
				C[i][j]=1;
			else
				C[i][j]=C[i-1][j]+C[i-1][j-1];
	for(int i=0;i<16;i++)
		for(int j=0;j<1<<16;j++)
			cache[i][j]=-2.0;
}

void makeG(vector<string> input)
{
	for(int i=0;i<16;i++)
		for(int j=0;j<16;j++)
			G[i][j]=0;
	gN=1;
	set<int> masaoFriends;
	for(int i=0;i<input.size();i++)
	{
		int a;
		stringstream ss(input[i]);
		vector<int> v;
		while(ss>>a)
			v.push_back(a);
		if(i==0)
		{
			masaoFriends.insert(1);
			for(int j=0;j<v.size();j++)
				masaoFriends.insert(v[j]);
			for(int j=0;j<v.size();j++)
			{
				int p=distance(masaoFriends.begin(),masaoFriends.find(i+1));
				int q=distance(masaoFriends.begin(),masaoFriends.find(v[j]));
				G[p][q]=G[q][p]=1;
			}
		}
		else if(masaoFriends.find(i+1)!=masaoFriends.end())
		{
			gN++;
			for(int j=0;j<v.size();j++)
				if(masaoFriends.find(v[j])!=masaoFriends.end())
				{
					int p=distance(masaoFriends.begin(),masaoFriends.find(i+1));
					int q=distance(masaoFriends.begin(),masaoFriends.find(v[j]));
					G[p][q]=G[p][q]=1;
				}
		}
		E[distance(masaoFriends.begin(),masaoFriends.find(i+1))]=v.size();
	}

	/*printf("masaoFriends =");
	for(set<int>::iterator it=masaoFriends.begin();it!=masaoFriends.end();it++)
		printf("%d ",*it);
	printf("\n");
	for(int i=0;i<gN;i++)
	{
		for(int j=0;j<gN;j++)
			printf("%d ",G[i][j]);
		printf("\n");
	}*/
}

double rec(int pos,int left)
{
	//printf("rec(%d,%d)\n",pos,left);
	double &r=cache[pos][left];
	if(-1.0<r)
		;
	else if(left==0)
		r=1.0;
	else
	{
		r=0.0;
		vector<double> prob;
		for(int i=0;i<gN;i++)
			if(G[pos][i] && ((1<<i)&left))
				prob.push_back(rec(i,left^(1<<i)));
		sort(prob.begin(),prob.end());
		reverse(prob.begin(),prob.end());
		/*printf("prob=");
		for(int i=0;i<prob.size();i++)
			printf("%lf ",prob[i]);
		printf("\n");*/
		if(K<E[pos])
		{
			double p=1.0;
			for(int i=0;i<prob.size();i++)
			{
				r+=(p-C[E[pos]-i-1][K]/C[E[pos]][K])*prob[i];
				p=C[E[pos]-i-1][K]/C[E[pos]][K];
			}
		}
		else if(0<prob.size())
			r+=prob[0];
	}
	return r;
}

class FriendTour 
{
	public:
		double tourProbability(vector<string> input,int _K) 
		{
			K=_K;
			makeG(input);
			init();

			return rec(0,(1<<gN)-1-1);
		}

};

SRM331 div1 medium

2011-06-19 11:05:03 JST · 元記事

2^10は使わなかった。

#include <algorithm>
#include <vector>
using namespace std;

class Shopping 
{
	public:
	int minNumber(int X, vector <int> values) 
	{
		sort(values.begin(),values.end());
		if(values[0]!=1)
		{
			return -1;
		}
		else
		{
			int r=1;
			int s;
			for(s=1;s<X;)
			{
				int n=0;
				while(n+1<values.size() && values[n+1]<=s+1)
					n++;
				s+=values[n];
				r++;
			}
			return r;
		}
	}
};

SRM337 div1 medium

2011-06-17 21:06:29 JST · 元記事

ALGORITHM NOTE ヒストグラム中の最大の長方形の面積
または
Amazon.co.jp： プログラミングコンテストチャレンジブック: 秋葉 拓哉, 岩田 陽一, 北川 宜稔: 本
の280ページ、スタックの利用を参照。

#include <vector>
#include <stack>

using namespace std;
typedef long long LL;

class BuildingAdvertise 
{
	public:

		int n;
		vector<int> building;
		vector<int> L,R;

		void makeR(vector<int> h)
		{
			building=vector<int>(n,0);
			int j=0;
			int m=h.size();
			for(int i=0;i<n;i++)
			{
				building[i]=h[j];
				int s=(j+1)%m;
				h[j]=((h[j]^h[s])+13)%835454957;
				j=s;
			}
		}

		LL getMaxArea(vector <int> h, int _n) 
		{
			n=_n;
			makeR(h);

			L=vector<int>(n,0);
			R=vector<int>(n,0);

			stack<int> sta;
			for(int i=0;i<n;i++)
			{
				while(!sta.empty() && building[i]<=building[sta.top()])
					sta.pop();
				if(sta.empty())
					L[i]=0;
				else
					L[i]=sta.top()+1;
				sta.push(i);
			}
			while(!sta.empty())
				sta.pop();
			for(int i=n-1;0<=i;i--)
			{
				while(!sta.empty() && building[i]<=building[sta.top()])
					sta.pop();
				if(sta.empty())
					R[i]=n;
				else
					R[i]=sta.top();
				sta.push(i);
			}
			LL ans=0;
			for(int i=0;i<n;i++)
				ans=max(ans,(LL)(R[i]-L[i])*(LL)building[i]);
			return ans;

		}
};

SRM367 div1medium

2011-06-16 20:38:58 JST · 元記事

問題概要
あるデバイスに送らななければならないデータがその始まりのアドレスと大きさの形で渡される。
データはパケットに分けて送らなければならない。
そのパケットはメインデータとオーバーヘッド部分から成る。
このとき送らなければならない最小のバイト数を求めよ。


要はどこからどこまでのデータをまとめて送るかということが大切。
http://www.topcoder.com/tc?module=Static&d1;=match_editorials&d2;=srm367
↑の中にあるサンプルソースコードが非常にわかり易い。

#include <vector>

using namespace std;
typedef long long LL;

class DeviceProgramming 
{
	public:

		LL neededBytes(LL maxPacketSize,LL overhead,LL size)
		{
			LL maxInfo=maxPacketSize-overhead;
			LL packetsNeeded=((size%maxInfo==0)?size/maxInfo:size/maxInfo+1);
			return size+packetsNeeded*overhead;
		}

		LL minBytes(vector <int> offset, vector <int> size, int maxPacketSize, int overhead) 
		{
			int n=offset.size();
			LL DP[60];
			vector<int> st,fn;
			for(int i=0;i<n;i++)
			{
				st.push_back(offset[i]);
				fn.push_back(offset[i]+size[i]-1);
			}
			for(int i=0;i<n;i++)
				for(int j=i+1;j<n;j++)
					if(st[i]>st[j])
					{
						swap(st[i],st[j]);
						swap(fn[i],fn[j]);
					}
			for(int i=0;i<n;i++)
			{
				DP[i]=neededBytes(maxPacketSize,overhead,fn[i]-st[0]+1);
				for(int j=0;j<i;j++)
					DP[i]=min(DP[i],DP[j]+neededBytes(maxPacketSize,overhead,fn[i]-st[j+1]+1));
			}
			return DP[n-1];
		}
};

SRM384 div1 medium

2011-06-12 20:58:18 JST · 元記事

問題要約
ーーーーーーーーーーーー
n人の人が円形に並んでいて、
時計回りの順にボールを任意の他の人に向かって投げる。
ボールがあたった人は円から離脱し、
番号iの人が投げたボールが当たる確率はprobably[i]である。
この時円形に並んでいる人が１人になるまでのボールが投げられる回数の最小の期待値を求めよ。
ーーーーーーーーーーーー
まず、制限の
"Each element of probability will be between 10 and 100, inclusive."
に解法のヒントを見つけた。
10%からということは10000回くらいボールを投げ続けたら
失敗する確率は最悪でも　0.9^10000≒0　に収束する。
それならこの問題は rec(誰が残っている(2^6),誰の番(6),何回目(10000))でメモ化再帰。
つまずいたのはボールを投げるときに誰に向かって投げるかということ。
様々なケースが考えられるので投げられる全員に投げてみて一番期待値が小さいものを選ぶと良い。

#include <algorithm>
#include <string>
#include <vector>
#include <iostream>
#include <sstream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <set>
#include <map>
#include <climits>
#include <cfloat>

using namespace std;
double cache[(1<<6)+1][6+1][10000+10];

class SchoolTrip 
{
	public:

		int n;
		vector<double> p;

		int popCount(int _mask)
		{
			int r=0;
			for(int i=0;i<28;i++)
				if(_mask & (1<<i))
					r++;
			return r;
		}

		double rec(int left,int turn,int time)
		{
			double &r=cache[left][turn][time];
			if(-1.0<r)
				;
			else if(time==10000 || popCount(left)==1)
				r=0.0;
			else
			{
				r=1.0;
				int ln=popCount(left)-1;
				double minThrow=FLT_MAX;
				for(int hit=turn+1;;hit++)
				{
					hit%=n;
					if(hit==turn)
						break;
					int sleft=left^(1<<hit);
					for(int nTurn=turn+1;;nTurn++)
					{
						nTurn%=n;
						if(sleft&(1<<nTurn))
						{
							minThrow=min(minThrow,p[turn]*rec(sleft,nTurn,time+1));
							break;
						}
					}
				}
				r+=minThrow;
				for(int nTurn=turn+1;;nTurn++)
				{
					nTurn%=n;
					if(left & (1<<nTurn))
					{
						r+=(1.0-p[turn])*rec(left,nTurn,time+1);
						break;
					}
				}
			}
			return r;
		}

		double duration(vector <int> probability) 
		{
			n=probability.size();
			p.clear();
			for(int i=0;i<n;i++)
				p.push_back(probability[i]/100.0);
			for(int i=0;i<(1<<6);i++)
				for(int j=0;j<6;j++)
					for(int k=0;k<10000;k++)
						cache[i][j][k]=-2.0;
			return rec((1<<n)-1,0,0);
		}
};

SRM397 div1 medium

2011-06-12 17:54:26 JST · 元記事

k乗和の公式を使う。
求め方は数学Bの教科書に載っている。

#include <algorithm>
#include <string>
#include <vector>
#include <iostream>
#include <sstream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <set>
#include <map>
#include <climits>

using namespace std;
typedef long long LL;

static const LL MOD=1000000007;
LL ADD(LL x,LL y){	return (x+y)%MOD;	}
LL SUB(LL x,LL y){	return (x-y+MOD)%MOD;	}
LL MUL(LL x,LL y){	return ((x*y)%MOD+MOD)%MOD;	}
LL POW(LL x,LL e)
{
	LL r=1;
	for(;e;x=MUL(x,x),e>>=1)
		if(e&1)
			r=MUL(r,x);
	return r;
}
LL DIV(LL x,LL y){	return MUL(x,POW(y,MOD-2));	}
LL COMB(LL n,LL k)
{
	LL r=1;
	for(LL i=1;i<=k;i++)
		r=DIV(MUL(r,n-i+1),i);
	return r;
}

class SumOfPowers 
{
	public:

		LL n,k;
		LL bin[60][60];
		LL cache[60];

		LL rec(int p)
		{
			LL &r=cache[p];
			if(r!=-1)
				;
			else if(p==0)
				r=n;
			else if(p==1)
				r=n*(n+1)/2%MOD;
			else
			{
				r=(POW(n+1,p+1)-1)%MOD;
				for(int i=0;i<p;i++)
				{
					r-=MUL(bin[p+1][i],rec(i));
					r=(r%MOD+MOD)%MOD;
				}
				r=DIV(r,bin[p+1][p]);
			}
			return r;
		}

		int value(int _n, int _k) 
		{
			n=_n,k=_k;
			for(int i=0;i<60;i++)
				for(int j=0;j<60;j++)
					if(i<j)
						bin[i][j]=0;
					else if(j==0)
						bin[i][j]=1;
					else
						bin[i][j]=(bin[i-1][j]+bin[i-1][j-1])%MOD;
			/*for(int i=1;i<10;i++)
				for(int j=1;j<=i;j++)
					printf("bin[%d][%d]=%lld\n",i,j,bin[i][j]);*/
			memset(cache,-1,sizeof(cache));
			return rec(k);
		}
};

SRM508div1

2011-06-02 02:29:20 JST · 元記事

なにそれオイシーの？

SRM388div1 medium

2011-05-31 22:12:15 JST · 元記事

使うギターを全探索するだけ。

#include <algorithm>
#include <string>
#include <vector>
#include <set>

using namespace std;

class GuitarConcert 
{
	public:

		vector<string> comp(vector<string> vs1,vector<string> vs2)
		{
			if(vs1.size()==vs2.size())
				return (vs1<vs2)?vs1:vs2;
			else
				return (vs1.size()<vs2.size())?vs1:vs2;
		}

		vector <string> buyGuitars(vector <string> guitarNames, vector <string> guitarSongs) 
		{
			int nG=guitarNames.size(),nS=guitarSongs[0].size(),maxPlay=0;

			vector<string> ans;
			for(int use=0;use<(1<<nG);use++)
			{
				set<int> se;
				vector<string> vs;
				for(int i=0;i<nG;i++)
					if(use & (1<<i))
					{
						vs.push_back(guitarNames[i]);
						for(int j=0;j<nS;j++)
							if(guitarSongs[i][j]=='Y')
								se.insert(j);
					}
				sort(vs.begin(),vs.end());

				if(maxPlay<se.size())
				{
					maxPlay=se.size();
					ans=vs;
				}
				else if(maxPlay==se.size())
					ans=comp(ans,vs);
				else
					;
			}
			return ans;
		}
};

SRM388div1 medium

2011-05-30 22:43:15 JST · 元記事

分からなかった。
ビットを使って部分集合を求めると結構速いらしい。
このDPの計算量はそれぞれのビットについて
(mask:0,sub:0),(mask:1,sub:0),(mask:1,sub:1)
の場合があるので3^nなんだって。
2^n*2^nしか思いつかなかった。

#include <vector>

using namespace std;

int OK[1<<15];
int DP[1<<15];

class InformFriends 
{
	public:

		int maximumGroups(vector<string> _friends) 
		{
			int n;
			vector<int> friends;
			n=_friends.size();
			for(int i=0;i<n;i++)
			{
				int a=(1<<i);
				for(int j=0;j<n;j++)
					if(_friends[i][j]=='Y')
						a+=(1<<j);
				friends.push_back(a);
			}
			for(int i=0;i<(1<<n);i++)
			{
				int a=0;
				for(int j=0;j<n;j++)
					if(i&(1<<j))
						a|=friends[j];
				OK[i]=(a==(1<<n)-1);
			}
			for(int mask=0;mask<(1<<n);mask++)
			{
				DP[mask]=0;
				int m=mask;
				while(m!=0)
				{
					if(OK[m])
						DP[mask]=max(DP[mask],DP[mask^m]+1);
					m--;
					m&=mask;
				}
			}
			return DP[(1<<n)-1];
		}
};

SRM386div1 medium

2011-05-29 15:40:21 JST · 元記事

凸な図形＝三角形の集まりなので
rec(どの座標を使ったか(=mask)、今どの座標か(=k))は
kを含む三角形を(i,j,k)の面積をsとすると
s+rec(maskにi,j,kを追加した物,k+1)
の内最小化のものを返せば良い。

あとはそれをメモ化再帰して終わり。
三角形の面積を求めるにはヘロンの公式が便利。

#include <vector>
#include <cmath>
#include <cfloat>

using namespace std;
double cache[(1<<15)][20];

class PolygonCover 
{
	public:

		int n;
		vector<double> x;
		vector<double> y;

		double S(int v1,int v2,int v3)
		{
			double e1,e2,e3,s,x1,y1,x2,y2,x3,y3;
			x1=x[v1]-x[v2];
			y1=y[v1]-y[v2];
			x2=x[v2]-x[v3];
			y2=y[v2]-y[v3];
			x3=x[v3]-x[v1];
			y3=y[v3]-y[v1];

			e1=sqrt(x1*x1+y1*y1);
			e2=sqrt(x2*x2+y2*y2);
			e3=sqrt(x3*x3+y3*y3);
			s=(e1+e2+e3)/2;

			return sqrt(s*(s-e1)*(s-e2)*(s-e3));
		}

		double rec(int mask,int k)
		{
			double &r=cache[mask][k];

			if(-1<r)
				return r;
			else if(k==n)
				r=0.0;
			else if(mask & (1<<k))
				r=rec(mask,k+1);
			else
			{
				r=DBL_MAX;
				for(int i=0;i<n;i++)
					for(int j=i+1;j<n;j++)
						if(i!=k && j!=k)
						{
							int m=(mask|(1<<i)|(1<<j)|(1<<k));
							r=min(r,S(i,j,k)+rec(m,k+1));
						}
			}
			return r;
		}

		double getArea(vector<int> _x, vector<int> _y) 
		{
			x.clear();
			y.clear();
			n=_x.size();
			for(int i=0;i<n;i++)
			{
				x.push_back(_x[i]);
				y.push_back(_y[i]);
			}
			for(int i=0;i<(1<<n);i++)
				for(int j=0;j<n;j++)
					cache[i][j]=-2;

			return rec(0,0);
		}
};

SRM329div1 medium

2011-05-29 12:48:13 JST · 元記事

最大値の最小化なのでセオリー通りに二分探索。

#include <algorithm>
#include <string>
#include <vector>
#include <sstream>

using namespace std;

class LogCutter 
{
	public:

		int n;
		vector<int> canCut;

		bool check(int L,int C,int mid)
		{
			int beforeCut=0;
			while(mid<L-canCut[beforeCut] && 0<C)
			{
				if(mid<canCut[beforeCut+1]-canCut[beforeCut])
					return false;
				int nextCut=beforeCut+1;
				while(nextCut+1<n && canCut[nextCut+1]-canCut[beforeCut]<=mid)
					nextCut++;
				C--;
				beforeCut=nextCut;
			}

			return (L-canCut[beforeCut]<=mid);
		}

		string bestCut(int L, int A, int K, int C) 
		{
			long long ll=L,aa=A;
			canCut.clear();
			for(long long i=1;i<=K;i++)
				canCut.push_back((aa*i)%(ll-1)+1);
			canCut.push_back(0);
			canCut.push_back(L);
			sort(canCut.begin(),canCut.end());
			canCut.erase(unique(canCut.begin(),canCut.end()),canCut.end());
			n=canCut.size();

			int ceil=L,bottom=1;
			while(1<ceil-bottom)
			{
				int mid=(ceil+bottom)/2;
				if(check(L,C,mid))
					ceil=mid;
				else
					bottom=mid;
			}

			int opt=(check(L,C,bottom)?bottom:ceil);
			int bc=n-1,fc,fisrtCut;
			while(1)
			{
				fc=bc-1;
				while(canCut[bc]-canCut[fc-1]<=opt && 0<=fc-1)
					fc--;
				if(fc==0)
				{
					if(0<C)
						fisrtCut=1;
					else
						fisrtCut=bc;
					break;
				}
				bc=fc;
				C--;
			}
			stringstream ans;
			ans << opt << " " << canCut[fisrtCut];
			return ans.str();
		}
};

SRM410 div1 medium

2011-05-27 23:14:26 JST · 元記事

#include <algorithm>
#include <vector>
#include <cstdio>
#include <cstring>
#include <climits>

using namespace std;
typedef long long LL;

int n;
LL k;
int range[55][55][2][2];
LL cache[55][55][2];

class ContiguousCache 
{
	public:

		LL intersect(int a1,int a2,int b1,int b2)
		{
			if(a1<=b1 && b1<=a2 && a2<=b2)
				return a2-b1+1;
			else if(a1<=b1 && b2<=a2)
				return b2-b1+1;
			else if(b1<=a1 && a2<=b2)
				return a2-a1+1;
			else if(b1<=a1 && a1<=b2 && b2<=a2)
				return b2-a1+1;
			else
				return 0;
		}

		LL rec(int a1,int a2,int slide)
		{
			int b1=a2+1,b2,nextSlide;
			LL &r=cache[a1][a2][slide];

			if(r!=-1)
				;
			else if(b1==n)
				r=0;
			else	
			{
				r=LLONG_MAX;
				for(b2=b1;b2<n && range[b1][b2][0][0]!=-1;b2++)
					for(nextSlide=0;nextSlide<2;nextSlide++)
						r=min(r,k-intersect(range[a1][a2][slide][0],range[a1][a2][slide][1],range[b1][b2][nextSlide][0],range[b1][b2][nextSlide][1])+rec(b1,b2,nextSlide));
			}

			return r;
		}

		LL minimumReads(int _n,int _k,vector<int> addresses) 
		{
			reverse(addresses.begin(),addresses.end());
			addresses.push_back(INT_MIN/2);
			reverse(addresses.begin(),addresses.end());

			n=addresses.size(),k=_k;
			for(int i=0;i<n;i++)
				for(int j=0;j<n;j++)
				{
					LL mi=INT_MAX/2,ma=INT_MIN/2;
					for(int l=min(i,j);l<=max(i,j);l++)
					{
						mi=min(mi,(LL)addresses[l]);
						ma=max(ma,(LL)addresses[l]);
					}
					range[i][j][0][0]=-1,range[i][j][0][1]=-1;
					range[i][j][1][0]=-1,range[i][j][1][1]=-1;
					if(ma-mi<k)
					{
						range[i][j][0][0]=mi,range[i][j][0][1]=mi+k-1;
						range[i][j][1][0]=ma-k+1,range[i][j][1][1]=ma;
					}
				}

			memset(cache,-1,sizeof(cache));
			return rec(0,0,0);
		}
};

SRM439 div1 medium

2011-05-08 22:33:16 JST · 元記事

メモ化再帰。
難しくて解説みた。

#include <string>
#include <vector>
#include <cstdio>
#include <cstring>

using namespace std;
typedef long long LL;

LL DP[(1<<13)][2][15][15];

class PalindromePhrases 
{
	public:

		int n;
		vector<string> words;

		bool isPalindome(string s)
		{
			for(int i=0;i<s.size();i++)
				if(s[i]!=s[s.size()-1-i])
					return false;
			return true;
		}

		LL rec(int left,int isLeft,int id,int len)
		{
			LL &r=DP[left][isLeft][id][len];
			if(r!=-1)
				;
			else if(isLeft)
			{
				r=0;
				int pos=words[id].size()-len;
				string s=words[id].substr(pos);
				if(isPalindome(s))
					r++;
				for(int next=0;next<n;next++)
					if(left & (1<<next))
						if(len<words[next].size())
						{
							int i;
							for(i=0;i<len;i++)
								if(words[id][pos+i]!=words[next][words[next].size()-1-i])
									break;
							if(i==len)
								r+=rec(left^(1<<next),0,next,words[next].size()-i);
						}
						else
						{
							int i;
							for(i=0;i<words[next].size();i++)
								if(words[id][pos+i]!=words[next][words[next].size()-1-i])
									break;
							if(i==words[next].size())
								r+=rec(left^(1<<next),1,id,len-i);
						}
			}
			else
			{
				r=0;
				int pos=len-1;
				string s=words[id].substr(0,len);
				if(isPalindome(s))
					r++;
				for(int next=0;next<n;next++)
					if(left & (1<<next))
						if(len<words[next].size())
						{
							int i;
							for(i=0;i<len;i++)
								if(words[id][i]!=words[next][len-1-i])
									break;
							if(i==len)
								r+=rec(left^(1<<next),1,next,words[next].size()-i);
						}
						else
						{
							int i;
							for(i=0;i<words[next].size();i++)
								if(words[id][pos-i]!=words[next][i])
									break;
							if(i==words[next].size())
								r+=rec(left^(1<<next),0,id,len-i);
						}
			}
			return r;
		}

		LL getAmount(vector <string> _words) 
		{
			words=_words;
			n=words.size();

			memset(DP,-1,sizeof(DP));
			LL ans=0;
			for(int i=0;i<n;i++)
				ans+=rec(((1<<n)-1)^(1<<i),1,i,words[i].size());
			return ans;
		}
};

Google Code Jam Qualification Round

2011-05-08 22:33:15 JST · 元記事

満点で予選通過した。
Round1は中間テストの真っ最中だけど、頑張って出るかもしれない。

SRM479 div1 medium

2011-05-07 19:49:20 JST · 元記事

二部探索＋ダイクストラ。

#include <algorithm>
#include <string>
#include <vector>
#include <sstream>
#include <queue>
#include <utility>

using namespace std;
typedef long long LL;

int n;
LL limit;
LL INF=((LL)1<<LL(50));
LL visited[500];
vector<int> F[500];
vector<int> T[500];
vector<int> P[500];
vector<int> D[500];

class TheAirTripDivOne 
{
	public:

		bool check(LL wait)
		{
			for(int i=0;i<n;i++)
				visited[i]=INF;
			priority_queue<pair<LL,int>,vector<pair<LL,int> >,greater<pair<LL,int> > > que;
			que.push(make_pair(0,0));

			bool succes=false;
			while(!que.empty())
			{
				LL now=que.top().first;
				int pos=que.top().second;
				que.pop();

				if(visited[pos]<now || limit<now)
					continue;
				else if(pos==n-1)
				{
					succes=true;
					break;
				}
				else
				{
					visited[pos]=now;
					if(pos!=0)
						now+=wait;
					else
						;

					for(int i=0;i<F[pos].size();i++)
						if(F[pos][i]!=INF)
						{
							LL ne=F[pos][i]+((now<F[pos][i])?0:
								(now-F[pos][i]+P[pos][i]-1)/P[pos][i]*P[pos][i])+T[pos][i];
							que.push(make_pair(ne,D[pos][i]));
						}
				}
			}
			return succes;
		}

		int find(int _n, vector <string> flights, int time) 
		{
			limit=time;
			n=_n;

			for(int i=0;i<500;i++)
				F[i].clear(),T[i].clear(),P[i].clear(),D[i].clear();

			string input;
			for(int i=0;i<flights.size();i++)
			{
				for(int j=0;j<flights[i].size();j++)
					if(flights[i][j]==',')
						flights[i][j]=' ';
				input+=flights[i];
			}
			stringstream ss(input);
			int a,b,f,t,p;
			while(ss >> a >> b >> f >> t >> p)
			{
				a--,b--;
				F[a].push_back(f);
				T[a].push_back(t);
				P[a].push_back(p);
				D[a].push_back(b);
			}

			int ceil=time,bottom=1;
			while(5<ceil-bottom)
			{
				int wait=(ceil+bottom)/2;

				if(check(wait))
					bottom=wait;
				else
					ceil=wait;
			}
			for(int i=ceil;bottom<=i;i--)
				if(check(i))
					return i;

			return -1;
		}
};

SRM464 div1 medium

2011-05-05 12:22:11 JST · 元記事

二部探索＋２SATらしい。

#include <vector>
#include <cstdio>
#include <cstring>

using namespace std;

class ColorfulDecoration 
{
	public:

		int n;
		int len;
		int dist[110][110];
		int G[110][110];

		int ABS(int n)
		{
			return (n<0)?-n:n;
		}

		bool check()
		{
			memset(G,0,sizeof(G));
			for(int i=0;i<2*n;i++)
				for(int j=0;j<2*n;j++)
					if(i%n!=j%n && dist[i][j]<len)
					{
						int ii=(i<n)?i+n:i-n;
						int jj=(j<n)?j+n:j-n;
						G[i][jj]=1;
					}

			for(int i=0;i<2*n;i++)
				G[i][i]=1;

			for(int k=0;k<2*n;k++)
				for(int i=0;i<2*n;i++)
					for(int j=0;j<2*n;j++)
						if(G[i][k]==1 && G[k][j]==1)
							G[i][j]=1;

			for(int i=0;i<n;i++)
				if(G[i][i+n] && G[i+n][i])
					return false;
			return true;
		}

		int getMaximum(vector <int> xa, vector <int> ya, vector <int> xb, vector <int> yb) 
		{
			n=xa.size();
			for(int i=0;i<2*n;i++)
				for(int j=0;j<2*n;j++)
				{
					int x1=(i<n)?xa[i]:xb[i-n];
					int y1=(i<n)?ya[i]:yb[i-n];
					int x2=(j<n)?xa[j]:xb[j-n];
					int y2=(j<n)?ya[j]:yb[j-n];
					dist[i][j]=max(ABS(x1-x2),ABS(y1-y2));

					//printf("dist[%d][%d]=%d\n",i,j,dist[i][j]);
				}
			int ceil=1000000000,bottom=0;
			while(5<ceil-bottom)
			{
				len=(ceil+bottom)/2;
				if(check())
					bottom=len;
				else
					ceil=len;
			}
			for(len=ceil;bottom<=len;len--)
			{
				if(check())
					return len;
			}
			return 0;
		}
};