tsujimotterのノートブック

日曜数学者 tsujimotter の「趣味で数学」実践ノート

Wieferich 商を計算する

こんばんは。今日は簡単な記事を書きたいと思います。

経緯

インテジャーズさんの以下の記事:
integers.hatenablog.com

でこんなことが書かれていました。

 1093 3511 が実際にWieferich素数であることを証明しようと思うとき、愚直にコンピュータで  2^{p−1}−1 を計算して実際に  p^2 で割り切れることを確認するという方法があるでしょう。

誰かに  \frac{2^{1092}−1}{1093^2} および \frac{2^{3510}−1}{3511^2}
の数値を教えていただいて、ここに掲載したいという願望がありますが、とりあえずコンピュータの知識がない私には計算できません。

「愚直にコンピュータの計算」

これはやるしかないですね。笑

ということで,2つの Wieferich 素数に対応する Wieferich 商

 \displaystyle \frac{2^{1092}−1}{1093^2}

 \displaystyle \frac{2^{3510}−1}{3511^2}

を計算したいと思います。

結論

Ruby のインタプリタの irb を起動します。そして,以下のコマンドを打ち込みます。以上です。

$ irb
irb(main):001:0> 2**1092
=> 53058536290991634737396540170283858539198978395771271474551549598742698935712572215646062825029533563666321339663466341699370021653261445199826360714064955944866263669556212233526813063143284104557957658305559283253118889734882642786544086063293282737405311454375777285526890640894984855797452638426498665888535738081213096656896
irb(main):002:0> (2**1092 - 1) / (1093*1093)
=> 44413494081518198849533662331181676408048705850648409260420047728448020243362336732919931147165011282532627859449483774480512704278211797105113184470137216826755192252750567098391923538330743259784219179278230914061886704575890192672947523551514530826548476962166943834989934818423641467742786909315203600294760836095968855

というわけで以下が得られました。

(2**1092 - 1) / (1093*1093) =
444134940815181988495336623311816764080487058506484092604200
477284480202433623367329199311471650112825326278594494837744
805127042782117971051131844701372168267551922527505670983919
235383307432597842191792782309140618867045758901926729475235
515145308265484769621669438349899348184236414677427869093152
03600294760836095968855

もう1つの方もこのように計算できます。一瞬です。

$ irb
irb(main):001:0> 2**3510
=> 4123678330405087685498374221366984954089317057900104262462784875123790425105688284686404828183633224942349909964891086086676925759584441340507088115762723035704344160664907942351632995309234873230720724730672756917176305067668035542444896101177345052548810184785328926774725935717835798651444434276670015032998880347558303319286200881282911325381616533722737115450145930916961537131913757215277371800137275864121652932255403611649883620607541869985422611380925136858353400304606133084255357430301930613665490583193975084451066213800348446463658847351044167257072998067401824974454641569423021718445571724956556610925411248121271018867103767845126877315661086297321467468234450388231735471475672432407097826847947552577946732477743671085802532103177658853957960953937628759295649081969101068141306491454515657017075167637032474324430969643933568914915807923425951587598028215460799776130796600200956964482108313216473786634870222110938349411949520063457665149056813193039601852026593109271569256781038420076239187040954693149511676801395787300912507360641024
irb(main):002:0> (2**3510 - 1) / (3511*3511)
=> 334520796088972249521877348439021970668521632739721161369535098675821420517060576000381989288791212882744471313690446137964973797173276821125312886582578611478247366977650981307933376764066392568931604121568430853982556435332145725059800751625407510200379324968525004887574798342438254532542061871273107080963907172449942149451295309041171196857856472222730442529942387270877079662957292072924194692348462862019578856430094554247490847263326276263972959410467791859782458556592908683564910041063272650091249253024609321548078112788894377402773838867245982841985001856264883339301580764026168131102596601830756476790112731766101023821142322513515270703975493247557273711212411266850689262438137212444584410816438611463126445540507282364292727564139076663071447173588839499449680836423127595497870629440119526450423839243326359360343016803674886367621102114875480786438133300992243020582891706847118395648270858476725732361584689735010985080129376523801272425983067189252024203544898529776057950334148453647549917538811754435566234549121063004160704463

したがって,以下の結果が得られました。

(2**3510 - 1) / (3511*3511) =
334520796088972249521877348439021970668521632739721161369535
098675821420517060576000381989288791212882744471313690446137
964973797173276821125312886582578611478247366977650981307933
376764066392568931604121568430853982556435332145725059800751
625407510200379324968525004887574798342438254532542061871273
107080963907172449942149451295309041171196857856472222730442
529942387270877079662957292072924194692348462862019578856430
094554247490847263326276263972959410467791859782458556592908
683564910041063272650091249253024609321548078112788894377402
773838867245982841985001856264883339301580764026168131102596
601830756476790112731766101023821142322513515270703975493247
557273711212411266850689262438137212444584410816438611463126
445540507282364292727564139076663071447173588839499449680836
423127595497870629440119526450423839243326359360343016803674
886367621102114875480786438133300992243020582891706847118395
648270858476725732361584689735010985080129376523801272425983
067189252024203544898529776057950334148453647549917538811754
435566234549121063004160704463

やったね!ひと仕事おしまい!

自作プログラムの供養

今回の私の記事「どうしてこんなに淡々としているのか?」と不思議に思った方いませんか?(いないですかね・・・?)

実は,1つ残念なことがあったのです。

本当は,この記事の計算は,自作のプログラムで挑もうと思ったのです。Ruby のべき乗で,まさかこんなに大きな値の計算ができるとは思わなかったのです。

べき乗の計算をそのままやると,時間がかかりすぎるだろうから「繰り返し自乗法」を使わないとダメだな「よしがんばるぞ」と。繰り返し自乗法については,私の旧ブログの記事が参考になります。

繰り返し自乗法

せっかくがんばってプログラムを作ったので,供養として公開させてください。

# 数値クラスを拡張
class Numeric

   # 繰り返し自乗法により,べき乗を効率的に計算
   def power e
      v = 1
      p = self
   
      while e > 0
         if (e % 2) == 1
             v = (v * p) 
         end 
         
         p = (p * p)
         e = e >> 1     
      end
      v
   end
end

# 1093 を確認
base = 2
p = 1093

# 商を計算
puts "(#{base}^#{p-1} - 1) / (#{p}*#{p}) = "
puts (base.power(p-1) - 1) / (p*p)
puts ""

# あまりを確認
puts "(#{base}^#{p-1} - 1) % (#{p}*#{p}) = "
puts (base.power(p-1) - 1) % (p*p)
puts ""


# 3511 を確認
p = 3511

# 商を計算
puts "(#{base}^#{p-1} - 1) / (#{p}*#{p}) = "
puts (base.power(p-1) - 1) / (p*p)
puts ""

# あまりを確認
puts "(#{base}^#{p-1} - 1) % (#{p}*#{p}) = "
puts (base.power(p-1) - 1) % (p*p)
puts ""


一応,プログラムは問題なく動くことを確認しています。

しかしながら,策士策に溺れる,とはこのことですね。

それでは今日はこの辺で。