ピタゴラ数の逆襲

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

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

ピタゴラ数とは a²+b²=c² となる自然数のこと。有名どころでは 3²+4²=5² とか 5²+12²=13²とか。

３つの数の最大公約数が1の場合、ちょっと特別で原始ピタゴラ数と呼ぶ。 6²+8²=10² とか 30²+40²=50² はどっちも 3²+4²=5² に何かをかけて出来るものだから、やや負けている。

とある整数 m と n が与えられたら、 (m²+n²) と (m²-n²) と 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 になるように何倍かすればよい。それを効率よく計算するプログラム。

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

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

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

☆☆☆

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)
}

コメントの投稿