Talk:Magnanimous numbers: Difference between revisions

From Rosetta Code
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)

Revision as of 17:31, 3 October 2021

How to find bigger numbers faster

Just a simple observation.
Numbers with more than 2 digits that end in an even number must have all other digits odd and vice versa, so the sum of of partial numbers can't get even and therefor not prime.The exception ist 101, which sums to 2, the only even prime.

 557: 882428665
 558: 995955112

So one can construct those numbers as base 5 numbers converted to decimal and append the last digit.
dgt[x] = 0..4 goes into even decimal digits 2*dgt[x] -> /0/,2,4,6,8 and than odd digits 2*dgt[x]+1 -> 1,3,5,7,9
/0/ ist tested one digit before.
lets go with total 3 digits, so 2 digits in base 5 in front once even, once odd interpretated

Dgt[0,0] = 00 even :  001, 003, 005, 007, 009 odd 11 :  110, 112, 114, 116, 118
Dgt[0,1] = 02 even :  021, 023, 025, 027, 029 odd 13 :  130, 132, 134, 136, 138
Dgt[0,2] = 04 even :  041, 043, 045, 047, 049 odd 15 :  150, 152, 154, 156, 158
Dgt[0,3] = 06 even :  061, 063, 065, 067, 069 odd 17 :  170, 172, 174, 176, 178
Dgt[0,4] = 08 even :  081, 083, 085, 087, 089 odd 19 :  190, 192, 194, 196, 198
Dgt[1,0] = 20 even :  201, 203, 205, 207, 209 odd 31 :  310, 312, 314, 316, 318
Dgt[1,1] = 22 even :  221, 223, 225, 227, 229 odd 33 :  330, 332, 334, 336, 338
Dgt[1,2] = 24 even :  241, 243, 245, 247, 249 odd 35 :  350, 352, 354, 356, 358
Dgt[1,3] = 26 even :  261, 263, 265, 267, 269 odd 37 :  370, 372, 374, 376, 378
Dgt[1,4] = 28 even :  281, 283, 285, 287, 289 odd 39 :  390, 392, 394, 396, 398
Dgt[2,0] = 40 even :  401, 403, 405, 407, 409 odd 51 :  510, 512, 514, 516, 518
Dgt[2,1] = 42 even :  421, 423, 425, 427, 429 odd 53 :  530, 532, 534, 536, 538
Dgt[2,2] = 44 even :  441, 443, 445, 447, 449 odd 55 :  550, 552, 554, 556, 558
Dgt[2,3] = 46 even :  461, 463, 465, 467, 469 odd 57 :  570, 572, 574, 576, 578
Dgt[2,4] = 48 even :  481, 483, 485, 487, 489 odd 59 :  590, 592, 594, 596, 598
Dgt[3,0] = 60 even :  601, 603, 605, 607, 609 odd 71 :  710, 712, 714, 716, 718
Dgt[3,1] = 62 even :  621, 623, 625, 627, 629 odd 73 :  730, 732, 734, 736, 738
Dgt[3,2] = 64 even :  641, 643, 645, 647, 649 odd 75 :  750, 752, 754, 756, 758
Dgt[3,3] = 66 even :  661, 663, 665, 667, 669 odd 77 :  770, 772, 774, 776, 778
Dgt[3,4] = 68 even :  681, 683, 685, 687, 689 odd 79 :  790, 792, 794, 796, 798
Dgt[4,0] = 80 even :  801, 803, 805, 807, 809 odd 91 :  910, 912, 914, 916, 918
Dgt[4,1] = 82 even :  821, 823, 825, 827, 829 odd 93 :  930, 932, 934, 936, 938
Dgt[4,2] = 84 even :  841, 843, 845, 847, 849 odd 95 :  950, 952, 954, 956, 958
Dgt[4,3] = 86 even :  861, 863, 865, 867, 869 odd 97 :  970, 972, 974, 976, 978
Dgt[4,4] = 88 even :  881, 883, 885, 887, 889 odd 99 :  990, 992, 994, 996, 998

Now one can see, that 101 isn't in there
No big win for 3 digits 20*10 versus 900.
But for 9 digits 5^8*10 vs 1e9-1e8 : 3,906,250 vs 900,000,000 = 1/230.4 reducing a week (168 h ) to 45 min
To speed up prime test of the partial numbers of digits, start with the one of same size in the middle.
Show with 995,955,112

first 
9959     +     55112 =    65071
99595    +      5112 =   104707

995955   +       112 =   996067
995      +    955112 =   956107

9959551  +       12 =   
99       +  5955112 =   

9        + 95955112 =   
99595511 +        2 =   

Especially for large numbers it is important to test first the smallest values, which can be checked in a sieve of erathostenes, to reject them as soon as possible.Notice, to check 1e9-1 = 999,999,999 , you need a sieve of size 1e8+8. I can not imagine, how to get such a large number like the 571.th 97,393,713,331,910 in "acceptable" time.

-Horsth (talk) 17:49, 1 October 2021 (UTC)

It is done :-)
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>

Output:
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

--Horst.h (talk) 17:30, 3 October 2021 (UTC)