Help us understand the problem. What is going on with this article?

IchigoJam で円周率をファインマン・ポイントまで計算表示したい

More than 3 years have passed since last update.

はじめに

IchigoJam は32桁×24行 = 768文字 の表示画面を持ちます。
円周率を 10進数で表現すると小数点以下 762桁から767桁に 9 が 6つ並ぶ箇所があり、俗に『ファインマン・ポイント』と呼ばれているようです。
IchigoJam の表示画面ではそれがギリ表示できるので挑戦してみました。

プログラム

pi.bas
1 ' 999999 in pi
2 Z=#C03:Y=#700
3 X=0
4 IFPEEK(Z)=39X=X-1:W=33-PEEK(Z-X):IFW<1V=V<<7-W+W/96*34:POKEY,V>>(X&7):Y=Y+(X&7<7):GOTO4
5 Z=Z+PEEK(Z-1)+4:IFPEEK(Z-3)+PEEK(Z-2):GOTO3
6 Z=USR(#700,0):WAIT600:GOTO6
10 'サNA3M7」"inJQ*(}9\!'!!$、*ソ"%1A!」!!]1pa=!&1m!!(bツサS!=エ|ソ!`."e4!!"ウ\チ(エrォ7"7b;:"%$<@dB1+タJ12aT2y(w9ィ65#9%ュ|!HIN?h!Gア/」#5エ.W
11 '!!&AGgAD!ゥDQ/FYHウ"ァ>9}7ang<R(MVTB!\}"/W!!)SU1ェ[+4aT-A」Ei1'!(d、'b-%J);)(}i)!!*1QBY9)」ii(1"ゥ!#a,[$m!H0h=#\エシ」4Mag"m'」>9」」R
12 '!my6Fj&hォa%!Nd、^ツツツコツソ1!!`」A!!BV!)」3B"%$ォ7ヲ*;)"tngi!a4Aj!CiPG・3;シZヲ!!cス!2q!&y#Ac!キウ"q%ゥ"=a&@WeQ!Y2チ,J:6(yAヲ.y!*4O'!#9/A&
13 'ird」-nゥ=5U!Zy!b@.e%q*QRqsqa>:UXセカaEQ5m#D0qa~!+SC!+!0ツ「"ツ"E0ウa<AJ!+!#1$ツコ}`」q(Xo#m)a"B(ツコ、ツ」:y%ソ<-!SQbカJェ#ゥaj}r-ヲC(!!"9IF
14 '9..A"E"C.!E0ツ。「ツ#0G、a!|クォaT`!」0QサPY%Xツツツ0a

※ あれ? よく見たら 5行目の GOTO3の前の ':' 要らんな? これ削除すると動画から撮り直しになってメンドいのでとりあえず放っておこう。

実行中の動画

IchigoJam で円周率をファインマン・ポイントまで計算表示 pic.twitter.com/0egsGyTpka

— fujita nozomu (@fujitanozomu) 2016年11月25日

入力方法

プログラムは ASCIIコード + 半角カナ のみで構成されているため IchigoJam に接続したキーボードから普通に入力可能ですが、大変なので、USB シリアル変換装置を介して接続した PC からコピペするのが現実的でしょう。Teraterm 等汎用のターミナルプログラムを使用するのであれば、半角カナのコードが IchigoJam のそれと等しくなるようターミナルプログラムの文字コードの設定を SJIS か JIS に設定する必要があります。
文字の取りこぼしが発生することがあるので、入力したプログラムは LIST コマンドで表示し、正しく入力できたか確認するべきです。10行目以下は一見ランダムなコメント文に見えますが、一文字でも間違えるとプログラムは正しく実行されずエラーも表示されずに暴走する可能性があるため、実行の前には SAVE コマンドで保存をした方が良いでしょう。

プログラムの説明

パソコン等で円周率を計算するプログラムは数多く発表されていますが、多くのものは任意精度演算を使用しているために多くのメモリを使用します。IchigoJam は標準的に使用できる言語は BASIC ですが、その仕様では使用できるメモリが BASIC プログラムを格納する領域の 1024バイトと -32768から32767の値を格納できる A から Z までの 26個の変数と 102個の要素数の配列変数だけであり、今回の、768桁の円周率の計算を任意精度演算で行うには不足しています。他の方法を検討する必要があります。

IchigoJam の RAM

IchigoJam は、以下の領域が RAM として公開されており、ユーザーが本来の機能とぶつからない様気を付ければデータとして読み書きしたり機械語プログラムを配置したり自由に使用できるようなっています。

アドレス領域 本来の機能 注釈
#0700~#07FF PCG 領域 PCG 機能を使用しなければ自由に使用できる
#0800~#08FF 変数領域 配列変数 0~101, 変数 A~Zの順で各変数はアドレス固定で格納される。使用していない変数の領域は自由に使用できる
#0900~#0BFF VRAM 領域 画面左上隅から右方向、次行と続き、画面右下隅まで
#0C00~#0FFF BASIC プログラム格納領域 BASIC プログラムとぶつからない様気を付ければ使用可能

今回のプログラムでは隣接する PCG 領域と変数領域を跨いで機械語プログラムを配置し、VRAM 領域を作業領域兼表示用として使用することとしました。

BBP

円周率を計算する多くのプログラムは任意精度演算を使用するものであり、求める桁数の数倍のメモリを要します。今回は IchigoJam で使用できるメモリに限りがあるという理由により、BBP(Bailey–Borwein–Plouffe formula)というアルゴリズムを使用しています。内容の詳細は Wikipedia 等に譲りますが、円周率を 2進数で任意の桁の値を僅かな記憶領域の使用で求められるというものであり、今回の目的に合っています。

プログラム詳細

今回のプログラムは 4000 Stellen von Pi mit ATtiny2313 を下敷きにしています。上記にて公開されているプログラムを整数演算化し、2進→10進の変換を加えています。
10進数で 768桁の数値は 2進数でおよそ 2552ビット(319バイト)の情報量となりますが、2進→10進の変換の際に捨てられる値分の誤差を考慮して 2560ビット(320バイト)の計算を行っています。出力先として VRAM を使用し 320バイト(=32桁×10行)の書き込みを行った後、VRAM 上で 10/256 の計算を桁分行うことで 10進数へと変換しています。2進数での円周率の計算と比べて 2進→10進 の変換は軽いものですが、視覚的にカッコイイのでウェイト処理を追加し変換の様子が見えるようしています。
ソースコードは以下の内容となります。

pi.c
#include <stdint.h>

#define Digits (1 + 767)
#define ComputeSize (1 + ((int)(/*log(10)/log(2)*/3.32192809489 * (Digits - 1)) + 7) / 8)
#define DecimalPointPosition 20
#define Wait 100

static unsigned pi_n(unsigned n);

int main(int16_t arg, uint8_t* ram, const uint8_t* font)
{
    uint8_t* vram = ram + 0x900;

    for (int i = 1; i < ComputeSize; i++) {
        unsigned x = pi_n(2 * (i - 1)) >> (DecimalPointPosition - 8);
        if (i == 1) {
            vram[0] = x >> 8;
        }
        vram[i] = x;
    }
    for (int i = 0, j = ComputeSize; i < Digits; i++) {
        int c;
        if (i == 0) {
            c = vram[i];
        } else {
            c = 0;
            for (int k = j - i / 8, l = (k < Digits) ? k : Digits; l > i; l--) {
                c += 10 * vram[l - 1];
                if (l < Digits) {
                    vram[l] = c;
                }
                c >>= 8;
#if Wait > 0
                for (volatile int w = Wait; w > 0; w--) {
                    ;
                }
#endif
            }
            j++;
        }
        vram[i] = '0' + c;
    }
    return 0;
}

static unsigned mod_mul(unsigned a, unsigned b, unsigned n)
{
    unsigned ab = 0;

    while (1)
    {
        if (a % 2 == 1)
        {
            ab += b;
            if (ab >= n)
                ab -= n;
        }

        a = a / 2;
        if (a == 0)
            return ab;

        b = b * 2;
        if (b >= n)
            b -= n;
    }
}

static unsigned mod_pow16(unsigned k, unsigned n)
{
    unsigned p = 1;
    unsigned _16 = 16;

    if (n == 1)
        return 0;

    while (_16 >= n)
        _16 -= n;

    while (1)
    {
        if (k % 2 == 1)
            p = mod_mul(_16, p, n);

        k = k / 2;
        if (k == 0)
            break;

        _16 = mod_mul(_16, _16, n);
    }

    return p;
}

static __attribute__((noinline)) unsigned divide(unsigned a, unsigned b, unsigned c)
{
    unsigned r = 0;
    do {
        a <<= 1;
        r <<= 1;
        if (a >= c) {
            a -= c;
            r++;
        }
    } while (--b);
    return r;
}

static unsigned sigma_a(unsigned n, uint8_t j)
{
    unsigned s = 0;
    for (unsigned k = n-1; k+1 != 0; k--)
    {
        unsigned j_8k = j + 8*k;

        s += divide(mod_pow16(n-k, j_8k), DecimalPointPosition, j_8k);
    }

    return s;
}

static unsigned sigma_b(unsigned n, uint8_t j)
{
    unsigned s = 0;
    for (unsigned i = DecimalPointPosition; i > 0; i -= 4) {
        s += divide(1, i, j + 8*n);
        n++;
    }

    return s;
}

static unsigned sigma(unsigned n, uint8_t j)
{
    return sigma_a(n, j) + sigma_b(n, j);
}

static unsigned pi_n(unsigned n)
{
    return 4 * sigma(n, 1) - 2 * sigma(n, 4) - sigma(n, 5) - sigma(n, 6);
}

これを https://launchpad.net/gcc-arm-embedded で公開されている gcc-arm-none-eabi-5_4-2016q3-20160926-win32.exe の gcc と ld、objcopy を使用して機械語プログラムに変換しています。
変換した機械語プログラムは 456バイトのサイズとなりました。これをメモリに書き込む BASIC プログラムは単純には以下のような内容となりますが

10 POKE#700,240,181,1,37,133,176,1,145,53,75,1,33,238,24,118,0
11 POKE#710,48,0,0,240,137,248,4,33,4,0,48,0,0,240,132,248
12 POKE#720,7,0,5,33,48,0,0,240,127,248,100,0,231,27,127,0
13 POKE#730,63,26,6,33,48,0,0,240,119,248,63,26,58,11,1,45
14 POKE#740,4,209,144,33,1,155,63,13,9,1,95,84,144,35,1,153
15 POKE#750,27,1,235,24,202,84,160,35,1,53,91,0,157,66,211,209
16 POKE#760,192,38,220,59,0,34,156,70,182,0,0,42,4,209,144,35
17 POKE#770,1,153,27,1,201,92,39,224,7,33,211,23,11,64,155,24
18 POKE#780,219,16,235,26,179,66,1,221,192,35,155,0,0,33,147,66
19 POKE#790,25,221,10,39,1,152,196,24,18,72,32,24,0,120,120,67
20 POKE#7a0,65,24,17,72,131,66,3,220,144,32,0,1,36,24,33,112
21 POKE#7b0,96,70,9,18,3,144,3,152,0,40,2,221,3,152,1,56
22 POKE#7c0,248,231,1,59,227,231,1,53,144,35,1,152,27,1,211,24
23 POKE#7d0,48,49,1,50,193,84,178,66,199,209,0,32,5,176,240,189
24 POKE#7e0,255,255,255,127,255,8,0,0,255,2,0,0,16,181,0,35
25 POKE#7f0,1,36,32,66,3,208,91,24,147,66,0,211,155,26,64,8
26 POKE#800,4,208,73,0,138,66,244,216,137,26,242,231,24,0,16,189
27 POKE#810,0,35,64,0,91,0,144,66,1,211,128,26,1,51,1,57
28 POKE#820,0,41,246,209,24,0,112,71,240,181,38,74,135,176,131,24
29 POKE#830,219,0,4,147,92,24,0,35,3,144,5,145,70,30,1,147
30 POKE#840,115,28,40,208,3,155,0,32,159,27,16,37,1,44,24,208
31 POKE#850,165,66,1,211,45,27,251,231,1,35,2,147,1,35,31,66
32 POKE#860,5,208,2,153,34,0,40,0,255,247,192,255,2,144,127,8
33 POKE#870,6,208,41,0,40,0,34,0,255,247,184,255,5,0,237,231
34 POKE#880,2,152,34,0,20,33,255,247,195,255,1,155,1,62,27,24
35 POKE#890,1,147,8,60,212,231,5,154,4,155,148,70,99,68,28,0
36 POKE#8a0,0,38,20,37,48,52,106,0,41,0,162,26,1,32,255,247
37 POKE#8b0,175,255,4,61,54,24,0,45,245,209,1,155,240,24,7,176
38 POKE#8c0,240,189,192,70,255,255,255,31

約 1.7kB のサイズとなってしまい IchigoJam の BASIC のフリーエリアである 1024バイトを超えてしまうので一工夫しています。
次の Makefile と perl スクリプトにより機械語プログラムを BASIC のコメント文に変換しています。

Makefile(※行頭の空白をタブに変換してください)
pi.apo: pi.c
    arm-none-eabi-gcc -mcpu=cortex-m0 -mthumb -fpic -Os -g -c $< -o $(<:.c=.o)
    arm-none-eabi-ld --entry 0x700 -Ttext 0x700 $(<:.c=.o) -o $(@:.apo=.x) -Map $(@:.apo=.map)
    arm-none-eabi-objcopy -O srec $(@:.apo=.x) $(@:.apo=.mot)
    perl hex2apo.pl $(@:.apo=.mot) > $@
hex2apo.pl
#!/usr/bin/perl
$start = 0xffffffff;
$end = 0x00000000;
while(<>) {
    $S = substr($_, 0, 1);
    $type = substr($_, 1, 1);
    die if ($S != "S" || !($type >= 0 && $type <= 9));
    $size = hex(substr($_, 1 + 1, 2));
    for ($sum = 0, $i = 0; $i <= $size; $i++) {
        $sum += hex(substr($_, 1 + 1 + 2 * $i, 2));
    }
    die if (($sum & 0xff) != 0xff);
    if ($type >= 1 && $type <= 3) {
        $size = $size - 1 - ($type + 1);
        $address = hex(substr($_, 1 + 1 + 2, 2 * ($type + 1)));
        if ($start > $address) {
            $start = $address;
        }
        if ($end < $address + $size - 1) {
            $end = $address + $size - 1;
        }
        for ($i = 0; $i < $size; $i++) {
            $mem[$address + $i] = hex(substr($_, 1 + 1 + 2 + 2 * ($type + 1) + 2 * $i, 2));
        }
    }
}
$line = 10;
for ($i = $start; $i <= $end;) {
    printf("%d '", $line);
    for ($j = 0; $j < 15; $j++) {
        $n = 0;
        for ($k = 0; $k < 8; $k++) {
            $n = ($n << 8) & 0xff00;
            if ($k < 7) {
                $n += $mem[$i++];
            }
            &putc((($n >> ($k + 1)) & 0x7f));
            if ($i > $end) {
                if ($k < 7) {
                    &putc((($n << (6 - $k)) & 0x7f));
                }
                last;
            }
        }
        if ($i > $end) {
            last;
        }
    }
    printf("\n");
    $line += 1;
}

sub putc {
    my ($c) = @_;
    $c += 0x21;
    if ($c >= 0x7f) {
        $c += (0xa1 - 0x7f);
    }
    printf("%c", $c);
}
pi.apo(変換後のBASICのコメント文)
10 'サNA3M7」"inJQ*(}9\!'!!$、*ソ"%1A!」!!]1pa=!&1m!!(bツサS!=エ|ソ!`."e4!!"ウ\チ(エrォ7"7b;:"%$<@dB1+タJ12aT2y(w9ィ65#9%ュ|!HIN?h!Gア/」#5エ.W
11 '!!&AGgAD!ゥDQ/FYHウ"ァ>9}7ang<R(MVTB!\}"/W!!)SU1ェ[+4aT-A」Ei1'!(d、'b-%J);)(}i)!!*1QBY9)」ii(1"ゥ!#a,[$m!H0h=#\エシ」4Mag"m'」>9」」R
12 '!my6Fj&hォa%!Nd、^ツツツコツソ1!!`」A!!BV!)」3B"%$ォ7ヲ*;)"tngi!a4Aj!CiPG・3;シZヲ!!cス!2q!&y#Ac!キウ"q%ゥ"=a&@WeQ!Y2チ,J:6(yAヲ.y!*4O'!#9/A&
13 'ird」-nゥ=5U!Zy!b@.e%q*QRqsqa>:UXセカaEQ5m#D0qa~!+SC!+!0ツ「"ツ"E0ウa<AJ!+!#1$ツコ}`」q(Xo#m)a"B(ツコ、ツ」:y%ソ<-!SQbカJェ#ゥaj}r-ヲC(!!"9IF
14 '9..A"E"C.!E0ツ。「ツ#0G、a!|クォaT`!」0QサPY%Xツツツ0a

BASIC のコメント文である '(アポストロフィ)以降に 1文字7ビットの情報として、ASCII コードの '!'(#21)~'~'(#7E) を`0000000~`1011101、半角カナの '。'(#A1)~'ツ'(#C2)を`1011110~`1111111 の値を表すよう使用しており、8文字で 7バイトの情報を記録しています。これを以下の BASIC のプログラムで #00~#FF の値に変換しメモリに書き込んでいます。

2 Z=#C03:Y=#700
3 X=0
4 IFPEEK(Z)=39X=X-1:W=33-PEEK(Z-X):IFW<1V=V<<7-W+W/96*34:POKEY,V>>(X&7):Y=Y+(X&7<7):GOTO4
5 Z=Z+PEEK(Z-1)+4:IFPEEK(Z-3)+PEEK(Z-2):GOTO3

BASIC プログラムのフリーエリアをなるべく消費したいくないのでコードゴルフ的な書き方をしており少々見づらいものとなっています。
BASIC のプログラム中で使用している変数は V~Z の 5つであり、変数領域の #08F6~#08FF に該当する為この領域への書き込みはできませんが、今回のプログラムで使用する領域は #0700~#08C7 の範囲なので問題とはなりません。
コメント文の '(アポストロフィ)以降に '!'(#21) より小さい文字コードがあればその文字より行末までは無視されるので、

1 ' 999999 in pi

空白(#20)を使用すれば上記のようなコメントを書くことも可能となっています。

おわりに

私はまだ IchigoJam を触り始めて日も浅く、購入前は 1024バイトの BASIC のフリーエリアは小さすぎるだろうと考えていたのですが、工夫次第で意外と使い出があり驚いています。
今回はひとつの BASIC プログラムに詰め込むことに労力を払い、1024バイトのフリーエリア中、3/4 のサイズに収めることができました。まだメモリ領域には余裕があり、また複数のプログラムに分割することも考慮すればまだまだ色々なことができそうで夢が広がる感じですね。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした