Next: , Previous:   [Contents][Index]

29 Number Theory


Previous: , Up: Number Theory   [Contents][Index]

29.1 Functions and Variables for Number Theory

関数: bern (n)

整数nについてn番目のBernoulli数を返します。 もしzerobernfalseなら ゼロに等しいBernoulli数は抑制されます。

burnも参照してください。

(%i1) zerobern: true$
(%i2) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
                  1  1       1      1        1
(%o2)       [1, - -, -, 0, - --, 0, --, 0, - --]
                  2  6       30     42       30
(%i3) zerobern: false$
(%i4) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
            1  1    1   5     691   7    3617  43867
(%o4) [1, - -, -, - --, --, - ----, -, - ----, -----]
            2  6    30  66    2730  6    510    798
関数: bernpoly (x, n)

変数xに関するn番目のBernoulli多項式を返します。

関数: bfzeta (s, n)

引数sに関するRiemannのゼータ関数を返します。 戻り値は多倍長浮動小数点(bfloat)です; nは戻り値の小数点以下の桁数です。

関数: bfhzeta (s, h, n)

引数shに関するHurwitzのゼータ関数を返します。 戻り値は多倍長浮動小数点(bfloat)です; nは戻り値の小数点以下の桁数です。

Hurwitzゼータ関数は以下のように定義されます。

                        inf
                        ====
                        \        1
         zeta (s,h)  =   >    --------
                        /            s
                        ====  (k + h)
                        k = 0

load ("bffac")はこの関数をロードします。

関数: burn (n)

n番目のBernoulli数の近似の有理数をを返します。 burnは、(有理)Bernoulli数が まあまあの効率で(超越的)ゼータによって近似できるという観察を利用します。

                   n - 1  1 - 2 n
              (- 1)      2        zeta(2 n) (2 n)!
     B(2 n) = ------------------------------------
                                2 n
                             %pi

bernは、返す前にインデックスnまでのBernoulli数すべてを計算するので、 burnは、大きな、孤立したn(たぶん105以上のn) に対してbernより効率的かもしれません。 burnは、255よりおおきな偶数nに対して近似を呼び出します。 奇数と255以下のnに対しては、関数bernが呼び出されます。

load ("bffac")はこの関数をロードします。bernも参照してください。

関数: cf (expr)

exprを連分数に変換します。 exprは、 連分数と整数の平方根から成る式です。 式の中のオペランドは 代数演算子を組み合わせられます。 連分数と平方根は別にして、式の中の因子は整数か有理数でなければいけません。 Maximaは、 cfの外側で連分数に関する演算について知りません。

listarithfalseにバインドした後、 cfは、引数を評価します。 cfは、リストとして表現された連分数を返します。

連分数a + 1/(b + 1/(c + ...))は、 リスト[a, b, c, ...]で表現されます。 リストの要素a, b, c, ...は、 整数に評価されなければいけません。 exprは、 may also contain sqrt (n)も含むかもしれません。nは整数です。 この場合、cfは、 変数cflengthの値掛ける周期と同じ数の連分数の項を与えます。

連分数は、 cfdisrepが返す代数表現を評価することで、 数に評価することができます。 連分数を評価する別の方法に関しては、 cfexpandも参照してください。

cfdisrep, cfexpand, cflengthも参照してください。

例:

  • exprは、連分数と整数の平方根から成る式です。
    (%i1) cf ([5, 3, 1]*[11, 9, 7] + [3, 7]/[4, 3, 2]);
    (%o1)               [59, 17, 2, 1, 1, 1, 27]
    (%i2) cf ((3/17)*[1, -2, 5]/sqrt(11) + (8/13));
    (%o2)        [0, 1, 1, 1, 3, 2, 1, 4, 1, 9, 1, 9, 2]
    
  • cflengthは、 連分数の何周期を代数的無理数のために計算するかを制御します。
    (%i1) cflength: 1$
    (%i2) cf ((1 + sqrt(5))/2);
    (%o2)                    [1, 1, 1, 1, 2]
    (%i3) cflength: 2$
    (%i4) cf ((1 + sqrt(5))/2);
    (%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
    (%i5) cflength: 3$
    (%i6) cf ((1 + sqrt(5))/2);
    (%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
    
  • 連分数は、cfdisrepが返す代数的表現を評価することによって評価されることができます。
    (%i1) cflength: 3$
    (%i2) cfdisrep (cf (sqrt (3)))$
    (%i3) ev (%, numer);
    (%o3)                   1.731707317073171
    
  • Maximaは、 cfの外側で連分数に関する演算について知りません。
    (%i1) cf ([1,1,1,1,1,2] * 3);
    (%o1)                     [4, 1, 5, 2]
    (%i2) cf ([1,1,1,1,1,2]) * 3;
    (%o2)                  [3, 3, 3, 3, 3, 6]
    
関数: cfdisrep (list)

連分数[a, b, c, ...]のリスト表現から、 形式a + 1/(b + 1/(c + ...))の通常の代数式を構成し返します。

(%i1) cf ([1, 2, -3] + [1, -2, 1]);
(%o1)                     [1, 1, 1, 2]
(%i2) cfdisrep (%);
                                  1
(%o2)                     1 + ---------
                                    1
                              1 + -----
                                      1
                                  1 + -
                                      2
関数: cfexpand (x)

連分数xのコンバージェントの 最後(列1)とその1つ前(列2)の 分子と分母の行列を返します。

(%i1) cf (rat (ev (%pi, numer)));

`rat' replaced 3.141592653589793 by 103993/33102 =3.141592653011902
(%o1)                  [3, 7, 15, 1, 292]
(%i2) cfexpand (%); 
                         [ 103993  355 ]
(%o2)                    [             ]
                         [ 33102   113 ]
(%i3) %[1,1]/%[2,1], numer;
(%o3)                   3.141592653011902
オプション変数: cflength

デフォルト値: 1

cflengthは、 値cflength掛ける周期として 関数cfが与える連分数の項の数を制御します。 従って、デフォルトは1周期を与えます。

(%i1) cflength: 1$
(%i2) cf ((1 + sqrt(5))/2);
(%o2)                    [1, 1, 1, 1, 2]
(%i3) cflength: 2$
(%i4) cf ((1 + sqrt(5))/2);
(%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
(%i5) cflength: 3$
(%i6) cf ((1 + sqrt(5))/2);
(%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]
関数: divsum (n, k)
関数: divsum (n)

divsum (n, k)は、 nの約数のk乗した和を返します。

divsum (n)nの約数の和を返します。

(%i1) divsum (12);
(%o1)                          28
(%i2) 1 + 2 + 3 + 4 + 6 + 12;
(%o2)                          28
(%i3) divsum (12, 2);
(%o3)                          210
(%i4) 1^2 + 2^2 + 3^2 + 4^2 + 6^2 + 12^2;
(%o4)                          210
関数: euler (n)

非負の整数nに対して n番目のEuler数を返します。

Euler-Mascheroni定数に関しては、%gammaを参照してください。

(%i1) map (euler, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]);
(%o1)    [1, 0, - 1, 0, 5, 0, - 61, 0, 1385, 0, - 50521]
関数: fib (n)

n項のFibonacci数を返します。 fib(0)は0に等しく、fib(1)は1に等しく、 fib (-n)(-1)^(n + 1) * fib(n)に等しい。

fibをコールした後, prevfibfib (x - 1)、 計算された最後の1つ前のFibonacci数に等しい。

(%i1) map (fib, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]);
(%o1)         [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
関数: fibtophi (expr)

exprに関するFibonacci数を 定数%phiを使って表現します。 %phiは、(1 + sqrt(5))/2, 近似的に1.61803399です。

例:

(%i1) fibtophi (fib (n));
                           n             n
                       %phi  - (1 - %phi)
(%o1)                  -------------------
                           2 %phi - 1
(%i2) fib (n-1) + fib (n) - fib (n+1);
(%o2)          - fib(n + 1) + fib(n) + fib(n - 1)
(%i3) fibtophi (%);
            n + 1             n + 1       n             n
        %phi      - (1 - %phi)        %phi  - (1 - %phi)
(%o3) - --------------------------- + -------------------
                2 %phi - 1                2 %phi - 1
                                          n - 1             n - 1
                                      %phi      - (1 - %phi)
                                    + ---------------------------
                                              2 %phi - 1
(%i4) ratsimp (%);
(%o4)                           0
関数: ifactors (n)

正の整数nに対して、 nの素因数分解を返します。 もしn=p1^e1..pk^nknの素因数への分解なら、 ifactorsは[[p1, e1], ... , [pk, ek]]を返します。

使われる素因数分解法は9973までの素数による試行除算と、 Pollardのロー法、楕円曲線法です。

(%i1) ifactors(51575319651600);
(%o1)     [[2, 4], [3, 2], [5, 2], [1583, 1], [9050207, 1]]
(%i2) apply("*", map(lambda([u], u[1]^u[2]), %));
(%o2)                        51575319651600
関数: igcdex (n, k)

リスト [a, b, u]を返します。 ここで、 unkの最大公約数で、 ua n + b kに等しいです。 引数 nkは整数でなければいけません。

igcdexはユークリッドのアルゴリズムを実装します。 gcdex.も参照してください。

コマンド load("gcdex")はこの関数をロードします。

例:

(%i1) load("gcdex")$

(%i2) igcdex(30,18);
(%o2)                      [- 1, 2, 6]
(%i3) igcdex(1526757668, 7835626735736);
(%o3)            [845922341123, - 164826435, 4]
(%i4) igcdex(fib(20), fib(21));
(%o4)                   [4181, - 2584, 1]
関数: inrt (x, n)

xの絶対値の整数n乗根を返します。

(%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
(%i2) map (lambda ([a], inrt (10^a, 3)), l);
(%o2) [2, 4, 10, 21, 46, 100, 215, 464, 1000, 2154, 4641, 10000]
関数: inv_mod (n, m)

mを法とするnの逆元を計算します。 もしnmを法とするゼロ因子なら、 inv_mod (n,m)falseを返します。

(%i1) inv_mod(3, 41);
(%o1)                           14
(%i2) ratsimp(3^-1), modulus=41;
(%o2)                           14
(%i3) inv_mod(3, 42);
(%o3)                          false
関数: isqrt (x)

整数 xの絶対値の「整数平方根」を返します。

関数: jacobi (p, q)

pqのJacobi記号を返します。

(%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
(%i2) map (lambda ([a], jacobi (a, 9)), l);
(%o2)         [1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]
関数: lcm (expr_1, ..., expr_n)

引数の最小公倍数を返します。 引数は、整数はもちろん一般式を取り得ます。

load ("functs")はこの関数をロードします。

関数: mod (x, y)

もしxyが実数で、yがゼロでないなら、 x - y * floor(x / y)を返します。 さらにすべての実数xに関して、mod (x, 0) = xが成り立ちます。 定義mod (x, 0) = xの議論に関しては、 Graham, Knuth, Patashnik著の「コンピュータの数学」の3.4節を参照してください。 関数mod (x, 1) は、周期が1でmod (1, 1) = 0mod (0, 1) = 0ののこぎり波関数です。

複素数の偏角の主値(区間(-%pi, %pi]での数)を見つけるためには、 関数x |-> %pi - mod (%pi - x, 2*%pi)を使います。 xは引数です。

xyが定数式(例えば、10 * %pi)の時、 modは、floorceilingが使うのと同じ多倍長浮動小数点評価スキームを 使います。 再び同様に、まれですが、modは間違った値を返すことがありえます。

数値でない引数xyに関して,modは、いくつかの式整理規則を知っています:

(%i1) mod (x, 0);
(%o1)                           x
(%i2) mod (a*x, a*y);
(%o2)                      a mod(x, y)
(%i3) mod (0, x);
(%o3)                           0
関数: next_prime (n)

nよりも大きな最も小さな素数を返します。

(%i1) next_prime(27);
(%o1)                       29
関数: partfrac (expr, var)

主変数varに関する部分分数式exprを展開します。 partfracは、完全な部分分数分解を行います。 利用したアルゴリズムは、 部分分数展開(元の分母の因子)の分母は互いに素であるという事実に基づいています。 分子は分母の線形結合として書け、結果が展開ということになります。

(%i1) 1/(1+x)^2 - 2/(1+x) + 2/(2+x);
                      2       2        1
(%o1)               ----- - ----- + --------
                    x + 2   x + 1          2
                                    (x + 1)
(%i2) ratsimp (%);
                                 x
(%o2)                 - -------------------
                         3      2
                        x  + 4 x  + 5 x + 2
(%i3) partfrac (%, x);
                      2       2        1
(%o3)               ----- - ----- + --------
                    x + 2   x + 1          2
                                    (x + 1)
関数: power_mod (a, n, m)

a^n mod mを計算するために 剰余アルゴリズムを使います。 ここで、anは整数で、mは正の整数です。 もしnが負なら、inv_modが剰余逆元を見つけるために使われます。

(%i1) power_mod(3, 15, 5);
(%o1)                          2
<(%i2) mod(3^15,5);
(%o2)                          2
(%i3) power_mod(2, -1, 5);
(%o3)                          3
(%i4) inv_mod(2,5);
(%o4)                          3
関数: primep (n)

素数テスト。 もしprimep (n)falseを返すなら、 nは合成数であり、 もしtrueを返すなら、nは非常に高い確立で素数です。

3317044064679887385961981より小さなnに対して、 Miller-Rabinのテストの決定的バージョンが使われます。 もしprimep (n)trueを返すなら、 nは素数です。

3317044064679887385961981よりの大きなnに対して、 primepは、 primep_number_of_tests個のMiller-Rabinの疑似素数テストと 1つのLucasの疑似素数テストを使います。 nがMiller-Rabinのテスト1つを通過する確率は 1/4より小さいです。 primep_number_of_testsに関してデフォルト値25を使うと、 通過したnが合成である確率は 10^-15よりもはるかに小さいです。

オプション変数: primep_number_of_tests

デフォルト値: 25

primepの中で使われるMiller-Rabinのテストの回数。

関数: prev_prime (n)

nよりも小さな最大の素数を返します。

(%i1) prev_prime(27);
(%o1)                       23
関数: qunit (n)

実二次数体sqrt (n)の基本単数、 すなわち、ノルムが1の要素を返します。 ここで、nは整数です。 これは、結果的にペル方程式a^2 - n b^2 = 1を解くことになります。

(%i1) qunit (17);
(%o1)                     sqrt(17) + 4
(%i2) expand (% * (sqrt(17) - 4));
(%o2)                           1
関数: totient (n)

n以下の、 nと互いに素な整数の数を返します。

オプション変数: zerobern

デフォルト値: true

zerobernfalseの時、 bernはBernoulli数を除外し、eulerはゼロに等しいEuler数を除外します。 berneulerを参照してください。

関数: zeta (n)

Riemannのゼータ関数を返します。 もしxが負の整数、0, 1,または、正の偶数なら、 Reimannのゼータ関数は厳密な値に整理されます。 正の偶数に対しては、 加えて、オプション変数zeta%pitrueでなければいけません。 (zeta%piを参照してください。) 浮動小数点または多倍長浮動小数点数に対して、Reimannゼータ関数は数値的に評価されます。 Maximaは、 有理非整数、浮動小数点数、複素数の引数を含む 他の引数すべてに対して、また、もしzeta%piが値falseなら偶数に対して、 名詞形zeta (n)を返します。

zeta(1)は未定義ですが、 Maximaは上からと下からの極限limit(zeta(x), x, ,1)を知っています。

bfzetazeta%piも参照してください。

例:

(%i1) zeta([-2, -1, 0, 0.5, 2, 3, 1+%i]);
                                             2
            1     1                       %pi
(%o1) [0, - --, - -, - 1.460354508809586, ----, zeta(3), 
            12    2                        6
                                                    zeta(%i + 1)]
(%i2) limit(zeta(x),x,1,plus);
(%o2)                          inf
(%i3) limit(zeta(x),x,1,minus);
(%o3)                         minf
オプション変数: zeta%pi

デフォルト値: true

zeta%pitrueの時、 偶数nに対して、zeta%pi^nに比例する式を返します。 そうでないなら、 偶数nに対して、zetaは名詞形zeta (n)を返します。

例:

(%i1) zeta%pi: true$
(%i2) zeta (4);
                                 4
                              %pi
(%o2)                         ----
                               90
(%i3) zeta%pi: false$
(%i4) zeta (4);
(%o4)                        zeta(4)

Next: , Previous:   [Contents][Index]