条件付きの倍数

三角関数表, 対数表をはじめとして, 数表の存在感の薄い昨今である. なにしろ手元の計算機で関数値が一発で計算出来る世の中だからだ. だが今でも, 理科年表には, 10.0から100.9までの四桁の対数表があり, なんとなく空しい.

計算機が手元にない, 古き良き時代には, 数表は重要であった. Charles Babbageが階差機関を考案したのも, 数表を自動作成したいからであった.

江戸時代に歩いて日本地図を作成した伊能忠敬に関連した展示会にいってみると, 漢数字で書かれた各種の数表も並べられていて, そんな昔から数表が使われていたことに驚く.

ところで対数表の話だ. 高木先生の解析概論には気になる脚注がある. 「Wolframノ表ニハ1000以下ノ素数ノ自然対数ノ50桁ノ表ガ掲ゲラレテヰル. コノ表ハ既ニ少年がうすガ愛用シタモノデアル. 」だ.

なるほど, ある数の対数は, その素因数の対数の和なので, 素数の対数表だけあれば充分なのであった. しかし50桁も足すのは大変であったろう. また常用対数に変換するには, ln 10で割らなければならず, これも一仕事であるが, 昔の人は計算は得意であったからだとまずは納得する.

常用対数表は, スコットランド人のJohn Napierが1614年に作成, 出版したといわれている. そこで, 2014年に, 出版400年を記念してエディンバラで会合を開こうという話が進んでいる.

Napierの後, 何人もが常用対数表を計算し, 刊行したが, 先頃, Edward Sangの対数表の話をウェブで見つけた. Sangは1875年に, 素数について28桁の常用対数表を完成させた. その計算方法は次のようだ.

小さい素数の方から順に計算していく. さて, 素数1619の常用対数を計算するには, その1619を何倍かし, 下の方の桁を...999999のようにする. 9は何桁でもよいが, 多い方がよい. 例えば

714021*1619=1155999999

を作る. この計算が, このブログのタイトルの, ...9999で終るという「条件付きの倍数」である.

この倍数の作り方は次の通り.

まず1619の1から9までの倍数表を用意する.

(map (lambda (x) (display (list x (* x 1619))) (newline))
 (a2b 1 10))
(1 1619)
(2 3238)
(3 4857)
(4 6476)
(5 8095)
(6 9714)
(7 11333)
(8 12952)
(9 14571)

被乗数が2と5以外の素数だから, 最後の桁には1から9が出揃う.

次にこの表を見ながら, 下の桁から積の各桁が9になるように1619の倍数を足していく. まず1倍の1619が9で終るから, 1619と書く. 下から2桁目は, 1だから, 最後が8で終る倍数, 3238をそこに足す. すると33999になる. 4桁目は3だから, 最後が6で終る倍数, 6476を足す. 今度は6509になる. 0のところを9にしなければならにから, 9で終る1619をもう一度足す. 2269になる. 6のところを9にするので, 11333を足す. こうしてここまでで1155999999が得られる.

      1619       1
     3238       2
----------
     33999
   6476       40
   ----
   6509
  1619       1
  -----
  2269
11333       7
-----
11559

この数を素因数分解すると, 1155999999=3×7×11×11×281×1619 である. いま常用対数を求めようとしている1619以外の素因数は, すべてこれより小さいから, それらの常用対数は分かっているはずだ.

従って, 1155999999の常用対数が得られれば, 1619の常用対数は,
log 1619=log 1155999999 - log 3 - log 7 - log 11 - log 11 - log 281
で得られるわけである.

ところで, log 1155999999=log(1156×10⁶-1).
1156=17×17×2×2だから, この常用対数は既知である.

一方, x=1/(1156×10⁶)とおいて,

log(1156×10⁶-1)=log(1156×10⁶×(1-x))=log(1156×10⁶)+log(1-x).

自然対数をlnと書くと, ln(1+x)のTaylor展開は

ln(1+x)=x-x²/2+x³/3-x⁴/4+...    (|x|≤1, x≠-1)

だったから,

log(1-x)=ln(1-x)/ln 10=(-x-x²/2-x³/3-x⁴/4-...)/ln 10

Schemeで計算してみよう. Schemeのlog関数は自然対数lnのことである.

(define x (/ 1. (* 1156 1000000)))

x=>8.650519031141868e-10

(/ (* x x) 2)=>3.741573975407383e-19

(/ (* x x x) 3)=>2.157770458712447e-28

(log 10)=>2.302585092994046

常用対数のlog 1156×10⁶ は

(/ (log (* 1156 (expt 10 6))) (log 10))=>9.06295783408451

従って log 11559999として

(- 9.06295783408451
(/ (+ x (/ (* x x) 2) (/ (* x x x) 3)) 2.302585092994046))
=> 9.062957833708824

が得られる. そこでlog 1619は

(- 9.062957833708824
(/ (+ (log 3) (log 7) (log 11) (log 11) (log 281)) (log 10)))
=>3.2092468487533754

log 1619をSchemeで直接計算すると

(/ (log 1619) (log 10)) => 3.2092468487533736

最後の2桁が違うがまぁよかろう.

こんな数値計算をするもの久しぶりだ.

条件付きの倍数

「とっておきの数学パズル」に0,1,2という問題があった. 2.2

「十進表記で0と1のみからなる自然数nの0でない倍数が存在する」というのである.

計算してみると0から9のnについて, 次のようになる. つまり3は37倍すると1だけの十進数111となる.

n       倍  十進表記
0        1         0
1        1         1
2        5        10
3       37       111
4       25       100
5        2        10
6      185      1110
7      143      1001
8      125      1000
9 12345679 111111111

この9の場合はタイガー計算器でよく遊んだ計算である.

十進表記の求め方を考えてみた. %で剰余を表すと, 3の場合は, 1%3=1, 10%3=1, 100%3=1だから, 111ではその和が3, つまり剰余が0, よって111は3の倍数である.

7の場合は, 1%7=1, 10%7=3, 100%7=2, 1000%7=6だから, 和が7の倍数になるように1と1000を選び, 1001が7の倍数となるわけだ.

そこで早速プログラムを考えてみる.

10のべきをnで割った剰余のべき集合の和のnによる剰余が0のものを求めたい.

まず(0)なるリストを用意する. これは空のべき集合である. 次に10の0乗, 1をnで割った剰余1を, これまでのリストの各要素に足し, (1)が出来, これを(0)の前にappendする. (1 0)が出来る.

次は10の1乗, 10をnで割った剰余を, これまでのリストの各要素に足す. nが7なら, 剰余3を足し, (4 3)が出来る. appendすると(4 3 1 0)になる. それぞれ(11 10 1 0)を7で割った剰余である.

同様に100を7で割った剰余2を足すと(6 5 3 2)が出来, appendして(6 5 3 2 4 3 1 0)になる. この段階では, 最後の0以外に7の倍数はないので, まだ続ける.

1000でやると剰余は10の時の剰余3と100の時の剰余2の積で6. よって新しいリストは(12 11 9 8 10 9 7 6)となり, べき集合に7の倍数7が現れた.

7は新しいリストにあるから1000番台は1, リストの右半分にあるから, 100番台は0, そのまた右半分にあるから, 10番台も0, その(7 0)では左にあるから, 1番台は1で, 結局1001が7の0と1による倍数(の最小のもの)であった.

Schemeのプログラムは次のとおり.

(define (zeroone n e rs)
 (let* ((r (modulo e n))
  (rs1 (map (lambda (x) (modulo (+ x r) n)) rs)))
  (if (member 0 rs1)
   (- (length (member 0 (append rs1 rs))) 1)
   (zeroone n (* e 10) (append rs1 rs)))))

引数のeは10の順次のべきで, 最初の起動では1, rsはリストで最初は(0). let*で用意するrはこの桁の剰余, rs1は新しいリストである. let*の後の本体のifは, 新しいリストに0があるかを見, あればその位置を返す. なければ次の桁へ進む.

起動は次のようにする.

(zeroone 13 1 '(0))

9が結果の値だが, これを二進で1001としたのが求める十進表記になる. n=2から20までの値は

(map (lambda (x) (zeroone x 1 '(0))) (a2b 2 21)) =>
 (2 7 4 2 14 9 8 511 2 3 28 9 18 14 16 29 1022 25 4)

二進で示すと

(map (lambda (x) (number->string (zeroone x 1 '(0)) 2))
 (a2b 2 21)) =>
 ("10" "111" "100" "10" "1110" "1001" "1000" "111111111"
  "10" "11" "11100" "1001" "10010" "1110" "10000" "11101"
  "1111111110" "11001" "100")

ところでこういう倍数が存在する理由はなにか. 10のべきをnで割った剰余を表示してみると:

(define (foo n)
 (map (lambda (x) (modulo (expt 10 x) n)) (a2b 0 20)))

(foo 2) =>
 (1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 3) =>
 (1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)
(foo 4) =>
 (1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 5) =>
 (1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 6) =>
 (1 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4)
(foo 7) =>
 (1 3 2 6 4 5 1 3 2 6 4 5 1 3 2 6 4 5 1 3)
(foo 8) =>
 (1 2 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 9) =>
 (1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)

のようで, 0があるものはその桁を1つとれば完成. 0がないものは剰余だから当然同じ数が何回も現れる. その桁をn個足せば当然nの倍数が得られるわけだ.

面白い問題であった. ところでこういう条件付き倍数の生成を対数表の計算で利用した人がいた. Edward Sangのその話は次回に.