Talk:Magnanimous numbers: Difference between revisions
Content added Content deleted
(→How to find bigger numbers faster: show/explain how it works) |
(→How to find bigger numbers faster: get up to 571 : 97,393,713,331,910 in 9.5 min) |
||
Line 57: | Line 57: | ||
I can not imagine, how to get such a large number like the 571.th 97,393,713,331,910 in "acceptable" time. |
I can not imagine, how to get such a large number like the 571.th 97,393,713,331,910 in "acceptable" time. |
||
-[[User:Horsth|Horsth]] ([[User talk:Horsth|talk]]) 17:49, 1 October 2021 (UTC) |
-[[User:Horsth|Horsth]] ([[User talk:Horsth|talk]]) 17:49, 1 October 2021 (UTC) |
||
It is done :-)<br> |
|||
modified using gmp |
|||
<lang pascal> |
|||
uses |
|||
strUtils,SysUtils,gmp; |
|||
const |
|||
MaxLimit = 10*1000*1000*1000+8; |
|||
MAXHIGHIDX = 13;//trunc(ln(MaxLimit)/ln(10)); |
|||
type |
|||
tprimes = array of byte; |
|||
tBaseType = word; |
|||
tpBaseType = pWord; |
|||
tBase =array[0..15] of tBaseType; |
|||
tNumType = NativeUint; |
|||
tSplitNum =array[0..15] of tNumType; |
|||
tMagList = array[0..1000] of Uint64; |
|||
var |
|||
{$ALIGN 32} |
|||
dgtBase5, |
|||
dgtEvenBase10, |
|||
dgtOddBase10: tbase; |
|||
MagList : tMagList; |
|||
primes : tprimes; |
|||
z : mpz_t;// mpz for testing probably prime |
|||
pPrimes0 : pByte; |
|||
T0: int64; |
|||
HighIdx,lstIdx, cnt,num,MagIdx: NativeUint; |
|||
.... |
|||
function isMagn(var dgtBase10: tBase):boolean; |
|||
//split number into sum of all "partitions" of digits |
|||
//check if sum is always prime |
|||
//1234 -> 1+234,12+34 ; 123+4 |
|||
var |
|||
LowSplitNum{,HighSplitNum} :tSplitNum; |
|||
i,fac,n: NativeInt; |
|||
Begin |
|||
n := 0; |
|||
fac := 1; |
|||
For i := 0 to HighIdx-1 do |
|||
begin |
|||
n := fac*dgtBase10[i]+n; |
|||
fac *=10; |
|||
LowSplitNum[HighIdx-1-i] := n; |
|||
end; |
|||
n := 0; |
|||
fac := HighIdx; |
|||
For i := 0 to fac-1 do |
|||
begin |
|||
n := n*10+dgtBase10[fac-i]; |
|||
// HighSplitNum[i]:= n; |
|||
LowSplitNum[i] += n; |
|||
end; |
|||
// check numbers within prime sieve |
|||
result := true; |
|||
For i := 0 to fac-1 do |
|||
begin |
|||
n := LowSplitNum[i]; |
|||
if n <=MAXLIMIT then |
|||
result := result AND (pPrimes0[n] = 0); |
|||
if NOT(result) then |
|||
EXIT; |
|||
end; |
|||
//only check remaining numbers out of prime sieve if needed |
|||
For i := 0 to fac-1 do |
|||
begin |
|||
n := LowSplitNum[i]; |
|||
if n >MAXLIMIT then |
|||
Begin |
|||
mpz_set_ui(z,n); |
|||
result := result AND (mpz_probab_prime_p(z,1) >0); |
|||
if NOT(result) then |
|||
EXIT; |
|||
end; |
|||
end; |
|||
end; |
|||
... |
|||
BEGIN |
|||
mpz_init_set_ui(z,0); |
|||
T0 := Gettickcount64; |
|||
InitPrimes; |
|||
... |
|||
InsertSort(@MagList[0],0,MagIdx-1); |
|||
mpz_clear(z); |
|||
//Output format of the list of https://oeis.org/A252996/b252996.txt |
|||
For cnt := 0 to MagIdx-1 do |
|||
writeln(cnt+1:3,' ',MagList[cnt]); |
|||
T0 -= Gettickcount64; |
|||
writeln(-T0 / 1000: 0: 3, ' s'); |
|||
{$IFDEF WINDOWS} |
|||
readln; |
|||
{$ENDIF} |
|||
end. |
|||
</lang> |
|||
{{out}} |
|||
<pre> |
|||
getting primes 27.061 s .. |
|||
,,, |
|||
559 1777137770 |
|||
560 2240064227 |
|||
561 2444402809 |
|||
562 5753779594 |
|||
563 6464886245 |
|||
564 9151995592 |
|||
565 22226226625 |
|||
566 31993717930 |
|||
567 39393115598 |
|||
568 46884486265 |
|||
569 84448000009 |
|||
570 5391391551358 |
|||
571 97393713331910 |
|||
536.366 s |
|||
real 9m23,751s user 9m21,133s sys 0m2,293s |
|||
</pre> |
|||
--[[User:Horst.h|Horst.h]] ([[User talk:Horst.h|talk]]) 17:30, 3 October 2021 (UTC) |