飽きるまでProject Eulerに挑戦中

はじめに

ネタバレ注意

タイトルの通り。 まずは一番簡単なレベルの問題を制覇したい。

Project Euler とは

About - Project Euler

一言で言うと数学絡みのプログラミングの問題を出しているサイト。 どの問題も1分以内で答えが出せる設定らしい。 コンピュータで力任せに計算して答えを出すこともできるが、数学的にエレガントな解き方があったりするとかしないとか。

問題例

ネタバレは推奨されていないが、そこら中に解法が転がっている。 はじめから解く気がない人や既に解いた人が見るのはいいと思う。 そもそもサイト上に解いた人が読める解説PDF付の問題が存在している。

Wikipediaで解説されているProblem 1と似たようなProblem 28 Number spiral diagonalsを問題例として紹介する。

Problem 28 - Project Euler

 1001\times1001の螺旋状に並んだ数字の対角線部分の数字の和を求めよという問題。

21 22 23 24 25
20  7  8  9 10
19  6  1  2 11
18  5  4  3 12
17 16 15 14 13

5\times5の例では 1 + 3 + 5 + 7 + 9 + 13 + 17 + 21 + 25 = 101になる。

実際に螺旋状に辿って角に来たときに足して行くのもよいが、正直自分にはそれは実装できない(できるが)。

足すべき数字を列挙すると次のようになる。

 1 | 3, 5, 7, 9 | 13, 17, 21, 25 | \dots

懐かしき群数列。

今考えるのは 1001\times1001なのでとりあえず1は置いておく。 まず螺旋の大きさは(奇数)x(奇数)になる。 数字は1から始まり1ずつ増えていくので最大の数字はその(奇数)の2乗。

|で区切られた(群という)四つの項は外側二つの和と内側二つの和が等しいので群の最初の項と最後の項がわかれば、それらを足して2倍するとその群の項の和が得られる。

例:  | 3, 5, 7, 9|に注目すると、 3 + 9 = 5 + 7 = 12。この群の項の和は2 \times 12 = 24

1を第0群とする。 i \geq 1のとき、第 i群の最後の項は  { (2i + 1)^{2} } であり、最初の項は前の群の最後の項より 2i大きい。

したがって、第i群の項の和は \displaystyle 2\{(2i - 1)^{2} + 2i + (2i + 1)^{2}\} = 16i^{2} + 4i + 4 \, (i \geq 1)である。

よって、第0群から第n群までの項の和は

 {\displaystyle
1 + \sum_{i=1}^{n} (16i^{2} + 4i + 4) = 1 + \frac{16}{3}n^{3} + 10n^{2} + \frac{26}{3}n = 1 + \frac{2n(n(8n + 15) + 13)}{3}
}

である(n=0のときも成立つ)。

 1001\times1001の螺旋で最大の数は 1001^{2}だから、第500群まである。 n = 500で計算すれば答えが出る。

#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>

uint64_t sum(uint64_t n)
{
    return 1 + 2 * n * (n * (8 * n + 15) + 13) / 3;
}

int main(void)
{
    uint64_t ans = sum(500);
    printf("%"PRIu64"\n", ans);
    return 0;
}

Rust テンプレート

自分用

Rust勉強中

まだCっぽい書き方しかできない

入力

1行標準入力

fn readline() -> String {
    let mut buf = String::new();
    std::io::stdin().read_line(&mut buf).unwrap();
    buf.trim_end()
}

ファイル入力

use std::io::{BufRead, BufReader};
use std::fs::File;

fn main() {
    let reader = BufReader::new(File::open("file.txt").unwrap());
    for line in reader.lines() {
        let buf = line.unwrap();
    }
}

素数

素数判定

TODO: AKS素数判定法を理解・実装

fn is_prime(n: usize) -> bool {
    if n <= 3 {
        return n > 1;
    } else if n % 2 == 0 || n % 3 == 0 {
        return false;
    }
    let mut i = 5;
    while i * i <= n {
        if n % i == 0 || n % (i + 2) == 0 {
            return false;
        }
        i += 6;
    }
    true
}

エラトステネスのふるい

そのまんま

// All the elements of `is_prime` must be initialized to true
// prior to calling this function.
fn sieve(is_prime: &mut [bool]) {
    is_prime[0] = false;
    is_prime[1] = false;
    let n = is_prime.len();
    let mut i = 2;
    while i * i <= n {
        if is_prime[i] {
            for j in (i * i..n).step_by(i) {
                is_prime[j] = false;
            }
        }
        i += 1;
    }
}

反転

fn reverse_num(mut n: usize) -> usize {
    let mut r = 0;
    while n > 0 {
        r *= 10;
        r += n % 10;
        n /= 10;
    }
    r
}

桁分解

fn digits_vec(mut n: usize, base: usize) -> Vec<usize> {
    let mut digits = vec![];
    while n > 0 {
        digits.push(n % base);
        n /= base;
    }
    digits.into_iter().rev().collect()
}

最大公約数と最小公倍数

fn gcd(mut a: usize, mut b: usize) -> usize {
    if a < b {
        std::mem::swap(&mut a, &mut b);
    }
    while b != 0 {
        let tmp = a;
        a = b;
        b = tmp % b;
    }
    a
}

fn lcm(a: usize, b: usize) -> usize {
    a / gcd(a, b) * b
}

約数の和

TODO: bを消す?(もっと簡単にできる?)

proper divisors(日本語でなんと言うか知らず)の和も載せとく

fn sum_of_divisors(n: usize) -> usize {
    if n == 0 {
        return 0;
    } else if n == 1 {
        return 1;
    }
    let mut sum = n + 1;
    let mut a = 2;
    let mut b = n;
    while a < b {
        if n % a == 0 {
            b = n / a;
            sum += if a == b { a } else { a + b };
        }
        a += 1;
    }
    sum
}

fn sum_of_proper_divisors(n: usize) -> usize {
    sum_of_divisors(n) - n
}

多倍長整数

num を使う

Cargo.toml

[package]
name = "problem_0xxx"
version = "0.1.0"

[dependencies]
num-bigint = "0.2"
num-traits = "0.2"

ソースコード

extern crate num_bigint;
extern crate num_traits;

use num_bigint::BigUint;
use num_traits::{pow, Zero, One};

fn main() {
}

階乗

このまま使うことはないだろうけど

fn fac(n: usize) -> BigUint {
    let mut r: BigUint = One::one();
    for i in 2..=n {
        r *= i;
    }
    r
}

フィボナッチ数

同じく

fn fib(n: usize) -> BigUint {
    if n == 0 {
        return Zero::zero();
    }
    let mut pre = Zero::zero();
    let mut cur = One::one();
    for i in 0..n {
        let next = &pre + &cur;
        pre = cur;
        cur = next;
    }
    cur
}

その他

next_permutation

stableにはこれ系がない

Next lexicographical permutation algorithm

ツェラーの公式

使うのか?

fn zeller(mut y: i32, mut m: i32, d: i32) -> i32 {
    let c = y / 100;
    let g = 5 * c + c / 4;
    if m == 1 || m == 2 {
        m += 12;
        y -= 1;
    }
    y %= 100;
    (d + 26 * (m + 1) / 10 + y + y / 4 + g) % 7
}

やっとLevel 1になった

A*の実装

はじめに

A*を実装してみた。 実装したのはA* だが、ヒューリスティック関数が常に0を返すため、動作はダイクストラ法と変わらない。

A*とは

A*は経路探索でよく用いられるアルゴリズムダイクストラ法に「現在の点から終点までの推定コスト」を追加することで、より効率的に最短経路を見つけることが可能になる。 その推定コストを計算する関数をヒューリスティック関数というが、その設計がよくないと効率的にはならない。 詳しくは調べたほうが早い。 気が向いたら書く。

ソースコード

ほとんど出席しなかった大学の授業で札幌駅-新千歳空港駅間の最安乗り換えの話がでていたので、 最安乗り換えの一つを求めるプログラムを自分で実装してみた。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# A* search algorithm
# reference: https://en.wikipedia.org/wiki/A*_search_algorithm

import sys
from collections import defaultdict
from dataclasses import dataclass, field
from typing import Any, Dict, List
from heapq import heappush, heappop

@dataclass(order=True)
class PrioritizedItem(object):
    priority: int
    item: Any=field(compare=False)

@dataclass
class PriorityQueue(object):
    que: List['PrioritizedItem']=field(default_factory=list, init=False)

    def push(self, priority, item):
        heappush(self.que, PrioritizedItem(priority, item))

    def pop(self):
        return heappop(self.que).item

    def top(self):
        return self.que[0].item

    def __len__(self):
        return len(self.que)

INF = 10000000

NULL = 0
OPEN = 1
CLOSE = 2

@dataclass
class Vertex(object):
    name: str=field(default=None)
    costs: Dict[int, int]=field(default_factory=dict, init=False)
    # other fields for the heuristic function

def estimate_cost(u, v):
    return 0 # Dijkstra's algorithm

def reconstruct_path(paths, current_id):
    path_ids = [current_id]
    while current_id in paths:
        current_id = paths[current_id]
        path_ids.append(current_id)
    path_ids.reverse()
    return path_ids

def a_star(vertices, start_id, goal_id):
    paths = {}
    goal = vertices[goal_id]
    start = vertices[start_id]
    states = [NULL] * len(vertices)
    states[start_id] = OPEN
    start_priority = estimate_cost(start, goal)
    priorities = { start_id:start_priority }
    costs = defaultdict(lambda: INF)
    costs[start_id] = 0
    que = PriorityQueue()
    que.push(start_priority, start_id)
    while len(que):
        current_id = que.pop()
        if states[current_id] == CLOSE:
            continue
        current_cost = costs[current_id]
        if current_id == goal_id:
            return current_cost, reconstruct_path(paths, current_id)
        states[current_id] = CLOSE
        current = vertices[current_id]
        for neighbor_id, cur_nbr_cost in current.costs.items():
            if states[neighbor_id] == CLOSE:
                continue
            neighbor = vertices[neighbor_id]
            neighbor_cost = current_cost + cur_nbr_cost
            if states[neighbor_id] != OPEN:
                states[neighbor_id] = OPEN
            elif neighbor_cost >= costs[neighbor_id]:
                continue
            paths[neighbor_id] = current_id
            costs[neighbor_id] = neighbor_cost
            priority = neighbor_cost + estimate_cost(neighbor, goal)
            priorities[neighbor_id] = priority
            que.push(priority, neighbor_id)
    return 0, list()

def main():
    n, start_id, goal_id = map(int, sys.stdin.readline().split())
    vertices = [Vertex(sys.stdin.readline().rstrip()) for i in range(n)]
    while True:
        line_split = sys.stdin.readline().split()
        if len(line_split) == 1 and int(line_split[0]) == -1:
            break
        i, j, cost = map(int, line_split)
        vertices[i].costs[j] = cost
        vertices[j].costs[i] = cost

    total_cost, path_ids = a_star(vertices, start_id, goal_id)
    print('from:', vertices[start_id].name)
    print('to:', vertices[goal_id].name)
    print('cost:', total_cost)
    print('path:')
    for i in path_ids:
        v = vertices[i]
        print(v.name)

if __name__ == '__main__':
    main()

実行例

2019-02-27現在の運賃を使って計算した。

入力

8 0 7
新千歳空港
南千歳
千歳
恵庭
北広島
新札幌
白石
札幌
0 1 310
0 2 350
0 3 400
0 4 590
0 5 880
0 6 980
0 7 1070
1 2 170
1 3 260
1 4 450
1 5 640
1 6 740
1 7 840
2 3 220
2 4 360
2 5 640
2 6 740
2 7 840
3 4 260
3 5 450
3 6 540
3 7 640
4 5 260
4 6 360
4 7 450
5 6 210
5 7 260
6 7 210
-1

実行結果

from: 新千歳空港
to: 札幌
cost: 1040
path:
新千歳空港
恵庭
札幌

30円安くなる。時は金なり。

参考サイト

A* search algorithm - Wikipedia

1から100までFizz Buzz

現実逃避。1から100までFizz Buzzを表示するというよくあるやつ。

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

char *x[] = {
    "Fizz",
    "Buzz",
    "Fizz Buzz"
};
int y[][5]={
    {2, -2, -2, -2, -2},
    {-1, 0, -1,-1, -2},
    {-2, 1, 0, -1, -2},
    {2, -2, 1, 0, -1}
};
char z[100][4];

char *f(int a, int b, int c)
{
    if (a == 1)
        return "1";
    if (!c || y[abs(y[3][b])][abs(c) - 1] == -2) {
        b = (b + 1) % 5;
        c = y[3][b] < 0 ? -4 : 1;
    }
    while (y[abs(y[3][b])][abs(c) - 1] == -2)
        c++;
    puts(f(a - 1, b, c + 1));
    return y[abs(y[3][b])][abs(c) - 1] == -1 ?
        sprintf(z[a - 1], "%d", a), z[a - 1] :
        x[y[abs(y[3][b])][abs(c) - 1]];
}

int main(void)
{
    puts(f(100, -1, 0));
    return 0;
}

某VPNの設定を半自動化する小汚い小スクリプト群

何をしているのかわからない人はやらないほうがいいです

俺もわからん.

VPNはセットアップが楽ちん

さらに楽にするために半自動化します. その場しのぎで書いただけで何も考えてません. 使用, 参考にする場合は要注意.

筆者の環境はArch Linuxだが, 執筆時にちょっと書き換えて動作確認はしていない.

はい

連番振りたいが行をあけると無理?

1.こちらを導入:

GitHub - flamusdiu/python-pia: Commandline tool to auto configure PIA services

ポート番号でアルゴリズムを設定してくれるのでこことか参考(特に気にしないのであればUDP 1197でよいと思う):

What's the difference between the OVPN files? - Knowledgebase / Technical / OVPN Files - PIA Support Portal

2. vpn-hosts.txtのリストにプレフィクスを付けるスクリプト:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

def main():
    path = '/etc/private-internet-access/'
    hosts_file = path + 'vpn-hosts.txt'
    with open(hosts_file) as f:
        lines = f.readlines()
    newlines = ['PIA ' + line for line in lines]
    print(*newlines, sep='')
    yn = input('write to file?[y/N]: ')
    if yn == 'y':
        with open(hosts_file, 'w') as f:
            f.writelines(newlines)

if __name__ == '__main__':
    main()

3. pia.confvpn-hosts.txtを全部突っ込むスクリプト:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

def main():
    path = '/etc/private-internet-access/'
    hosts_file = path + 'vpn-hosts.txt'
    conf_file = path + 'pia.conf'

    with open(hosts_file, 'r') as f:
        countries = [line.strip().split(',')[0] for line in f.readlines()]
    with open(conf_file, 'r') as f:
        conf_lines = f.readlines()

    for i, line in enumerate(conf_lines):
        if line.startswith('hosts'):
            conf_lines[i] = 'hosts = {}\n'.format(', '.join(countries))
            break

    print(*conf_lines, sep='')

    yn = input('write to file?[y/N]: ')
    if yn == 'y':
        with open(conf_file, 'w') as f:
            f.writelines(conf_lines)

if __name__ == '__main__':
    main()

4. /etc/openvpn/client/以下にOpenVPN設定ファイルを展開:

# pia -a

5. 公開鍵証明書ファイル名にプレフィクス付ける(俺の環境では動いてるが初回は自分で調整して):

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import os

def main():
    path = '/etc/openvpn/client/'
    pia_confs = [fname for fname in os.listdir(path) if os.path.isfile(path + fname) and fname.startswith('PIA')]

    conf = pia_confs[0]
    with open(path + conf) as f:
        lines = f.readlines()
    names = dict()
    targets = ('ca', 'crl-verify')
    for i, line in enumerate(lines):
        l = line.split(' ')
        if l[0] in targets:
            names[l[0]] = l[1].rstrip()
            if len(names) >= len(targets):
                break
    for tgt, name in names.items():
        pathname, basename = os.path.split(name)
        newname = '{}/pia_{}'.format(pathname, basename)
        if os.path.exists(name):
            os.rename(name, newname)
        names[tgt] = newname

    for cnf in pia_confs:
        with open(path + cnf) as f:
            lines = f.readlines()
        for i, line in enumerate(lines):
            l = line.split(' ')
            for tgt, name in names.items():
                if l[0] == tgt:
                    lines[i] = '{} {}\n'.format(tgt, name)
                    break
        with open(path + cnf, 'w') as f:
            f.writelines(lines)

if __name__ == '__main__':
    main()

6. キルスイッチ

キルスイッチについてはこちらを参照(固定化のデメリットがある):

Internet Kill Switch by Firewall (OpenVPN + iptables) - vanaestea’s blog

DNS

そのままだと/etc/resolv.confVPNサービス側のDNSに書き換えられたと記憶しているが, それが嫌なら, たとえば, resolv.conf

nameserver 1.1.1.1
nameserver 1.0.0.1

のように好きに設定したあと,

# chattr +i /etc/resolv.conf

でイミュータブル属性を付加するとよい(参考: Private Internet Access - ArchWiki).

SOLVED: なぜセグフォ?

最近現実逃避が捗りすぎて困っている.

最近は低レベルのコードとかよくわからないので勉強している(遊んでいるだけ).

このC言語ソースコードは, x86_64のLinuxなら動く環境があると思う.

const unsigned char main[] = { 72, 199, 192, 1, 0, 0, 0, 72, 199, 195, 1, 0, 0, 0, 72, 141, 53, 21, 0, 0, 0, 72, 199, 194, 21, 0, 0, 0, 15, 5, 72, 199, 192, 60, 0, 0, 0, 72, 49, 255, 15, 5, 97, 108, 119, 97, 121, 115, 32, 109, 105, 115, 115, 32, 116, 104, 101, 32, 109, 97, 114, 107, 10 };

Wandboxでは動作した: https://wandbox.org/permlink/KeSJ6KkTDQGUXc3z

しかし自分の環境(Arch Linux)ではセグフォで落ちてしまう. main()に入った瞬間に落ちているようだ.

$ gdb ./a.out 
GNU gdb (GDB) 8.2
Copyright (C) 2018 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
...
(gdb) r
Starting program: /tmp/tmp.lX2FcrqB9C/a.out 

Program received signal SIGSEGV, Segmentation fault.
0x0000555555556020 in main ()
(gdb) disass 0x0000555555556020
Dump of assembler code for function main:
=> 0x0000555555556020 <+0>:     mov    $0x1,%rax
   0x0000555555556027 <+7>:     mov    $0x1,%rbx
...
End of assembler dump.
(gdb) 

なんもわからん.完全に理解した.

おそらくExec Shieldが働いている(?)のでコンパイル時に-z execstackを渡す必要がある.