while True:
    p = next_prime(random.randint(0, 10**500))
    if len(str(p)) != 500:
        continue
    q = Integer(int(str(p)[250:] + str(p)[:250]))
    if q.is_prime():
        break
        

>> p * q
6146024643941503757217715363256725297474582575057128830681803952150464985329239705861504172069973746764596350359462277397739134788481500502387716062571912861345331755396960400668616401300689786263797654804338789112750913548642482662809784602704174564885963722422299918304645125966515910080631257020529794610856299507980828520629245187681653190311198219403188372517508164871722474627810848320169613689716990022730088459821267951447201867517626158744944551445617408339432658443496118067189012595726036261168251749186085493288311314941584653172141498507582033165337666796171940245572657593635107816849481870784366174740265906662098222589242955869775789843661127411493630943226776741646463845546396213149027737171200372484413863565567390083316799725434855960709541328144058411807356607316377373917707720258565704707770352508576366053160404360862976120784192082599228536166245480722359263166146184992593735550019325337524138545418186493193366973466749752806880403086988489013389009843734224502284325825989
>> pow(m, 65537, p * q)
3572030904528013180691184031825875018560018830056027446538585108046374607199842488138228426133620939067295245642162497675548656988031367698701161407333098336631469820625758165691216722102954230039803062571915807926805842311530808555825502457067483266045370081698397234434007948071948000301674260889742505705689105049976374758307610890478956315615270346544731420764623411884522772647227485422185741972880247913540403503772495257290866993158120540920089734332219140638231258380844037266185237491107152677366121632644100162619601924591268704611229987050199163281293994502948372872259033482851597923104208041748275169138684724529347356731689014177146308752441720676090362823472528200449780703866597108548404590800249980122989260948630061847682889941399385098680402067366390334436739269305750501804725143228482932118740926602413362231953728010397307348540059759689560081517028515279382023371274623802620886821099991568528927696544505357451279263250695311793770159474896431625763008081110926072287874375257

If we write pp as:

p=10250x+yp = 10^{250}x+y

with 0≤x,y<102500\leq x,y<10^{250}, then qq is obtained by swapping xx and yy:

q=10250y+xq = 10^{250}y+x

This means that for n=pqn=pq we have:

n=pq=(10250x+y)(10250y+x)=n = pq = (10^{250}x+y)(10^{250}y+x) = =10500xy+10250x2+10250y2+xy=(10500+1)xy+10250(x2+y2)=10^{500}xy + 10^{250}x^2 + 10^{250}y^2 + xy = (10^{500}+1)xy+10^{250}(x^2+y^2)

But now n≡0xy+10250(x2+y2)mod  10500+1n \equiv 0xy + 10^{250}(x^2+y^2) \mod 10^{500}+1 or x2+y2≡n10250mod  10500+1 .x^2+y^2 \equiv \frac{n}{10^{250}} \mod 10^{500}+1\,. As 1025010^{250} and 10500+110^{500}+1 are coprime, this is well-defined.

In sage:

sage: R=Zmod(10^500+1)
sage: s = Integer(R(n)*R(10^250)^-1)

But on the other hand, as xx and yy are less than 1025010^{250}, the sum of their squares must be less than (10250)2+(10250)2=2⋅10500(10^{250})^2+(10^{250})^2=2\cdot 10^{500}. As we already know the residue of x2+y2x^2+y^2 mod 10500+110^{500}+1, this means that we only have two possibilities left: ss and s+10500+1s+10^{500}+1. It is possible to quite easily exclude the first one: ss is divisible by 33 but not 99, which would be in contradiction to the Sum of two squares theorem. Of course, one can also just try both values.

This means that x2+y2x^2+y^2 must be s+10500+1s+10^{500}+1. Using n=(10500+1)xy+10250(x2+y2)n=(10^{500}+1)xy+10^{250}(x^2+y^2) we now have two equations for two variables. Those can be solved by substituting and solving the resulting biquadratic equation, or just using sage:

sage: x,y=var("x,y")
sage: solve([n==(10^500+1)*x*y+10^250*(x^2+y^2),x^2+y^2==s+10^500+1],(x,y))
[[x == 6704087713997099865507815769768764401477034704508108261256143985284139367993820913412851771729239540206881848078387952673084087851136091482870695733338884807660355569782830813482768223955545908153884900037741101021721778447622913463893844919116113159, y == 9167577910876006891858257597237629440949705918335724259091862313200766026122707397405770653228405765712936180484223756343702067850288957263434298651591281476391443835610092615155677011406076889438184185284410734629391841138106293296764467469014716371], ...]

(symmetric and negative solutions omitted). Now it’s just a matter of recovering pp and qq and deciphering:

sage: p=10^250*x+y
sage: q=10^250*y+x
sage: p*q==n
True
sage: d=ZZ.quo((p-1)*(q-1))(65537)**-1
sage: m=int(pow(c,d,n))
sage: binascii.a2b_hex(hex(m)[2:-1])
'midnight{w3ll_wh47_d0_y0u_kn0w_7h15_15_4c7u4lly_7h3_w0rld5_l0n6357_fl46_4nd_y0u_f0und_17_50_y0u_5h0uld_b3_pr0ud_0f_y0ur53lf_50_uhmm_60_pr1n7_7h15_0n_4_75h1r7_0r_50m37h1n6_4nd_4m4z3_7h3_p30pl3_4r0und_y0u}\x99\xd3\x84\x00\xdd\x98\xbf\xedv\xe8|\xd3#\x01sR\x83\x9f\xce\x9fHg\xef\xfb\x07\x05\xc7\xd1R4*\xbc\x9e`\x9aW\x10\x0b5\xc0\xc0\xf9\xcc2dD\x00\xd5\x89\xfd2\xd2l\xe3N3\x8bU6}[\x92\xd5\xf5\x0fO\xde-\xf1\xb0\xbe\xaf\xc5\xcfM\xadyo\xd9\xbf\xff\xeci\xd5$\xf99\xde\xad\xaaP{\xfc\xf7\x91 o2\x99M\x9cE\xfe@,\xc0\x8d\x88wU\xb0\x82X\xa2;r\xeaq\x8eV\x05t\x94\xb8\x8c\xba\x90\xf5\xa8\xf9\x17$\x94\xf4:\x11\x9e\xfc\x0b]\x97\xbbMv9\x865f*$\x93r\xf65j\xc0jk\xf0\xab\xd5\t\xb3\xda\x17\xb0~\x05}\xc4@(\xa8\r\x16\x01V\xcdm\x901\x17\x18\xf3\xfd\xd6L\xa7\x13|\xaa\x9d\x1e_\xb4%g/Q+\xff\xc1\xbe\xf1fB#g\xa8\xda\xdd4'