Yuki Nakata's Blog

One color just reflects another


 フィボナッチ数列を計算してみよう1 C言語 の続きです。

前回はC言語だけでフィボナッチ数列を計算していました。

それだけだとunsigned long型の上限があり、大きな数字が計算できません。

今回はgmpライブラリを使ってフィボナッチ数列 1000001番目を計算してみます。

fibo.cというファイル名で保存してください。 

//gcc fibo.c -o fibo -lgmp
#include <stdio.h>
#include <string.h>
#include <stdint.h>
#include <gmp.h>

void omg_i_love_leonardo_of_pisa(uint32_t num, mpz_t * result) 
mpz_t retval, last, tmp;
mpz_init(retval);
  mpz_init(last);
mpz_init(tmp);

uint32_t i = 1;
if(num == 0) {
return;
}
  mpz_set_ui(retval, 1U);
mpz_set_ui(last, 0U);

for(; i < num; i++) {
mpz_set(tmp, retval);
mpz_add(retval, retval, last);

mpz_set(last, tmp);
}

mpz_set(*result, retval);
 }

int main() 
{
uint32_t num;
mpz_t fibo;
mpz_init(fibo);

//1000001
omg_i_love_leonardo_of_pisa(1000001, &fibo);
mpz_out_str(stdout, 10, fibo);
printf("\n");
return 1;
}

$ gcc fibo.c -o fibo -lgmp
$ ./fibo

すぐに計算結果が出てきます。
科学技術すげー。

http://www.wolframalpha.com/input/?i=fibonacci+1000001


実はgmpにはmpz_fib_uiというフィボナッチ関数が用意されています。
mpz_fib_ui関数を実行するだけフィボナッチ数列を求めることができます。

fibo3.cというファイル名で保存してください。
#include <stdio.h>
#include <string.h>
#include <stdint.h>
#include <gmp.h>

int main()
{
  int n = 1000000;
  mpz_t fm2;
  mpz_inits(fm2, NULL);
  mpz_fib_ui(fm2, n);
  gmp_printf("fib(%d) = %Zd\n", n, fm2);
  return 1;
}
 
$ gcc fibo3.c -o fibo3 -lgmp
$ ./fibo3
 
参考サイト
1,000,000th Fibonacci Number One-Liner in C

フィボナッチ数列を計算してみよう C言語

フィボナッチ数列とは直前の2つの数字の和になっている数列のことです。

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89

数列を見たら分かりますね。C言語でプログラミングしていきましょう。

3通りの方法で書いています。

fibo1, fibo2, fibo3 と進むにつれて処理が早くなっています。

fibo1が最も単純で何度も再起呼び出しをしていて遅い方法です。

fibo2はデータを保存して再起呼び出しの回数を少なくしています。

fibo3が一番シンプルです。for文で回しているだけです。
初めからfibo3で書けよという感じですね。

関数の再起呼び出しは難しい。プログラムを読むだけで頭を使います。

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

#define MAX_VAL 47
unsigned long fibo2Arr[MAX_VAL];
unsigned long fibo3Arr[MAX_VAL];

unsigned long fibo1(int n)
{
unsigned long f;
switch (n) {
case 1:
case 2:
f = 1L;
break;
default:
f = fibo1(n - 1) + fibo1(n - 2);
break;
}
return (f);
}

unsigned long fibo2(int n)
{
unsigned long f;
if (fibo2Arr[n]) {
f = fibo2Arr[n];
return (f);
}
switch (n) {
case 1:
case 2:
f = 1L;
break;
default:
f = fibo2(n  - 1) + fibo2(n - 2);
break;
}
fibo2Arr[n] = f;
return (f);
}

unsigned long fibo3(int n)
{
int i;
fibo3Arr[1] = 1;
fibo3Arr[2] = 1;
for (i = 3; i <= n; i++) {
fibo3Arr[i] = fibo3Arr[i - 1] + fibo3Arr[i - 2];
}
return (fibo3Arr[n]);
}

int main(void)
{
int n = MAX_VAL;
clock_t start, end;
unsigned long result1, result2, result3;
start = clock();
result1 = fibo1(n);
end = clock();
printf("fibo1 %lu %f\n", result1, (double)end - start);
start = clock();
result2 = fibo2(n);
end = clock();
printf("fibo2 %lu %f\n", result2, (double)end - start);
start = clock();
result3 = fibo3(n);
end = clock();
printf("fibo3 %lu %f\n", result3, (double)end - start);
return 0;


フィボナッチ数列の47まで計算して、それぞれ掛かった時間を計測しています。

私の環境ではこんな感じ。fibo1だけがやたら遅いです。

fibo1 2971215073 15680165.000000

fibo2 2971215073 2.000000

fibo3 2971215073 1.000000

フィボナッチ数列の47は
2971215073です。

たった47番目で 約30億になりますよ、と。

Wolfram Alphaで検索して合っているか確認してみましょう。

http://www.wolframalpha.com/input/?i=fibonacci+47 

同様に100番目でも検索してみてください。 

http://www.wolframalpha.com/input/?i=fibonacci+100 


今回はラズベリーパイにしゃべらせます。

AquesTalk Piを使うとラズベリーパイで簡単に音声を出力させることができます。

もちろんブラウザからリモートコントロールできるようにしますよ。

21世紀的ですね。

以下のリンクからAquesTalk Piをダウンロードしてください。

http://www.a-quest.com/products/aquestalkpi.html 

$ sudo chmod +x  AquesTalk

$ sudo cp AquesTalk /usr/bin/

私の環境ではウェブサーバーとしてnginxを利用しています。

/usr/share/nginx/www/ ここがウェブルートとなっています。

ここにindex.phpを置くとブラウザで表示されるわけですね。

ここにtalk サブディレクトリを作成してindex.phpファイルを置いています。

http://192.168.0.101/talk/

ここにアクセスするとラズベリーパイが話してくれるというわけ。

talkディレクトリの中にAquesTalkに付いてきた aq_dicを入れてください。

talk/aq_dic/aq_user.dic
talk/aq_dic/aq_dic.bin
talk/aq_dic/readme.txt
という構成になっています。

index.phpという名前で/usr/share/nginx/www/talk/に保存してください。


<?php
//AquesTalk Pi Web

if($_POST["talk"]) {
exec("AquesTalkPi '".$_POST["talk"]."' | aplay");
}
?>

<!DOCTYPE html>
<html lang="ja">
<head>
<meta charset="UTF-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>AquesTalk Pi Web</title>
<link rel="stylesheet"
  href="http://code.jquery.com/mobile/1.4.2/jquery.mobile-1.4.2.min.css" />
<script src="http://code.jquery.com/jquery-1.11.1.min.js"></script>
<script src="http://code.jquery.com/mobile/1.4.2/jquery.mobile-1.4.2.min.js">
</script>

<style type="text/css">  
<!-- 
    #header, #footer {
          background-color: #478384;
color: #fef4f4;
}
-->  
</style>
</head>
<body>
<div data-role="page">
  <div data-role="header" id="header">
    <h1>AquesTalk Pi Web</h1>
  </div>
  <div role="main" class="ui-content">
    

<form action="index.php" method="post">
  
  しゃべる内容:<br />
  <textarea name="talk" cols="30" rows="5"></textarea><br />
  <br />
  <input type="submit" value="しゃべる" />
</form>


  </div>
  <div data-role="footer" id="footer">
    <h3><a href="http://www.nakatayuki.com/" >nakata yuki</a></h3>
  </div>
</div>
</body>
</html>

 
参考にさせていただいたサイト
ラズベリーパイがWEB経由でしゃべるよ
 

radiko2


ラズベリーパイでラジオのradikoを再生しよう! の続きになります。

前回作成したindex.phpにcssデザインを施して見栄えよくしましょう。

スマホからアクセスするので、JQUERY MOBILEを使って
CSSデザインしてみました。

格好良くなりましたね。

index.phpで保存してください。


<?php
//Radiko Player
//http://www.nakatayuki.com


if(isset($_GET['id'])){

$id = $_GET['id'];
if($id=="stop"){
exec("killall mplayer");
}else if($id=="CCL"){
exec("killall mplayer");
exec("play_radiko.sh CCL", $e);
}else if($id=="802"){
exec("killall mplayer");
exec("play_radiko.sh 802", $e);
}else if($id=="ABC"){
exec("killall mplayer");
exec("play_radiko.sh ABC", $e);
}else if($id=="MBS"){
exec("killall mplayer");
exec("play_radiko.sh MBS", $e);
}else if($id=="OBC"){
exec("killall mplayer");
exec("play_radiko.sh OBC", $e);
}else if($id=="RN1"){
exec("killall mplayer");
exec("play_radiko.sh RN1", $e);
}else if($id=="RN2"){
exec("killall mplayer");
exec("play_radiko.sh RN2", $e);
}else if($id=="CRK"){
exec("killall mplayer");
exec("play_radiko.sh CRK", $e);
}else if($id=="KISSFMKOBE"){
exec("killall mplayer");
exec("play_radiko.sh KISSFMKOBE", $e);
}else if($id=="HOUSOU-DAIGAKU"){
exec("killall mplayer");
exec("play_radiko.sh HOUSOU-DAIGAKU", $e);
}
}
?>
<!DOCTYPE html>
<html lang="ja">
<head>
<meta charset="UTF-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Radiko Player</title>
<link rel="stylesheet"
  href="http://code.jquery.com/mobile/1.4.2/jquery.mobile-1.4.2.min.css" />
<script src="http://code.jquery.com/jquery-1.11.1.min.js"></script>
<script src="http://code.jquery.com/mobile/1.4.2/jquery.mobile-1.4.2.min.js">
</script>

<style type="text/css">  
<!-- 
    #header, #footer {
          background-color: #ffea00;
color: #493759;
}
-->  
</style>
</head>
<body>
<div data-role="page">
  <div data-role="header" id="header">
    <h1>Radiko Player</h1>
  </div>
  <div role="main" class="ui-content">
    <!--リストを準備-->

<ul data-role="listview" data-inset="true">
        <li><a href="index.php?id=stop">stop</a></li>
    </ul>

    <ul data-role="listview" data-inset="true">
        <li><a href="index.php?id=802">FM 802</a></li>
<li><a href="index.php?id=CCL">FM COCOLO</a></li>
<li><a href="index.php?id=ABC">ABC Radio</a></li>
<li><a href="index.php?id=MBS">MBS Radio</a></li>
<li><a href="index.php?id=OBC">OBC Radio Osaka</a></li>
<li><a href="index.php?id=RN1">Radio NIKKEI 1</a></li>
<li><a href="index.php?id=RN2">Radio NIKKEI 2</a></li>
<li><a href="index.php?id=CRK">CRK Radio Kansai</a></li>
<li><a href="index.php?id=KISSFMKOBE" >Kiss FM KOBE</a></li>
<li><a href="index.php?id=HOUSOU-DAIGAKU" >HOUSOU DAIGAKU</a></li>
    </ul>
  </div>
  <div data-role="footer" id="footer">
    <h3><a href="http://www.nakatayuki.com/" >nakata yuki</a></h3>
  </div>
</div>
</body>
</html>

 

arimalobby

有馬グランドホテルに行ってまいりました。

日帰りプランです。

アクセスは南方面の夙川から車で30分。意外と近いんですね。

http://www.arima-gh.jp/dayplan/ticket.html

入場券付館内共通券というものです。

平日大人で3650円。これで温泉に入れて2000円分までレストランで食事ができます。

2000円を超えた場合は超えた分だけ有料になります。

ホテルに着くと和服のお姉さんが出迎えてくれます。
arimagarden


ホテル、庭園、温泉どれも見事です。

贅沢にお金が掛かっています。

ホテルの温泉というのも良いですね。

平日の午前中だったので客はほとんどいません。温泉をほぼ独占できました。

露天風呂も一人。人がいない温泉は本当に素晴らしいと実感しました。

11時から温泉につかり、12時に食事。

お寿司の和食と洋食を選べます。
arimalunch

今回は洋食にしました。

贅沢な気分に浸れました。

日帰りで良いのでたまにはこういうのもいいですね。大満足です。



radiko
ラズベリーパイでFM COCOLOや FM 802のラジオを聞きたいですよね。

ラズベリーパイを使えば1日中radikoを聞くことができますよ。1か月聞きっぱなしでもコストは電気代の約50円だけ。

しかも、スマホやPCからラズベリーパイのラジオをリモートコントロールできます。

21世紀的ですね。さっそくやってみましょう。

. mplayerとrtmpdumpをインストールします。

$ sudo apt-get install mplayer

$ sudo apt-get install rtmpdump swftools libxml2-utils

mplayerで持っている音声ファイルmp3の音が出るか試してみてください。

$ mplayer 音声ファイル.mp3

音が流れたら成功です。


2. play_radiko.shをダウンロード

play_radiko.sh シェルスクリプトを下のリンクからダウンロードしてください。
 https://gist.github.com/ihsoy-s/5292735

play_radiko.sh を編集します。
作成するディレクトリを指定します。

私はoutdirを下の変更しました。
outdir="/home/pi/radio/"

これで /home/pi/radio/ディレクトリが作成されるようになります。

$ sudo chmod +x play_radiko.sh

$ sudo cp play_radiko.sh /usr/bin/play_radiko.sh

3. ラジオを再生してみよう

試しにFM COCOLOを流してみましょう。
$ play_radiko.sh CCL 

FM 802を再生してみましょう。
$ play_radiko.sh 802

CCLがFM COCOLOです。802がFM 802。
ラジオ曲の放送局IDは以下のリンクで確認できます。
http://www.dcc-jpl.com/foltia/wiki/radikomemo
 
4. ブラウザからリモートコントロールしよう。
ラズベリーパイからラジオを再生することができたら、ブラウザからリモートコントロールできるようにしましょう。

ここまでは簡単でした。
実はここからが本題で、難しいかもしれません。

今回ここは解説を省略します。

IPアドレスを固定してウェブサーバーを立ててください。

私は軽量のnginxとphpでウェブサーバーを立てています。定番のapacheでもいいです。

IPアドレスは192.168.0.101にしています。

IPアドレスは毎回起動時に変動すると面倒なので固定IPアドレスにします。

5. ウェブサーバーから音を鳴らせるようにしよう。
ブラウザからコマンドを実行させて音を鳴るようにする必要があります。

$ usermod -G www-data, audio www-data
要するにウェブサーバーをaudioグループに追加しています。

apacheやnginx, ともにデフォルトでwww-dataになっていますので、これで良いと思います。

更新の仕方が分からないので、ここで一度再起動しました。

6. index.phpを作成しよう。

ウェブサーバーを無事に立てれたとします。

以下のphpをindex.phpで保存して、ウェブルートに置いてください。


<?php
//radiko player
//http://www.nakatayuki.com/


if(isset($_GET['id'])){

$id = $_GET['id'];
if($id=="stop"){
exec("killall mplayer");
}else if($id=="CCL"){
exec("killall mplayer");
exec("play_radiko.sh CCL", $e);
}else if($id=="802"){
exec("killall mplayer");
exec("play_radiko.sh 802", $e);
}else if($id=="ABC"){
exec("killall mplayer");
exec("play_radiko.sh ABC", $e);
}else if($id=="MBS"){
exec("killall mplayer");
exec("play_radiko.sh MBS", $e);
}else if($id=="OBC"){
exec("killall mplayer");
exec("play_radiko.sh OBC", $e);
}else if($id=="RN1"){
exec("killall mplayer");
exec("play_radiko.sh RN1", $e);
}else if($id=="RN2"){
exec("killall mplayer");
exec("play_radiko.sh RN2", $e);
}else if($id=="CRK"){
exec("killall mplayer");
exec("play_radiko.sh CRK", $e);
}else if($id=="KISSFMKOBE"){
exec("killall mplayer");
exec("play_radiko.sh KISSFMKOBE", $e);
}else if($id=="HOUSOU-DAIGAKU"){
exec("killall mplayer");
exec("play_radiko.sh HOUSOU-DAIGAKU", $e);
}
}
?>
<html>
<head>
<meta name="viewport" content="width=device-width">
</head>
<body>
<p><a href="index.php?id=stop">stop</a></p>
<ul>
<li><a href="index.php?id=802">FM 802</a></li>
<li><a href="index.php?id=CCL">FM COCOLO</a></li>
<li><a href="index.php?id=ABC">ABC Radio</a></li>
<li><a href="index.php?id=MBS">MBS Radio</a></li>
<li><a href="index.php?id=OBC">OBC Radio Osaka</a></li>
<li><a href="index.php?id=RN1">Radio NIKKEI 1</a></li>
<li><a href="index.php?id=RN2" >Radio NIKKEI2</a></li>
<li><a href="index.php?id=CRK" >CRK Radio Kansai</a></li>
<li><a href="index.php?id=KISSFMKOBE" >Kiss FM KOBE</a></li>
<li><a href="index.php?id=HOUSOU-DAIGAKU" >HOUSOU DAIGAKU</a></li>
</ul>
</body>
</html>
 
 radiko player
パソコンからラズベリーパイにアクセスしたところ。

IPアドレスにアクセスします。私の環境ではhttp://192.168.0.101 です。
ラジオ曲をクリックするとラジオの再生が始まります。
STOPで停止ができます。

ラズベリーパイでラジオライフを楽しんでくださあい。でわ。 

前回のlinuxで科学技術計算をしようの続きです。

円周率を1億桁計算しました!からのプログラムです。

pai.c


/* Calculate pi based on Chudnovsky algorithm (& Binary Splitting method), using GMP */
/* [1] Computation of 2700 billion decimal digits of Pi using a Desktop Computer,
       Fabrice Bellard, Feb 11 2010 (4th revision),
       http://bellard.org/pi/pi2700e9/pipcrecord.pdf */
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>

mpz_t A, B, C, C3over24;

void computePQT(mpz_t P, mpz_t Q, mpz_t T, int n1, int n2) {
  int m;
  if (n1 + 1 == n2) {
    mpz_t a;
    mpz_init(a);
    /* P(n2-1, n2) = - (2 * n2 - 1) * (6 * n2 - 5) * (6 * n2 - 1) */
    mpz_set_ui(P, (2 * n2 - 1));
    mpz_mul_ui(P, P, (6 * n2 - 5));
    mpz_mul_ui(P, P, (6 * n2 - 1));
    mpz_neg(P, P);
    /* Q(n2-1, n2) = C^3 * n2^3 / 24 */
    mpz_mul_ui(Q, C3over24, n2);
    mpz_mul_ui(Q, Q, n2);
    mpz_mul_ui(Q, Q, n2);
    /* a_n2 = (A + B * n2) */
    mpz_mul_ui(a, B, n2);
    mpz_add(a, a, A);
    /* T(n2-1, n2) = a_n2 * P(n2-1, n2) */
    mpz_mul(T, a, P);
    mpz_clear(a);
  } else {
    mpz_t P1, Q1, T1, P2, Q2, T2;
    mpz_init(P1);
    mpz_init(Q1);
    mpz_init(T1);
    mpz_init(P2);
    mpz_init(Q2);
    mpz_init(T2);
    m = (n1 + n2) / 2;
    computePQT(P1, Q1, T1, n1, m);
    computePQT(P2, Q2, T2, m, n2);
    /* P = P1 * P2 */
    mpz_mul(P, P1, P2);
    /* Q = Q1 * Q2 */
    mpz_mul(Q, Q1, Q2);
    /* T = T1 * Q2 + P1 * T2 */
    mpz_mul(T, T1, Q2);
    mpz_mul(P1, P1, T2);
    mpz_add(T, T, P1);
    mpz_clear(P1);
    mpz_clear(Q1);
    mpz_clear(T1);
    mpz_clear(P2);
    mpz_clear(Q2);
    mpz_clear(T2);
  }
}

int main (int argc, char* argv[]) {

  clock_t start, end;
  FILE* fp;

  int digits = 10000;
  int prec = digits * log2(10);
  int n = digits / 14;

  start = clock();

  mpf_t pi, temp;
  mpz_t P, Q, T;

  /* Initialize GMP numbers */
  mpf_set_default_prec(prec);
  mpf_init(pi);
  mpf_init(temp);
  mpz_init(P);
  mpz_init(Q);
  mpz_init(T);
  mpz_init(A);
  mpz_init(B);
  mpz_init(C);
  mpz_init(C3over24);

  /* Assignment */
  mpz_set_str(A, "13591409", 10);
  mpz_set_str(B, "545140134", 10);
  mpz_set_str(C, "640320", 10);
  mpz_mul(C3over24, C, C);
  mpz_mul(C3over24, C3over24, C);
  mpz_div_ui(C3over24, C3over24, 24);

  computePQT(P, Q, T, 0, n);

  /* pi = C ^ (1 / 2)     */
  mpf_set_z(temp, C);
  mpf_sqrt(pi, temp);
  /*      * C             */
  mpf_mul(pi, pi, temp);
  /*      * Q             */
  mpf_set_z(temp, Q);
  mpf_mul(pi, pi, temp);
  /*      / (T + A * Q)   */
  mpz_mul(Q, A, Q);
  mpz_add(Q, Q, T);
  mpf_set_z(temp, Q);
  mpf_div(pi, pi, temp);
  /*      / 12            */
  mpf_div_ui(pi, pi, 12);

  /* mpf_out_str(stdout, 10, digits, pi); */
  fp = fopen("output.txt", "w");
  if (fp == NULL) {
    printf("couldn't open output.txt");
    exit(EXIT_FAILURE);
  }
  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);
  mpf_out_str(fp, 10, digits, pi);

  mpf_clear(pi);
  mpf_clear(temp);
  mpz_clear(P);
  mpz_clear(Q);
  mpz_clear(T);
  mpz_clear(A);
  mpz_clear(B);
  mpz_clear(C);
  mpz_clear(C3over24);

  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

  return 0;
}

桁を変更するには
int digits = 10000;

この箇所で変更してください。
1億桁にするには
int digits = 100000000;

10万桁からやってみましょう。
$ #10万桁
$ gcc pai.c -o pai -lgmp
$ ./pai
0.001 s
0.002 s

$#1億桁
$ gcc pai.c -o pai -lgmp
$ ./pai
147.206 s
180.558 s

1億桁がたった3分で計算できてしまいました。
MacBook Pro 2015 メモリ8Gです。
凄いパソコンを全く有効利用できていませんね。
プログラムの中身は後で読みたいと思います。ちゃんと読みますよ。 

↑このページのトップヘ