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 $p$ as:

$p = 10^{250}x+y$

with $0\leq x,y<10^{250}$, then $q$ is obtained by swapping $x$ and $y$:

$q = 10^{250}y+x$

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

$n = pq = (10^{250}x+y)(10^{250}y+x) =$ $=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 \equiv 0xy + 10^{250}(x^2+y^2) \mod 10^{500}+1$ or $x^2+y^2 \equiv \frac{n}{10^{250}} \mod 10^{500}+1\,.$ As $10^{250}$ and $10^{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 $x$ and $y$ are less than $10^{250}$, the sum of their squares must be less than $(10^{250})^2+(10^{250})^2=2\cdot 10^{500}$. As we already know the residue of $x^2+y^2$ mod $10^{500}+1$, this means that we only have two possibilities left: $s$ and $s+10^{500}+1$. It is possible to quite easily exclude the first one: $s$ is divisible by $3$ but not $9$, which would be in contradiction to the Sum of two squares theorem. Of course, one can also just try both values.

This means that $x^2+y^2$ must be $s+10^{500}+1$. Using $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 $p$ and $q$ 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'
`