Magnanimous numbers: Difference between revisions

=={{header|Pascal}}==

{{works with|Free Pascal}}

Version nearly like on Talk. Generating the sieve for primes takes most of the time.<br>

found all til 564 : 9151995592 in 0.826 s.

<lang>program Magnanimous;

//Magnanimous Numbers

//algorithm find only numbers where all digits are even except the last

//or where all digits are odd except the last

//Magnanimous Numbers that can not be found by this algorithm

//0,1,11,20,101,1001

{$IFDEF FPC}

{$MODE DELPHI}

{$Optimization ON,ALL}

{$CODEALIGN proc=16}

{$ENDIF}

uses

strUtils,SysUtils;

const

MaxLimit = 1000*1000*1000;

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..571] of Uint64;

var

{$ALIGN 32}

dgtBase5,

dgtEvenBase10,

dgtOddBase10: tbase;

MagList : tMagList;

primes : tprimes;

pPrimes0 : pByte;

T0: int64;

HighIdx,lstIdx, cnt,num,MagIdx: NativeUint;

procedure InsertSort(pMag:pUint64; Left, Right : NativeInt );

var

I, J: NativeInt;

Pivot : Uint64;

begin

for i:= 1 + Left to Right do

begin

Pivot:= pMag[i];

j:= i - 1;

while (j >= Left) and (pMag[j] > Pivot) do

begin

pMag[j+1]:=pMag[j];

Dec(j);

end;

pMag[j+1]:= pivot;

end;

end;

procedure InitPrimes;

const

smallprimes :array[0..5] of byte = (2,3,5,7,11,13);

var

pPrimes : pByte;

p,i,j,l : NativeUint;

begin

l := 1;

for j := 0 to High(smallprimes) do

l*= smallprimes[j];

//scale primelimit to be multiple of l

i :=((MaxLimit-1) DIV l+1)*l+1;//+1 should suffice

setlength(primes,i);

pPrimes := @primes[0];

for j := 0 to High(smallprimes) do

begin

p := smallprimes[j];

i := p;

if j <> 0 then

p +=p;

while i <= l do

begin

pPrimes[i] := 1;

inc(i,p)

end;

end;

//turn the rime wheel

for p := length(primes) div l -1 downto 1 do

move(pPrimes[1],pPrimes[p*l+1],l);

l := High(primes);

//reinsert smallprimes

for j := 0 to High(smallprimes) do

pPrimes[smallprimes[j]] := 0;

pPrimes[1]:=1;

pPrimes[0]:=1;

p := smallprimes[High(smallprimes)];

repeat

repeat

inc(p)

until pPrimes[p] = 0;

i := p*p;

j := 2*p;

if i > l then

BREAK;

while i<= l do

begin

pPrimes[i] := 1;

inc(i,p)

end;

until false;

pPrimes0 := pPrimes;

end;

procedure OutBase5;

var

pb: tpBaseType;

i : NativeUint;

begin

pb:= @dgtBase5[0];

for i := HighIdx downto 0 do

write(pB[i]:2);

writeln;

end;

function IncDgtBase5:nativeUint;

var

pb: tpBaseType;

num: nativeint;

begin

pb:= @dgtBase5[0];

result := 0;

repeat

num := pb[result] + 1;

if num < 5 then

begin

pb[result] := num;

break;

end;

pb[result] := 0;

Inc(result);

until False;

if HighIdx < result then

begin

HighIdx := result;

pb[result] := 0;

end;

end;

procedure CnvEvenBase10(lastIdx:NativeInt);

var

pdgt : tpBaseType;

idx: nativeint;

begin

pDgt := @dgtEvenBase10[0];

for idx := lastIdx downto 1 do

pDgt[idx] := 2 * dgtBase5[idx];

pDgt[0] := 2 * dgtBase5[0]+1;

end;

procedure CnvOddBase10(lastIdx:NativeInt);

var

pdgt : tpBaseType;

idx: nativeint;

begin

pDgt := @dgtOddBase10[0];

for idx := lastIdx downto 1 do

pDgt[idx] := 2 * dgtBase5[idx] + 1;

pDgt[0] := 2 * dgtBase5[0];

end;

function Base10toNum(var dgtBase10: tBase):NativeUint;

var

i : NativeInt;

begin

Result := 0;

for i := HighIdx downto 0 do

Result := Result * 10 + dgtBase10[i];

end;

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 //not needed for small numbers

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

result := true;

For i := 0 to fac-1 do

begin

n := n*10+dgtBase10[fac-i];

//HighSplitNum[i]:= n;

result := result AND ( pPrimes0[LowSplitNum[i] +n] = 0);

IF not(result) then

break;

end;

end;

begin

T0 := Gettickcount64;

InitPrimes;

T0 -= Gettickcount64;

writeln('getting primes ',-T0 / 1000: 0: 3, ' s');

T0 := Gettickcount64;

fillchar(dgtBase5,SizeOf(dgtBase5),#0);

fillchar(dgtEvenBase10,SizeOf(dgtEvenBase10),#0);

fillchar(dgtOddBase10,SizeOf(dgtOddBase10),#0);

//Magnanimous Numbers that can not be found by this algorithm

MagIdx := 0;

MagList[MagIdx] := 1;inc(MagIdx);

MagList[MagIdx] := 11;inc(MagIdx);

MagList[MagIdx] := 20;inc(MagIdx);

MagList[MagIdx] := 101;inc(MagIdx);

MagList[MagIdx] := 1001;inc(MagIdx);

HighIdx := 0;

lstIdx := 0;

repeat

if dgtBase5[highIdx] <> 0 then

begin

CnvEvenBase10(lstIdx);

num := Base10toNum(dgtEvenBase10);

if isMagn(dgtEvenBase10)then

Begin

MagList[MagIdx] := num;

inc(MagIdx);

// write(num:10,' ');OutBase5;

end;

end;

CnvOddBase10(lstIdx);

num := Base10toNum(dgtOddBase10);

if isMagn(dgtOddBase10) then

Begin

MagList[MagIdx] := num;

inc(MagIdx);

// write(num:10,' ');OutBase5;

end;

lstIdx := IncDgtBase5;

until HighIdx > trunc(ln(MaxLimit)/ln(10));

InsertSort(@MagList[0],0,MagIdx-1);

For cnt := 0 to MagIdx-1 do

writeln(cnt+1:3,' ',MagList[cnt]);

T0 -= Gettickcount64;

writeln(-T0 / 1000: 0: 3, ' s');

readln;

end.

</lang>

{{out}}

<pre style="height:180px">

TIO.RUN

getting primes 22.435 s

// copied last lines 564 9,151,995,592--0.826 s--Real time: 23.476 s User time: 20.033 s Sys. time: 2.904 s CPU share: 97.70 %

1 0

2 1

3 2

4 3

5 4

6 5

7 6

8 7

9 8

10 9

11 11

12 12

13 14

14 16

15 20

16 21

17 23

18 25

19 29

20 30

21 32

22 34

23 38

24 41

25 43

26 47

27 49

28 50

29 52

30 56

31 58

32 61

33 65

34 67

35 70

36 74

37 76

38 83

39 85

40 89

41 92

42 94

43 98

44 101

45 110

46 112

47 116

48 118

49 130

50 136

51 152

52 158

53 170

54 172

55 203

56 209

57 221

58 227

59 229

60 245

61 265

62 281

63 310

64 316

65 334

66 338

67 356

68 358

69 370

70 376

71 394

72 398

73 401

74 403

75 407

76 425

77 443

78 449

79 467

80 485

81 512

82 518

83 536

84 538

85 554

86 556

87 574

88 592

89 598

90 601

91 607

92 625

93 647

94 661

95 665

96 667

97 683

98 710

99 712

100 730

101 736

102 754

103 772

104 776

105 790

106 794

107 803

108 809

109 821

110 845

111 863

112 881

113 889

114 934

115 938

116 952

117 958

118 970

119 974

120 992

121 994

122 998

123 1001

124 1112

125 1130

126 1198

127 1310

128 1316

129 1598

130 1756

131 1772

132 1910

133 1918

134 1952

135 1970

136 1990

137 2209

138 2221

139 2225

140 2249

141 2261

142 2267

143 2281

144 2429

145 2447

146 2465

147 2489

148 2645

149 2681

150 2885

151 3110

152 3170

153 3310

154 3334

155 3370

156 3398

157 3518

158 3554

159 3730

160 3736

161 3794

162 3934

163 3974

164 4001

165 4027

166 4063

167 4229

168 4247

169 4265

170 4267

171 4427

172 4445

173 4463

174 4643

175 4825

176 4883

177 5158

178 5176

179 5374

180 5516

181 5552

182 5558

183 5594

184 5752

185 5972

186 5992

187 6001

188 6007

189 6067

190 6265

191 6403

192 6425

193 6443

194 6485

195 6601

196 6685

197 6803

198 6821

199 7330

200 7376

201 7390

202 7394

203 7534

204 7556

205 7592

206 7712

207 7934

208 7970

209 8009

210 8029

211 8221

212 8225

213 8801

214 8821

215 9118

216 9172

217 9190

218 9338

219 9370

220 9374

221 9512

222 9598

223 9710

224 9734

225 9752

226 9910

227 11116

228 11152

229 11170

230 11558

231 11930

232 13118

233 13136

234 13556

235 15572

236 15736

237 15938

238 15952

239 17716

240 17752

241 17992

242 19972

243 20209

244 20261

245 20861

246 22061

247 22201

248 22801

249 22885

250 24407

251 26201

252 26285

253 26881

254 28285

255 28429

256 31370

257 31756

258 33118

259 33538

260 33554

261 35116

262 35776

263 37190

264 37556

265 37790

266 37930

267 39158

268 39394

269 40001

270 40043

271 40049

272 40067

273 40427

274 40463

275 40483

276 42209

277 42265

278 44009

279 44443

280 44447

281 46445

282 48089

283 48265

284 51112

285 53176

286 53756

287 53918

288 55516

289 55552

290 55558

291 55576

292 55774

293 57116

294 57754

295 60007

296 60047

297 60403

298 60443

299 60667

300 62021

301 62665

302 64645

303 66667

304 66685

305 68003

306 68683

307 71536

308 71572

309 71716

310 71752

311 73156

312 75374

313 75556

314 77152

315 77554

316 79330

317 79370

318 80009

319 80029

320 80801

321 80849

322 82265

323 82285

324 82825

325 82829

326 84265

327 86081

328 86221

329 88061

330 88229

331 88265

332 88621

333 91792

334 93338

335 93958

336 93994

337 99712

338 99998

339 111112

340 111118

341 111170

342 111310

343 113170

344 115136

345 115198

346 115772

347 117116

348 119792

349 135158

350 139138

351 151156

352 151592

353 159118

354 177556

355 193910

356 199190

357 200209

358 200809

359 220021

360 220661

361 222245

362 224027

363 226447

364 226681

365 228601

366 282809

367 282881

368 282889

369 311156

370 319910

371 331118

372 333770

373 333994

374 335156

375 339370

376 351938

377 359794

378 371116

379 373130

380 393554

381 399710

382 400049

383 404249

384 408049

385 408889

386 424607

387 440843

388 464447

389 484063

390 484445

391 486685

392 488489

393 515116

394 533176

395 551558

396 559952

397 595592

398 595598

399 600881

400 602081

401 626261

402 628601

403 644485

404 684425

405 686285

406 711512

407 719710

408 753316

409 755156

410 773554

411 777712

412 777776

413 799394

414 799712

415 800483

416 802061

417 802081

418 804863

419 806021

420 806483

421 806681

422 822265

423 864883

424 888485

425 888601

426 888643

427 911390

428 911518

429 915752

430 931130

431 975772

432 979592

433 991118

434 999994

435 1115756

436 1137770

437 1191518

438 1197370

439 1353136

440 1379930

441 1533736

442 1593538

443 1711576

444 1791110

445 1795912

446 1915972

447 1951958

448 2000221

449 2008829

450 2442485

451 2604067

452 2606647

453 2664425

454 2666021

455 2828809

456 2862445

457 3155116

458 3171710

459 3193198

460 3195338

461 3195398

462 3315358

463 3373336

464 3573716

465 3737534

466 3751576

467 3939118

468 4000483

469 4408603

470 4468865

471 4488245

472 4644407

473 5115736

474 5357776

475 5551376

476 5579774

477 5731136

478 5759594

479 5959774

480 6462667

481 6600227

482 6600443

483 6608081

484 6640063

485 6640643

486 6824665

487 6864485

488 6866683

489 7113710

490 7133110

491 7139390

492 7153336

493 7159172

494 7311170

495 7351376

496 7719370

497 7959934

498 7979534

499 8044009

500 8068201

501 8608081

502 8844449

503 9171170

504 9777910

505 9959374

506 11771992

507 13913170

508 15177112

509 17115116

510 19337170

511 19713130

512 20266681

513 22086821

514 22600601

515 22862885

516 26428645

517 28862465

518 33939518

519 37959994

520 40866083

521 44866043

522 48606043

523 48804809

524 51137776

525 51513118

526 53151376

527 53775934

528 59593574

529 60402247

530 60860603

531 62202281

532 64622665

533 66864625

534 66886483

535 71553536

536 77917592

537 82486825

538 86842265

539 91959398

540 95559998

541 117711170

542 222866845

543 228440489

544 244064027

545 280422829

546 331111958

547 400044049

548 460040803

549 511151552

550 593559374

551 606202627

552 608844043

553 622622801

554 622888465

555 773719910

556 844460063

557 882428665

558 995955112

559 1777137770

560 2240064227

561 2444402809

562 5753779594

563 6464886245

564 9151995592

0.826 s

Real time: 23.476 s User time: 20.033 s Sys. time: 2.904 s CPU share: 97.70 %</pre>