Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 03-23-2012

ピタゴラ数をみつけるプログラムを書いてから3年になる。当時は結果だけ出してプログラムは載せていない。ま、大したプログラムではないんだけど。言語はRだ。プログラムと呼べるようなものはRしか出来ない。

ピタゴラ数とは a2+b2=c2 となる自然数のこと。有名どころでは 32+42=52 とか 52+122=132とか。

3つの数の最大公約数が1の場合、ちょっと特別で原始ピタゴラ数と呼ぶ。 62+82=102 とか 302+402=502 はどっちも 32+42=52 に何かをかけて出来るものだから、やや負けている。

とある整数 m と n が与えられたら、 (m2+n2) と (m2-n2) と 2mn がピタゴラ数になっている。m=2, n=1 で (3,4,5) が出来るし、m=3, n=2 で (5,12,13) が出来る。

課題:ある数字を与えられた時にその数字を含むピタゴラ数を全て求めよう。

役に立つ事実:
1. 原始ピタゴラ数でないものは全て原始ピタゴラ数の何倍かになっている。
2. 原始ピタゴラ数を作るには m と n の一つが偶数でもう一つが奇数で、さらに互いに素でなくてはいけない。
3. 4の倍数でない偶数を含む原始ピタゴラ数は存在しない。

ある数、例えば 724 を含むピタゴラ数を全部みつけて!という問題だったら、724 の約数である 2, 4, 181, 362, 724 を含む原始ピタゴラ数をみつけて後で 724 になるように何倍かすればよい。それを効率よく計算するプログラム。

ちなみに素因数分解に役立つ、2つの数の最大公約数をみつけるには R では 一行でできる

gcd <- function(a,b) ifelse(b==0, a, Recall(b,a%%b))

この ピタゴラ数計算ページ は下記のプログラムを使っている。1億までの数ならサクサクっとピタゴラ数をみつけることができる。

☆☆☆

pt <- function(x){ gcd.simple <- function(a, b) ifelse(b==0, a, Recall(b, a%%b)) # 最大公約数を求める。 pt1 <- function(x){ # lapply で使えるように。 if(x%%2 == 0){ # 偶数の場合 n <- rev( (1:sqrt(x/2-1))[(x/2) %% (1:sqrt(x/2-1))==0] ) m <- (x/2) / n } if(x%%2 == 1){ # 奇数の場合 ## m^2 + n^2 = x # m と n を力技で探す。 m1 <- seq(1, sqrt((x-1)/2)) n1 <- sqrt( x-m1^2 ) ok <- n1==round(n1,0) ## m^2 - n^2 = x k2 <- rev( (1:sqrt(x-1))[ x %% (1:sqrt(x-1))==0 ] ) j2 <- x/k2 k2 <- k2[(k2+j2)%%2 == 0] j2 <- j2[(k2+j2)%%2 == 0] m <- c(m1[ok], (k2+j2)/2) n <- c(n1[ok], (j2-k2)/2) } q <- cbind(abs(m^2-n^2), 2*m*n, m^2+n^2, x) c( t(q[gcd.simple(abs(m^2-n^2), 2*m*n) == 1, ]) ) } # pt1 の終わり。 prim.fac <- unique( gcd.simple(x, 1:sqrt(x)) ) y <- unique( c(prim.fac, x/rev(prim.fac)) ) y <- y[y>2 & y%%4!=2] u <- matrix( unlist(lapply(y, pt1)), ncol=4, byrow=TRUE) z <- u * matrix(rep(x/u[, 4], each=4), ncol=4, byrow=TRUE) z[, 4] <- x / u[, 4] o <- z[,2] < z[,1] z[o, 1:2] <- z[o, 2:1] z <- z[order(z[, 1], z[, 2]), ] if(is.null(nrow(z))){ z <- t(matrix(z)) } z1 <- formatC(z[,1], format='f', digit=0) z2 <- formatC(z[,2], format='f', digit=0) z3 <- formatC(z[,3], format='f', digit=0) z4 <- formatC(z[,4], format='f', digit=0) data.frame(A=z1, B=z2, C=z3, gcd=z4) }

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 金曜日|どうしても分類不可コメント(0)

コメント

コメントの投稿

ブログ検索

カレンダー

JanFebMar AprMay JunJulAugSepOctNovDec 2005200620072008200920102011
- - - 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 -

最近のコメント

履歴

2018年07月 (1)

2018年06月 (1)

2018年05月 (1)

2018年04月 (2)

2018年03月 (1)

2018年01月 (1)

2017年12月 (1)

2017年11月 (1)

2017年10月 (1)

2017年09月 (1)

2017年08月 (1)

2017年07月 (1)

2017年06月 (1)

2017年05月 (1)

2017年04月 (1)

2017年03月 (1)

2017年02月 (1)

2017年01月 (1)

2016年12月 (1)

2016年11月 (1)

2016年09月 (1)

2016年08月 (1)

2016年07月 (1)

2016年05月 (1)

2016年04月 (1)

2016年03月 (1)

2016年02月 (1)

2016年01月 (1)

2015年12月 (1)

2015年07月 (1)

2015年06月 (2)

2015年05月 (2)

2015年04月 (3)

2015年03月 (2)

2015年02月 (1)

2015年01月 (3)

2014年12月 (1)

2014年11月 (1)

2014年10月 (2)

2014年09月 (2)

2014年08月 (2)

2014年07月 (1)

2014年06月 (1)

2014年05月 (2)

2014年04月 (1)

2014年03月 (2)

2014年02月 (1)

2014年01月 (2)

2013年12月 (1)

2013年11月 (3)

2013年10月 (2)

2013年09月 (1)

2013年08月 (4)

2013年07月 (1)

2013年06月 (2)

2013年05月 (2)

2013年04月 (3)

2013年03月 (1)

2013年02月 (1)

2013年01月 (4)

2012年12月 (1)

2012年11月 (3)

2012年10月 (1)

2012年09月 (1)

2012年08月 (3)

2012年07月 (3)

2012年06月 (2)

2012年05月 (6)

2012年04月 (2)

2012年03月 (8)

2012年02月 (2)

2012年01月 (1)

2011年12月 (6)

2011年11月 (5)

2011年10月 (4)

2011年09月 (6)

2011年08月 (9)

2011年07月 (5)

2011年06月 (5)

2011年05月 (5)

2011年04月 (6)

2011年03月 (17)

2011年02月 (6)

2011年01月 (10)

2010年12月 (10)

2010年11月 (4)

2010年10月 (6)

2010年09月 (5)

2010年08月 (11)

2010年07月 (8)

2010年06月 (8)

2010年05月 (3)

2010年04月 (8)

2010年03月 (11)

2010年02月 (4)

2010年01月 (8)

2009年12月 (6)

2009年11月 (6)

2009年10月 (6)

2009年09月 (7)

2009年08月 (6)

2009年07月 (10)

2009年06月 (10)

2009年05月 (10)

2009年04月 (6)

2009年03月 (7)

2009年02月 (9)

2009年01月 (12)

2008年12月 (6)

2008年11月 (10)

2008年10月 (8)

2008年09月 (9)

2008年08月 (12)

2008年07月 (8)

2008年06月 (12)

2008年05月 (12)

2008年04月 (12)

2008年03月 (11)

2008年02月 (10)

2008年01月 (10)

2007年12月 (12)

2007年11月 (14)

2007年10月 (13)

2007年09月 (11)

2007年08月 (16)

2007年07月 (10)

2007年06月 (10)

2007年05月 (6)

2007年04月 (10)

2007年03月 (13)

2007年02月 (10)

2007年01月 (8)

2006年12月 (13)

2006年11月 (15)

2006年10月 (9)

2006年09月 (8)

2006年08月 (18)

2006年07月 (14)

2006年06月 (16)

2006年05月 (23)

2006年04月 (20)

2006年03月 (12)

2006年02月 (14)

2006年01月 (21)

2005年12月 (20)

2005年11月 (17)

2005年10月 (18)

2005年09月 (16)

2005年08月 (10)

2005年07月 (6)