experimental-projects/oeis-research/César Aguilera/Perfect numbers pseudo-relation/prog.sf

#!/usr/bin/ruby

# https://oeis.org/A305936 (recycled)

var M = 5000.primes.grep{.is_mersenne_prime}

M.shift

for p,q in (M.cons(2)) {
    var M1 = (2**p - 1)
    var M2 = (2**q - 1)

    var p1 = M1*(2**(p-1))
    var p2 = M2*(2**(q-1))

    #say [p1, p2, p, q]

    #say (sqrt(p1 * M1 * p2 * M2) / (p1 * M2))
    #(sqrt(496*31*137438691328*524287))/(496*524287)

    #say ("", 2**((p + q - 2) / 2))

    say ("2^((#{q} - #{p})/2) = ", 2**((q - p)/2))

    assert_eq((isqrt(p1 * M1 * p2 * M2) / (p1 * M2)),  2**((q - p)/2))
    assert_eq(isqrt((2**p - 1)**2 * 2**(p-1) * (2**q - 1)**2 * 2**(q-1)) / ((2**p - 1) * 2**(p-1) * (2**q - 1)),  2**((q - p)/2))
}


#sqrt((2^p - 1)^2 * 2^(p-1) * (2^q - 1)^2 * 2^(q-1)) / ((2^p - 1) * (2^q - 1) * 2^(q-1))

__END__
2^((5 - 3)/2) = 2
2^((7 - 5)/2) = 2
2^((13 - 7)/2) = 8
2^((17 - 13)/2) = 4
2^((19 - 17)/2) = 2
2^((31 - 19)/2) = 64
2^((61 - 31)/2) = 32768
2^((89 - 61)/2) = 16384
2^((107 - 89)/2) = 512
2^((127 - 107)/2) = 1024
2^((521 - 127)/2) = 200867255532373784442745261542645325315275374222849104412672
2^((607 - 521)/2) = 8796093022208
2^((1279 - 607)/2) = 139984046386112763159840142535527767382602843577165595931249318810236991948760059086304843329475444736
2^((2203 - 1279)/2) = 11908525658859223294760121268437066290850060053501019099651935423375594096449911575776314174894302258147533153997065059263030913083222523904
2^((2281 - 2203)/2) = 549755813888
2^((3217 - 2281)/2) = 762145642166990290864647761179972242614403843424065222377723867096038022172794340849684107193235344521442121855812163792833978437326241529856
2^((4253 - 3217)/2) = 858099707516326214372737599885174152158679412517913176174307932398192897924707006515319955082681819372162038923935107254640248499964580476571753536389382144
2^((4423 - 4253)/2) = 38685626227668133590597632