2015年7月31日金曜日

150731(2)

Sympy


多項式の展開・因数分解

以下のサイトにしたがって、インストール。
(http://sango.lab.tkikuchi.net/Members/tkikuchi/courses/ip2009/sympy)

from sympy import *

x = Symbol("x")      
f = (x + 1) ** 3          
f1 = expand(f)    
f2 = factor(f1)      
print("f = " + str(f))  
print("f1 = " + str(f1))
print("f2 = " + str(f2))

出力結果
f = (1 + x)**3
f1 = 1 + 3*x + 3*x**2 + x**3
f2 = (1 + x)**3

150731

Ruby


ラマヌジャン予想(3)

ラマヌジャン予想(2)を少し高速化してみた。

require 'prime'

# m次以下を取り出す
def mul(f_ary, b_ary, m)
  s1, s2 = f_ary.size, b_ary.size
  ary = Array.new(s1 + s2 - 1, 0)
  s10 = [s1 - 1, m].min
  (0..s10).each{|i|
    s20 = [s2 - 1, m - i].min
    (0..s20).each{|j|
      ary[i + j] += f_ary[i] * b_ary[j]
    }
  }
  ary
end

# m次以下を取り出す
def power(ary, n, m)
  return [1] if n == 0
  k = power(ary, n >> 1, m)
  k = mul(k, k, m)
  return k if n & 1 == 0
  return mul(k, ary, m)
end

def A000594(n)
  ary = Array.new(n + 1, 0)
  # ヤコビの公式の必要なところだけ取り出す
  i = 0
  j, k = 2 * i + 1, i * (i + 1) / 2
  while k <= n
    i & 1 == 1? ary[k] = -j : ary[k] = j
    i += 1
    j, k = 2 * i + 1, i * (i + 1) / 2
  end
  # 8乗してx倍
  power(ary, 8, n).unshift(0)[1..n]
end

# |t(p)| / p^5.5 の最大値
def ramanujan_conjecture(ary)
  max = 0
  r_max_ary = []
  p_ary = Prime.each(ary.size).to_a
  p_ary.each{|p|
    t = ary[p - 1]
    u = t.abs / p ** 5.5
    if max < u
      r_max_ary.push([p, t, u])
      max = u
    end
  }
  r_max_ary
end

ary = A000594(1000)
p ary
p ramanujan_conjecture(ary)

2015年7月30日木曜日

150730

Ruby


ラマヌジャン予想(2)

Ramanujan's tau function(1)および(2)、ラマヌジャン予想(1)
のコードではいちいち24乗していたので遅かった。
ヤコビの公式を使った後、8乗した方が速い。
ラマヌジャン予想(1)で行った計算もパッと求まる。

require 'prime'

# m次以下を取り出す
def mul(f_ary, b_ary, m)
  s1, s2 = f_ary.size, b_ary.size
  ary = Array.new(s1 + s2 - 1, 0)
  (0..s1 - 1).each{|i|
    (0..s2 - 1).each{|j|
      ary[i + j] += f_ary[i] * b_ary[j]
    }
  }
  ary[0..m]
end

# m次以下を取り出す
def power(ary, n, m)
  return [1] if n == 0
  k = power(ary, n >> 1, m)
  k = mul(k, k, m)
  return k if n & 1 == 0
  return mul(k, ary, m)
end

def A000594(n)
  ary = Array.new(n + 1, 0)
  # ヤコビの公式の必要なところだけ取り出す
  i = 0
  j, k = 2 * i + 1, i * (i + 1) / 2
  while k <= n
    i & 1 == 1? ary[k] = -j : ary[k] = j
    i += 1
    j, k = 2 * i + 1, i * (i + 1) / 2
  end
  # 8乗してx倍
  power(ary, 8, n).unshift(0)[1..n]
end

# |t(p)| / p^5.5 の最大値
def ramanujan_conjecture(ary)
  max = 0
  r_max_ary = []
  p_ary = Prime.each(ary.size).to_a
  p_ary.each{|p|
    t = ary[p - 1]
    u = t.abs / p ** 5.5
    if max < u
      r_max_ary.push([p, t, u])
      max = u
    end
  }
  r_max_ary
end

ary = A000594(1000)
p ary
p ramanujan_conjecture(ary)

2015年7月28日火曜日

150728(2)

Ruby


ラマヌジャン予想(1)

先ほどのコードを使って、|τ(p)| / p^5.5 の最大値を記録してみる。

require 'prime'

# m次以下を取り出す
def mul(f_ary, b_ary, m)
  s1, s2 = f_ary.size, b_ary.size
  ary = Array.new(s1 + s2 - 1, 0)
  (0..s1 - 1).each{|i|
    (0..s2 - 1).each{|j|
      ary[i + j] += f_ary[i] * b_ary[j]
    }
  }
  ary[0..m]
end

# m次以下を取り出す
def power(ary, n, m)
  return [1] if n == 0
  mul(ary, power(ary, n - 1, m), m)
end

def A000594(n)
  ary = [1]
  # 無限積のうち必要なところだけ取り出す
  (1..n - 1).each{|i|
    b_ary = Array.new(i + 1, 0)
    b_ary[0], b_ary[-1] = 1, -1
    ary = mul(ary, b_ary, n)
  }
  # 24乗してx倍
  power(ary, 24, n).unshift(0)[1..n]
end

# |t(p)| / p^5.5 の最大値
def ramanujan_conjecture(ary)
  max = 0
  r_max_ary = []
  p_ary = Prime.each(ary.size).to_a
  p_ary.each{|p|
    t = ary[p - 1]
    u = t.abs / p ** 5.5
    if max < u
      r_max_ary.push([p, t, u])
      max = u
    end
  }
  r_max_ary
end

ary = A000594(1000)
p ary
p ramanujan_conjecture(ary)

実行結果
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920, 534612, -3709
44, -577738, 401856, 1217160, 987136, -6905934, 2727432, 10661420, -7109760, -42
19488, -12830688, 18643272, 21288960, -25499225, 13865712, -73279080, 24647168,
128406630, -29211840, -52843168, -196706304, 134722224, 165742416, -80873520, 16
7282496, -182213314, -255874080, -145589976, 408038400, 308120442, 101267712, -1
7125708, -786948864, -548895690, -447438528, 2687348496, 248758272, -1696965207,
 611981400, -1740295368, 850430336, -1596055698, 1758697920, 2582175960, -141453
3120, 2686677840, -3081759120, -5189203740, -1791659520, 6956478662, 1268236032,
 1902838392, 2699296768, -2790474540, -3233333376, -15481826884, 10165534848, 46
98104544, 1940964480, 9791485272, -9600560640, 1463791322, 4373119536, -64258047
00, -15693610240, -8951543328, 3494159424, 38116845680, 4767866880, 1665188361,
-7394890608, -29335099668, 6211086336, -33355661220, 411016992, 32358470760, 451
64021760, -24992917110, 13173496560, 9673645072, -27442896384, -13316478336, -64
496363904, 51494658600, -49569988608, 75013568546, 40727164968, -60754911516, 37
534859200, 81742959102, 41767088832, -225755128648, -48807306240, -20380127040,
38305336752, 90241258356, 107866805760, 73482676310, -61972223040, -45917755128,
 -16528605184, -85146862638, -64480268160, 90047003760, -189014559360, 656558795
34, 124540889760, 115632958896, 102825676800, 498319933, -166955487888, 77646351
384, 77785143296, -359001100500, -45668121408, -262717201024, 338071388160, -431
5678416, 66971388960, 631528759932, -198311113728, -178514816480, 371563845216,
-353937956400, -583413304320, -297198746214, -112754509056, 596793577940, 119045
821440, 677211820992, -234995646528, -308865667656, -112181096448, 620204022900,
 -35130991728, -427635232164, 268217998208, -1115433620850, 154219312800, -82444
7297848, 900676761600, 784811057562, 214837039872, -255232501440, 214308444672,
1315116754406, -914804296320, -402206035896, -950091448320, -312162946368, -3996
4520664, -357832759588, -453553290624, 650708341920, 704042392032, 2754833892216
, -356462346240, -1458379197393, 800535869280, -1211595753060, 25209042176, -950
387449578, -776603298240, 426959023400, 527734751232, -1307679342480, 5998300106
40, 1681384224780, 807974455680, -996774496018, -232167481728, 1753032622824, 15
74983618560, -880090306620, 319595480064, -3691995187608, -3955776986112, 122698
4915520, -1235871806400, 2762403350592, 680222785536, 5442387685442, -1800325645
104, -703199584080, 2497932784704, -2876091504354, 1458117876384, 728391402200,
-2154174528000, -3901420374768, -1961831018448, -2150040612720, 2561714781696, 1
488221734860, 5418123087552, -2118677359896, -570305978368, 5699723069040, 48912
3048960, -6793168439188, 2349393987456, 2467454288544, -2165790200544, -82717169
640, -6190616678400, 884806004992, -1763584231440, 368875413144, -3800963013120,
 3989820497292, 1102026123072, 7334863021472, 3293650354176, 2897808426675, 2043
524703312, -1359839565924, -3954789780480, -11824411223170, -2161128090240, -225
5788918656, 10847792102400, -17563353448518, -1575741108816, 12979893235680, 763
8507905280, 9605445111360, -2775191013504, -7139577462960, 1201502453760, -23130
6909358, -11959678392, 13400796651732, -10239936590464, -8196341949810, -1863512
433216, -6159507467960, -4464190832640, -7392445116336, 8616026412000, 129830535
45252, -2800978113024, 9966916930464, 6305212824576, -8405626627440, -1364187309
6704, 23961192565506, 103576281984, 3050979729616, 4107578522880, -1459251465309
0, -15156690238368, -24273728464488, 11381333483520, -7708949021340, 42843555955
20, -6298215111720, 22789249173248, 25837706543670, 8494510953600, -376793236052
8, -6817096065024, 2437758558144, 7132769909136, -13632191675700, -6915609888768
, -16418932005874, -14323045870560, 6005256141024, -6832194969600, 2103572290708
2, -16253083703808, 16713176326532, -14413066320384, 12976653967200, 74127760237
44, -5159168680848, 22354294505472, 13420028104723, -14884896549600, 18903419273
592, -2154700825984, -23926858987458, 10263245571936, -25063854064200, -15393380
766720, -39175875516960, 26770406900400, -10770926678736, 9458784518400, 2867528
54752, 19786735148352, 20599225693704, 10524271493120, 33599791937460, -18835465
381488, 15311092828556, 13176671778816, -56890292419296, 6125580034560, 49875160
575912, -12299441172480, -99480832756438, -31562802105744, 9190709433360, -56107
996840960, 83369248359366, 9652944861504, 68647725277560, 13037603389440, 227407
97105712, 7491910712832, -73627062866280, -2451157267392, 14731871253050, 858798
6230112, 18517634430120, 26030014940160, -44996963217024, -15617000206080, -6358
4021925868, 43181266711296, 20707267642902, -66116013413184, -74777223849720, -4
165208506368, 121001428335986, 35001100737432, -21457009384776, 49099533315840,
-28250591730816, 29078298073440, 61522344410800, -1446779811840, 22691844947520,
 22809298789872, -155661561078204, -47631668958720, -25643022194650, -1024701656
1600, 42336109121040, -105161550594048, 24909815245602, 31384304219520, 47292873
863760, 36789573985920, 29139505641792, -40353221394720, 157584150853560, -46370
707891200, -2824382481819, 23922587904432, 125576623116, -14239605545984, 707011
2085260, -42072782947776, -177901220129584, 18403444948992, -35015731390206, 211
22167358880, 26724356607312, 19601856110592, -55161734023378, 88607884502592, -9
0468277326000, 227027200942080, -74185389602940, -29447637972480, 14646311632298
0, -75800137459200, -66204734658048, -66297680414208, 231449571733632, 851939898
16320, -43235954274240, -130617304450608, 1946216834244, -110419972899712, -1498
71571611810, 16876790017920, -128749205976048, -143359620687360, 159145247502864
, 69026196104496, 184104364634400, 89431229751552, 208110680273846, -17481393652
800, -44985733752960, -25171202969600, -133407937691598, 93634088994432, 3052950
6193984, -120325635798144, 8042859783630, 51600974705280, -97413424224168, -1470
20152688640, -206167580638390, -35717321636640, -74894084045928, 332311549369856
, 86888027422560, 50848256637504, -141688531396440, 113644706660352, 15039198164
0880, -136793353656960, 73403515193820, 29999547002880, 171111932338622, 1630360
42540512, -305398345130928, -134834785367040, 176095964901150, -59218902925056,
-116479278716528, -132835132300032, -77834148249312, 1985212071360, -71775829446
768, -72336417914880, 99881248225682, -21235344119808, 156291413770800, -1081664
99528320, 198763752966240, -8853009915456, -29031220908760, 218142225100800, 192
848217019101, -95755691935008, 328369848718692, 67590935548416, -120715789641300
, -176036712515328, -281089272454200, -45197025083392, -612368143631550, -695474
02240200, 164724885738504, 125336181803136, -207760719057696, 32636149582176, 46
723705697760, 226970543923200, 303483032911706, 283785869356080, 506060490060720
, -132549189534720, -729307946668938, 54138934047744, 122188164073712, 126754807
111680, -64318590362880, 421520482764432, -617380683662484, -96645454674048, 259
227709345696, -311517437656320, 331409422110312, -438383931955200, -915560900529
6, -230530682672640, -271857947399500, -170211715494912, 181380557687814, 171349
859111040, 1050837984850080, -239423044976640, 105271555603732, 5551365824592, -
78665062484736, -733526941376, 362315536077180, -321619119641568, -2199099717618
64, 587683317365760, -90173855416176, 196712206795440, -483863128068108, -114295
429237248, -886767711942420, 147828179231040, -293446222622280, -52163393486848,
 -163948629394368, 177418682792064, -108877719272500, 528449619936000, 694218140
838432, -311593285086048, 506588355787752, 160751787356160, 394818492462660, -23
9206006331136, -367511557743036, 386719719907328, 85753393288710, 20173503905856
0, -24509721895568, -364965248630784, -781259049093600, -575068621572144, -10903
97271369840, 6352678628352, 1436688754143552, -73223513510784, -239497637293656,
 -235739289139200, 927574652509722, 350220351674160, -21818651341228, -929610334
619904, 107593673896800, 582569483147712, 364931430558912, 132989157310464, -605
238167047943, 185014776512160, 589716680624820, 262773809858560, -17801288792019
6, 151157162681280, 435865277859480, -1307904735160320, 423708824644560, -620104
957048080, -907217963244684, 520996671820800, -1695266465052058, 90430376652672,
 -251187172996536, 1358440752807936, 354921326577300, -58506205395456, 752144751
323996, 437476554427008, -790555104585666, 327172600216800, 1368997013214600, 39
6895871877120, -638228464065920, 394054368140976, -221782757268240, -87848014672
7680, 187489032934806, -144126147384576, 9894172288504, -79833163038720, -930382
787277216, -504857349769968, 244971285916812, -996855800500224, -411259346541540
, -401116231836768, -27881913916584, 827184675778560, 1352425875836970, -3114396
95212800, 1432228477996772, 454650262789632, 696125644349184, 123820048340352, -
475388987464200, -306756182605824, -877659324192574, -322080674513352, 137148169
6731384, -912940321708800, 491186908840992, -453682062566208, -853270528819176,
123661090882560, 317117898149220, 574244615698992, -2434250501782764, 6294790617
45408, -563383208178560, 601532497540800, -724775059097208, -179869321928704, -3
03318212278158, 940221012407040, 558507191467680, 1641918289891200, 183554633354
400, 258502240289664, -1701984261329400, -542851981056000, 2339216023326602, -68
82068514048, 1759401252578412, 1213586422432256, 2406885276390, -494381416648896
, -2496072486898144, -2097168523591680, -541810234405440, -806395006499040, -155
2583345382048, -1155241876731264, 2473009012516862, -367466227885344, 3750318771
84720, -756226380349440, 24336777156666, 1365367018063104, 4225450805148020, 375
702242119680, -1366161820349760, -1197003853821888, 418481404089840, -1437171065
48736, -488895969711875, 2387539986154512, 1436330213398080, -1935851862485632,
1258353120405276, -220577026400640, -4263261111420568, 3220111123046400, -171187
8446675376, -2000861960624784, -1268924080945920, 592047284838912, 9804012847617
66, -1647545406661440, -1112733760765896, 1632884804812800, 1008295928583042, -5
45779130537088, 303982416160892, 459503857053696, -20844726749280, 1767049508790
720, 3435828424131096, 140675112737280, -2774210589848880, -353564910073200, 222
971113257984, 526729822113536, -1185387156769098, -444423226322880, 305028391047
1560, 304156780634112, -166349637206046, 1079927117208576, -2265102374509140, -9
57842679306240, -5330121520762738, 1526016526220832, 1005434765317584, -24782292
19952640, -862226563598400, -496974423429648, 2393919729693360, -405511548934195
2, 1848385481410944, 1794653372393280, 3719016970449144, 829999889252352, 474119
8635421922, -2904034280063664, 1868559748713000, 2146734178562496, -141307163864
1474, 514968225234624, -1256027191734224, -2817886259865600, -342679570612848, 6
78014201539584, -3031157627077068, 1783468948504320, -1435469944213620, -1476536
265859200, -2979751628238840, -16905402892288, 922102026851124, -544604278740480
, -2747313442193908, 1398970325778816, 1017280238423904, 3735877465876896, 28825
12981450200, 2733643609804800, -2127859436502828, 615432532671600, -442596506902
6536, -628483682444800, 5727465442629702, -1016066618904960, -1942652670145880,
1443076443734016, 3270933095391360, -597835565894448, -1368704107203888, 1924903
992130560, 698326457818910, -1135028972730240, -4331712693612240, -2111401637452
800, -985169554365696, -699348135403008, -1491821174778480, -2474997578876160, -
1799173520665920, -3782019620485440, 9709787093595120, -541834695843840, 3780043
874082112, 67785179563656, -58289341158216, 1467252058138496, -3274269549861750,
 -3013838954784, 2464686886685576, 817229535682560, 3082017633650397, -169682690
046240, 118269009151272, -2580464020796928, 7912852398754982, 4269629283110016,
-2065478171352120, -3667249129586688, -8276770434109008, 840377553364944, -84069
36179151460, 1295492931344640, -1552195881925920, -641384558575488, 136286979447
1992, -1124976089825280, -5387544388705500, 1323881616561072, 3333728731570524,
5434616916158976, -1510999629912864, 2171238655824000, 6817220120892752, 2652778
444947456, 3271729493403504, 1780449350470560, -3982080448605840, -1806121795645
440, -667048975014994, -3515114791751520, 2511663066476928, 4350268758528000, -7
744079729695638, 1588913631793152, -1230393932134640, -4066257732071424, 3790637
408024460, -5554789721607168, 2998000190340120, -3437752020369408, 2524111226919
170, 1037662902581760, 6038220526507512, -8011194672970624, -11145268112466978,
-46709204021856, 1347459830544800, 6337146270766080, 768846891863232, 3596917718
683440, 3285001442747640, 1035109787765760, 5234645524234464, 3089980943425152,
-9409519712300400, -1675135446577152, 6352013923780980, -3819485940068736, 13227
104255899436, 4233606694409088, -6116979573050976, -4418504751225600, 1425699068
010672, -5132574924871680, -4019022069226556, -4994656326572304, -19426551533776
80, -1072192144038400, 2302477758246246, 1079657610071040, -18558651348375264, 5
015858304614400, 2840270079131730, 3201790504598352, 782560406237064, 5742890791
658496, -1507747030957440, -732708148655616, 6511102049004840, 6905645184936960,
 5604717670694010, -193028634807120, -5085164185024588, 3164859781923840, -94951
8954853056, 2337922181380032, -1728332228810040, -1717908208386048, -18258436578
5360, 4948021935321360, -1099342046917296, -2190662393713920, 279110773799022, 1
797458017102272, -13526517861603928, -19071793268183040, -3435312302276400, -208
5312658141440, 272543929352676, 3118693073766912, 18045917610367430, 34005247535
14560, -4137570865480248, -1559486316150784, 11719129719838338, -360940755938112
0, 13305847699403280, -8389992357626880, 3872298735325440, -1761684364651680, -7
961834707368360, -1721713132339200, 4287752862251071, -4106686376126928, 5301002
172584664, 9999543942484736, -7043971523408190, 7329560283142272, -8343868958152
, -1575524037500928, 4211720434286064, -4226303157627600, -3397052374923408, -36
32092712736768, -14982644001390898, 2795502689196672, -5852007487279800, 7623581
505914880, -22256130175851894, 1868019557983488, 5442366024533060, 1217596737100
80, -1300110507573696, 1722619906722432, 10810987426382112, 14414456987320320, -
4590371381461740, -2397149957416368, 3381847082390196, -1302434439348224, 203777
23102676160, -3750993930499200, 8944439700308392, 6207816494668800, -85247669702
73078, -4770330071189760, 6011114426772000, -542984608147968, -28102437759471274
, 696749301810240, -6029568464839416, 2548958848450560, 4222089494633682, -46283
57208458424, 516091995437132, -5873015772013824, -6316091224178400, -78808763692
48608, 5719060466578536, -3879131953213440, 4398936813945856, 2897178951391200,
890229680050932, -10796918367606784, 28650951002224320, 6746142538900800, 812108
5805687400, -5660667323351040, -2714273523041472, 14696835447157200, -6785413121
403840, -4265574004065600, 11022255310711932, -3953397257724096, 72261719397504,
 -7193206955658240, -4814420815766940, 4986257257384704, -84377784300844, 200168
3841040128, -9289515101228586, -1121368936746240, -11009058285505488, 2652116416
266240, -15682896303708816, -7283592789880944, 8467147568239920, 174055333205062
40, -10574317556301408, -12145451761457280, -4863513801242680, 7607170877644800,
 3858395392796112, 17503390720054512, -5656913118074736, 3320521288261632, 46462
98291681650, -2932515937769088, 25655490084944664, -25258393596395520, 357534397
6258530, 1543646168709120, -18092058797213940, 25853256276218496, 12568540465129
824, 14817136407899616, -17832336756146640, 5546608703032320, 38637331724206586,
 -6221465024296704, -25069169854622376, -19106402842920960, -34899727878407658,
-7953826130647488, 5744373208966224, -5122449823088640, 5926337141961600, 219734
616127104, -28512291118415604, -14139215203921920, -845687870789636, 65245907375
88000, 21009050586560232, 9768672367534080, 40033434235820202, -4353133384507536
, 13342408183359360, 10509458025477120, 17299226769945120, -25220111636401920, 4
976295806607216, 3285476054138880, -22616076492128607, -2526517334489568, -10255
287323350908, 340483770574976, 26286732520684860, 1887961499633664, 184953043083
97016, 42098067939840, -18554019842302560, -8695572865852320, -21487665696030828
, -19725972671349504, -9992711669027360, 5277839322284736, 3712431555768600, 686
6990520492032, -8738799717564174, 2164172529988224, -13361513402011320, 12065015
350120320, -8350791783897330, 11612715073634592, 11892426157508232, 655956376492
0320, -13891521966029820, 21282425086618080, -11339234730690048, 906679499283712
0, -319279232436576, 7042709342934720, 23440857788319392, 10394584268931072, -16
023173525318736, 3934767105464832, 3518130472626000, 10881679211246592, -2140041
5987399554, 2613065262540000, 13352424013671120, -30328412970240000]
[[2, -24, 0.5303300858899106], [3, 252, 0.5987336124929452], [5, 4830, 0.6912133
33204735], [11, 534612, 1.0008729094970374], [17, -6905934, 1.1796504994162533],
 [47, 2687348496, 1.7091720053051764], [103, -225755128648, 1.9188140482390927]]

150728

Ruby


Ramanujan's tau function(2)

次の工夫を思いつく。

# m次以下を取り出す
def mul(f_ary, b_ary, m)
  s1, s2 = f_ary.size, b_ary.size
  ary = Array.new(s1 + s2 - 1, 0)
  (0..s1 - 1).each{|i|
    (0..s2 - 1).each{|j|
      ary[i + j] += f_ary[i] * b_ary[j]
    }
  }
  ary[0..m]
end

# m次以下を取り出す
def power(ary, n, m)
  return [1] if n == 0
  mul(ary, power(ary, n - 1, m), m)
end

p mul([1, 1], [1, 3, 3, 1], 0)
p mul([1, 1], [1, 3, 3, 1], 3)
p mul([1, 1], [1, 3, 3, 1], 6)
p power([1, 1], 4, 0)
p power([1, 1], 4, 3)
p power([1, 1], 4, 6)

出力結果
[1]
[1, 4, 6, 4]
[1, 4, 6, 4, 1]
[1]
[1, 4, 6, 4]
[1, 4, 6, 4, 1]

これを用いると昨日のコードは次のようになる。

# m次以下を取り出す
def mul(f_ary, b_ary, m)
  s1, s2 = f_ary.size, b_ary.size
  ary = Array.new(s1 + s2 - 1, 0)
  (0..s1 - 1).each{|i|
    (0..s2 - 1).each{|j|
      ary[i + j] += f_ary[i] * b_ary[j]
    }
  }
  ary[0..m]
end

# m次以下を取り出す
def power(ary, n, m)
  return [1] if n == 0
  mul(ary, power(ary, n - 1, m), m)
end

def A000594(n)
  ary = [1]
  # 無限積のうち必要なところだけ取り出す
  (1..n - 1).each{|i|
    b_ary = Array.new(i + 1, 0)
    b_ary[0], b_ary[-1] = 1, -1
    ary = mul(ary, b_ary, n)
  }
  # 24乗してx倍
  power(ary, 24, n).unshift(0)[1..n]
end
ary = A000594(28)

# OEIS A000594のデータ
ary0 =
[1,-24,252,-1472,4830,-6048,-16744,84480,-113643,
 -115920,534612,-370944,-577738,401856,1217160,
 987136,-6905934,2727432,10661420,-7109760,
 -4219488,-12830688,18643272,21288960,-25499225,
 13865712,-73279080,24647168]
# 一致の確認
p ary == ary0

2015年7月27日月曜日

150727(3)

Ruby


Ramanujan's tau function(1)

オンライン整数列大辞典の
A000594(http://oeis.org/A000594/list)
と比較し、答え合わせしてみる。

def mul(f_ary, b_ary)
  s1, s2 = f_ary.size, b_ary.size
  ary = Array.new(s1 + s2 - 1, 0)
  (0..s1 - 1).each{|i|
    (0..s2 - 1).each{|j|
      ary[i + j] += f_ary[i] * b_ary[j]
    }
  }
  ary
end

def power(ary, n)
  return [1] if n == 0
  mul(ary, power(ary, n - 1))
end

def A000594(n)
  ary = [1]
  # 無限積のうち必要なところだけ取り出す
  (1..n - 1).each{|i|
    b_ary = Array.new(i + 1, 0)
    b_ary[0], b_ary[-1] = 1, -1
    ary = mul(ary, b_ary)
  }
  # 24乗してx倍
  power(ary, 24).unshift(0)[1..n]
end
ary = A000594(28)

# OEIS A000594のデータ
ary0 =
[1,-24,252,-1472,4830,-6048,-16744,84480,-113643,
 -115920,534612,-370944,-577738,401856,1217160,
 987136,-6905934,2727432,10661420,-7109760,
 -4219488,-12830688,18643272,21288960,-25499225,
 13865712,-73279080,24647168]
# 一致の確認
p ary == ary0

150727(2)

Ruby


Reverse and Add(2)

この作業により、数字は大きくなっていくが、
その最終の値を並べるコードを書いてみた。
オンライン整数列大辞典の
A033865(http://oeis.org/A033865/list)
と比較し、答え合わせしてみる。

def reverse_and_add2(n)
  m = n.to_s.reverse.to_i
  while !(n == m)
    n += m
    m = n.to_s.reverse.to_i
  end
  n
end

def A033865(n)
  (0..n).map{|i| reverse_and_add2(i)}
end
ary = A033865(80) # n < 196なら大丈夫

# OEIS A033865のデータ
ary0 =
[0,1,2,3,4,5,6,7,8,9,11,11,33,44,55,66,77,88,99,
 121,22,33,22,55,66,77,88,99,121,121,33,44,55,33,
 77,88,99,121,121,363,44,55,66,77,44,99,121,121,
 363,484,55,66,77,88,99,55,121,363,484,1111,66,77,
 88,99,121,121,66,484,1111,4884,77,88,99,121,121,
 363,484,77,4884,44044,88]
# 一致の確認
p ary == ary0

150727

Ruby


Reverse and Add(1)

オンライン整数列大辞典の
A033665(http://oeis.org/A033665/list)
と比較し、答え合わせしてみる。

def reverse_and_add(n)
  cnt = 0
  m = n.to_s.reverse.to_i
  while !(n == m)
    n += m
    cnt += 1
    m = n.to_s.reverse.to_i
  end
  cnt
end

def A033665(n)
  (0..n).map{|i| reverse_and_add(i)}
end
ary = A033665(103)

# OEIS A033665のデータ
ary0 =
[0,0,0,0,0,0,0,0,0,0,1,0,1,1,1,1,1,1,1,2,1,1,0,1,
 1,1,1,1,2,1,1,1,1,0,1,1,1,2,1,2,1,1,1,1,0,1,2,1,2,
 2,1,1,1,1,1,0,1,2,2,3,1,1,1,1,2,1,0,2,3,4,1,1,1,2,
 1,2,2,0,4,6,1,1,2,1,2,2,3,4,0,24,1,2,1,2,2,3,4,6,
 24,0,1,0,1,1]
# 一致の確認
p ary == ary0

2015年7月26日日曜日

150726(5)

Ruby


回文数式

1〜9 を用いて左右対称な式を作ってみた。

ary = (1..9).to_a
ary.permutation(9){|p_ary|
  a = p_ary[0..2].join.to_i
  b = p_ary[3..5].join.to_i
  c = p_ary[6..8].join.to_i
  p_ary.reverse!
  d = p_ary[0..2].join.to_i
  e = p_ary[3..5].join.to_i
  f = p_ary[6..8].join.to_i
  m = a + b * c
  if m == d * e + f
    puts "#{a} + #{b} * #{c} = #{d} * #{e} + #{f} (= #{m})"
  end
}

出力結果
152 + 483 * 769 = 967 * 384 + 251 (= 371579)
152 + 769 * 483 = 384 * 967 + 251 (= 371579)
173 + 482 * 569 = 965 * 284 + 371 (= 274431)
173 + 569 * 482 = 284 * 965 + 371 (= 274431)
182 + 459 * 763 = 367 * 954 + 281 (= 350399)
182 + 763 * 459 = 954 * 367 + 281 (= 350399)
251 + 384 * 967 = 769 * 483 + 152 (= 371579)
251 + 967 * 384 = 483 * 769 + 152 (= 371579)
275 + 369 * 481 = 184 * 963 + 572 (= 177764)
275 + 481 * 369 = 963 * 184 + 572 (= 177764)
281 + 367 * 954 = 459 * 763 + 182 (= 350399)
281 + 954 * 367 = 763 * 459 + 182 (= 350399)
371 + 284 * 965 = 569 * 482 + 173 (= 274431)
371 + 965 * 284 = 482 * 569 + 173 (= 274431)
374 + 289 * 561 = 165 * 982 + 473 (= 162503)
374 + 561 * 289 = 982 * 165 + 473 (= 162503)
436 + 291 * 578 = 875 * 192 + 634 (= 168634)
436 + 578 * 291 = 192 * 875 + 634 (= 168634)
473 + 165 * 982 = 289 * 561 + 374 (= 162503)
473 + 982 * 165 = 561 * 289 + 374 (= 162503)
572 + 184 * 963 = 369 * 481 + 275 (= 177764)
572 + 963 * 184 = 481 * 369 + 275 (= 177764)
578 + 329 * 461 = 164 * 923 + 875 (= 152247)
578 + 461 * 329 = 923 * 164 + 875 (= 152247)
634 + 192 * 875 = 578 * 291 + 436 (= 168634)
634 + 875 * 192 = 291 * 578 + 436 (= 168634)
647 + 291 * 385 = 583 * 192 + 746 (= 112682)
647 + 385 * 291 = 192 * 583 + 746 (= 112682)
738 + 241 * 569 = 965 * 142 + 837 (= 137867)
738 + 569 * 241 = 142 * 965 + 837 (= 137867)
739 + 248 * 561 = 165 * 842 + 937 (= 139867)
739 + 561 * 248 = 842 * 165 + 937 (= 139867)
746 + 192 * 583 = 385 * 291 + 647 (= 112682)
746 + 583 * 192 = 291 * 385 + 647 (= 112682)
798 + 245 * 361 = 163 * 542 + 897 (= 89243)
798 + 361 * 245 = 542 * 163 + 897 (= 89243)
837 + 142 * 965 = 569 * 241 + 738 (= 137867)
837 + 965 * 142 = 241 * 569 + 738 (= 137867)
869 + 271 * 345 = 543 * 172 + 968 (= 94364)
869 + 345 * 271 = 172 * 543 + 968 (= 94364)
875 + 164 * 923 = 329 * 461 + 578 (= 152247)
875 + 923 * 164 = 461 * 329 + 578 (= 152247)
879 + 245 * 361 = 163 * 542 + 978 (= 89324)
879 + 361 * 245 = 542 * 163 + 978 (= 89324)
897 + 163 * 542 = 245 * 361 + 798 (= 89243)
897 + 542 * 163 = 361 * 245 + 798 (= 89243)
937 + 165 * 842 = 248 * 561 + 739 (= 139867)
937 + 842 * 165 = 561 * 248 + 739 (= 139867)
968 + 172 * 543 = 345 * 271 + 869 (= 94364)
968 + 543 * 172 = 271 * 345 + 869 (= 94364)
978 + 163 * 542 = 245 * 361 + 879 (= 89324)
978 + 542 * 163 = 361 * 245 + 879 (= 89324)

150726(4)

Ruby


Bell number(1)

オンライン整数列大辞典の
A000110(http://oeis.org/A000110/list)
と比較し、答え合わせしてみる。
(余談になるが、ベル数を求めていると、
「争いのない状態」というものが奇跡に近い
ことのように思えてきた。)

def bell(n)
  return 1 if n == 0
  i = 1
  bell = [1]
  while i < n
    next_bell = [bell[-1]]
    # Bell triangle
    i.times{|j| next_bell[j + 1] = next_bell[j] + bell[j]}
    bell = next_bell
    i += 1
  end
  bell[-1]
end

def A000110(n)
  (0..n).map{|i| bell(i)}
end
ary = A000110(26)

# OEIS A000110のデータ
ary0 =
[1,1,2,5,15,52,203,877,4140,21147,115975,678570,
 4213597,27644437,190899322,1382958545,10480142147,
 82864869804,682076806159,5832742205057,
 51724158235372,474869816156751,4506715738447323,
 44152005855084346,445958869294805289,
 4638590332229999353,49631246523618756274]
# 一致の確認
p ary == ary0

150726(3)

Ruby


1 / n の性質について(1)

require 'bigdecimal'

p BigDecimal(1r / 9, 100).to_s("F")
p BigDecimal(1r / 99, 200).to_s("F")
p BigDecimal(1r / 999, 300).to_s("F")
p BigDecimal(1r / 9999, 400).to_s("F")

# 次の性質はオンライン整数列大辞典のA059988のCOMMENTSに載っている
p BigDecimal(1r / (9 ** 2), 100).to_s("F")
p BigDecimal(1r / (99 ** 2), 200).to_s("F")
p BigDecimal(1r / (999 ** 2), 300).to_s("F")
p BigDecimal(1r / (9999 ** 2), 400).to_s("F")

p BigDecimal(1r / (9 ** 3), 100).to_s("F")
p BigDecimal(1r / (99 ** 3), 200).to_s("F")
p BigDecimal(1r / (999 ** 3), 300).to_s("F")
p BigDecimal(1r / (9999 ** 3), 400).to_s("F")

p BigDecimal(1r / (9 ** 4), 100).to_s("F")
p BigDecimal(1r / (99 ** 4), 200).to_s("F")
p BigDecimal(1r / (999 ** 4), 300).to_s("F")
p BigDecimal(1r / (9999 ** 4), 400).to_s("F")

出力結果
"0.11111111111111111111111111111111111111111111111111111111111111111111111111111
11111111111111111111111"
"0.01010101010101010101010101010101010101010101010101010101010101010101010101010
10101010101010101010101010101010101010101010101010101010101010101010101010101010
1010101010101010101010101010101010101010101"
"0.00100100100100100100100100100100100100100100100100100100100100100100100100100
10010010010010010010010010010010010010010010010010010010010010010010010010010010
01001001001001001001001001001001001001001001001001001001001001001001001001001001
001001001001001001001001001001001001001001001001001001001001001"
"0.00010001000100010001000100010001000100010001000100010001000100010001000100010
00100010001000100010001000100010001000100010001000100010001000100010001000100010
00100010001000100010001000100010001000100010001000100010001000100010001000100010
00100010001000100010001000100010001000100010001000100010001000100010001000100010
00100010001000100010001000100010001000100010001000100010001000100010001000100010
001"
"0.01234567901234567901234567901234567901234567901234567901234567901234567901234
567901234567901234567901"
"0.00010203040506070809101112131415161718192021222324252627282930313233343536373
83940414243444546474849505152535455565758596061626364656667686970717273747576777
879808182838485868788899091929394959697990001"
"0.00000100200300400500600700800901001101201301401501601701801902002102202302402
50260270280290300310320330340350360370380390400410420430440450460470480490500510
52053054055056057058059060061062063064065066067068069070071072073074075076077078
0790800810820830840850860870880890900910920930940950960970980991001"
"0.00000001000200030004000500060007000800090010001100120013001400150016001700180
01900200021002200230024002500260027002800290030003100320033003400350036003700380
03900400041004200430044004500460047004800490050005100520053005400550056005700580
05900600061006200630064006500660067006800690070007100720073007400750076007700780
07900800081008200830084008500860087008800890090009100920093009400950096009700980
09901"
"0.00137174211248285322359396433470507544581618655692729766803840877914951989026
0631001371742112482853224"
"0.00000103061015212836455566789206213754729212335579032854821039700133670136731
04888287012560045923988378839924600561270288849107337016734017040108255290379563
41292735538220692796756463729221611070402010102"
"0.00000000100300601001502102803604505506607809110512013615317119021023125327630
03253513784064354654965285615956306667037417808208619039469910360821291772262763
27379432486541597654712771831892955018082147213280348417487558630703777852929006
084163243324406489573658744831920009099190282375469564660757855955055"
"0.00000000000100030006001000150021002800360045005500660078009101050120013601530
17101900210023102530276030003250351037804060435046504960528056105950630066607030
74107800820086109030946099010351081112811761225127513261378143114851540159616531
71117701830189119532016208021452211227823462415248525562628270127752850292630033
08131603240332134033486357036553741382839164005409541864278437144654560465647534
85149505050515"
"0.00015241579027587258039932937052278616064624295076969974089315653101661332114
00701112635269013869836915"
"0.00000001041020355685216722896860668825805345579147262958129505451516501718552
83887760474874343898222104836766915157184548372213205424411454717576851063539187
50922148843820411038351075391204258697680828384869"
"0.00000000000100401002003505608412016522028636445556068081697014133154177302630
26029282796580644999654619905521477784451488906714923542582051962323144436208471
23450830262749290888543256027859752707725807955168448797214702260891595373226155
16124540865197638287244610484968160161070990018255802759225301086682187703"
"0.00000000000000010004001000200035005600840120016502200286036404550560068008160
96911401330154017712024230026002925327636544060449549605456598465457140777084369
13998810661148123423245419151816216729784259602082721023428480662377722926308592
51242235993782397151668368457647910012123994745716096462202483075320307315760839
08721685328856918905293878023506005974735757492149755929779405784292894745521116
86517166667171768"

150726(2)

Ruby


Primes of the form identical odd digits followed by a 1

オンライン整数列大辞典の
A089346(http://oeis.org/A089346/list)
と比較し、答え合わせしてみる。
(実行時間の関係で、333333333333333331 以下のものとの比較に限定。)

require 'prime'

def A089346(n)
  odds = [1, 3, 5, 7, 9]
  ary = []
  i = 1
  while
    odds.each{|j|
      m = ([j] * i).push(1).join.to_i
      break if m > n
      ary.push(m) if m.prime?
    }
    i += 1
  end
  ary
end
ary = A089346(333333333333333331)

# OEIS A089346のデータの一部
ary0 =
[11,31,71,331,991,3331,33331,99991,333331,3333331,
 9999991,33333331,555555555551,5555555555551,
 7777777777771,333333333333333331]
# 一致の確認
p ary == ary0

150726

Ruby


Prime numbers of the form 33…331

掲題について興味深いサイトを見つけた。
(http://stdkmd.com/nrr/3/33331.htm)
なお、オンライン整数列大辞典にも
A051200(http://oeis.org/A051200/list)
があるが、初項が異なる。

最初の8項なら、数分で求まる。

require 'prime'

m = 18
(0..m).each{|i|
  n = (10 ** i - 7) / 3
  p [i, n] if n.prime?
}

出力結果
[2, 31]
[3, 331]
[4, 3331]
[5, 33331]
[6, 333331]
[7, 3333331]
[8, 33333331]
[18, 333333333333333331]

2015年7月20日月曜日

150720(3)

Ruby


Collatz conjecture(2)

次の予想について実験してみた。
(http://www.numbertheory.org/gnubc/challenge)

結果の先頭に
3の倍数で終わるなら、0 を
-1で終わるなら、-1 を
-2、-4の繰り返しになるなら、-2 をそれぞれつけた。
また、仮に予想が正しくないとすると、これら以外が生じることになるので、
その場合は ? をつけることにする。

def matthews(n, start = n)
  return [-2, -4] if start == -2
  ary = [n]
  return ary if n == -1 || n == -2
  return ary if n % 3 == 0
  return matthews((7 * n + 2) / 3, start).unshift(n) if n % 3 == 1
  return matthews((n - 2) / 3, start).unshift(n)
end

n = 41
i = -n * 6
while i <= n
  ary = matthews(i)
  last = ary[-1]
  if last % 3 == 0
    p [0, ary]
  elsif last == -1
    p [-1, ary]
  elsif last == -2 || last == -4
    p [-2, ary]
  else # 予想が正しければこのようなことはない
    p ['?', ary]
  end
  i += 1
end

出力結果
[0, [-246]]
[0, [-245, -571, -191, -445, -149, -347, -809, -1887]]
[-2, [-244, -82, -28, -10, -4, -2]]
[0, [-243]]
[0, [-242, -564]]
[0, [-241, -81]]
[0, [-240]]
[0, [-239, -557, -1299]]
[0, [-238, -80, -186]]
[0, [-237]]
[0, [-236, -550, -184, -62, -144]]
[0, [-235, -79, -27]]
[0, [-234]]
[0, [-233, -543]]
[0, [-232, -78]]
[0, [-231]]
[0, [-230, -536, -1250, -2916]]
[0, [-229, -77, -179, -417]]
[0, [-228]]
[0, [-227, -529, -177]]
[0, [-226, -76, -26, -60]]
[0, [-225]]
[0, [-224, -522]]
[0, [-223, -75]]
[0, [-222]]
[0, [-221, -515, -1201, -401, -935, -2181]]
[0, [-220, -74, -172, -58, -20, -46, -16, -6]]
[0, [-219]]
[0, [-218, -508, -170, -396]]
[0, [-217, -73, -25, -9]]
[0, [-216]]
[0, [-215, -501]]
[0, [-214, -72]]
[0, [-213]]
[0, [-212, -494, -1152]]
[0, [-211, -71, -165]]
[0, [-210]]
[0, [-209, -487, -163, -55, -19, -7, -3]]
[0, [-208, -70, -24]]
[0, [-207]]
[0, [-206, -480]]
[0, [-205, -69]]
[0, [-204]]
[0, [-203, -473, -1103, -2573, -6003]]
[0, [-202, -68, -158, -368, -858]]
[0, [-201]]
[0, [-200, -466, -156]]
[0, [-199, -67, -23, -53, -123]]
[0, [-198]]
[0, [-197, -459]]
[0, [-196, -66]]
[0, [-195]]
[0, [-194, -452, -1054, -352, -118, -40, -14, -32, -74, -172, -58, -20, -46, -16
, -6]]
[0, [-193, -65, -151, -51]]
[0, [-192]]
[0, [-191, -445, -149, -347, -809, -1887]]
[0, [-190, -64, -22, -8, -18]]
[0, [-189]]
[0, [-188, -438]]
[0, [-187, -63]]
[0, [-186]]
[0, [-185, -431, -1005]]
[0, [-184, -62, -144]]
[0, [-183]]
[0, [-182, -424, -142, -48]]
[0, [-181, -61, -21]]
[0, [-180]]
[0, [-179, -417]]
[0, [-178, -60]]
[0, [-177]]
[0, [-176, -410, -956, -2230, -744]]
[0, [-175, -59, -137, -319, -107, -249]]
[0, [-174]]
[0, [-173, -403, -135]]
[0, [-172, -58, -20, -46, -16, -6]]
[0, [-171]]
[0, [-170, -396]]
[0, [-169, -57]]
[0, [-168]]
[0, [-167, -389, -907, -303]]
[0, [-166, -56, -130, -44, -102]]
[0, [-165]]
[0, [-164, -382, -128, -298, -100, -34, -12]]
[0, [-163, -55, -19, -7, -3]]
[0, [-162]]
[0, [-161, -375]]
[0, [-160, -54]]
[0, [-159]]
[0, [-158, -368, -858]]
[0, [-157, -53, -123]]
[0, [-156]]
[0, [-155, -361, -121, -41, -95, -221, -515, -1201, -401, -935, -2181]]
[0, [-154, -52, -18]]
[0, [-153]]
[0, [-152, -354]]
[0, [-151, -51]]
[0, [-150]]
[0, [-149, -347, -809, -1887]]
[0, [-148, -50, -116, -270]]
[0, [-147]]
[0, [-146, -340, -114]]
[0, [-145, -49, -17, -39]]
[0, [-144]]
[0, [-143, -333]]
[0, [-142, -48]]
[0, [-141]]
[0, [-140, -326, -760, -254, -592, -198]]
[0, [-139, -47, -109, -37, -13, -5, -11, -25, -9]]
[0, [-138]]
[0, [-137, -319, -107, -249]]
[0, [-136, -46, -16, -6]]
[0, [-135]]
[0, [-134, -312]]
[0, [-133, -45]]
[0, [-132]]
[0, [-131, -305, -711]]
[0, [-130, -44, -102]]
[0, [-129]]
[0, [-128, -298, -100, -34, -12]]
[0, [-127, -43, -15]]
[0, [-126]]
[0, [-125, -291]]
[0, [-124, -42]]
[0, [-123]]
[0, [-122, -284, -662, -1544, -3602, -8404, -2802]]
[0, [-121, -41, -95, -221, -515, -1201, -401, -935, -2181]]
[0, [-120]]
[0, [-119, -277, -93]]
[0, [-118, -40, -14, -32, -74, -172, -58, -20, -46, -16, -6]]
[0, [-117]]
[0, [-116, -270]]
[0, [-115, -39]]
[0, [-114]]
[0, [-113, -263, -613, -205, -69]]
[0, [-112, -38, -88, -30]]
[0, [-111]]
[0, [-110, -256, -86, -200, -466, -156]]
[0, [-109, -37, -13, -5, -11, -25, -9]]
[0, [-108]]
[0, [-107, -249]]
[0, [-106, -36]]
[0, [-105]]
[0, [-104, -242, -564]]
[0, [-103, -35, -81]]
[0, [-102]]
[0, [-101, -235, -79, -27]]
[0, [-100, -34, -12]]
[0, [-99]]
[0, [-98, -228]]
[0, [-97, -33]]
[0, [-96]]
[0, [-95, -221, -515, -1201, -401, -935, -2181]]
[0, [-94, -32, -74, -172, -58, -20, -46, -16, -6]]
[0, [-93]]
[0, [-92, -214, -72]]
[0, [-91, -31, -11, -25, -9]]
[0, [-90]]
[0, [-89, -207]]
[0, [-88, -30]]
[0, [-87]]
[0, [-86, -200, -466, -156]]
[0, [-85, -29, -67, -23, -53, -123]]
[0, [-84]]
[0, [-83, -193, -65, -151, -51]]
[-2, [-82, -28, -10, -4, -2]]
[0, [-81]]
[0, [-80, -186]]
[0, [-79, -27]]
[0, [-78]]
[0, [-77, -179, -417]]
[0, [-76, -26, -60]]
[0, [-75]]
[0, [-74, -172, -58, -20, -46, -16, -6]]
[0, [-73, -25, -9]]
[0, [-72]]
[0, [-71, -165]]
[0, [-70, -24]]
[0, [-69]]
[0, [-68, -158, -368, -858]]
[0, [-67, -23, -53, -123]]
[0, [-66]]
[0, [-65, -151, -51]]
[0, [-64, -22, -8, -18]]
[0, [-63]]
[0, [-62, -144]]
[0, [-61, -21]]
[0, [-60]]
[0, [-59, -137, -319, -107, -249]]
[0, [-58, -20, -46, -16, -6]]
[0, [-57]]
[0, [-56, -130, -44, -102]]
[0, [-55, -19, -7, -3]]
[0, [-54]]
[0, [-53, -123]]
[0, [-52, -18]]
[0, [-51]]
[0, [-50, -116, -270]]
[0, [-49, -17, -39]]
[0, [-48]]
[0, [-47, -109, -37, -13, -5, -11, -25, -9]]
[0, [-46, -16, -6]]
[0, [-45]]
[0, [-44, -102]]
[0, [-43, -15]]
[0, [-42]]
[0, [-41, -95, -221, -515, -1201, -401, -935, -2181]]
[0, [-40, -14, -32, -74, -172, -58, -20, -46, -16, -6]]
[0, [-39]]
[0, [-38, -88, -30]]
[0, [-37, -13, -5, -11, -25, -9]]
[0, [-36]]
[0, [-35, -81]]
[0, [-34, -12]]
[0, [-33]]
[0, [-32, -74, -172, -58, -20, -46, -16, -6]]
[0, [-31, -11, -25, -9]]
[0, [-30]]
[0, [-29, -67, -23, -53, -123]]
[-2, [-28, -10, -4, -2]]
[0, [-27]]
[0, [-26, -60]]
[0, [-25, -9]]
[0, [-24]]
[0, [-23, -53, -123]]
[0, [-22, -8, -18]]
[0, [-21]]
[0, [-20, -46, -16, -6]]
[0, [-19, -7, -3]]
[0, [-18]]
[0, [-17, -39]]
[0, [-16, -6]]
[0, [-15]]
[0, [-14, -32, -74, -172, -58, -20, -46, -16, -6]]
[0, [-13, -5, -11, -25, -9]]
[0, [-12]]
[0, [-11, -25, -9]]
[-2, [-10, -4, -2]]
[0, [-9]]
[0, [-8, -18]]
[0, [-7, -3]]
[0, [-6]]
[0, [-5, -11, -25, -9]]
[-2, [-4, -2]]
[0, [-3]]
[-2, [-2, -4]]
[-1, [-1]]
[0, [0]]
[0, [1, 3]]
[0, [2, 0]]
[0, [3]]
[0, [4, 10, 24]]
[0, [5, 1, 3]]
[0, [6]]
[0, [7, 17, 5, 1, 3]]
[0, [8, 2, 0]]
[0, [9]]
[0, [10, 24]]
[0, [11, 3]]
[0, [12]]
[0, [13, 31, 73, 171]]
[0, [14, 4, 10, 24]]
[0, [15]]
[0, [16, 38, 12]]
[0, [17, 5, 1, 3]]
[0, [18]]
[0, [19, 45]]
[0, [20, 6]]
[0, [21]]
[0, [22, 52, 122, 40, 94, 220, 514, 1200]]
[0, [23, 7, 17, 5, 1, 3]]
[0, [24]]
[0, [25, 59, 19, 45]]
[0, [26, 8, 2, 0]]
[0, [27]]
[0, [28, 66]]
[0, [29, 9]]
[0, [30]]
[0, [31, 73, 171]]
[0, [32, 10, 24]]
[0, [33]]
[0, [34, 80, 26, 8, 2, 0]]
[0, [35, 11, 3]]
[0, [36]]
[0, [37, 87]]
[0, [38, 12]]
[0, [39]]
[0, [40, 94, 220, 514, 1200]]
[0, [41, 13, 31, 73, 171]]

150720(2)

Ruby


Collatz conjecture(1)

1に到達するまでの数列の長さを求めるコードを書いてみた。
オンライン整数列大辞典の
A006577(http://oeis.org/A006577/list)
と比較し、答え合わせしてみる。

def collatz(n)
  return 0 if n == 1
  return 1 + collatz(n / 2) if n % 2 == 0
  return 1 + collatz(3 * n + 1)
end

def A006577(n)
  (1..n).map{|i| collatz(i)}
end
ary = A006577(72)

# OEIS A006577のデータ
ary0 =
[0,1,7,2,5,8,16,3,19,6,14,9,9,17,17,4,12,20,20,7,
 7,15,15,10,23,10,111,18,18,18,106,5,26,13,13,21,
 21,21,34,8,109,8,29,16,16,16,104,11,24,24,24,11,
 11,112,112,19,32,19,32,19,19,107,107,6,27,27,27,
 14,14,14,102,22]
# 一致の確認
p ary == ary0

150720

Ruby


隣接素数の和で表す表し方の数

オンライン整数列大辞典の
A054845(http://oeis.org/A054845/list)
と比較し、答え合わせしてみる。
A054845(http://oeis.org/A054845)によれば、
Moser shows that the average order of a(n) is log 2, that is, sum(i=1..n, a(i)) ~ n log 2
ということなのだが、以下のコードではこのことが正しいと確信できなかった。

require 'prime'

def A054845(n)
  ary = Array.new(n + 1, 0)
  p_ary = Prime.each(n).to_a
  size = p_ary.size
  # 先頭をp_ary[i]に固定する
  (0..size - 1).each{|i|
    j = i
    sum = p_ary[j]
    while sum <= n
      ary[sum] += 1
      j += 1
      break if j > size - 1
      sum += p_ary[j]
    end
  }
  ary
end
ary = A054845(101)

# OEIS A054845のデータ
ary0 =
[0,0,1,1,0,2,0,1,1,0,1,1,1,1,0,1,0,2,1,1,0,0,0,2,
 1,0,1,0,1,1,1,2,0,0,0,0,2,1,0,1,0,3,1,1,0,0,0,1,1,
 1,0,0,1,2,0,0,1,0,1,2,2,1,0,0,0,0,0,2,1,0,0,2,2,1,
 0,1,0,1,1,1,0,0,0,3,1,0,0,0,1,1,2,0,0,0,0,1,0,2,1,
 0,2,2]
# 一致の確認
p ary == ary0

def accumulation(ary)
  i = 1
  while i <= ary.size - 1
    ary[i] += ary[i - 1]
    i += 1
  end
  ary
end
# 元々
p ary
# 累積
p accumulation(ary)

M = 8
ary = A054845(10 ** M)
a_ary = accumulation(ary)
p ''
p Math.log(2) # a_ary[i]はi * log2くらい
(0..M).each{|m| p a_ary[10 ** m]}

出力結果
true
[0, 0, 1, 1, 0, 2, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 2, 1, 1, 0, 0, 0, 2, 1, 0, 1
, 0, 1, 1, 1, 2, 0, 0, 0, 0, 2, 1, 0, 1, 0, 3, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1,
2, 0, 0, 1, 0, 1, 2, 2, 1, 0, 0, 0, 0, 0, 2, 1, 0, 0, 2, 2, 1, 0, 1, 0, 1, 1, 1,
 0, 0, 0, 3, 1, 0, 0, 0, 1, 1, 2, 0, 0, 0, 0, 1, 0, 2, 1, 0, 2, 2]
[0, 0, 1, 2, 2, 4, 4, 5, 6, 6, 7, 8, 9, 10, 10, 11, 11, 13, 14, 15, 15, 15, 15,
17, 18, 18, 19, 19, 20, 21, 22, 24, 24, 24, 24, 24, 26, 27, 27, 28, 28, 31, 32,
33, 33, 33, 33, 34, 35, 36, 36, 36, 37, 39, 39, 39, 40, 40, 41, 43, 45, 46, 46,
46, 46, 46, 46, 48, 49, 49, 49, 51, 53, 54, 54, 55, 55, 56, 57, 58, 58, 58, 58,
61, 62, 62, 62, 62, 63, 64, 66, 66, 66, 66, 66, 67, 67, 69, 70, 70, 72, 74]
""
0.6931471805599453
0
7
72
684
6667
65716
650999
6475810
64581536

2015年7月19日日曜日

150719(2)

Ruby


Aliquot sequence(2)

もう少し調べてみた。
(実行時間は30分ほどかかる。)

require 'prime'

def power(a, n)
  return 1 if n == 0
  k = power(a, n >> 1)
  k *= k
  return k if n & 1 == 0
  return k * a
end

def sigma(i)
  sum = 1
  pq = i.prime_division
  pq.each{|a, n| sum *= (power(a, n + 1) - 1) / (a - 1)}
  sum
end

def aliquot(i)
  # Lehmer Five がらみは調べない
  return '?' if [276, 306, 396, 552, 564, 660, 696, 780, 828, 888, 966, 996].include?(i)
  # iが840のときは http://www.aliquot.de/sequenzen/840.elf を見た方がいい
  return 'Hard!' if i == 840
  a_ary = []
  j = sigma(i) - i
  while !a_ary.include?(j) && !(j == 0)
    a_ary.push(j)
    j = sigma(j) - j
  end
  return a_ary.size if a_ary.include?(i) || j == 0
  return a_ary.size + 1
end

p (1..1000).map{|i| aliquot(i)}

出力結果
[0, 1, 1, 2, 1, 1, 1, 2, 3, 3, 1, 6, 1, 4, 4, 5, 1, 3, 1, 6, 2, 5, 1, 4, 2, 6, 2
, 1, 1, 14, 1, 2, 5, 7, 2, 3, 1, 6, 2, 3, 1, 13, 1, 4, 6, 7, 1, 5, 3, 2, 3, 8, 1
, 12, 2, 4, 2, 3, 1, 10, 1, 8, 2, 3, 2, 11, 1, 4, 3, 5, 1, 8, 1, 4, 4, 4, 2, 10,
 1, 6, 4, 5, 1, 5, 2, 8, 6, 6, 1, 9, 3, 5, 3, 3, 3, 8, 1, 2, 3, 4, 1, 17, 1, 6,
7, 5, 1, 9, 1, 6, 2, 8, 1, 16, 2, 4, 3, 9, 3, 11, 7, 4, 7, 5, 2, 15, 1, 2, 2, 5,
 1, 10, 3, 6, 8, 7, 1, 177, 1, 4, 4, 5, 3, 8, 3, 5, 5, 10, 1, 176, 1, 11, 5, 7,
2, 7, 1, 6, 3, 10, 2, 3, 1, 6, 8, 9, 1, 174, 5, 8, 2, 8, 1, 15, 2, 4, 3, 6, 1, 5
1, 1, 8, 3, 5, 2, 14, 2, 11, 2, 9, 1, 12, 1, 5, 5, 3, 1, 13, 1, 3, 2, 7, 2, 9, 2
, 6, 8, 6, 2, 51, 1, 10, 5, 7, 4, 9, 3, 9, 3, 2, 2, 175, 1, 16, 7, 5, 1, 7, 1, 8
, 6, 10, 1, 174, 2, 6, 2, 6, 1, 10, 1, 2, 8, 10, 2, 14, 6, 11, 7, 10, 1, 16, 3,
6, 4, 5, 1, 13, 7, 10, 3, 7, 1, 30, 2, 8, 4, 7, 1, 12, 1, 11, 3, 5, 2, "?", 1, 6
, 2, 15, 1, 16, 1, 2, 6, 10, 4, 9, 4, 11, 2, 6, 1, 15, 3, 6, 4, 12, 2, 8, 4, 8,
8, 12, 2, "?", 1, 12, 2, 9, 1, 173, 1, 11, 3, 11, 1, 34, 2, 13, 3, 7, 2, 2, 2, 1
0, 2, 9, 3, 33, 1, 6, 3, 9, 2, 13, 1, 2, 4, 9, 2, 13, 3, 12, 7, 5, 1, 15, 1, 5,
3, 11, 1, 51, 3, 6, 4, 9, 1, 29, 7, 6, 6, 11, 2, 50, 1, 14, 4, 12, 2, 12, 1, 6,
8, 13, 2, 49, 1, 12, 2, 6, 1, 8, 2, 4, 3, 13, 1, 14, 2, 2, 9, 4, 3, "?", 1, 8, 2
, 10, 1, 14, 7, 10, 4, 12, 4, 27, 1, 6, 5, 12, 2, 13, 2, 8, 4, 9, 1, 10, 1, 8, 3
, 5, 4, 13, 4, 10, 9, 7, 1, 4, 1, 10, 7, 10, 2, 12, 1, 14, 9, 12, 1, 15, 4, 6, 6
, 8, 1, 11, 2, 6, 3, 9, 4, 29, 1, 11, 4, 11, 1, 13, 1, 8, 9, 7, 1, 14, 5, 5, 3,
11, 3, 16, 4, 15, 8, 3, 1, 65, 4, 11, 7, 7, 2, 15, 1, 13, 2, 9, 1, 14, 2, 6, 10,
 1, 2, 33, 1, 9, 3, 7, 1, 9, 2, 10, 8, 14, 1, 32, 5, 6, 5, 11, 2, 13, 2, 5, 4, 2
5, 1, 12, 1, 11, 2, 9, 4, 172, 5, 13, 8, 14, 3, 32, 2, 8, 4, 12, 4, 22, 1, 6, 3,
 11, 3, 31, 1, 10, 2, 12, 4, "?", 7, 16, 5, 6, 1, 13, 3, 17, 9, 3, 1, "?", 4, 11
, 2, 7, 1, 172, 1, 13, 6, 11, 6, 6, 1, 4, 2, 21, 4, 48, 3, 10, 9, 7, 1, 13, 4, 1
0, 3, 8, 1, 47, 2, 10, 3, 7, 1, 69, 1, 10, 2, 12, 2, 14, 1, 3, 4, 11, 2, 46, 1,
10, 10, 10, 1, 13, 1, 14, 5, 12, 2, 18, 8, 12, 4, 4, 3, 12, 1, 8, 5, 14, 4, 7, 3
, 13, 5, 13, 1, 30, 1, 16, 6, 11, 1, 11, 2, 3, 2, 2, 1, 29, 2, 12, 3, 7, 1, "?",
 1, 10, 8, 11, 4, 28, 2, 15, 2, 17, 2, 26, 1, 10, 5, 3, 1, 13, 8, 20, 7, 6, 1, 1
3, 4, 12, 2, 4, 2, 12, 1, 10, 6, 6, 4, "?", 2, 12, 3, 15, 1, 300, 3, 10, 7, 7, 2
, 16, 1, 8, 4, 14, 3, 10, 2, 12, 9, 7, 1, 194, 3, 2, 3, 13, 3, 172, 1, 18, 12, 1
1, 2, 21, 1, 13, 6, 8, 2, 4, 1, 24, 2, 17, 1, 28, 3, 14, 8, 23, 3, 30, 1, 9, 5,
11, 2, 20, 1, 7, 10, 19, 1, 8, 4, 22, 6, 5, 2, 22, 1, 13, 4, 8, 1, 7, 4, 7, 2, 3
, 2, "?", 2, 12, 5, 2, 2, 31, 1, 8, 5, 4, 8, 17, 5, 11, 4, 13, 1, 30, 3, 2, 5, 1
1, 3, 25, 2, 13, 4, 3, 1, 28, 1, 11, 3, 17, 6, 27, 3, 13, 4, 9, 1, 48, 1, 9, 9,
11, 1, "?", 1, 7, 2, 15, 2, 47, 2, 15, 2, 9, 1, "Hard!", 15, 6, 7, 14, 4, 46, 4,
 12, 5, 10, 2, 20, 1, 15, 8, 23, 1, 167, 1, 17, 8, 11, 1, 59, 2, 11, 8, 10, 4, 1
71, 5, 4, 2, 12, 2, 6, 1, 13, 5, 11, 1, 6, 1, 12, 6, 7, 1, "?", 9, 12, 10, 4, 2,
 14, 3, 15, 7, 7, 2, 3, 2, 12, 3, 14, 3, 13, 1, 5, 5, 19, 1, 4, 4, 12, 7, 7, 2,
12, 1, 21, 2, 9, 3, 9, 4, 8, 5, 16, 1, 45, 3, 8, 4, 6, 8, 184, 1, 7, 2, 19, 1, 6
, 3, 8, 3, 14, 1, 19, 7, 20, 4, 17, 1, 5, 2, 14, 8, 12, 4, 171, 3, 15, 10, 12, 2
, "?", 1, 5, 4, 12, 1, 15, 6, 10, 2, 15, 1, 299, 2, 9, 2, 7, 1, 17, 3, 15, 3, 16
, 2, 298, 1, 8, 3, 14, 3, "?", 1, 8, 2, 17]

150719

Ruby


Aliquot sequence(1)

オンライン整数列大辞典の
A000203(http://oeis.org/A000203/list)、
A003023(http://oeis.org/A003023/list)
と比較し、答え合わせしてみる。

# Catalanの予想が正しいと仮定
require 'prime'

def power(a, n)
  return 1 if n == 0
  k = power(a, n >> 1)
  k *= k
  return k if n & 1 == 0
  return k * a
end

def sigma(i)
  sum = 1
  pq = i.prime_division
  pq.each{|a, n| sum *= (power(a, n + 1) - 1) / (a - 1)}
  sum
end

def aliquot(i)
  a_ary = []
  j = sigma(i) - i
  while !a_ary.include?(j) && !(j == 0)
    a_ary.push(j)
    j = sigma(j) - j
  end
  return a_ary.size if a_ary.include?(i) || j == 0
  return a_ary.size + 1
end

def A000203(n)
  (1..n).map{|i| sigma(i)}
end
ary = A000203(70)

# OEIS A000203のデータ
ary0 =
[1,3,4,7,6,12,8,15,13,18,12,28,14,24,24,31,18,39,
 20,42,32,36,24,60,31,42,40,56,30,72,32,63,48,54,
 48,91,38,60,56,90,42,96,44,84,78,72,48,124,57,93,
 72,98,54,120,72,120,80,90,60,168,62,96,104,127,84,
 144,68,126,96,144]
# 一致の確認
p ary == ary0

def A003023(n)
  (1..n).map{|i| aliquot(i)}
end
ary = A003023(102)

# OEIS A003023のデータ
ary0 =
[0,1,1,2,1,1,1,2,3,3,1,6,1,4,4,5,1,3,1,6,2,5,1,4,
 2,6,2,1,1,14,1,2,5,7,2,3,1,6,2,3,1,13,1,4,6,7,1,5,
 3,2,3,8,1,12,2,4,2,3,1,10,1,8,2,3,2,11,1,4,3,5,1,
 8,1,4,4,4,2,10,1,6,4,5,1,5,2,8,6,6,1,9,3,5,3,3,3,
 8,1,2,3,4,1,17]
# 一致の確認
p ary == ary0

# 完全数のとき
p aliquot(6)     # [6]より1
# 友愛数のとき
p aliquot(220)   # [284, 220]より2
# 社交数のとき
p aliquot(12496) # [14288, 15472, 14536, 14264, 12496]より5

2015年7月12日日曜日

150712

OSQL


1> BEGIN
2>   DECLARE @a VARCHAR(4), @b INT
3>   SET @a = '文字'
4>   SET @b = 5
5>   PRINT @a
6>   PRINT @b
7> END
8> GO
文字
5
1>

2015年7月10日金曜日

150710

SQLite


次のぺージを参考にさわってみた。
(https://www.db.soc.i.kyoto-u.ac.jp/lec/le4db/index.php?SQLite)

C:\pg\sqlite>sqlite3 sample.db
SQLite version 3.8.10.2 2015-05-20 18:17:19
Enter ".help" for usage hints.
sqlite> CREATE TABLE products (pid INTEGER, name VARCHAR(20), price INTEGER, PRI
MARY KEY (pid));
sqlite> INSERT INTO products VALUES(1, 'AAA', 100);
sqlite> SELECT * FROM products;
1|AAA|100
sqlite> INSERT INTO products VALUES(2, 'BBB', 80);
sqlite> INSERT INTO products VALUES(3, 'CCC', 220);
sqlite> SELECT * FROM products;
1|AAA|100
2|BBB|80
3|CCC|220
sqlite> DROP TABLE products;
sqlite> .exit

2015年7月9日木曜日

150709(2)

Ruby


素数の個数(5)

e^n (n ≦ 33) 未満の素数の個数を求めてみた。

0
1
4
8
16
34
79
183
429
1019
2466
6048
14912
37128
93117
234855
595341
1516233
3877186
9950346
25614562
66124777
171141897
443963543
1154106844
3005936865
7842921261
20496470801
53645077679
140599114669
368973074565
969455391690
2550043255883
6714639721322

150709

Ruby


素数の個数(4)

少し改良してみた。
念のため、
オンライン整数列大辞典の
A000720(http://oeis.org/A000720/list)、
A040014(http://oeis.org/A040014/list)
と比較し、答え合わせしてみる。
(実行時間は30分ほどかかる。)

def pi(n)
  m = Math.sqrt(n).to_i
  keys = (1..m).map{|i| n / i}
  keys += (1..keys[-1] - 1).to_a.reverse
  h = {}
  keys.each{|i| h[i] = i - 1}
  (2..m).each{|i|
    if h[i] > h[i - 1] # このときiは素数
      hp = h[i - 1]
      i2 = i * i
      h[n] -= h[n / i] - hp
      # h[n / j]は「i = jのときのh[n]に呼び出されるまでは計算しておく」ことにし、
      # i = jのとき、h[n / (j + 1)]から計算していくことにする
      keys[i..-1].each{|j|
        break if j < i2
        h[j] -= h[j / i] - hp
      }
    end
  }
  h[n]
end

def A000720(n)
  (1..n).map{|i| pi(i)}
end
ary = A000720(78)

# OEIS A000720のデータ
ary0 =
[0,1,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,8,8,8,8,9,9,
 9,9,9,9,10,10,11,11,11,11,11,11,12,12,12,12,13,13,
 14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16,
 17,17,18,18,18,18,18,18,19,19,19,19,20,20,21,21,
 21,21,21,21]
# 一致の確認
p ary == ary0

def A040014(n)
  (1..n).map{|i| pi(Math.exp(i).to_i)}.unshift(0)
end
ary = A040014(29)

# OEIS A040014のデータ
ary0 =
[0,1,4,8,16,34,79,183,429,1019,2466,6048,14912,
 37128,93117,234855,595341,1516233,3877186,9950346,
 25614562,66124777,171141897,443963543,1154106844,
 3005936865,7842921261,20496470801,53645077679,
 140599114669]
# 一致の確認
p ary == ary0

2015年7月8日水曜日

150708

Ruby


素数の個数(3)

(確認する順序が逆になってしまったが、)
オンライン整数列大辞典の
A000720(http://oeis.org/A000720/list)
と比較し、答え合わせしてみる。

def pi(n)
  m = Math.sqrt(n).to_i
  keys = (1..m).map{|i| n / i}
  keys += (1..keys[-1] - 1).to_a.reverse
  h = {}
  keys.each{|i| h[i] = i - 1}
  (2..m).each{|i|
    if h[i] > h[i - 1] # このときiは素数
      hp = h[i - 1]
      i2 = i * i
      keys.each{|j|
        break if j < i2
        h[j] -= h[j / i] - hp
      }
    end
  }
  h[n]
end

def A000720(n)
  (1..n).map{|i| pi(i)}
end
ary = A000720(78)

# OEIS A000720のデータ
ary0 =
[0,1,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,8,8,8,8,9,9,
 9,9,9,9,10,10,11,11,11,11,11,11,12,12,12,12,13,13,
 14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16,
 17,17,18,18,18,18,18,18,19,19,19,19,20,20,21,21,
 21,21,21,21]
# 一致の確認
p ary == ary0

2015年7月7日火曜日

150707(3)

Ruby


素数の個数(2)

e^n (n ≦ 31) 未満の素数の個数を求めてみた。
(実行時間は2時間半ほどかかる。)

0
1
4
8
16
34
79
183
429
1019
2466
6048
14912
37128
93117
234855
595341
1516233
3877186
9950346
25614562
66124777
171141897
443963543
1154106844
3005936865
7842921261
20496470801
53645077679
140599114669
368973074565
969455391690

150707(2)

Ruby


素数の個数(1)

e^n 未満の素数の個数を求めてみた。
オンライン整数列大辞典の
A040014(http://oeis.org/A040014/list)
と比較し、答え合わせしてみる。
(実行時間は35分ほどかかる。)

def pi(n)
  m = Math.sqrt(n).to_i
  keys = (1..m).map{|i| n / i}
  keys += (1..keys[-1] - 1).to_a.reverse
  h = {}
  keys.each{|i| h[i] = i - 1}
  (2..m).each{|i|
    if h[i] > h[i - 1] # このときiは素数
      hp = h[i - 1]
      i2 = i * i
      keys.each{|j|
        break if j < i2
        h[j] -= h[j / i] - hp
      }
    end
  }
  h[n]
end

def A040014(n)
  (1..n).map{|i| pi(Math.exp(i).to_i)}.unshift(0)
end
ary = A040014(29)

# OEIS A040014のデータ
ary0 =
[0,1,4,8,16,34,79,183,429,1019,2466,6048,14912,
 37128,93117,234855,595341,1516233,3877186,9950346,
 25614562,66124777,171141897,443963543,1154106844,
 3005936865,7842921261,20496470801,53645077679,
 140599114669]
# 一致の確認
p ary == ary0

p ary0
# 大雑把に求めた近似値
p (1..29).map{|i| (Math.exp(i) / i).to_i}.unshift('?')

出力結果
true
[0, 1, 4, 8, 16, 34, 79, 183, 429, 1019, 2466, 6048, 14912, 37128, 93117, 234855
, 595341, 1516233, 3877186, 9950346, 25614562, 66124777, 171141897, 443963543, 1
154106844, 3005936865, 7842921261, 20496470801, 53645077679, 140599114669]
["?", 2, 3, 6, 13, 29, 67, 156, 372, 900, 2202, 5443, 13562, 34031, 85900, 21793
4, 555381, 1420879, 3647776, 9393805, 24258259, 62800749, 162950583, 423687106,
1103713422, 2880195973, 7528061901, 19705490392, 51652038010, 135563251625]

150707

Ruby


素数の和(1)

10^n 未満の素数の和を求めてみた。
オンライン整数列大辞典の
A046731(http://oeis.org/A046731/list)
と比較し、答え合わせしてみる。
(実行時間の関係で、n = 10^11, … , 10^15 のときの比較は行なっていません。)

def sum_of_primes(n)
  m = Math.sqrt(n).to_i
  keys = (1..m).map{|i| n / i}
  keys += (1..keys[-1] - 1).to_a.reverse
  h = {}
  keys.each{|i| h[i] = i * (i + 1) / 2 - 1}
  (2..m).each{|i|
    if h[i] > h[i - 1] # このときiは素数
      hp = h[i - 1]
      i2 = i * i
        keys.each{|j|
          break if j < i2
          h[j] -= i * (h[j / i] - hp)
        }
    end
  }
  h[n]
end

def A046731(n)
  (1..n).map{|i| sum_of_primes(10 ** i - 1)}.unshift(0)
end
ary = A046731(10)

# OEIS A046731のデータの一部
ary0 =
[0,17,1060,76127,5736396,454396537,37550402023,
 3203324994356,279209790387276,24739512092254535,
 2220822432581729238]
# 一致の確認
p ary == ary0

2015年7月6日月曜日

150706

Ruby


Carmichael number

オンライン整数列大辞典の
A002997(http://oeis.org/A002997/list)
と比較し、答え合わせしてみる。

require 'prime'

def pow(a, m, mod)
  return 1 % mod if m == 0
  k = pow(a, m >> 1, mod)
  k *= k
  return k % mod if m & 1 == 0
  return k * a % mod
end

def is_carmichael?(n)
  return false if n.prime?
  i = 2
  while pow(i, n, n) == i && i < n
    i += 1
  end
  return true if i == n # nはカーマイケル数
  return false
end

def A002997(n)
  (2..n).select{|i| is_carmichael?(i)}
end
ary = A002997(512461)

# OEIS A002997のデータ
ary0 =
[561,1105,1729,2465,2821,6601,8911,10585,15841,
 29341,41041,46657,52633,62745,63973,75361,101101,
 115921,126217,162401,172081,188461,252601,278545,
 294409,314821,334153,340561,399001,410041,449065,
 488881,512461]
# 一致の確認
p ary == ary0