Av: Didrik Pemmer Aalen
Et Harshadtall er definert som et positivt tall som er delelig med summen av alle sifferne i tallet. I noen tilfeller er summen av alle sifferne i tallet et primtall. Disse kaller vi Harshadprimtall.
1729 er et Harshadtall fordi 1 + 7 + 2 + 9 = 19 og 1729 % 19 = 0. Dette er også et Harshadprimtall, fordi 19 er et primtall.
1730 er ikke et Harshadtall fordi 1 + 7 + 3 + 0 = 11 og 1730 % 11 = 3.
Hvor mange tall fra 1 til og med 98765432 er Harshadprimtall?
Høres ut som vi har en ny rekord? Veldig kult! Liker den getCrosssumChunked_fastmod_unrolled, det ble en veldig effektiv og konsis implementasjon. Perfekt også at den kan ta 5 sifre om gangen, for da kan den jo egentlig skalere opp vilkårlig langt, bare til kosten av en fastmod +addisjon for hver gruppe med 5 siffer. Har lyst til å tilbake til noen av de forrige tverrsumoppgavene og gjøre samme optimaliseringer på dem.
Jeg er bare på telefonen nå men skal se om jeg får noe tid til å se på det etterpå.
Jeg har fortsatt én idé til jeg vil teste ut, har håp om at vi kan barbere bort enda noen prosent...
Godt nyttår folkens!
Det må, i rettferdighetens navn, sies at vi har spekulert en drøy uke på dette mens @DarioSucic sin CUDA-løsningn var mer spontan(?), så vi skal ikke bli for kjepphøye. 😄
@michaelo, jeg prøvde å få denne til å virke på en Raspberry Pi 3, men den hadde en litt gammel clang (3.8.1) og i tillegg finnes det ikke noen 128-bits integer. Det er, så vidt jeg ser, ingen ting i koden som bruker int128, men jeg må endre en del i fastmod.h (eller google en Raspberry-versjon) når vi finner dette:
// fastmod computes (a % d) given precomputed M
FASTMOD_API uint32_t fastmod_u32(uint32_t a, uint64_t M, uint32_t d) {
uint64_t lowbits = M * a;
return (uint32_t)(mul128_u32(lowbits, d));
}og
FASTMOD_API uint64_t mul128_u32(uint64_t lowbits, uint32_t d) {
return ((__uint128_t)lowbits * d) >> 64;
}Ettersom vi begrenser oss til 32-bits (signed) integer selv på det store datasettet skulle dette være grei skuring å skrive om til ren 64-bits integer, men det ble for mye god drikke i går kveld. Man kjenner sin begrensning etter mange års erfaring. 😉
Edit: Det er visst ikke så rett fram likevel...
Det man gjør her er vel egentlig en 64bit * 32bit multiplikasjon og henter ut de øverste 64 bitene. På en del plattformer finnes det en egen instruksjon som gjør 64bit * 64bit og returnerer de øverste bitene samtidig som de laveste bitene havner i et register, og jeg tror det er den de prøver å trigge her. Men det kan virke som Armv8 mangler denne:
A 64 × 64 to 128-bit multiply requires a sequence of two instructions to generate a pair of 64-bit result registers:
± (64 × 64) gives the lower 64 bits of the result [63:0].
(64 × 64) gives the higher 64 bits of the result [127:64].
Men vi trenger vel her bare de øverste 64 bitene, så da tror jeg instruksjonen UMULH gjør det du trenger uten behov for noe uint128. Så spørs det om det finnes noen compiler intrinsic for denne eller om du må ty til inline assembly?
EDIT: Ser ut som __uint128_t gjør nettopp dette og burde vært støttet på både GCC og clang (men kanskje din er for gammel?)
Se mer diskusjon rundt temaet her:
https://stackoverflow.com/questions/28868367/getting-the-high-part-of-64-bit-integer-multiplication/51587262
Godt nyttår folkens!
Det må, i rettferdighetens navn, sies at vi har spekulert en drøy uke på dette mens @DarioSucic sin CUDA-løsningn var mer spontan(?), så vi skal ikke bli for kjepphøye. 😄
Godt nyttår ja! 🍾
Enig, på ingen måte meningen å være kjepphøye, det har tatt mye jobb og mange iterasjoner å komme hit, mens @DarioSucic kom opp med sin løsning i Cuda på presumptivt noen minutter, uten å egentlig kunne cuda 😄
Likevel synes jeg det er spennende at vi nærmer oss den løsningen i kjøretid til tross for at spådommen tidlig var at det var nyttesløst å konkurrere mot slike resultater! Om jeg hadde hatt en moderne entusiast-nivå desktop-cpu med 12 eller 16 fysiske kjerner ville vi antakeligvis være forbi cuda-løsningen på ytelse nå, det synes jeg er gøy å tenke på 😄
Det man gjør her er vel egentlig en 64bit * 32bit multiplikasjon og henter ut de øverste 64 bitene.
Yes, jeg kom til samme konklusjon. En 64 bits left-shit var en "dead giveaway".
Men vi trenger vel her bare de øverste 64 bitene, så da tror jeg instruksjonen UMULH gjør det du trenger uten behov for noe uint128. Så spørs det om det finnes noen compiler intrinsic for denne eller om du må ty til inline assembly?
Dette kommer godt fram i fastmod.h:
#ifdef _MSC_VER
// __umulh is only available in x64 mode under Visual Studio: don't compile to 32-bit!
FASTMOD_API uint64_t mul128_u32(uint64_t lowbits, uint32_t d) {
return __umulh(lowbits, d);
}
#else // _MSC_VER NOT definedMin Raspberry har en armv7 (rev 4, altså v7l) og den er ikke 64-bits, og den mangler UMULH. Det kan jo fikses i software (shift og ta vare på carry), men det blir en omfattende operasjon og da spørs det om vinningen går opp i spinningen.
En liten test-funksjon:
uint64_t test(uint64_t a , uint64_t b) {
return a % b;
}Blir til:
(gcc)
0: e92d4010 push {r4, lr}
4: ebfffffe bl 0 <__aeabi_uldivmod>
8: e1a00002 mov r0, r2
c: e1a01003 mov r1, r3
10: e8bd8010 pop {r4, pc}(clang)
0: e92d4800 push {fp, lr}
4: e1a0b00d mov fp, sp
8: ebfffffe bl 0 <__umoddi3>
c: e8bd8800 pop {fp, pc}Som vi ser så har ikke ARM v71 støtte for mod (eller div) for 64-bits heltall, noe som ikke bør komme som et sjokk når vi er på en 32-bits platform. Ser ikke ut som ARM støtter divisjon av mindre heltall heller (når jeg tester tilsvarende med clang), men det kan være fordi ARM før Cortext faktisk manglet divisjon i hardware, og at min clang er så gammel at den tror det fremdeles er tilfellet (eller mener den er smartere enn ARM Cortex). Det er egentlig ganske forvirrende om ARM har hardware-divisjon eller ikke. Division on ARM Cores
__umoddi3.c
COMPILER_RT_ABI du_int
--
__umoddi3(du_int a, du_int b)
{
du_int r;
__udivmoddi4(a, b, &r);
return r;
}Kildekoden til __udivmoddi4
EDIT: Ser ut som __uint128_t gjør nettopp dette og burde vært støttet på både GCC og clang (men kanskje din er for gammel?)
Yes, men den glimrer med sitt fravær dessverre. Lurer på om dette faktisk krever 64-bits platform. Det nok en god grunn for at ARMv7 ikke støtter 128-bits integer. Spørs om det faktisk blir for få registre (ARM har registre i bøtter og spann) når et register i ARMv7 er 32 bits.
Det man gjør her er vel egentlig en 64bit * 32bit multiplikasjon og henter ut de øverste 64 bitene.
Yes, jeg kom til samme konklusjon. En 64 bits left-shit var en "dead giveaway".
Right shift var det vel strengt tatt 🔪 🌲 🤣
EDIT: Ser ut som __uint128_t gjør nettopp dette og burde vært støttet på både GCC og clang (men kanskje din er for gammel?)
Yes, men den glimrer med sitt fravær dessverre. Lurer på om dette faktisk krever 64-bits platform. Det nok en god grunn for at ARMv7 ikke støtter 128-bits integer. Spørs om det faktisk blir for få registre (ARM har registre i bøtter og spann) når et register i ARMv7 er 32 bits.
Jeg guglet bare raskt Raspberry Pi 3 og fikk opp en variant med 64-bit Cortex v8 (Pi 3 modell B+) så jeg antok at den CPU-en var standard for alle modellene. Det er nok bare 64-bit CPU-er som har 64->128 multiplikasjon i hardware ja.
Right shift var det vel strengt tatt
😄 Right!
(selv om det spiller liten rolle hvilken vei man skyver 64 bits, de forsvinner uansett)
Jeg guglet bare raskt Raspberry Pi 3 og fikk opp en variant med 64-bit Cortex v8 (Pi 3 modell B+) så jeg antok at den CPU-en var standard for alle modellene. Det er nok bare 64-bit CPU-er som har 64->128 multiplikasjon i hardware ja.
Hmm... Jeg var nesten 100% sikker på at jeg hadde en sånn en. Og jeg måtte ta bilde av CPU-en (gammel, svaksynt mann bruker telefonens kamera som forstørrelsesglass). Der står det BROADCOM BCM2837RIFBG. Det er en "Raspberry Pi 3 with 64-bit quad-core SoC, built-in Wi-Fi and Bluetooth".
Har jeg prestert å installere 32-bits Linux på den???
...etter litt googling... Ja, det ser ut som om Raspbian er 32-bits. "By default, the RPi3 chooses the 32-bit variation" Og videre: "The RPi3 constrains itself so much that it acts as if the provided 32-bit mode is based on ARMv7 architecture rather than the built-in ARMv8"
Det var da som svarte!
😄
(selv om det spiller liten rolle hvilken vei man skyver 64 bits, de forsvinner uansett)
Når du har 128 til å begynne med så vil jeg påstå det spiller ganske stor rolle for hva du sitter igjen med etter skyvingen ;)
(det å caste en uint64_t til uint128_t og så bitshifte 64 steg for så å caste til uint64_t igjen høres jo unektelig ut som en NOP, men her er jo nettopp det å caste til uint128_t og RIGHT shifte 64 bits en måte å signalisere til kompilatoren at her vil jeg at multiplikasjonen skal gjøres som 64->128 og jeg vil ha de øverste 64 bittene tilbake av de 128 som genereres i multiplikasjonen)
Har jeg prestert å installere 32-bits Linux på den???
...etter litt googling... Ja, det ser ut som om Raspbian er 32-bits. "By default, the RPi3 chooses the 32-bit variation" Og videre: "The RPi3 constrains itself so much that it acts as if the provided 32-bit mode is based on ARMv7 architecture rather than the built-in ARMv8"
Det var da som svarte!
😄
Høres ut som du har flere prosjekter å ta tak i de neste dagene 😄
Jeg har fått inn nok en optimalisering her, noe som kutter kjøretiden på min versjon fra ca 9.2s til 6.6s på 1000 kjøringer (bygget med cl) med lite datasett, en besparelse på ca 25%.
Det er en variant av unrollingen for 10 og 10 steg jeg gjorde opprinnelig, men denne unroller flere steg og hopper samtidig også over tall som er delelige med 3.
Jeg vil ikke poste koden her ennå, for jeg har ikke fått inn din optimalisering med constexpr ennå, men det burde kunne gi en ytterligere 1-2% besparelse. Det ser ut til å være noen problemer med constexpr og fastmod i visual studio, jeg må finne ut hva som er greia der.
(selv om det spiller liten rolle hvilken vei man skyver 64 bits, de forsvinner uansett)
Når du har 128 til å begynne med så vil jeg påstå det spiller ganske stor rolle for hva du sitter igjen med etter skyvingen ;)
Det har du helt rett i. Nå var det vel helst det at jeg ikke har 128-bits, og så en >> 64 som var åpenbaringen her. For at >> 64 skal gi noe annet enn 0 må man ha noe som er større enn uint64_t. Jeg så på signaturene og trodde først at det ikke var behov for int128, i tillegg til at datasettet (selv det store) bare har skarve 31 bits, men det ble naturligvis helt feil. Vi trenger ikke 128 bits, men vi trenger mer enn 64 bits. 😄
(det å caste en uint64_t til uint128_t og så bitshifte 64 steg for så å caste til uint64_t igjen høres jo unektelig ut som en NOP,
det å caste en uint64_t til uint128_t og så bitshifte 64 steg for så å caste til uint64_t er en NOP! 😄
return ((__uint128_t)lowbits * d) >> 64;
Men det å caste produktet av en uint64_t og en uint32_t til uint128_t og så bitshifte 64 steg for så å caste til uint64_t er ikke en NOP. 😉
Right-shift er divisjon. Det var derfor jeg opprinnelig sa left-shift fordi jeg trodd det var snakk om multiplikasjon. Og her er det både multiplikasjon og divisjon. Gikk litt fort i svingene der.
Høres ut som du har flere prosjekter å ta tak i de neste dagene
Hmm... Ja. Og jeg som har blitt så glad i Raspbian. Nå må det bli noe annet.
Godt nyttår ja folkens!
Hvor lenge klarer vi å holde denne tråden i live? 😆
Det må, i rettferdighetens navn, sies at vi har spekulert en drøy uke på dette mens @DarioSucic sin CUDA-løsningn var mer spontan(?), så vi skal ikke bli for kjepphøye. 😄
Haha, helt sant! Og skal vi snakke verdi/sekund (både utviklersekunder og CPU/GPU-sekunder) så er det spesielt vanskelig å overgå den!
@michaelo
Har endelig fått tid til å teste din løsning her, og du har gått forbi meg! 👏 🏅 😄
Nåvel. Det er vel en ganske kollektiv effort dette her! Lett skumming av fløta (i en liten periode, ref nedenfor) 😉
@michaelo
Jeg har fått inn nok en optimalisering her, noe som kutter kjøretiden på min versjon fra ca 9.2s til 6.6s på 1000 kjøringer (bygget med cl) med lite datasett, en besparelse på ca 25%.
Pent!! 🙌
Jeg testet litt med pregenerering av tverrsummer (droppet constexpr for å gøre det litt mer "fair", pregenereringen tar jo uansett bare noen mikrosekunder). Med min siste versjon av loopingen ser det faktisk ut som det er mindre effektivt å bruke pregenererte tverrsummer enn å summere kun de sifrene som endrer seg.
Jeg har gjort det konfigurerbart fra kommandolinjen, så dere kan jo teste på deres systemer og se om dere observerer det samme.
Mine resultater er omtrent som følger:
./knowit23 --small --count 1000
cl: 6600ms
gcc: 6870ms
clang: 6750ms
./knowit23 --small --count 1000 --precalc
cl: 6800ms
gcc: 7100ms
clang: 6800ms
#include <iostream>
#include <thread>
#include <algorithm>
#include <inttypes.h>
#include <climits>
#include <math.h>
#include "cxxopts.hpp"
#include "fastmod.h"
typedef uint64_t digit_t;
bool* buildPrimeSieve(size_t maxPrime, bool inverted = true)
{
auto maxSquareRoot = (int)sqrt(maxPrime);
bool* eliminated = new bool[maxPrime + 1]{ false };
eliminated[0] = true;
eliminated[1] = true;
for (int j = 4; j <= maxPrime; j += 2)
{
eliminated[j] = true;
}
for (int i = 3; i <= maxPrime; i += 2)
{
if (!eliminated[i])
{
if (i <= maxSquareRoot)
{
for (int j = i * i; j <= maxPrime; j += i)
{
eliminated[j] = true;
}
}
}
}
if (!inverted) {
for (int i = 0; i < maxPrime; i++) eliminated[i] = !eliminated[i];
}
return eliminated;
}
inline void check(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
if (prime[digitsum] && fastmod::is_divisible(i, M[digitsum]))
{
count++;
}
}
inline void check2(uint32_t& count, bool* prime, uint64_t* M, uint32_t i, uint32_t digitsum)
{
if (prime[digitsum] && fastmod::is_divisible(i, M[digitsum]))
{
count++;
}
}
inline void check100rest0(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
check2(count, prime, M, i + 1, digitsum + 1);
check2(count, prime, M, i + 2, digitsum + 2);
check2(count, prime, M, i + 3, digitsum + 3);
check2(count, prime, M, i + 5, digitsum + 5);
check2(count, prime, M, i + 7, digitsum + 7);
check2(count, prime, M, i + 10, digitsum + 1);
check2(count, prime, M, i + 11, digitsum + 2);
check2(count, prime, M, i + 12, digitsum + 3);
check2(count, prime, M, i + 14, digitsum + 5);
check2(count, prime, M, i + 16, digitsum + 7);
check2(count, prime, M, i + 20, digitsum + 2);
check2(count, prime, M, i + 21, digitsum + 3);
check2(count, prime, M, i + 23, digitsum + 5);
check2(count, prime, M, i + 25, digitsum + 7);
check2(count, prime, M, i + 29, digitsum + 11);
check2(count, prime, M, i + 30, digitsum + 3);
check2(count, prime, M, i + 32, digitsum + 5);
check2(count, prime, M, i + 34, digitsum + 7);
check2(count, prime, M, i + 38, digitsum + 11);
check2(count, prime, M, i + 41, digitsum + 5);
check2(count, prime, M, i + 43, digitsum + 7);
check2(count, prime, M, i + 47, digitsum + 11);
check2(count, prime, M, i + 49, digitsum + 13);
check2(count, prime, M, i + 50, digitsum + 5);
check2(count, prime, M, i + 52, digitsum + 7);
check2(count, prime, M, i + 56, digitsum + 11);
check2(count, prime, M, i + 58, digitsum + 13);
check2(count, prime, M, i + 61, digitsum + 7);
check2(count, prime, M, i + 65, digitsum + 11);
check2(count, prime, M, i + 67, digitsum + 13);
check2(count, prime, M, i + 70, digitsum + 7);
check2(count, prime, M, i + 74, digitsum + 11);
check2(count, prime, M, i + 76, digitsum + 13);
check2(count, prime, M, i + 83, digitsum + 11);
check2(count, prime, M, i + 85, digitsum + 13);
check2(count, prime, M, i + 89, digitsum + 17);
check2(count, prime, M, i + 92, digitsum + 11);
check2(count, prime, M, i + 94, digitsum + 13);
check2(count, prime, M, i + 98, digitsum + 17);
}
inline void check100rest1(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
check2(count, prime, M, i + 0, digitsum + 0);
check2(count, prime, M, i + 1, digitsum + 1);
check2(count, prime, M, i + 2, digitsum + 2);
check2(count, prime, M, i + 4, digitsum + 4);
check2(count, prime, M, i + 6, digitsum + 6);
check2(count, prime, M, i + 10, digitsum + 1);
check2(count, prime, M, i + 11, digitsum + 2);
check2(count, prime, M, i + 13, digitsum + 4);
check2(count, prime, M, i + 15, digitsum + 6);
check2(count, prime, M, i + 19, digitsum + 10);
check2(count, prime, M, i + 20, digitsum + 2);
check2(count, prime, M, i + 22, digitsum + 4);
check2(count, prime, M, i + 24, digitsum + 6);
check2(count, prime, M, i + 28, digitsum + 10);
check2(count, prime, M, i + 31, digitsum + 4);
check2(count, prime, M, i + 33, digitsum + 6);
check2(count, prime, M, i + 37, digitsum + 10);
check2(count, prime, M, i + 39, digitsum + 12);
check2(count, prime, M, i + 40, digitsum + 4);
check2(count, prime, M, i + 42, digitsum + 6);
check2(count, prime, M, i + 46, digitsum + 10);
check2(count, prime, M, i + 48, digitsum + 12);
check2(count, prime, M, i + 51, digitsum + 6);
check2(count, prime, M, i + 55, digitsum + 10);
check2(count, prime, M, i + 57, digitsum + 12);
check2(count, prime, M, i + 60, digitsum + 6);
check2(count, prime, M, i + 64, digitsum + 10);
check2(count, prime, M, i + 66, digitsum + 12);
check2(count, prime, M, i + 73, digitsum + 10);
check2(count, prime, M, i + 75, digitsum + 12);
check2(count, prime, M, i + 79, digitsum + 16);
check2(count, prime, M, i + 82, digitsum + 10);
check2(count, prime, M, i + 84, digitsum + 12);
check2(count, prime, M, i + 88, digitsum + 16);
check2(count, prime, M, i + 91, digitsum + 10);
check2(count, prime, M, i + 93, digitsum + 12);
check2(count, prime, M, i + 97, digitsum + 16);
check2(count, prime, M, i + 99, digitsum + 18);
}
inline void check100rest2(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
check2(count, prime, M, i + 0, digitsum + 0);
check2(count, prime, M, i + 1, digitsum + 1);
check2(count, prime, M, i + 3, digitsum + 3);
check2(count, prime, M, i + 5, digitsum + 5);
check2(count, prime, M, i + 9, digitsum + 9);
check2(count, prime, M, i + 10, digitsum + 1);
check2(count, prime, M, i + 12, digitsum + 3);
check2(count, prime, M, i + 14, digitsum + 5);
check2(count, prime, M, i + 18, digitsum + 9);
check2(count, prime, M, i + 21, digitsum + 3);
check2(count, prime, M, i + 23, digitsum + 5);
check2(count, prime, M, i + 27, digitsum + 9);
check2(count, prime, M, i + 29, digitsum + 11);
check2(count, prime, M, i + 30, digitsum + 3);
check2(count, prime, M, i + 32, digitsum + 5);
check2(count, prime, M, i + 36, digitsum + 9);
check2(count, prime, M, i + 38, digitsum + 11);
check2(count, prime, M, i + 41, digitsum + 5);
check2(count, prime, M, i + 45, digitsum + 9);
check2(count, prime, M, i + 47, digitsum + 11);
check2(count, prime, M, i + 50, digitsum + 5);
check2(count, prime, M, i + 54, digitsum + 9);
check2(count, prime, M, i + 56, digitsum + 11);
check2(count, prime, M, i + 63, digitsum + 9);
check2(count, prime, M, i + 65, digitsum + 11);
check2(count, prime, M, i + 69, digitsum + 15);
check2(count, prime, M, i + 72, digitsum + 9);
check2(count, prime, M, i + 74, digitsum + 11);
check2(count, prime, M, i + 78, digitsum + 15);
check2(count, prime, M, i + 81, digitsum + 9);
check2(count, prime, M, i + 83, digitsum + 11);
check2(count, prime, M, i + 87, digitsum + 15);
check2(count, prime, M, i + 89, digitsum + 17);
check2(count, prime, M, i + 90, digitsum + 9);
check2(count, prime, M, i + 92, digitsum + 11);
check2(count, prime, M, i + 96, digitsum + 15);
check2(count, prime, M, i + 98, digitsum + 17);
}
inline void check100rest3(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
check2(count, prime, M, i + 0, digitsum + 0);
check2(count, prime, M, i + 2, digitsum + 2);
check2(count, prime, M, i + 4, digitsum + 4);
check2(count, prime, M, i + 8, digitsum + 8);
check2(count, prime, M, i + 11, digitsum + 2);
check2(count, prime, M, i + 13, digitsum + 4);
check2(count, prime, M, i + 17, digitsum + 8);
check2(count, prime, M, i + 19, digitsum + 10);
check2(count, prime, M, i + 20, digitsum + 2);
check2(count, prime, M, i + 22, digitsum + 4);
check2(count, prime, M, i + 26, digitsum + 8);
check2(count, prime, M, i + 28, digitsum + 10);
check2(count, prime, M, i + 31, digitsum + 4);
check2(count, prime, M, i + 35, digitsum + 8);
check2(count, prime, M, i + 37, digitsum + 10);
check2(count, prime, M, i + 40, digitsum + 4);
check2(count, prime, M, i + 44, digitsum + 8);
check2(count, prime, M, i + 46, digitsum + 10);
check2(count, prime, M, i + 53, digitsum + 8);
check2(count, prime, M, i + 55, digitsum + 10);
check2(count, prime, M, i + 59, digitsum + 14);
check2(count, prime, M, i + 62, digitsum + 8);
check2(count, prime, M, i + 64, digitsum + 10);
check2(count, prime, M, i + 68, digitsum + 14);
check2(count, prime, M, i + 71, digitsum + 8);
check2(count, prime, M, i + 73, digitsum + 10);
check2(count, prime, M, i + 77, digitsum + 14);
check2(count, prime, M, i + 79, digitsum + 16);
check2(count, prime, M, i + 80, digitsum + 8);
check2(count, prime, M, i + 82, digitsum + 10);
check2(count, prime, M, i + 86, digitsum + 14);
check2(count, prime, M, i + 88, digitsum + 16);
check2(count, prime, M, i + 91, digitsum + 10);
check2(count, prime, M, i + 95, digitsum + 14);
check2(count, prime, M, i + 97, digitsum + 16);
}
inline void check100rest4(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
check2(count, prime, M, i + 1, digitsum + 1);
check2(count, prime, M, i + 3, digitsum + 3);
check2(count, prime, M, i + 7, digitsum + 7);
check2(count, prime, M, i + 9, digitsum + 9);
check2(count, prime, M, i + 10, digitsum + 1);
check2(count, prime, M, i + 12, digitsum + 3);
check2(count, prime, M, i + 16, digitsum + 7);
check2(count, prime, M, i + 18, digitsum + 9);
check2(count, prime, M, i + 21, digitsum + 3);
check2(count, prime, M, i + 25, digitsum + 7);
check2(count, prime, M, i + 27, digitsum + 9);
check2(count, prime, M, i + 30, digitsum + 3);
check2(count, prime, M, i + 34, digitsum + 7);
check2(count, prime, M, i + 36, digitsum + 9);
check2(count, prime, M, i + 43, digitsum + 7);
check2(count, prime, M, i + 45, digitsum + 9);
check2(count, prime, M, i + 49, digitsum + 13);
check2(count, prime, M, i + 52, digitsum + 7);
check2(count, prime, M, i + 54, digitsum + 9);
check2(count, prime, M, i + 58, digitsum + 13);
check2(count, prime, M, i + 61, digitsum + 7);
check2(count, prime, M, i + 63, digitsum + 9);
check2(count, prime, M, i + 67, digitsum + 13);
check2(count, prime, M, i + 69, digitsum + 15);
check2(count, prime, M, i + 70, digitsum + 7);
check2(count, prime, M, i + 72, digitsum + 9);
check2(count, prime, M, i + 76, digitsum + 13);
check2(count, prime, M, i + 78, digitsum + 15);
check2(count, prime, M, i + 81, digitsum + 9);
check2(count, prime, M, i + 85, digitsum + 13);
check2(count, prime, M, i + 87, digitsum + 15);
check2(count, prime, M, i + 90, digitsum + 9);
check2(count, prime, M, i + 94, digitsum + 13);
check2(count, prime, M, i + 96, digitsum + 15);
}
inline void check100rest5(uint32_t digitsum, uint32_t& count, bool* prime, uint32_t i, uint64_t* M)
{
check2(count, prime, M, i + 0, digitsum + 0);
check2(count, prime, M, i + 2, digitsum + 2);
check2(count, prime, M, i + 6, digitsum + 6);
check2(count, prime, M, i + 8, digitsum + 8);
check2(count, prime, M, i + 11, digitsum + 2);
check2(count, prime, M, i + 15, digitsum + 6);
check2(count, prime, M, i + 17, digitsum + 8);
check2(count, prime, M, i + 20, digitsum + 2);
check2(count, prime, M, i + 24, digitsum + 6);
check2(count, prime, M, i + 26, digitsum + 8);
check2(count, prime, M, i + 33, digitsum + 6);
check2(count, prime, M, i + 35, digitsum + 8);
check2(count, prime, M, i + 39, digitsum + 12);
check2(count, prime, M, i + 42, digitsum + 6);
check2(count, prime, M, i + 44, digitsum + 8);
check2(count, prime, M, i + 48, digitsum + 12);
check2(count, prime, M, i + 51, digitsum + 6);
check2(count, prime, M, i + 53, digitsum + 8);
check2(count, prime, M, i + 57, digitsum + 12);
check2(count, prime, M, i + 59, digitsum + 14);
check2(count, prime, M, i + 60, digitsum + 6);
check2(count, prime, M, i + 62, digitsum + 8);
check2(count, prime, M, i + 66, digitsum + 12);
check2(count, prime, M, i + 68, digitsum + 14);
check2(count, prime, M, i + 71, digitsum + 8);
check2(count, prime, M, i + 75, digitsum + 12);
check2(count, prime, M, i + 77, digitsum + 14);
check2(count, prime, M, i + 80, digitsum + 8);
check2(count, prime, M, i + 84, digitsum + 12);
check2(count, prime, M, i + 86, digitsum + 14);
check2(count, prime, M, i + 93, digitsum + 12);
check2(count, prime, M, i + 95, digitsum + 14);
check2(count, prime, M, i + 99, digitsum + 18);
}
inline void addToDigit(digit_t* digits, uint32_t& digitsum, uint32_t digitno, digit_t amount)
{
digits[digitno] += amount;
digit_t diff = 0;
diff += amount;
if (digits[digitno] >= 10) {
digits[digitno] -= 10;
diff -= 9;
uint32_t i = digitno + 1;
while (true) {
digits[i]++;
if (digits[i] == 10) {
diff -= 9;
digits[i] = 0;
i++;
}
else {
break;
}
}
}
digitsum += (uint32_t)diff;
}
inline void loadDigits(digit_t* digits, uint32_t number)
{
int rest = number;
for (uint32_t i = 0; i < 10; i++) {
digits[i] = rest % 10;
rest /= 10;
}
}
inline void sumDigits(digit_t* digits, uint32_t& sum)
{
digit_t sumt = digits[0];
sumt += digits[1];
sumt += digits[2];
sumt += digits[3];
sumt += digits[4];
sumt += digits[5];
sumt += digits[6];
sumt += digits[7];
sumt += digits[8];
sumt += digits[9];
sum = (uint32_t)sumt;
}
bool* prime;
uint64_t * M;
uint32_t * digitsums;
uint64_t digitsums_M;
uint32_t countHarshadPrimes10(uint32_t start, uint32_t end)
{
uint32_t count = 0;
uint32_t digitsum = 0;
digit_t digits[10];
//head:
uint32_t i = start;
loadDigits(digits, start);
sumDigits(digits, digitsum);
uint32_t headEnd = start + (10 - start % 10);
for (; i < headEnd; i++) {
check(digitsum++, count, prime, i, M);
}
//body:
loadDigits(digits, i);
sumDigits(digits, digitsum);
uint32_t bodyEnd = (end / 10) * 10;
for (; i < bodyEnd; i += 10)
{
check(digitsum, count, prime, i, M);
if (digitsum % 2 != 0) {
check(digitsum + 2, count, prime, i + 2, M);
check(digitsum + 4, count, prime, i + 4, M);
check(digitsum + 6, count, prime, i + 6, M);
check(digitsum + 8, count, prime, i + 8, M);
}
else {
check(digitsum + 1, count, prime, i + 1, M);
check(digitsum + 3, count, prime, i + 3, M);
check(digitsum + 5, count, prime, i + 5, M);
check(digitsum + 7, count, prime, i + 7, M);
check(digitsum + 9, count, prime, i + 9, M);
}
addToDigit(digits, digitsum, 1, 1);
}
//tail:
while (i < end) {
check(digitsum++, count, prime, i++, M);
}
return count;
}
inline uint32_t getDigitSum(uint32_t num, const uint32_t * precalc, const uint64_t precalc_M, const uint32_t count)
{
uint32_t upper = fastmod::fastdiv_u32(num, precalc_M);
uint32_t lower = num - (upper * count);
return precalc[lower] + precalc[upper];
}
uint32_t countHarshadPrimes100(uint32_t start, uint32_t end)
{
uint32_t count = 0;
uint32_t digitsum = 0;
digit_t digits[10];
//head:
uint32_t i = std::min(start + (100 - start % 100), end);
count = countHarshadPrimes10(start, i);
//body:
uint32_t bodyEnd = (end / 100) * 100;
loadDigits(digits, i - 100);
sumDigits(digits, digitsum);
for (; i < bodyEnd; i += 100)
{
addToDigit(digits, digitsum, 2, 1);
uint32_t mod = fastmod::fastmod_u32(digitsum, M[6], 6);
switch (mod)
{
case 0:
check100rest0(digitsum, count, prime, i, M);
break;
case 1:
check100rest1(digitsum, count, prime, i, M);
break;
case 2:
check100rest2(digitsum, count, prime, i, M);
break;
case 3:
check100rest3(digitsum, count, prime, i, M);
break;
case 4:
check100rest4(digitsum, count, prime, i, M);
break;
case 5:
check100rest5(digitsum, count, prime, i, M);
break;
}
}
//tail:
if (i < end) count += countHarshadPrimes10(i, end);
return count;
}
uint32_t countHarshadPrimes100Precalc(uint32_t start, uint32_t end)
{
uint32_t count = 0;
uint32_t digitsum = 0;
const uint32_t* precalc = digitsums;
const uint64_t precalc_M = digitsums_M;
const uint32_t precalc_count = 100000;
//head:
uint32_t i = std::min(start + (100 - start % 100), end);
count = countHarshadPrimes10(start, i);
//body:
uint32_t bodyEnd = (end / 100) * 100;
for (; i < bodyEnd; i += 100)
{
digitsum = getDigitSum(i, precalc, precalc_M, precalc_count);
uint32_t mod = fastmod::fastmod_u32(digitsum, M[6], 6);
switch (mod)
{
case 0:
check100rest0(digitsum, count, prime, i, M);
break;
case 1:
check100rest1(digitsum, count, prime, i, M);
break;
case 2:
check100rest2(digitsum, count, prime, i, M);
break;
case 3:
check100rest3(digitsum, count, prime, i, M);
break;
case 4:
check100rest4(digitsum, count, prime, i, M);
break;
case 5:
check100rest5(digitsum, count, prime, i, M);
break;
}
}
//tail:
if (i < end) count += countHarshadPrimes10(i, end);
return count;
}
int countHarshadPrimesMultithreaded(uint32_t start, uint32_t end, uint32_t func(uint32_t,uint32_t), uint32_t numthreads)
{
if (numthreads == 0) {
unsigned concurentThreadsSupported = std::thread::hardware_concurrency();
numthreads = std::max(concurentThreadsSupported, 1U);
}
if (numthreads == 1) return func(start, end);
uint32_t partitionsize = (end + numthreads - start) / numthreads;
std::thread* threads = new std::thread[numthreads];
uint32_t* sums = new uint32_t[numthreads];
for (uint32_t threadNo = 0; threadNo < numthreads; threadNo++)
{
uint32_t threadBegin = start;
uint32_t threadEnd = std::min(start + partitionsize, end);
uint32_t currentThread = threadNo;
threads[threadNo] = std::thread(
[threadBegin, threadEnd, currentThread, &sums, func] {
sums[currentThread] = func(threadBegin, threadEnd);
}
);
start += partitionsize;
}
uint32_t sum = 0;
for (uint32_t threadNo = 0; threadNo < numthreads; threadNo++)
{
threads[threadNo].join();
sum += sums[threadNo];
}
delete[] threads;
delete[] sums;
return sum;
}
void precalcDigitSums(uint32_t count = 100000)
{
digitsums = new uint32_t[count];
uint32_t digitsum = 0;
digit_t digits[10];
loadDigits(digits, 0);
for (uint32_t i = 0; i < count; i++)
{
digitsums[i] = digitsum;
addToDigit(digits, digitsum, 0, 1);
}
digitsums_M = fastmod::computeM_u32(count);
std::cout << "Precalculated " << count << " digitsums" << std::endl;
}
int main(int argc, char** argv)
{
cxxopts::Options options("knowit23", "Count number of Harshad numbers whose digit sums are primes");
options.add_options()
("h,help", "Print this help message")
("s,small", "Run small set (2 to 98765432)")
("p,precalc", "Precalc array of digitsums and use in calculation")
("t,threads", "Thread count (0 to auto detect)", cxxopts::value<std::uint32_t>()->default_value("0"))
("c,count", "Run count", cxxopts::value<std::uint32_t>()->default_value("10"))
;
try {
auto opt = options.parse(argc, argv);
if (opt.count("help") > 0) {
std::cout << options.help();
return 0;
}
int count = opt["count"].as<uint32_t>();
bool small = opt["small"].as<bool>();
int numthreads = opt["threads"].as<uint32_t>();
uint32_t end = small ? 98765432 : INT_MAX - 256;
prime = buildPrimeSieve(90, false);
M = new uint64_t[91];
uint32_t (*func)(uint32_t, uint32_t) = countHarshadPrimes100;
if (opt["precalc"].as<bool>()) {
precalcDigitSums();
func = countHarshadPrimes100Precalc;
}
for (uint32_t d = 1; d < 91; d++) M[d] = fastmod::computeM_u32(d);
std::cout << "From 2 to " << end << std::endl;
int result = 0;
for (int i = 0; i < count; i++)
result = countHarshadPrimesMultithreaded(2, end, func, numthreads);
std::cout << result << std::endl;
}
catch (const cxxopts::OptionException & e)
{
std::cout << "error parsing options: " << e.what() << std::endl;
exit(1);
}
}Et lite screenshot fra codegen-verktøyet mitt:
-s -t 1 -c 100
3,466.97 msec task-clock # 1.000 CPUs utilized
3 context-switches # 0.001 K/sec
0 cpu-migrations # 0.000 K/sec
140 page-faults # 0.040 K/sec
13,495,975,038 cycles # 3.893 GHz
27,387,408,960 instructions # 2.03 insn per cycle
6,710,139,681 branches # 1935.450 M/sec
135,834,587 branch-misses # 2.02% of all branches
3.467150855 seconds time elapsed
3.467190000 seconds user
0.000000000 seconds sys
Samme jobb med @michaelo sin siste
(countHarshadPrimesOptimized(2, 98765432) i en loop fra 0-99)
4,965.99 msec task-clock # 1.000 CPUs utilized
63 context-switches # 0.013 K/sec
0 cpu-migrations # 0.000 K/sec
136 page-faults # 0.027 K/sec
19,330,579,736 cycles # 3.893 GHz
49,240,330,840 instructions # 2.55 insn per cycle
7,902,949,202 branches # 1591.416 M/sec
13,909,003 branch-misses # 0.18% of all branches
4.966442363 seconds time elapsed
4.964241000 seconds user
0.002000000 seconds sys
Litt tregere, men mange færre branch-misses.
Interessante tall. Du kan gjerne også sammenlikne med min med --precalc I argumentlisten.
terjew
Interessante tall. Du kan gjerne også sammenlikne med min med --precalc I argumentlisten.
Det glemte jeg å si ja. 😄
--precalc ga noe dårligere resultat. Ikke mye men målbart (mindre enn 2%)
Jeg klarer fortsatt ikke helt gi slipp på denne, og trengte litt motivasjon, så jeg laget følgende:
EDIT: Ser ut som oppløsningen ble strupet så jeg lastet det opp her i tillegg
Jeg klarer fortsatt ikke helt gi slipp på denne
Hehe. Kjenner følelsen!
så jeg laget følgende:
Den var fin! 👍
Jeg nevnte forøvrig i en tidligere post at jeg hadde noen spor på å utlede rekka av summer direkte, men uansett hvor lur og/eller enkel jeg har forsøkt å være så ender det opp i for mange løkker til å bli ordentlig effektivt (det ble fort noe ala kjøretidskompleksiteten til posten til @DarioSucic lenger opp her) selv om jeg jobbet direkte på summen og rullet ut det dypeste nivået. Fikk den vel ned mot 300ms st, så ikke noe å hoppe i taket av.
Jeg har vært inne på tanken selv. Rekka i seg selv er jo egentlig ganske rett frem, men 10 nøstede for-løkker er jo litt voldsomt. Kunne kanskje rullet ut 2 og 2 siffer, og nøste 5 slike rekker?
Min strategi er jo egentlig å droppe løkker og loop-variabler mest mulig slik at man kan slippe branchingen det fører til. Unntaket er addToDigit som jeg ikke fikk mer effektiv ved å unrolle, men denne kalles uansett bare 1 gang for hvert 100 tall.
@michaelo
Hmm, fikk akkurat en ny idé... Hvis jeg kan splitte arbeidet i litt penere aligned chunker før trådene setter i gang blir jo egentlig alt mye lettere. Da kan jeg kanskje bake inn noen antakelser om modulo 6 også... Fristende! Men det må bli en annen dag.
@michaelo @bobcat4 @DarioSucic
Jeg har tenkt litt... i rettferdighetens navn må man kunne argumentere for at når man sammenlikner benchmarks kjørt på et skjermkort til 9000 kroner så må det være rimelig å sammenlikne dem med tilsvarende benchmarks kjørt på en prosessor til samme pris.
Så, jeg fant en kollega med en 3950X og tigget meg til litt kjøretid på hans private boks. Her er resultatet:
Lite, ST: 34,5 ms
Lite, MT (precalc): 3,78ms
Lite, MT: 3,74ms
Stort, ST: 712ms
Stort. MT (precalc): 57ms
Stort. MT: 57ms
@terjew cc: @DarioSucic @bobcat4
Hehe. Bra bra. Gyldig argument det!
Så er da de neste spørsmålene:
?
Fikk ned CUDA fra 5.4ms til 4.1ms ved å bytte fra Atomic Add til Parallell blokk-reduksjon slik det er beskrevet her.
Tror det er mye mer å hente, men det er i alle fall en start. 😅
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <math.h>
#define u32 uint32_t
#define u8 uint8_t
#define f32 float
#define f64 double
#define BLOCK_SIZE 128
__device__ void warp_reduce(volatile u32 *block_data, u32 tid)
{
block_data[tid] += block_data[tid + 32];
block_data[tid] += block_data[tid + 16];
block_data[tid] += block_data[tid + 8];
block_data[tid] += block_data[tid + 4];
block_data[tid] += block_data[tid + 2];
block_data[tid] += block_data[tid + 1];
}
__device__ u32 compute_single(u32 k, u8 *prime)
{
const u32 pow_table[9] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000};
u32 digit_sum = 0;
for (u32 i = 0; i < 9; i++)
digit_sum += (k / pow_table[i]) % 10;
return (k % digit_sum == 0) & prime[digit_sum];
}
__global__ void kernel(u32 *block_counts, u8 *prime, u32 n)
{
extern __shared__ u32 block_data[];
u32 tid = threadIdx.x;
u32 idx = blockIdx.x * blockDim.x + tid;
if (idx > n) block_data[tid] = 0;
else block_data[tid] = compute_single(idx + 1, prime);
__syncthreads();
for(u32 s = blockDim.x / 2; s > 32; s >>= 1)
{
if(tid < s) block_data[tid] += block_data[tid + s];
__syncthreads();
}
if(tid < 32) warp_reduce(block_data, tid);
if(tid == 0) block_counts[blockIdx.x] = block_data[0];
}
u8 *sieve(u32 n)
{
u32 lim = n + 1;
u8 *prime = (u8 *)malloc(lim);
memset(prime, 1, lim);
prime[0] = 0;
prime[1] = 0;
for (u32 k = 2; k < lim; k++)
{
if (!prime[k])
continue;
for (u32 i = k * k; i < lim; i += k)
prime[i] = 0;
}
return prime;
}
int main()
{
const u32 n = 98765432;
const u32 prime_limit = 9 * ceil(log10(n));
cudaEvent_t st, et;
cudaEventCreate(&st);
cudaEventCreate(&et);
u8 *prime = sieve(prime_limit);
u8 *prime_device;
cudaMalloc(&prime_device, prime_limit);
cudaMemcpy(prime_device, prime, prime_limit, cudaMemcpyHostToDevice);
u32 num_blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
u32 shared_mem_size = BLOCK_SIZE * sizeof(u32);
u32 *counts = (u32 *) malloc(num_blocks * sizeof(u32));
u32 *counts_device;
cudaMalloc(&counts_device, num_blocks * sizeof(u32));
cudaEventRecord(st);
kernel<<<num_blocks, BLOCK_SIZE, shared_mem_size>>>(counts_device, prime_device, n);
cudaEventRecord(et);
cudaEventSynchronize(et);
float ms = 0;
cudaEventElapsedTime(&ms, st, et);
printf("Kernel execution: %.3f ms\n", ms);
cudaMemcpy((void *)counts, counts_device, num_blocks * sizeof(u32), cudaMemcpyDeviceToHost);
u32 total_count = 0;
for(u32 i = 0; i < num_blocks; i++)
total_count += counts[i];
printf("Count: %d\n", total_count);
cudaFree(counts_device);
cudaFree(prime_device);
free(prime);
return 0;
}
Fikk ned CUDA fra 5.4ms til 4.1ms ved å bytte fra Atomic Add til Parallell blokk-reduksjon slik det er beskrevet her.
Tror det er mye mer å hente, men det er i alle fall en start. 😅
Nais! 👏
Tar du en benchmark på stort sett også?
- Varer også kodekalenderjula helt til påske?
Minst! Og da må den vel avløses av noen kodepåskenøtter?
Tar du en benchmark på stort sett også?
Endte opp med 85ms på det store datasettet, men egentlig jukser jeg litt siden mesteparten av summeringen nå er flyttet over til CPU utenfor kernelen. Om jeg får tid kan jeg se om jeg ikke får slengt inn litt bedre timing slik at resultatet er sammenlignbart med resten av tidene deres.
Tar du en benchmark på stort sett også?
Endte opp med
85mspå det store datasettet, men egentlig jukser jeg litt siden mesteparten av summeringen nå er flyttet over til CPU utenfor kernelen. Om jeg får tid kan jeg se om jeg ikke får slengt inn litt bedre timing slik at resultatet er sammenlignbart med resten av tidene deres.
Tenkte på det samme når jeg så over koden. Du gjemmer bort et par ting i benchmarken, deriblant den siste summeringen men også memcpy til og fra CUDA-devicen. Tenker kanskje det hadde vært riktigere å la programmet loope f.eks 1000 ganger på lite datasett, og måle total tid brukt på dette inkludert de punktene, men uten å ta med opprettelse og teardown av CUDA-konteksten. Det utgjør vel ikke all verden, men når vi er nede på noen få millisekunder for hele kjøringen begynner alle sånne småting å utgjøre en relevant del av tiden.
Evt time hele programmet med 1000 interne kjøringer fra shellet slik vi gjør for de andre variantene. Da får du jo med de 200 ms på opprettelsen av CUDA-kontekst etc, men dette fordeles jo utover de 1000 kjøringene så det skal vel ikke utgjøre mer enn 0,2ms per kjøring.
@bobcat4
Det var bl.a. clang's (til sitt forsvar: velmente) brutalitet som slo til her også ja. Men hadde enkelte målinger som var ganske så erratiske også som jeg ikke har klart å gjenskape. Vel, det du ser nå er i tråd (pun initielt ikke intended) med hva jeg ser her også nå.
Takk for optimaliseringsmotgiften 👍 . Den var effektig. Jeg ble selv lurt her av at den ikke optimaliserte
countHarshadPrimesMultithreadedut av loopen (selv kjørt med numthreads=1 som argument), kun når jeg kjørtecountHarshadPrimesOptimizeddirekte derifra.Det er ganske så hyggelige tall, helt klart. 1.086s/0.270s for stort datasett 1x singletrådet/multitrådet er ikke grusomt. Og som du sier: vi holder det fortsatt lineært. Denne tilnærmingen vil riktignok få en bump i det øyeblikket vi overstiger 10 siffer og da må legge til enda en modulo i
getCrossumChunked_fastmod_unrolled.