Saturday, October 15, 2011

2平方数の和(2)

nがあたえられたときに、x^2+y^2=nを解くにはどうしたらよいか。nが素数pの場合については以前書いた。では合成数の場合はどのようにすればよいか。

n=pqのとき、以前の方法でp=u^2+v^2, q = a^2+b^2ということが分かっているものとして、つぎの恒等式をつかってnを平方数の和に分ける。

n = pq = (u^2+v^2)(a^2+b^2) = (ua+vb)^2 + (va-ub)^2

これでnはua+vbとabs(va-ub)に分けられた(どちらも正数だけ採用することにする)。これを繰り返してゆけばnが2個以上の素数で分解出来る場合も対処できる。

元々、pが4で割ったとき1余る場合、x^2+y^2=pの解はあり、pの余りが3の場合に解はない。nにそのようなpが含まれている場合、pの個数が奇数のとき、x^2+y^2=nに解がなくなる。偶数の場合には解ありである。例えば、n=9の時、p=3だが3が2個あるのでx^2+y^2=1の解 (1,0)を3倍した(3, 0)が解となる。また、n=3やn=27には解がない。

コードにしたのが以下。素因数分解は状況によって違うものを使いたい(試し割りのものと、フルイをつかったもの)ので、初期化で渡すようにした。以下ではfacs()がそれで、これに期待しているのは、n=p^r * q^sのとき、[p,r],[q,s]を渡してくれること。fscsの中で使っているfactors()が素因数分解のコア部分。

名前をprocess_primeに変えたが、これが以前書いた、素数のx^2+y^2=pを解くところ。comb()が上に書いた恒等式で展開するもの。

class square_sum(object):
  def __init__(self, func):
    self.func = func
  def square_mod(self, p):
    from random import randint
    while True:
      a = randint(2, p - 1)
      b = pow(a, (p - 1) / 4) % p
      if (b*b) % p == p - 1:
        return b

  def process_prime(self, p):
    if p == 2: return [1, 1]
    if p % 4 != 1: return [0, 0]

    a = self.square_mod(p)
    b = 1
    k = (a * a + b * b) / p

    while k > 1:
      u, v = a % k, b % k   
      # abs(u) <= k/2
      # abs(v) <= k/2
      if u > k / 2: u -= k
      if v > k / 2: v -= k
 
      # u * u + v * v = 0 (mod k)
      q = (u * u + v * v) / k

      # (u * u + v * v) * (a * a + b * b) = (u * a + v * b)^2 + (v * a - u * b)^2 = (q * k) * (k * p)
      c, d = (u * a + v * b) / k, (v * a - u * b) / k

      # next turn
      a, b = c, d
      k = q
    return [abs(a), abs(b)] 
  def comb(self, X, Y):
    u, v = X
    a, b = Y

    res = set()
    r1 = [u * a + v * b, abs(v * a - u * b)]
    r1.sort()
    res.add(tuple(r1))

    r2 = [u * b + v * a, abs(v * b - u * a)]
    r2.sort()
    res.add(tuple(r2))

    return res

  def times(self, M, seq):
    return map(lambda x: x*M, seq)

  def solv(self, n):
    fac = self.func(n)

    # p%4==3 prime
    M = 1
    for p,r in fac:
      if p % 4 == 3:
        if r % 2 == 0:
          M *= p**(r/2)
        else:
          return []

    # p==2 or p%4==1 prime
    seq = []
    for p,r in fac:
      if p == 2 or p % 4 == 1:
        sq = self.process_prime(p)
        seq += [sq]*r

    if len(seq) ==  0:
      return [self.times(M, [0,1])]

    res = set()
    res.add(tuple(seq[0]))
    for k in xrange(1, len(seq)):
      new = set()
      for e in res:
         new |= self.comb(e, seq[k])
      res = new

    return map(lambda x: self.times(M,x), res) 
# if n=p1^r1 x ... x pm^rm, returns [p1,r1]...[pm,rm]
def facs(n):
  seq = factors(n)
  s = set(seq)
  r = []
  for p in s:
    r.append([p, seq.count(p)])
  return r
使い方は、
sq = square_sum(facs)
sq.solv(1234)
のような感じ。
比べてみると、しらみつぶしでx^2+y^2=nを解くよりはかなりはやい。




Friday, October 14, 2011

禅について



私と禅について書いてみる。始まりは15年ほど前。禅というものにはちょっと興味があって鈴木大拙の岩波新書の禅の本などを読んでみたが何のことか全然わからなかった。茶道にも興味があったので表千家の先生について習った。茶道の掛け軸と言えば、禅僧のものが定番なのでここでも禅にたいする疑問がわいた。禅とは何なのか?しかしどうやら本を読んでも分からないだろうなと思っていた。

松門寺

松門寺に行ったのは1998年のこと。インターネットで見つけた。松門寺は曹洞宗のお寺。ここの坐禅会は月一回、2,30人の人が集まり、坐禅の時間、正法眼蔵の講話の時間があった。お寺は住職である老師が一人で取り仕切っている。坐禅堂はモダンな作りで新しく清潔だった。坐禅の時間には、独参(どくさん)がある。独参では老師の部屋に一人一人入り対座し、坐禅についての質問や見解を述べる。坐禅がはじまりしばらくすると、独参開始の合図が鐘の音で知らされ、坐禅をやめて独参の部屋の前に並ぶ。

坐禅とは仏教である。インドからやってきた達磨が唐の時代の中国に伝えたということになっている。老師の言い方では、坐禅とは「自分自身になる」ことであった。これは指導者の説明の仕方により、いろいろな言い方がある。ある人は「即心是仏」と言った。言葉で何と言っても当たらない。言葉にはとらわれないようにとのこと。

「自分自身」、それ以上のことはない。また、従うべき指針や善悪の基準があるわけではない。つまり禅には教義がないのである。「自分自身になる」にはどうすればよいのか?「何もしなくてよい」と老師は言う。何もしなくてよい?一体どういうことなんだ?「坐禅のときに、集中ができません」、それが悪いことだと誰が決めたのか?集中すべきと思ったのは自分。そこから離れてごらん。集中できない自分を自分として受け入れる。それが坐禅の一部。という調子。一筋縄ではないのです。

老師は東工大の化学を出た人。答えが分かるのは答えが分かったとき、という事の性質上、先生が信用出来るかどうかが鍵となる。こちらも理系だからということもあったのかもしれないが、この人が有ると言ったものは本当に有るのだろうし、出来ると言ったことは出来るのだろうと思った。

それにしても、ちょっとした疑問を解消しようとしているのにしては、何だか遠くにきてしまったなあというのが私の気持ちだった。


私の独参

私の独参は、最初のうち質問。松門寺に行けば毎回行っていた。3回くらい行ったところで、もう質問することがなくなった。「公案をやってみるか」と言われた。公案とは坐禅の先生から弟子への質問である。禅宗でも曹洞宗は公案をやらないとされていて、臨済宗では公案をやる。しかし曹洞宗も臨済宗も本は同じ、違いはない、曹洞宗ではあるが、うちでは公案をやりますとのこと。次に来るときにはこれに答えてくださいという意味のことを言われた。これができたら「悟ったことになるのかな」などと思った。一ヶ月後、何と答えたのか忘れたが、それは違うね、という感じのことを言われた。公案はテストであり、悟っているかどうかがたちどころにわかるのが公案だとのこと。何なのだろうか。その答えは、本の後ろの方に「解答」が書いてあるような類いのものでない。

「誰にでもでき」て、それをするのには「何もしなくてよく」、「いつ分かってもおかしくない」と言われればいてもたってもいられないだろう。絶対に分かってやる、と思っていた。

そういう独参が2年くらい続いた。一ヶ月ごとに答えを持っていって、否定されるまでに10秒とかからない。方向は良くなってきているけれどね、など。その後ちょっと雑談するなどして退席である。一体オレはなにをやっているのかなあというのが正直なところだった。老師が言うところによると「100メートルの棒の上を歩いてきた。その先っぽにいる。お前はそこから一歩を進めないでいる」とか、「柿は熟すれば自然に木から落ちる。それなのに、お前は自分の力で木にしがみついているだけだ。なんで自分の力を使おうとするのか」など。これらの言い回しは、古来の公案にある言い方である。

さて、およそのことは分かった。しかし一体なんと言えば通してくれるのだろうか。今にして思えば、そう思っていること自体に問題があったのだろうなあ。あるとき、「一人通った人がいる」、と言うことだった。うらやましいなあと思った。「他にも通りそうな人がここには何人かいる。」それ、私のことなんだろうなあ。

最後には、この他には何もない、棒の先にいるなら先にすすもう。通ろうという気持ちで独参に向かった。いくつかのやりとりのあと、老師は言った「これからもしっかり坐ってください。」

これが2001年の夏のこと。この最後の独参のあと、老師のところには行っていない。

肝心なところは、禅は紙に書いてある知識ではなく、自分自身との対話である点。本をいくら読んでも分からなかっただろうと思う。また、チェックしてくれる指導者が必要である点です。

Tuesday, August 30, 2011

食塩水の問題

子供の算数を見ていて、ちょっと困った問題。
6%の食塩水100gに、15%の食塩水を混ぜたところ、濃さが10%になりました。15%の食塩水を何g混ぜましたか?

中学生以上なら、
6+0.15*x = 0.1*(100+x)

として、
x=80

とします。私は小学生のときからこうやっていましたが。

小学生の解答は以下のように考えるようです。
薄い方の食塩水の水と、食塩を取り替えて、濃さが10%になるようにするには?
逆に濃い方の食塩と、水を取り替えて、濃さが10%にするには?

つまり
(0.1-0.06)*100=4


(0.15-0.1)*x = 0.05*x

が等しいので
x =80

と考えるのです。

少し前に、方程式をたてないのもいいなあと思ったばかりですが、今回のこれはやっぱり嫌だなあ。問題によってやり方を変える感じが気にいらない。これだったら最初から方程式を教えた方がいいんじゃないかと思います。

Friday, May 06, 2011

ババ抜きのスタートの枚数

トランプのババ抜きは、参加者にカードを配り、同じ数字のカードを捨てることからはじまります。
カードが切れていないと、同じカードが揃っていて有利なスタートとなります。カードがちゃんとシャッフルされているとするなら、何枚からスタートとなるのが妥当なのでしょう?

確率の計算をしてもいいのですが、場合分けが面倒くさかったのでシミュレーションでやってみました。そう、計算をしても答え合わせのためにどのみちシミュレーションはやるんです。

左側が配られた枚数。そのときの右はスタートの枚数です。100回づつ実験しました。
(配った枚数は人数により、ありえない数があるのでところどころ飛んでます。3人でやるなら17枚か18枚、4人なら13枚、5人なら10枚か11枚etc....)

2 1.86
3 2.62
4 3.22
5 3.56
6 4.24
7 4.38
8 4.78
9 5.34
10 5.46
11 5.32
13 5.84
17 6.12
18 6.22
26 5.62

これを見ると、大体5、6枚で安定してるようですね。
プログラムもつけとこ。


from random import randint

numcards = 52
kinds = 13

def f(n):
a = [0 for i in xrange(0, kinds)]

while n > 0:
c = randint(0, kinds-1)
if a[c] < 4:
a[c] += 1
n -= 1

r = 0
for x in a: r += x % 2
return r

def gen_num():
a = set([])
for p in xrange(2, 25):
if numcards % p != 0:
a.add(numcards/p + 1)
a.add(numcards/p)
b = list(a)
b.sort()
return b

def main():
N = 100
for n in gen_num():
r = 0.0
for i in xrange(N):
r += f(n)
print n, r/N

main()

Sunday, May 01, 2011

pythonのdictionaryのsetdefault, pop,getなど

自分用のメモ。

dictionaryでその要素がなければ何とかで、ある場合は何とかという処理をif文で書くことがあるけれど、if "a" in dのようなことを書かずに済む方法。

>>> d = {}
>>> d.setdefault("a", 1)
1
>>> d
{'a': 1}

まだ要素がない場合、1が登録される。
既にある場合は指定とちがって既にある値が返ってくるので分かるという寸法。

>>> d.setdefault("a", 2)
1
>>> d
{'a': 1}
>>>

popにも似た機能がある。

>>> d.pop("c", 1)
1

何もないときの振る舞いがこれ。

返り値を指定しない場合は怒られてしまう。

>>> d.pop("c")
Traceback (most recent call last):
File "", line 1, in
KeyError: 'c'


getにも似たことが言える。

>>> d.get("c")
>>> d.get("c", 1)
1

Saturday, February 05, 2011

スーパータッチボール

スーパータッチボールというものを子供が学校でやっていて、大会があるというので観戦にゆきました。スーパータッチがどんなゲームかというと、バスケットボールのコートで行うもので、ラグビーと同じボールを使い、タックルの代わりに相手にタッチするものです。タッチされた相手はボールを離さなければならないのです。タッチラグビーというものを見たことがあるのですが、何が違うのかはわかりません。横浜が発祥の地だそうです。はじめて見ましたがなかなか面白かった。