円周率10億桁を計算する #1


すずきひろのぶ氏が指摘している通り、gmp の precision のサイズの上限により、
円周率計算は6.4億桁が上限となる。
m(_ _)m 大変申し訳ございません。嘘つきました。勉強不足。
現在、gmp の precision のサイズは int (2^31-1=2,147,483,647)ではありません。従って、(2^31-1)/log210 ≑ 646,456,992 が円周率計算の上限ではありません。
nitobe@debian64:~/pi$ cat size.c
#include<stdio.h>
#include<gmp.h>

main()
{
printf("int: %d¥n", sizeof(int));
printf("long int: %d¥n", sizeof(long int));
printf("long long int: %d¥n", sizeof(long long int));
printf("mp_bitcnt_t: %d¥n", sizeof(mp_bitcnt_t));
}
nitobe@debian64:~/pi$ make size
gcc -static -O3 -m64 -mtune=amdfam10 -I/usr/local/gmp-5.0.5/include size.c -o size /usr/local/gmp-5.0.5/lib/libgmp.a
nitobe@debian64:~/pi$ ./size
int: 4
long int: 8
long long int: 8
mp_bitcnt_t: 8
nitobe@debian64:~/pi$ 
理論上の上限は、2^63-1 = 2進 9,223,372,036,854,775,807桁、(2^63-1)/log210 = 10進 2,776,511,644,261,678,565桁 ということになる。2.78E(エクサ10^18)Bit = 348P(ペタ10^15)Byteのメモリを実装したコンピュータなどそうそうあるもんじゃない。当面 gmp の上限は見られないということだ。

鋭意計算中。

12/15 23:40 追記:
pi2.c (diff)
nitobe@debian64:~/pi$ diff pi.c pi2.c
10,11c10,11
< #define LOG_TEN_TWO  3.32192809488736234789
< #define bprec(n) (int)(((n+10)*LOG_TEN_TWO)+2)
---
> #define LOG_TEN_TWO  3.3219280948873623478703194294894
> #define bprec(n) (long int)(((n+10)*LOG_TEN_TWO)+2)
16a17
>   mpz_t temp, temp2;/* add itchyny */
18c19
<   dprec = 1000000L;/* decimal precision */
---
>   dprec = 1000000000L;/* decimal precision */
21c22
<   loopmax = 21;
---
>   loopmax = 31;
22a24
>   fprintf(stderr, "%ld %d initial:", prec, loopmax); /* Add nitobe */
34,35c36,38
< 
<    mpf_set_ui (A, 1);
---
>   mpz_init (temp);/* add itchyny */
>   mpz_init (temp2);/* add itchyny */
>   mpf_set_ui (A, 1);
41a45
>   mpz_set_ui (temp2, 2);/* add itchyny */
43a48
>       fprintf(stderr, " %d:", k);
53c58,60
<       mpf_pow_ui (t2, c2, k);/* 2^k */
---
> //      mpf_pow_ui (t2, c2, k);/* 2^k */
>       mpz_pow_ui (temp, temp2, k);/* 2^k */
>       mpf_set_z (t2, temp);/* add itchyny */
59c66,67
<   mpf_pow_ui (t2, t1, 2);
---
> //  mpf_pow_ui (t2, t1, 2);
>   mpf_mul (t2, t1, t1);/* mod itchyny */
61a70
>   fprintf(stderr, " conversion:", prec, loopmax); /* Add nitobe */
64c73
<   exit(0);
---
>   return 0;/* mod nitobe */
nitobe@debian64:~/pi$ 
makefile
nitobe@debian64:~/pi$ cat makefile
CC=gcc
OPT= -static -O3 -m64 -mtune=amdfam10
GMP=gmp-5.0.5
GMP_INC=-I/usr/local/$(GMP)/include
GMP_LIB=/usr/local/$(GMP)/lib/libgmp.a
SRC= pi2.c
EXE= pi

.PHONY: clean

all: $(EXE)

pi: $(SRC)
        $(CC) $(OPT) $(GMP_INC) $(SRC) -o $(EXE) $(GMP_LIB)

size: size.c
        $(CC) $(OPT) $(GMP_INC) size.c -o size $(GMP_LIB)

clean:
rm -f pi
rm -f size
nitobe@debian64:~/pi$ 

nitobe@debian64:~/pi$ time ./pi >pi1G.txt
3321928130 31 initial: 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: conversion:
real    196m38.626s
user    192m55.963s
sys     3m43.246s
nitobe@debian64:~/pi$ time ./format.sed pi1G.txt >pi1g.txt

real2m54.283s
user2m39.994s
sys0m3.148s
nitobe@debian64:~/pi$
注:実行時間は全く当てにならない。このPC上では、seti@home Link on BOINCLink 及び、World Community GridLink on BOINC がリソースの75%を占有している。

pi1g.txt
       1 3.
       2 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
       3 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196

          :
          :
10000000 4219776753871319682188195635848934815504410194647387557034502943416861599324354199731814355060392734
10000001 6434543524276655356743570219396394581990548327874671398682093196353628204612755715171395115275045519
10000002 64388312
さて、この値が正しいかどうか?それが問題だ。つづく

 

— posted by nitobe at 09:35 pm   commentComment [0] 

円周率646456700桁弱を計算する #1

gmpの限界を見ておこうと思う。環境は以下の通り。
gmp-5.0.5 /gcc-4.4.5 /debian64-2.6.32-5-amd64 /VirtualBox 4.2.4 /Windows 7 /DELL XPS8300 /3.4G Intel Core i7-2600。

GNU/Linux上で円周率の計算をおこなうLink 」すずきひろのぶ氏の百万桁計算のソースを使用する。変更点は以下の通り。itchynyさんの修正も適用している。すずきひろのぶ氏が指摘している通り、gmp の precision のサイズの上限により、円周率計算は6.4億桁が上限となる。12/15追記:大嘘です円周率10億桁を計算する #1Link
nitobe@debian64:~/pi$ diff pi.c pi2.c
16a17
>   mpz_t temp, temp2;/* add itchyny */
18c19
<   dprec = 1000000L;/* decimal precision */
---
>   dprec = 800000000L;/* decimal precision */
21c22
<   loopmax = 21;
---
>   loopmax = 31;
22a24
>   fprintf(stderr, "%d initial:", loopmax); /* Add nitobe */
34,35c36,38
< 
<    mpf_set_ui (A, 1);
---
>   mpz_init (temp);/* add itchyny */
>   mpz_init (temp2);/* add itchyny */
>   mpf_set_ui (A, 1);
41a45
>   mpz_set_ui (temp2, 2);/* add itchyny */
43a48
>       fprintf(stderr, " %d:", k); /* Add nitobe */
53c58,60
<       mpf_pow_ui (t2, c2, k);/* 2^k */
---
> //      mpf_pow_ui (t2, c2, k);/* 2^k */
>       mpz_pow_ui (temp, temp2, k);/* 2^k */
>       mpf_set_z (t2, temp);/* add itchyny */
59c66,67
<   mpf_pow_ui (t2, t1, 2);
---
> //  mpf_pow_ui (t2, t1, 2);
>   mpf_mul (t2, t1, t1);/* mod itchyny */
64c72
<   exit(0);
---
>   return 0;/* mod nitobe */
nitobe@debian64:~/pi$ 
dprecで求める桁数を指定している。 loopmaxは演算繰り返し回数だが、求め得る桁数と、以下の関係がある。
loopmaxcorrect digits (about)
9<500
10<800
11<1500
12<2900
13<5700
14<11300
15<22500
16<44900
17<89600
18<179000
19<357800
20<715500
21<1430800
22<2861500
23<5722800
24<11445400
25<22890600
26<45781000
27<91561900
28<183123600
29<366247100
30<732494000

makefile
nitobe@debian64:~/pi$ cat makefile
CC=gcc
CFLAG=-static -O3
GMP=gmp-5.0.5
GMP_INC=-I/usr/local/$(GMP)/include
GMP_LIB=/usr/local/$(GMP)/lib/libgmp.a

.PHONY: clean

all: pi

pi:
        $(CC) $(CFLAG) $(GMP_INC) pi2.c -o pi $(GMP_LIB)

clean:
rm pi
nitobe@debian64:~/pi$

make して実行
nitobe@debian64:~/pi$ make
gcc -static -O3 -I/usr/local/gmp-5.0.5/include pi2.c -o pi /usr/local/gmp-5.0.5/lib/libgmp.a
nitobe@debian64:~/pi$ time ./pi > pi800M.txt
31 initial: 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31:
real68m26.877s
user66m50.763s
sys1m36.338s
nitobe@debian64:~/pi$

単一行巨大ファイルは扱い難いので、百桁/行に整形する。
format.sed
nitobe@debian64:~/pi$ cat format.sed
#!/bin/sed -f
/^$/d
s/0.3/3.¥n/
s/[0-9]¥{100¥}/&¥n/g
s/e1/¥n/
nitobe@debian64:~/pi$ time ./format.sed pi800M.txt >pi800m.txt

real1m1.508s
user1m0.404s
sys0m1.108s
nitobe@debian64:~/pi$ ls -al
合計 1271884
drwxr-xr-x  2 nitobe nitobe      4096 2012-12-13 21:15 .
drwxr-xr-x 37 nitobe nitobe      4096 2012-12-13 19:50 ..
-rwxr-xr-x  1 nitobe nitobe        62 2012-12-04 21:16 format.sed
-rw-r--r--  1 nitobe nitobe       208 2012-12-13 19:49 makefile
-rwxr-xr-x  1 nitobe nitobe    854298 2012-12-13 19:49 pi
-rw-r--r--  1 nitobe nitobe      1851 2012-12-13 17:05 pi.c
-rwxr-xr-x  1 nitobe nitobe    854298 2012-12-13 19:40 pi2
-rw-r--r--  1 nitobe nitobe      2264 2012-12-13 19:48 pi2.c
-rw-r--r--  1 nitobe nitobe 646456999 2012-12-13 20:58 pi800M.txt
-rw-r--r--  1 nitobe nitobe 652921567 2012-12-13 21:08 pi800m.txt
nitobe@debian64:~/pi

8億桁を指定しているが、6億4千万桁しかない。フォーマット変換後、改行「¥n」が6百万個増えている。百桁/行フォーマットなのでエディタで開ける。重いけど。これは vi で見てみた。
      1 3.
      2 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
      3 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196
    :
    :
6464569 9336326991629433891751838174109742989418332960270103331518724202832616454999878479262174537890538106
6464570 6455256718445892438033533034196466727204878408211469825653223223178441813002292106453186861390731397
6464571 2309985607826687727041129258325166030790611057861923604585439381176655074474471470989090334958

さて、この値が正しいかどうか?それが問題だ。つづく

 

— posted by nitobe at 05:34 pm   commentComment [0] 

2013カレンダーdebug


本宅の上の方に表示されるカレンダースクリプトの2013年分のテーブル作成とデバグ。
2013カレンダーdebug Link
誤謬の指摘を乞う。
2013カレンダー・ソースLink
テーブル作成スクリプト・ソースLink

— posted by nitobe at 11:19 am   commentComment [0] 

深夜食堂 10



著者: 安倍 夜郎
出版: 小学館
ISBN: 4091848257
定 価: ¥ 780
出版日: 2012-11-30

 

— posted by nitobe at 08:48 pm   commentComment [0] 

Repair Thinkpad x31

というわけで、Thikpad x31 TYPE 2672-CBJ を修理する。私の記憶が確かならば、2004年前半に購入しているので、8年選手。んが、まだまだ現役だ。

症状: 電源インジケータ点灯、CPU FAN 回転、bios - boot loader 起動せず。

方針: X31 の持病、 タンタルチップコンデンサ C683 の容量抜けと推察される。C683 交換。

施工:
DSC_0087
分解・施工前

DSC_0090
施工前

DSC_0101
施工後

DSC_0095
施工後


使用工具: +1ドライバ、ラジオペンチ、ANTEX 18W半田ごて、ソルダーウィック、スパークルハンダ、むしめがね(老眼)

施工時間: 2時間

結果: 起動確認 復活!

交換部品: タンタルチップコンデンサ 220uF4V KEMET T520D227M004ATE012
購入:    マルツパーツ館WebShopLink @370-

修理不能を想定して、タンタルチップコンデンサの発注と同時に、YahooオークションでThinkpad X31 2672-B2J (¥6,800-)の入札・落札を行っていた。最悪、HDDだけ移植すればよいという寸法だ。直っちゃったので、X31を2台所有することとなった。また、コンデンサの送料を無料にするため、radiation-watch Link の、ポケットガイガーType5 EMBEDDED(¥6,450-)Link も同時発注した。馬鹿者、俺。嫁にはご内密に。


— posted by nitobe at 11:30 pm   commentComment [1]  pingTrackBack [0]

T: Y: ALL: Online:
ThemeSwitch
  • Basic
Created in 1.2136 sec.
prev
2024.12
next
1 2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30 31        
 
strawberry-linux geigercounter Ver.2
Sibasaki, Cyofu City, Tokyo, JAPAN
blogBar