2012年12月25日火曜日

NZMATH 1.2.0

The NZMATH development group is pleased to announce the availability of NZMATH 1.2.0. Some of the highlights of new feature of NZMATH version 1.2 are here:

  • implementation of AKS (Cyclotomic congruence test) algorithm.
  • new curves added for ECM algorithms.
  • improvement multiple polynomial quadratic sieve.

For more details, Please read CHANGES.txt .

2012年6月27日水曜日

nzmath.arith1.log is faster than math.log

nzmath.arith1.log(n) is faster than int(math.log(n)) if base is 2:
$ python -m timeit -s 'import math' '[int(math.log(i, 2)) for i in range(1, 10000000, 1000)]'
10 loops, best of 3: 91.6 msec per loop
$ python -m timeit -s 'import nzmath.arith1' '[nzmath.arith1.log(i, 2) for i in range(1, 10000000, 1000)]'
10 loops, best of 3: 32.5 msec per loop
It's true also for base 10:
$ python -m timeit -s 'import math' '[int(math.log(i, 10)) for i in range(1, 1000000000, 100000)]'
10 loops, best of 3: 93.1 msec per loop
$ python -m timeit -s 'import nzmath.arith1' '[nzmath.arith1.log(i, 10) for i in range(1, 1000000000, 100000)]'
100 loops, best of 3: 16.5 msec per loop
It isn't fair if I stop here.
$ python -m timeit -s 'import math' '[int(math.log(i, 3)) for i in range(1, 1000000000, 100000)]'
10 loops, best of 3: 93.2 msec per loop
$ python -m timeit -s 'import nzmath.arith1' '[nzmath.arith1.log(i, 3) for i in range(1, 1000000000, 100000)]'
10 loops, best of 3: 121 msec per loop

2012年5月6日日曜日

Cornacchia

Iwao Kimura at Blogger の Cornacchia-Smithのアルゴリズム:mod 4で1の素数は2平方数の和という記事を読みました。 Cornacchia という名前に聞き覚えが(読み方は多分イタリア人なのでコルナッキア)あります。 そう確かあれは cubic_root モジュール。

あまりこのモジュールを使ったことがある、という人も多くないと思いますので、その紹介を簡単に。 平方剰余というのは有名ですが、3乗剰余というものもあるのです。 cubic_root モジュールはこの計算をします。 それ以外に Eisenstein 整域 \(\mathbb{Z}[\omega]\) での有理素数の因数分解や素数を法とする3乗根の計算なども提供しています。 この中に何故か、というかまあ実装の都合上ここに、cornacchia という関数があります。

cornacchia(d, p) で \(x^2 + d y^2 = p\) の解(の一つ)をタプルとして返します。 ということで早速使ってみましょう。

>>> import nzmath.cubic_root
>>> nzmath.cubic_root.cornacchia(1, 37)
(6, 1)

つまり \(6^2+1\times 1^2 = 37\) ということですね。

実装は H.Cohen の A Course in Computational Algebraic Number Theory を参考にしています。

2011年12月29日木曜日

tnt domain

I apologize you about this late announcement about domain change.

The domain of Tokyo Metropolitan University changed from metro-u.ac.jp to tmu.ac.jp. Accordingly tnt domain is moved from tnt.math.metro-u.ac.jp to tnt.math.se.tmu.ac.jp. So the official homepage of NZMATH is now http://tnt.math.se.tmu.ac.jp/nzmath/.

All the old links on this blog will remain unedited.

2010年6月7日月曜日

NZMATH 1.0 Release Announcement

Dear number theorists,

We are pleased to announce the release of NZMATH 1.0.

NZMATH is a Python based system for number theory developed for six years. Various modules useful for number theory are included in NZMATH package. This time, we release its first stable version 1.0 officially. It can be used easily even for programming beginners. Please try it!

NZMATH version 1.0 introduces many new features and many improvements of functionality in the earlier versions. Some of the highlights:

  • support of algebraic number fields.
  • TeX (PDF) documentations.
  • elliptic curve primality proving (ECPP).
  • refinement of polynomial factorization over the rational integers.
  • configuration for external modules.

For a complete list of features and more informations, access to
http://tnt.math.metro-u.ac.jp/nzmath/

Download

You can download NZMATH from
http://downloads.sourceforge.net/nzmath/
http://tnt.math.metro-u.ac.jp/nzmath/

Mailing List

There is a mailing list to discuss about NZMATH among the users. To join the mailing list, send a mail to
nzmath-user-request@tnt.math.metro-u.ac.jp
with only the next line in the message body.
subscribe

Project and Code Repository

NZMATH project is on SourceForge.net now.
http://sourceforge.net/projects/nzmath/
The latest source codes are provided by Mercurial. Details are in
http://sourceforge.net/scm/?type=hg&group_id=171032

Development Centre

For more information or questions about NZMATH, contact:

  NZMATH development Group
  Department of Mathematics and Information Sciences
  Tokyo Metropolitan University
  1-1 Minami-osawa
  Hachioji, Tokyo
  1920397, JAPAN

e-Mail:  nzmath-support at tnt.math.metro-u.ac.jp
Tel   :  +81 42 677 2463
Fax   :  +81 42 677 2481

Supervisor: NAKAMULA, Ken / Chief Developer: TANAKA, Satoru

(The chief developer of NZMATH has changed from former MATSUI, Tetsushi, to TANAKA, Satoru, in April, 2010.)

2010年5月20日木曜日

Partition HOWTO 2 (日本語版)

この HOWTO 第2部では、高度にカスタマイズされた分割の生成を扱う。 第1部と同じように "from nzmath.combinatorial import *" してあることを仮定する。 次期 NZMATH もまもなく登場予定だが、待ちきれない人は前回と同じく mercurial リポジトリから直接クローンして試すことができる。

1. PartitionDriver

第1部では、partition_generator をいかに使いこなすかという課題に取り組んだ。 しかし、この方法は何かしら全単射を与えることができないならば、 無駄に多く生成してフィルタリングするという効率の悪い方法を採らざるを得ない。 より柔軟に条件付けた分割の生成には、PartitionDriver というクラスを使う。 正確には、必要に応じて PartitionDriver から継承したクラスを作成する。

まずは基本形から見ていこう。 partition_generator の呼び出しと同等のことは次のようにすれば良い。

>>> driver = PartitionDriver()
>>> for partition in driver.generate(5):
...     print partition
...
(5,)
(4, 1)
(3, 2)
(3, 1, 1)
(2, 2, 1)
(2, 1, 1, 1)
(1, 1, 1, 1, 1)
>>> 

インスタンスを生成して、generate を呼び出すだけである。

generate は、拡張されることを考えて書かれたテンプレートメソッドである。 これからこの拡張を少しずつ学んでいくために、まずはそのコードを見ておこう。

def generate(self, *args):
    """
    Generate all partitions of integer

    Variable arguments *args will be used in inherited method.
    """
    self.set_parameters(*args)
    if not self.ispreconditioned():
        raise StopIteration

    while True:
        while self.rest:
            self.pull()

        if not self.rest:
            yield tuple(self.partition)

        self.backtrack()
        if self.shouldterminate():
            break
        self.push()

set_parameters, ispreconditioned などのメソッドを再定義することで 振る舞いを変えることができる。元々の振る舞いは以下の通りである。 (できれば実際のコードを参照して欲しい。)

set_parameters:
一つの自然数を受け取り、self.rest に代入する。 self.partition を空リスト([])に初期化する。
ispreconditioned:
常に True を返す。
pull:
self.rest から self.partition の最後の要素を超えない分だけ引き去り、 それを self.partition の最後に付け足す。self.partition が空のときは self.rest 全てを self.partition の1要素とする。
backtrack:
self.partition の最後に並ぶ 1 を全て取り去り、それらを合算したものを self.rest に代入する。
shouldterminate:
self.partition が空ならば True を返す。それ以外ならば False。
push:
self.partition の最後の要素から 1 を引き、self.rest に 1 を加える。

2. 奇数への分割, 再び

第1部の最後に、奇数への分割を扱った。いまそれをジェネレータ関数の形で書けば、

def partition_into_odd(n):
    for m in range(n//2 + 1, n + 1):
        for partition in partition_generator(n - m, 2*m - n):
            yield tuple(2*x - 1 for x in partition_conjugate((2*m - n,) + partition))

となる。これは n の奇数への分割を次のように構成している。 m を (n/2, n] の中にある整数を動くものとし、m の最大がちょうど 2m-n になる分割の共役の各部分 x を 2x-1 に変換する。 見通しが良いとは言い難い。

それでは、PartitionDriver を継承して、奇数への分割を生成するようにしてみよう。 全ての数を使う代わりに、奇数だけを self.partition に収めるようにすれば良い。 要素を選んでいるのは pull と push なので、そこだけ変更する。

pull は self.rest から self.partition に足される部分を送り込む動作である。 デフォルトでは送り込めるだけ送り込む方針だが、ここでは奇数しか送り込んではいけない。 そこで、self.partition の最後の要素(これは奇数しか入っていない)よりも self.rest が大きければ self.partition の最後の要素と同じだけ、そうでなければ self.rest から取れる最大の奇数を self.partition に送り込むようにする。

一方の push は、次の分割の生成に移るために 1 でない最後の要素を減らす動作である。 ここでも、デフォルトでは 1 ずつ減らすが、これではせっかく奇数だった部分が偶数 になってしまうので、2 減らして奇数性を保つようにする。

最終的に、奇数への分割を生成するクラスはこうなる。

class OddPartitionDriver(PartitionDriver):
    def __init__(self):
        PartitionDriver.__init__(self)

    def pull(self):
        if self.partition and self.rest >= self.partition[-1]:
            part = self.partition[-1]
        elif self.rest % 2:
            part = self.rest
        else:
            part = self.rest - 1

        self.partition.append(part)
        self.rest -= part

    def push(self):
        self.partition[-1] -= 2
        self.rest += 2

コードの量は増えるが、やっていることは奇数を使って分割を作っていくという目的そのものに対応しているので、解り難くはないだろう。

3. そこそこ相異なる

全ての部分が相異なる分割は良く知られているが、ここではもう少し制限を緩めて、同じ大きさの部分は高々二つという分割を考えてみよう。

まずは pull を考えよう。 同じものを二つまでしか使えない、という制約を与えるためには、self.partition の 最後二つを見なければならない。 つまり、self.partition の最後の要素と、最後から二番目の要素が等しくなければ、 self.partition と同じだけ self.partition に送り込んでも良い。 二つが等しければ、大きくとも最後の要素より 1 小さい数しか送り込んではいけない。

    def pull(self):
        if len(self.partition) >= 2 and self.partition[-1] == self.partition[-2]:
            maximum = self.partition[-1] - 1
        elif self.partition:
            maximum = self.partition[-1]
        else:
            maximum = 0

        if self.rest >= maximum:
            part = maximum
        else:
            part = self.rest

        self.partition.append(part)
        self.rest -= part

今度は backtrack も考えないといけない。 backtrack は要するに生成の様子を樹上に配置したとして分岐ノードまで巻き戻すという操作である。 デフォルトでは最後の 1 の並びだけを取り去るが、今考えているそこそこ相異なるものの生成では、たとえば (1, 1) の前が (2, 2) だったら、もはやそれ以上その並びから進むことはできないので、(2, 2) の部分も取り去らないといけない。 そうすると、単純に二つずつ並んでいたら取り除くということで良いように思ってしまうかもしれないが、残念ながら必要な条件はそれほど簡単ではない。 (4, 3, 2, 1, 1) という分割から進むことを考えてみると、まず (1, 1) を取り除くことになる。 ところが、次の 2 を崩すとすれば 1 が二つまでという制約に違反することになるから、2 も取り除かねばならない。 ところがところが、さらに次の 3 を崩すと、2 と 1 か 1 を三つ(この時点でアウト)にしないといけないので、結局 1 が二つまでという制約を満たすことができなくなる。 したがって、3 も取り除く。 この勢いだと 4 も取り除かなければいけなくなるかと思いきや、3 と 1 に分けた内の 1 を既にある 1 の一つと合わせて 2 にすれば、(3, 3, 2, 2, 1) にできるので、4 は取り除かない。 さて、この判断を実際に新しい分割を構成せずにしなければならないが、どのようにすれば良いのだろう。

1 から k までの和は k*(k+1)/2 である。 そこで、いくつかの数を分割から取り除いて残っている最後の(一番小さい)数を m としよう。 次の分割に移るには、m を小さくするしかないので、m と今までに取り除いた分の数 r の和を m-1 までの数のそれぞれの数が高々二つまでの和に直せないといけない。 つまり、m + r <= (m-1)*(m-1+1)/2*2 = m*(m-1) でなければここは分岐ノードではなく、m も取り除いてさらに上のノードまで巻き戻さなければならない。 逆に m + r <= m*(m-1) が成り立っていればどうだろう。 等号が成り立つならば 1 から m-1 まで二つずつ埋められることは明らかである。 不等号であれば、そこからいくつか抜けばいいので、少なくとも一通りは条件に合う分割が得られる。

まとめよう。 まず 1,2,2,3,3 または 1,1,2,2 のように下から順番に埋まっている部分を取り除いて足し合わせ r とする。 さらに残った部分の内最小のものを m として、m + r > m*(m-1) である間その m を取り除き続ける(取り除いたら r に加える)。 つまり、こうなる。

    def backtrack(self):
        if len(self.partition) >= 2 and self.partition[-1] == 1:
            if self.partition[-2] == 1:
                i = 1
            else:
                self.rest += self.partition.pop()
                i = 2
            while len(self.partition) >= 2 and self.partition[-2] == i:
                self.rest += self.partition.pop()
                self.rest += self.partition.pop()
                i += 1
            if self.partition and self.partition[-1] == i:
                self.rest += self.partition.pop()

        while len(self.partition) >= 2 and self.partition[-1] + self.rest > self.partition[-1] * (self.partition[-1] - 1):
            self.rest += self.partition.pop()

もう既に十分複雑化しているが、ここまでのものだけだと最後の分割を取り出した後、 さらに絞り出そうとして無限ループに陥ることがある。 その苦悩から救ってやるために shouldterminate も再定義する。 止まらなければならないケースでは True を返し、 続行できるケースでは False を返す。 だいたい先ほどの説明が解ってもらえていれば下のコードの内容は摑めると思う。

    def shouldterminate(self):
        if not self.partition:
            return True
        n = self.rest + sum(self.partition)
        if self.partition[0] >= self.rest + 1:
            return False
        return self.partition[-1] + self.rest > self.partition[0] * (self.partition[0] - 1)

思ったよりも長くなってしまったが、ともかく、分割の満たすべき性質をコードに落とせば、 毎回同じように使い回すループは書かなくてもいいようにクラスができている、という点さえ押さえてもらえれば十分である。

2010年4月5日月曜日

Partition HOWTO (English version)

This HOWTO explains how to use integer partitions in NZMATH. Note that a part of features explained here is not provided in the current released versions, but will be released in the next version. If you are interested in the features, please clone the mercurial repository available on sourceforge.net:

% hg clone http://nzmath.hg.sourceforge.net:8000/hgroot/nzmath/nzmath

1. Partition

A partition of a natural number n means a sum of natural numbers whose evaluated value is n. For example, 2+2+1 is a partition of 5. As sum does not change its value even if the order of terms are changed, we do not care about the order of summands. That is, we do not distinguish 2+2+1, 2+1+2 and 1+2+2. The partitions often appear when we consider a combinatorial structure. Thus NZMATH has the related functions in the combinatorial module. We assume to have imported the module as from nzmath.combinatorial import * in the following examples.

The first thing we have to know is how to obtain the partitions. We use the partition_generator.

>>> for partition in partition_generator(5):
...     print partition
...
(5,)
(4, 1)
(3, 2)
(3, 1, 1)
(2, 2, 1)
(2, 1, 1, 1)
(1, 1, 1, 1, 1)
>>>

These are the whole partitions of 5. Instead of an expression of addition, the generator yields tuples of summands.

There has been seven partitions of 5, you may wonder how the number of partitions increase along with increasing n from 5. In principle, you can get the numbers by len(list(partition_generator(n))), but it is inefficient; we have another function to calculate the number of partitions only:

>>> partition_number(6)
11
>>>

We shall not use the function in this HOWTO anymore. The number of partitions increases almost exponentially.

2. Partition with limited maximum

Let's think about the case when a limit of maximum of summands is set. That is for example, we would like to obtain the partitions of 5 whose summands are less than or equal to three. The simplest strategy discarding unnecessary partitions looks like this:

>>> for partition in partition_generator(5):
...     if all(d <= 3 for d in partition):
...         print partition
...
(3, 2)
(3, 1, 1)
(2, 2, 1)
(2, 1, 1, 1)
(1, 1, 1, 1, 1)
>>>

This kind of "filtering" is not a wrong approach. It can be the best way in general if there is no other means. However, this section is written to say: "partition_generator can take the second argument."

>>> for partition in partition_generator(5, 3):
...     print partition
...
(3, 2)
(3, 1, 1)
(2, 2, 1)
(2, 1, 1, 1)
(1, 1, 1, 1, 1)
>>>

3. Partition with limited number of summands

Next situation is that the number of summands is limited. In the partitions of 5, we only need the number of summands is less than or equal to 3, say. Again, we start with a general method of filtering.

>>> for partition in partition_generator(5):
...     if len(partition) <= 3:
...         print partition
...
(5, )
(4, 1)
(3, 2)
(3, 1, 1)
(2, 2, 1)
>>>

By the way, the number of partitions of 5 into parts less than or equal to 3 is five and into less than or equal to 3 parts is also five. This is not an accident, i.e., there is a one-to-one correspondence between them called "conjugate". The conjugation is easy to understand through diagrams (called the Young diagram or the Ferrers diagram). The left diagram is (4, 1), and right (2, 1, 1, 1).

. . . .    . .
.          .
           .
           .

Counting points in rows from top to bottom gives (4, 1) to the left diagram and (2, 1, 1, 1) right. The conjugation is an operation of these diagrams; fix the left top corner and turn over the plane. The partition (4, 1) can be seen as the conjugate partition of (2, 1, 1, 1).

NZMATH has the function partition_conjugate.

>>> for partition in partition_generator(5, 3):
...     print partition_conjugate(partition)
...
(2, 2, 1)
(3, 1, 1)
(3, 2)
(4, 1)
(5,)
>>>

The order is different (accidentally in the reverse order), but the partitions are the same.

4. Partition with maximum exactly equal to k

The previous sections merely introduce the pre-defined functions. Now, we start to construct some kinds of partitions as needed. We have already seen the second argument for partition_generator makes it possible to partition a number into parts "at most k". Then, how to construct a partition with the maximum part exactly equal to k? Of course, filtering strategy still works, but you can construct such partitions without any to-be-filtered partitions. Notice that the "part exactly equal to k" is common among the necessary partitions. Then, by concatenating it with partitions of n-k into parts at most k, the partitions of n with the maximum exactly equal to k can be obtained.

# n=6, k=3
>>> for partition in partition_generator(6-3, 3):
...     print (3,)+partition
...
(3, 3)
(3, 2, 1)
(3, 1, 1, 1)
>>>

Problem: How to construct the partitions of n into exactly k parts? (hint: use conjugate.)

5. Partition with limited minimum

We have already covered the way to restrict the maximum of summands or of the number of summands. Next problem is how to restrict the minimum of them. For example, how to obtain the partition of 7 whose summands are at least 2. The conjugation is useful for this situation, again. A partition into parts at least 2 is conjugate to a partition with at least 2 biggest parts. Thus, this is an extension of the previous section.

>>> n=7; k=2
>>> for m in range(1, n//k + 1):
...     for partition in partition_generator(n - m*k, m):
...         print (m,)*k + partition:
... 
(1, 1, 1, 1, 1, 1, 1)
(2, 2, 2, 1)
(2, 2, 1, 1, 1)
(3, 3, 1)
>>>

The partition we need are the conjugates of them.

>>> n=7; k=2
>>> for m in range(1, n//k + 1):
...     for partition in partition_generator(n - m*k, m):
...         print partition_conjugate((m,)*k + partition):
... 
(7,)
(4, 3)
(5, 2)
(3, 2, 2)
>>>

6. Partition into odd parts

Next target is, much different to the previous sections, partition into odd parts. The partitions of 9 into odd parts, for example, are:

(9,)
(7, 1, 1)
(5, 3, 1)
(5, 1, 1, 1, 1)
(3, 3, 3)
(3, 3, 1, 1, 1)
(3, 1, 1, 1, 1, 1, 1)
(1, 1, 1, 1, 1, 1, 1, 1, 1)

The number of partitions of 9 is thirty, but the number of those into odd parts is only eight. Considering the difference, filtering should be avoided. By the following correspondence:

9 = 2*5 - 1
7 = 2*4 - 1
...
1 = 2*1 - 1

odd numbers and natural numbers are mapped one-to-one. Then the partitions are mapped:

(9,)                        ⟷ (5,)
(7, 1, 1)                   ⟷ (4, 1, 1)
(5, 3, 1)                   ⟷ (3, 2, 1)
(5, 1, 1, 1, 1)             ⟷ (3, 1, 1, 1, 1)
(3, 3, 3)                   ⟷ (2, 2, 2)
(3, 3, 1, 1, 1)             ⟷ (2, 2, 1, 1, 1)
(3, 1, 1, 1, 1, 1, 1)       ⟷ (2, 1, 1, 1, 1, 1, 1)
(1, 1, 1, 1, 1, 1, 1, 1, 1) ⟷ (1, 1, 1, 1, 1, 1, 1, 1, 1)

You can see that the right hand sides consist of the partition of 5 into exactly 1 part, the partitions of 6 into exactly 3 parts, ... and the partition of 9 into exactly 9 parts. The partitions of n into exactly k parts is known by the problem of section 4, thus we can use it.

>>> n = 9
>>> for m in range(n//2 + 1, n + 1):
...     for partition in partition_generator(n - m, 2*m - n):
...         print tuple(2*x - 1 for x in partition_conjugate((2*m - n,) + partition))
... 
(9,)
(3, 3, 3)
(5, 3, 1)
(7, 1, 1)
(3, 3, 1, 1, 1)
(5, 1, 1, 1, 1)
(3, 1, 1, 1, 1, 1, 1)
(1, 1, 1, 1, 1, 1, 1, 1, 1)
>>> 

Problem: A partition into odd parts is mapped one-to-one to a partition into different parts by applying the rule "fuse two same parts." By this fact, obtain the partitions of n into different parts.

7. Conclusion

We have explained the way to construct the partitions, from plain partitions to partitions with some limitations. The last one or two seems rather easy if we write specialized generator, but anyway you can grasp how to construct the partitions through a certain one-to-one correspondence. We will continue the exploration in HOWTO II, where we will explain such specialized generators.

References:
G.E.Andrews and K.Eriksson. Integer Partitions. Cambridge University Press.

Partition HOWTO (日本語版)

この HOWTO は NZMATH で数の分割を使う方法を説明するものです。 注意として、一部の機能は次期リリースに含まれる予定のもので、現在のリリースバージョンには含まれていないということです。 もし、ここでの説明が気に入って次期リリースに先行して使ってみたいという方は、sourceforge.net で公開されている mercurial のリポジトリを clone して下さい。

% hg clone http://nzmath.hg.sourceforge.net:8000/hgroot/nzmath/nzmath

1. 数の分割

自然数 n の分割とは、n を自然数の和として表すことである。 たとえば、2+2+1 は 5 の分割の一つである。 足し算は順番を変えても結果が同じなので、現れる数の順番は気にしないことにする。 つまり 2+2+1 と 2+1+2 と 1+2+2 は区別しない。 このような数の分割は、組合せ構造を考えるときにしばしば現れる。 そこで NZMATH では combinatorial モジュールに関連する関数が収められている。 以下の例では from nzmath.combinatorial import * してあると仮定する。

まず分割そのものを得るには partition_generator を使う。

>>> for partition in partition_generator(5):
...     print partition
...
(5,)
(4, 1)
(3, 2)
(3, 1, 1)
(2, 2, 1)
(2, 1, 1, 1)
(1, 1, 1, 1, 1)
>>>

これが 5 の分割の全てである。 足し算の式を返すのではなく、足される数のタプルが返るようになっている。

ここで7個の分割が現れたが、n を 5 から増やすとどれぐらい分割の個数は増えるだろう、 と疑問に思うかもしれない。 原理的には len(list(partition_generator(n))) で求められるが、 これは非効率的であり、この個数(分割数という)だけ求める方法がある。

>>> partition_number(6)
11
>>>

この関数はこの HOWTO では以降使わない。 分割数の増え方はほぼ指数的である。

2. 上限付き分割

今度は足される数に上限がある場合を考えてみよう。 つまり 5 の分割の内、たとえば足される数が3以下のものだけ得たいという場合にどうするかである。 単純に要らないものを捨てるならばこう書けばよいだろう。

>>> for partition in partition_generator(5):
...     if all(d <= 3 for d in partition):
...         print partition
...
(3, 2)
(3, 1, 1)
(2, 2, 1)
(2, 1, 1, 1)
(1, 1, 1, 1, 1)
>>>

この手のフィルタリングは間違いではない。 というより、一般的に他に手段がないのならこれがベストだろう。 しかし、この節は次のことを言うためにある: 「partition_generator には第2引数が渡せます」

>>> for partition in partition_generator(5, 3):
...     print partition
...
(3, 2)
(3, 1, 1)
(2, 2, 1)
(2, 1, 1, 1)
(1, 1, 1, 1, 1)
>>>

3. 個数に上限のある分割

次は足される個数に上限がある場合を考えよう。 5 の分割の中で足される個数が3個以下のものが欲しい、と。 再び、まずは一般的な手法から。

>>> for partition in partition_generator(5):
...     if len(partition) <= 3:
...         print partition
...
(5, )
(4, 1)
(3, 2)
(3, 1, 1)
(2, 2, 1)
>>>

さて、先ほどの3以下のものの個数が5個、この足される個数3個以下のものも5個。 これは偶然ではない。 即ち、これらの間には共役と呼ばれる一対一対応がある。 この対応は図形的に見ると解り易い(ヤング図形またはフェラーズ図形と呼ばれる)。 下図の左が (4, 1) 右が (2, 1, 1, 1) である。

. . . .    . .
.          .
           .
           .

点を横に数えると左は上から4,1。 右は上から 2,1,1,1 である。 共役はこの図形を裏返す操作である。 左上の角を固定して、右のものを下に持ってくるようにひっくり返す。 (4, 1) はこの操作で (2, 1, 1, 1) と重なるのが見て取れると思う。

NZMATH には partition_conjugate という関数がある。

>>> for partition in partition_generator(5, 3):
...     print partition_conjugate(partition)
...
(2, 2, 1)
(3, 1, 1)
(3, 2)
(4, 1)
(5,)
>>>

先ほどと順番が異なる(たまたま逆になった)が、同じ分割が出てくる。

4. ちょうど k 以下への分割

今までの節は既に定義されている関数の紹介に過ぎなかったが、ここからは必要に応じて組み立てるやり方である。 partition_generator の第2引数を使うことで「高々 k」の部分に分割する方法は学んだが、では、一番大きい部分が「ちょうど k」となる分割を考えよう。フィルターする戦略はここでも有効であるが、少し考えると無駄になる分割を一切生成せずに済ませることができることが解る。 即ち「ちょうど k」の部分は全ての必要な分割に共通することに着目する。 すると、残る n-k を高々 k の部分に分割したものと一緒にする(接ぐ)ことで、 最大がちょうど k の分割が得られる。

# n=6, k=3
>>> for partition in partition_generator(6-3, 3):
...     print (3,)+partition
...
(3, 3)
(3, 2, 1)
(3, 1, 1, 1)
>>>

問題: n をちょうど k 個の部分に分割する方法を考えよ。 (ヒント: 共役を用いよ)

5. 最小を制限された分割

最大を制限する方法は既に学んだが、では最小を制限するにはどうすればよいだろう。 たとえば、「少なくとも 2」の部分への 7 の分割というような問題である。 ここでも共役を考えるのがよいだろう。 少なくとも 2 の部分への分割は、少なくとも二つの最大部分を持つ分割の共役である。 するとこれは前節の拡張である。

>>> n=7; k=2
>>> for m in range(1, n//k + 1):
...     for partition in partition_generator(n - m*k, m):
...         print (m,)*k + partition:
... 
(1, 1, 1, 1, 1, 1, 1)
(2, 2, 2, 1)
(2, 2, 1, 1, 1)
(3, 3, 1)
>>>

得たいのはこれの共役であったから、

>>> n=7; k=2
>>> for m in range(1, n//k + 1):
...     for partition in partition_generator(n - m*k, m):
...         print partition_conjugate((m,)*k + partition):
... 
(7,)
(4, 3)
(5, 2)
(3, 2, 2)
>>>

6. 奇数部分への分割

最大・最小の制限は大体理解できてきたと思うので、次は大きく方向を変えて、奇数部分への分割を考えよう。 たとえば 9 の奇数部分への分割は

(9,)
(7, 1, 1)
(5, 3, 1)
(5, 1, 1, 1, 1)
(3, 3, 3)
(3, 3, 1, 1, 1)
(3, 1, 1, 1, 1, 1, 1)
(1, 1, 1, 1, 1, 1, 1, 1, 1)

である。9 の分割は全部で30個あるが、その内奇数部分への分割は8個だけである。 この差を考えればフィルター的な作り方は避けたい。

9 = 2*5 - 1
7 = 2*4 - 1
...
1 = 2*1 - 1

という対応で 2N+1 (奇数) と N (自然数) を対応づけると、

(9,)                        ⟷ (5,)
(7, 1, 1)                   ⟷ (4, 1, 1)
(5, 3, 1)                   ⟷ (3, 2, 1)
(5, 1, 1, 1, 1)             ⟷ (3, 1, 1, 1, 1)
(3, 3, 3)                   ⟷ (2, 2, 2)
(3, 3, 1, 1, 1)             ⟷ (2, 2, 1, 1, 1)
(3, 1, 1, 1, 1, 1, 1)       ⟷ (2, 1, 1, 1, 1, 1, 1)
(1, 1, 1, 1, 1, 1, 1, 1, 1) ⟷ (1, 1, 1, 1, 1, 1, 1, 1, 1)

右側を眺めると、次のような対応が見て取れる。 5 のちょうど1個への分割、6 のちょうど3個への分割、…、9 のちょうど9個への分割を集めてきたものと対応がつく。 ちょうど k 個への分割は4節の問題として提示した。 これが利用できる。

>>> n = 9
>>> for m in range(n//2 + 1, n + 1):
...     for partition in partition_generator(n - m, 2*m - n):
...         print tuple(2*x - 1 for x in partition_conjugate((2*m - n,) + partition))
... 
(9,)
(3, 3, 3)
(5, 3, 1)
(7, 1, 1)
(3, 3, 1, 1, 1)
(5, 1, 1, 1, 1)
(3, 1, 1, 1, 1, 1, 1)
(1, 1, 1, 1, 1, 1, 1, 1, 1)
>>> 

問題: 奇数部分への分割は互いに相異なる部分への分割に「同じ部分があればその二つを融合した部分を作る」という規則を繰り返し適用することで一対一に対応づけられる。このことを利用して互いに相異なる部分への分割を実現せよ。

7. 結語

単純な分割から、いくつかの制限を付けた分割を作る方法を説明した。 最後の方は専用のジェネレータを作った方が楽そうであるが、適当な一対一対応により無駄なく生成できることが見て取れたと思う。 専用のジェネレータを作る方法は HOWTO 第2部で紹介する。

以下の文献を参照した: アンドリュース/エリクソン「整数の分割」数学書房

2009年9月5日土曜日

Osaka

I, the leader of NZMATH development, have moved to Osaka University.
Online communication should increase, but there are no other changes to the development.
Please wait for a while until release of NZMATH 1.0.

2009年6月1日月曜日

NZMATH 0.92.0

NZMATH 0.92.0 Released.

See the usual places for download.