Snabbaste sättet att beräkna ett 128-bitars heltal modulo ett 64-bitars heltal (2023)

Snabbaste sättet att beräkna ett 128-bitars heltal modulo ett 64-bitars heltal

Jag har ett 128-bitars osignerat heltal A och ett 64-bitars osignerat heltal B. Vad är det snabbaste sättet att beräknaA % B– det är (64-bitars) återstoden från att dividera A med B?

Jag funderar på att göra detta i antingen C- eller assemblerspråk, men jag måste rikta in mig på 32-bitars x86-plattformen. Detta betyder tyvärr att jag inte kan dra fördel av kompilatorstöd för 128-bitars heltal, inte heller av x64-arkitekturens förmåga att utföra den nödvändiga operationen i en enda instruktion.

Redigera:

Tack för svaren så här långt. Men det verkar för mig som att de föreslagna algoritmerna skulle vara ganska långsamma – skulle inte det snabbaste sättet att utföra en 128-bitars med 64-bitars division vara att utnyttja processorns inbyggda stöd för 64-bitars med 32-bitars division? Är det någon som vet om det finns ett sätt att utföra den större divisionen vad gäller några mindre divisioner?

Re: Hur ofta byter B?

I första hand är jag intresserad av en generell lösning – vilken beräkning skulle du göra om A och B sannolikt kommer att vara olika varje gång?

En andra möjlig situation är dock att B inte varierar lika ofta som A – det kan finnas så många som 200 som att dividera med varje B. Hur skulle ditt svar skilja sig åt i det här fallet?

Du kan använda divisionsversionen avRysk bondemultiplikation.

För att hitta resten, kör (i pseudokod):

X = B;medan(X <= A/2){ X <<=1;}medan(A >= B){om(A >= X) A -= X; X >>=1;}

Modulen är kvar i A.

Du måste implementera skiftningar, jämförelser och subtraktioner för att arbeta på värden som består av ett par 64-bitars tal, men det är ganska trivialt (sannolikt bör du implementera vänster-förskjutning-för-1 somX + X).

Detta kommer att loopa högst 255 gånger (med 128 bitars A). Naturligtvis måste du göra en förkontroll för en nolldelare.

Kanske letar du efter ett färdigt program, men de grundläggande algoritmerna för aritmetik med flera precision kan hittas i KnuthsKonsten att programmera, Volym 2. Du kan hitta divisionsalgoritmen som beskrivs onlinehär. Algoritmerna hanterar godtycklig aritmetik med flera precision, och är därför mer generella än du behöver, men du bör kunna förenkla dem för 128-bitars aritmetik gjord på 64- eller 32-bitars siffror. Var beredd på en rimlig mängd arbete (a) att förstå algoritmen och (b) att konvertera den till C eller assembler.

Du kanske också vill kolla inHackers Delight, som är full av mycket smart assembler och annat hackeri på låg nivå, inklusive en del aritmetik med flera precision.

Om ditt B är tillräckligt litet föruint64_t +operation för att inte linda:

GivenA = AH*2^64 + AL:

A % B == (((AH % B) * (2^64% B)) + (AL % B)) % B == (((AH % B) * ((2^64- B) % B)) + (AL % B)) % B

Om din kompilator stöder 64-bitars heltal är detta förmodligen den enklaste vägen att gå.
MSVC:s implementering av en 64-bitars modulo på 32-bitars x86 är en hårig slingfylld sammansättning (VCcrtsrcintelllrem.asmför de modiga), så jag skulle personligen gå med på det.

Detta är nästan oprövad, delvis hastighetsmodifierad Mod128by64 "Russian peasant"-algoritmfunktion. Tyvärr är jag en Delphi-användare så den här funktionen fungerar under Delphi. 🙂 Men montören är nästan densamma så...

fungeraMod128by64(Utdelning: PUInt128; Divisor: PUInt64): UInt64;//In : eax = @Dividend// : edx = @Divisor//Ut: eax:edx som återstodasm//Registrerar i rutinen//Divisor = edx:ebp//Dividend = bh:ebx:edx //Vi behöver 64 bitar + 1 bit i bh//Resultat = esi:edi//ecx = Slingräknare och utdelningsindextryck ebx//Lagra register att staplapush esi push edi push ebp mov ebp, [edx]//Divisor = edx:ebpmov edx, [edx +4] mov ecx, ebp//Div med 0 testeller ecx, edx jz @DivByZero xor edi, edi//Klart resultatxor esi, esi//Start av 64-bitars divisionsslingamov ecx,15 //Ladda byte loop shift räknare och utdelningsindex@SkipShift8Bits://Små utdelningssiffror förskjutningsoptimeringcmp [eax + ecx], kap//Noll testjnz @EndSkipShiftDividend loop @SkipShift8Bits//Hoppa över 8-bitars loop@EndSkipShiftDividend: test edx, $FF000000//Enorma divisionsnummerskifteoptimeringjz @Shift8Bits//Denna Divisor är > $00FFFFFF:FFFFFFFFmov ecx,8 //Ladda byte skifträknaremov front, [eax +12]//Gör snabb 56-bitars (7 byte) shift...shr esi, cl//esi = $00XXXXXXmov edi, [eax +9]//Ladda för en byte högerskiftat 32 bitars värde@Shift8Bits: mov bl, [eax + ecx]//Ladda 8 bitar av Dividend//Här kan vi ta bort partiell loop 8-bitars division för att öka exekveringshastigheten...mov ch,8 //Sätt partiell byte räknare värde@Do65BitsShift: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1setc bh//Spara 65:e bitensub edi, ebp//Jämför utdelning och divisorsbb esi, edx// Subtrahera divisornsbb bh,0 //Använd 65:e biten i bhjnc @NoCarryAtCmp//Testa...lägg till edi, ebp//Return privius dividendstatadc esi, edx@NoCarryAtCmp: dec ch//Minska räknarejnz @Do65BitsShift//Slutet av 8-bitars (byte) partiell divisionsslingadec cl//Minska byte loop skifträknarejns @Shift8Bits//Sista hoppet vid cl = 0!!!//Slutet av 64-bitars divisionsslingaflytta eax, redigera//Ladda in resultatet till eax:edxmov edx, esi@RestoreRegisters: pop ebp//Återställ registerpop edi pop esi pop ebx ret@DivByZero: xor eax, eax//Här kan du höja Div med 0 undantag, nu fungerar bara returnera 0.xor edx, edx jmp @RestoreRegistersend;

Åtminstone ytterligare en hastighetsoptimering är möjlig! Efter 'Huge Divisor Numbers Shift Optimization' kan vi testa divisorer med hög bit, om det är 0 behöver vi inte använda extra bh-register som 65:e bit för att lagra i det. Så utrullad del av slingan kan se ut så här:

shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1sub edi, ebp//Jämför utdelning och divisorsbb esi, edx// Subtrahera divisornjnc @NoCarryAtCmpX lägg till edi, ebp//Return privius dividendstatadc esi, edx@NoCarryAtCmpX:

Jag känner till frågan som specificerade 32-bitars kod, men svaret för 64-bitars kan vara användbart eller intressant för andra.

Och ja, division 64b/32b => 32b är en användbar byggsten för 128b % 64b => 64b. libgcc__umoddi3(källa länkad nedan) ger en uppfattning om hur man gör sånt, men den implementerar bara 2N % 2N => 2N ovanpå en 2N / N => N-division, inte 4N % 2N => 2N.

Bredare multiprecisionsbibliotek finns tillgängliga, t.ex.https://gmplib.org/manual/Integer-Division.html#Integer-Division.

GNU C på 64-bitars maskinerger en__du128typ, och libgcc-funktioner för att multiplicera och dividera så effektivt som möjligt på målarkitekturen.

x86-64div r/m64instruktionen gör 128b/64b => 64b division (producerar också återstoden som en andra utgång), men den misslyckas om kvoten svämmar över. Så du kan inte direkt använda den omA/B > 2^64-1, men du kan få gcc att använda den åt dig (eller till och med infoga samma kod som libgcc använder).

Detta sammanställer (Godbolt kompilator utforskare) till en eller tvådivinstruktioner (som händer i enlibgccfunktionsanrop). Om det fanns ett snabbare sätt skulle libgcc förmodligen använda det istället.

#omfatta uint64_t AmodB(osignerad__int128 A,uint64_tB){lämna tillbakaA % B;}

De__umodti3funktion som den anropar beräknar en hel 128b/128b modulo, men implementeringen av den funktionen kontrollerar det speciella fallet där divisorns övre halva är 0, som du kanse i libgcc-källan. (libgcc bygger si/di/ti-versionen av funktionen från den koden, som är lämpligt för målarkitekturen.udiv_qrnndär ett inline-asm-makro som gör osignerad 2N/N => N-division för målarkitekturen.

För x86-64(och andra arkitekturer med en hårdvarudelningsinstruktion),den snabba vägen(närhigh_half(A) < B; garanterardivkommer inte att fela)är bara två ej tagna grenar, lite ludd för processorer som inte fungerar som de ska att tugga igenom,och en singeldiv r64instruktion, som tar cirka 50-100 cykler1på moderna x86-processorer, enligtAgner Fogs insn tabeller. En del annat arbete kan pågå parallellt meddiv, men heltalsdelningsenheten är inte särskilt pipelined ochdivavkodar till många uops (till skillnad från FP division).

Reservvägen använder fortfarande bara två 64-bitarsdivinstruktioner för det fall därBär bara 64-bitars, menA/Bpassar inte i 64 bitar såA/Bdirekt skulle fel.

Observera att libgcc's__umodti3just inlines__udivmoddi4i ett omslag som bara returnerar resten.

Fotnot 1: 32-bitarsdivär över 2 gånger snabbare på Intel-processorer. På AMD-processorer beror prestanda bara på storleken på de faktiska ingångsvärdena, även om de är små värden i ett 64-bitars register. Om små värden är vanliga kan det vara värt att jämföra en gren till en enkel 32-bitars division innan du gör 64-bitars eller 128-bitars division.

För upprepad modulo av sammaB

Det kan vara värt att överväga att beräkna afixpunkt multiplikativ inversförB, om en sådan finns. Till exempel, med kompileringstidskonstanter, gör gcc optimeringen för typer som är smalare än 128b.

uint64_t modulo_by_constant64(uint64_tA){lämna tillbakaA %0x12345678ABULL; } movabs rdx,-2233785418547900415mov rax, rdi mul rdx mov rax, rdx# bortkastad instruktion, kunde ha fortsatt att använda RDX.movabs rdx,78187493547shr rax,36 # divisionsresultatimul rax, rdx# multiplicera och subtrahera för att få modulosub rdi, rax mov rax, rdi ret

x86:ormul r64instruktionen gör 64b*64b => 128b (rdx:rax) multiplikation, och kan användas som ett byggblock för att konstruera en 128b * 128b => 256b multiplikation för att implementera samma algoritm. Eftersom vi bara behöver den höga hälften av hela 256b-resultatet, sparar det några multiplikationer.

Moderna Intel-processorer har mycket hög prestandamul: 3c latens, en per klocka genomströmning. Den exakta kombinationen av skift och tillägg som krävs varierar dock med konstanten, så det allmänna fallet med att beräkna en multiplikativ invers vid körning är inte lika effektiv varje gång den används som en JIT-kompilerad eller statiskt kompilerad version (även ovanpå förberäkningskostnaderna).

IDK där brytpunkten skulle vara. För JIT-kompilering kommer det att vara fler än ~200 återanvändningar, såvida du inte cachelagrar genererad kod för vanligt användaBvärden. För det "normala" sättet kan det möjligen vara i intervallet 200 återanvändningar, men IDK hur dyrt det skulle vara att hitta en modulär multiplikativ invers för 128-bitars / 64-bitars division.

halkarkan göra detta åt dig, men bara för 32- och 64-bitarstyper. Ändå är det förmodligen en bra utgångspunkt.

Jag har gjort båda versionerna av Mod128by64 "Russian peasant" divisionsfunktion: klassisk och hastighetsoptimerad. Hastighetsoptimerad kan göra på min 3Ghz PC mer än 1000 000 slumpmässiga beräkningar per sekund och är mer än tre gånger snabbare än klassisk funktion.
Om vi ​​jämför exekveringstiden för att beräkna 128 gånger 64 och beräkna 64 gånger 64 bitars modulo så är denna funktion bara cirka 50% långsammare.

Klassisk rysk bonde:

fungeraMod128by64Classic(Utdelning: PUInt128; Divisor: PUInt64): UInt64;//In : eax = @Dividend// : edx = @Divisor//Ut: eax:edx som återstodasm//Registrerar i rutinen//edx:ebp = Divisor//ecx = Slingräknare//Resultat = esi:editryck ebx//Lagra register att staplapush esi push edi push ebp mov ebp, [edx]//Ladda divisor till edx:ebpmov edx, [edx +4] mov ecx, ebp//Div med 0 testeller ecx, edx jz @DivByZero push [eax]//Lagra Divisor till stackentryck [eax +4] tryck [eax +8] tryck [eax +12] refräng var, var//Klart resultatChor esi, esi mov ecx,128 //Ladda skifträknare@Do128BitsShift: shl [esp +12],1 //Skift utdelning från stack vänster för en bitrcl [esp +8],1rcl [esp +4],1rcl [esp],1rcl edi,1rcl esi,1setc bh//Spara 65:e bitensub edi, ebp//Jämför utdelning och divisorsbb esi, edx// Subtrahera divisornsbb bh,0 //Använd 65:e biten i bhjnc @NoCarryAtCmp//Testa...lägg till edi, ebp//Return privius dividendstatadc esi, edx@NoCarryAtCmp: loop @Do128BitsShift//Slutet av 128-bitars divisionsslingaflytta eax, redigera//Ladda in resultatet till eax:edxmov edx, esi@RestoreRegisters: lea esp, esp +16 //Återställ Divisors utrymme på stackpop ebp//Återställ registerpop edi pop esi pop ebx ret@DivByZero: xor eax, eax//Här kan du höja Div med 0 undantag, nu fungerar bara returnera 0.xor edx, edx jmp @RestoreRegistersend;

Hastighetsoptimerad rysk bonde:

fungeraMod128by64Oprimerad(Utdelning: PUInt128; Divisor: PUInt64): UInt64;//In : eax = @Dividend// : edx = @Divisor//Ut: eax:edx som återstodasm//Registrerar i rutinen//Divisor = edx:ebp//Dividend = ebx:edx //Vi behöver 64 bitar//Resultat = esi:edi//ecx = Slingräknare och utdelningsindextryck ebx//Lagra register att staplapush esi push edi push ebp mov ebp, [edx]//Divisor = edx:ebpmov edx, [edx +4] mov ecx, ebp//Div med 0 testeller ecx, edx jz @DivByZero xor edi, edi//Klart resultatxor esi, esi//Start av 64-bitars divisionsslingamov ecx,15 //Ladda byte loop shift räknare och utdelningsindex@SkipShift8Bits://Små utdelningssiffror förskjutningsoptimeringcmp [eax + ecx], kap//Noll testjnz @EndSkipShiftDividend loop @SkipShift8Bits//Skip Compute 8 Bits unrolled loop ?@EndSkipShiftDividend: test edx, $FF000000//Enorma divisionsnummerskifteoptimeringjz @Shift8Bits//Denna Divisor är > $00FFFFFF:FFFFFFFFmov ecx,8 //Ladda byte skifträknaremov front, [eax +12]//Gör snabb 56-bitars (7 byte) shift...shr esi, cl//esi = $00XXXXXXmov edi, [eax +9]//Ladda för en byte högerskiftat 32 bitars värde@Shift8Bits: mov bl, [eax + ecx]//Ladda 8 bitars del av Dividend// Beräkna 8 bitars orullad loopshl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividentAbove0//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow0 ja @DividentAbove0 cmp edi, ebp//utdelning är del större?jb @DividentBelow0@DividentAbove0: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividendBelow0: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividentAbove1//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow1 ja @DividentAbove1 cmp edi, ebp//utdelning är del större?jb @DividentBelow1@DividentAbove1: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividendBelow1: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividentAbove2//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow2 ja @DividentAbove2 cmp edi, ebp//utdelning är del större?jb @DividentBelow2@DividentAbove2: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividentBelow2: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividentAbove3//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow3 ja @DividentAbove3 cmp edi, ebp//utdelning är del större?jb @DividentBelow3@DividentAbove3: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividentBelow3: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividentAbove4//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow4 ja @DividentAbove4 cmp edi, ebp//utdelning är del större?jb @DividentBelow4@DividentAbove4: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividentBelow4: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividentAbove5//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow5 ja @DividentAbove5 cmp edi, ebp//utdelning är del större?jb @DividentBelow5@DividentAbove5: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividentBelow5: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividendAbove6//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow6 ja @DividentAbove6 cmp edi, ebp//utdelning är del större?jb @DividentBelow6@DividentAbove6: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividentBelow6: shl bl,1 //Skift utdelning kvar en bitrcl edi,1rcl esi,1jc @DividentAbove7//utdelning hi bit set?cmp esi, edx//utdelning hej del större?jb @DividentBelow7 ja @DividentAbove7 cmp edi, ebp//utdelning är del större?jb @DividentBelow7@DividentAbove7: sub edi, ebp//Return privius dividendstatsbb esi, edx@DividentBelow7://Beräkningsslut 8 bitar (unrolled loop)dec cl//Minska byte loop skifträknarejns @Shift8Bits//Sista hoppet vid cl = 0!!!//Slut på divisionsslinganflytta eax, redigera//Ladda in resultatet till eax:edxmov edx, esi@RestoreRegisters: pop ebp//Återställ registerpop edi pop esi pop ebx ret@DivByZero: xor eax, eax//Här kan du höja Div med 0 undantag, nu fungerar bara returnera 0.xor edx, edx jmp @RestoreRegistersend;

Jag skulle vilja dela med mig av några tankar.

Det är inte så enkelt som MSN föreslår är jag rädd.

I uttrycket:

(((AH % B) * ((2^64- B) % B)) + (AL % B)) % B

både multiplikation och addition kan svämma över. Jag tror att man kan ta hänsyn till det och fortfarande använda det allmänna konceptet med vissa modifieringar, men något säger mig att det kommer att bli riktigt läskigt.

Jag var nyfiken på hur 64-bitars modulo-operation implementerades i MSVC och jag försökte ta reda på något. Jag vet inte riktigt montering och allt jag hade tillgängligt var Express-utgåvan, utan källan till VCcrtsrcintelllrem.asm, men jag tror att jag lyckades få en uppfattning om vad som pågår, efter lite lekande med felsöknings- och demonteringsutgången. Jag försökte ta reda på hur resten beräknas vid positiva heltal och divisorn >=2^32. Det finns viss kod som handlar om negativa siffror så klart, men jag grävde inte i det.

Så här ser jag det:

Om divisor >= 2^32 skiftas både utdelningen och divisorn åt höger så mycket som behövs för att passa in divisorn i 32 bitar. Med andra ord: om det krävs n siffror för att skriva ner divisorn binärt och n > 32, slängs n-32 minst signifikanta siffror av både divisor och utdelning. Därefter utförs divisionen med hjälp av hårdvarustöd för att dividera 64-bitars heltal med 32-bitars ettor. Resultatet kan vara felaktigt, men jag tror att det går att bevisa att resultatet kan vara lägre med högst 1. Efter divisionen multipliceras divisorn (den ursprungliga) med resultatet och produkten subtraheras från utdelningen. Sedan korrigeras det genom att addera eller subtrahera divisorn om det behövs (om resultatet av divisionen var borta med ett).

Det är lätt att dela 128-bitars heltal med 32-bitars ett med hjälp av hårdvarustöd för 64-bitars med 32-bitars division. Om divisorn < 2^32, kan man beräkna resten genom att utföra bara 4 divisioner enligt följande:

Låt oss anta att utdelningen lagras i:

DWORD-utdelning[4] = ...

resten kommer att gå in på:

DWORD rest;1) Dela utdelning[3] efter divisor. Förvara resten i resten.2) Dela uppQWORD (resten:utdelning[2])efter divisor. Förvara resten i resten.3) DelaQWORD (resten:utdelning[1])efter divisor. Förvara resten i resten.4) DelaQWORD (resten:utdelning[0])efter divisor. Förvara resten i resten.

Efter dessa 4 steg kommer den variabla resten att innehålla det du letar efter.
(Snälla döda mig inte om jag missförstod mig. Jag är inte ens en programmerare)

Om divisorn är större än 2^32-1 har jag inga goda nyheter. Jag har inte ett fullständigt bevis på att resultatet efter skiftet inte är mer än 1, i proceduren jag beskrev tidigare, som jag tror att MSVC använder. Jag tror dock att det har något att göra med det faktum att den del som kasseras är minst 2^31 gånger mindre än divisorn, utdelningen är mindre än 2^64 och divisorn är större än 2^32-1 , så resultatet är mindre än 2^32.

Om utdelningen har 128 bitar fungerar inte tricket med att kassera bitar. Så generellt sett är förmodligen den bästa lösningen den som föreslagits av GJ eller caf. (Tja, det skulle förmodligen vara det bästa även om det fungerade att kassera bitar. Division, multiplikationssubtraktion och korrigering på 128 bitars heltal kan vara långsammare.)

Jag funderade också på att använda flyttals-hårdvaran. x87 flyttalsenhet använder 80 bitars precisionsformat med bråkdel 64 bitar långt. Jag tror att man kan få det exakta resultatet av 64 bitar för 64 bitars division. (Inte resten direkt, utan även resten med multiplikation och subtraktion som i "MSVC-proceduren"). OM utdelningen >=2^64 och < 2^128 att lagra den i flyttalsformatet verkar likna att kassera minst signifikanta bitar i "MSVC-proceduren". Kanske kan någon bevisa att felet i så fall är bundet och finner det användbart. Jag har ingen aning om det har en chans att vara snabbare än GJ:s lösning, men det kanske är värt det att prova.

Lösningen beror på exakt vad du försöker lösa.

T.ex. om du gör aritmetik i en ring modulo ett 64-bitars heltal sedan använda
Montgomerys minskningär mycket effektiv. Naturligtvis förutsätter detta att du har samma modul många gånger och att det lönar sig att omvandla elementen i ringen till en speciell representation.

För att bara ge en väldigt grov uppskattning av hastigheten på denna Montgomerys-reduktion: Jag har ett gammalt riktmärke som utför en modulär exponentiering med 64-bitars modul och exponent i 1600 ns på en 2,4Ghz Core 2. Denna exponentiering gör ungefär 96 modulära multiplikationer ( och modulära reduktioner) och behöver därför cirka 40 cykler per modulär multiplikation.

Det accepterade svaret från @caf var riktigt trevligt och högt betygsatt, men det innehåller en bugg som inte har setts på flera år.

För att hjälpa till att testa det och andra lösningar lägger jag upp en testsele och gör den till en community-wiki.

osignerad cafMod(osigneradA,osigneradB){ hävda(B);osigneradX = B;// while (X < A / 2) { Originalkod använd < medan(X <= A /2) { X <<=1; }medan(A >= B) {om(A >= X) A -= X; X >>=1; }lämna tillbakaA;}tomhet cafMod_test(osigneradnum,osigneradde){om(den ==0)lämna tillbaka;osignerady0 = num % då;osignerady1 = mod(antal, den);om(y0 != y1) {printf("FAIL num:%x den:%x %x %xn"num, the, y0, yl); flush(stdout);utgång(-1); }}osignerad rand_unsigned(){osigneradx = (osignerad) rand();lämna tillbakax *2^ (osignerad) rand();}tomhet cafMod_tests(tomhet){konst osigneradi[] = {0,1,2,3,0x7FFFFFFF,0x80000000, UINT_MAX -3, UINT_MAX -2, UINT_MAX -1, UINT_MAX };för(osigneradden =0; den <storlek avjag /storlek avjag[0]; den++) {om(i[den] ==0)Fortsätta;för(osigneradnum =0; num <storlek avjag /storlek avjag[0]; num++) { cafMod_test(i[num], i[den]); } } cafMod_test(0x8711dd11,0x4388ee88); cafMod_test(0xf64835a1,0xf64835a);tid_tt; tid(&t); srand((osignerad) t);printf("%och", (osignerad) t);fflush(stdout);för(lång långn =10 000 LL*1000 LL*1000 LL; n >0; n--) { cafMod_test(rand_unsigned(), rand_unsigned()); }sätter("Gjort");}int huvud(tomhet){ cafMod_tests();lämna tillbaka 0;}

Som en allmän regel är division långsam och multiplikation är snabbare, och bitskiftning är snabbare än. Från vad jag har sett av svaren hittills har de flesta av svaren använt en brute force-strategi med hjälp av bitskift. Det finns ett annat sätt. Om det är snabbare återstår att se (AKA profile it).

Istället för att dividera, multiplicera med det reciproka. För att upptäcka A % B, beräkna därför först den reciproka av B … 1/B. Detta kan göras med några slingor med hjälp av Newton-Raphson-metoden för konvergens. För att göra detta väl kommer att bero på en bra uppsättning initiala värden i en tabell.

För mer information om Newton-Raphson-metoden för att konvergera på det reciproka, sehttp://en.wikipedia.org/wiki/Division_(digital)

När du väl har ömsesidigt är kvoten Q = A * 1/B.

Resten R = A – Q*B.

För att avgöra om detta skulle vara snabbare än råkraften (eftersom det kommer att bli många fler multiplikationer eftersom vi kommer att använda 32-bitars register för att simulera 64-bitars och 128-bitars tal, profilera det.

Om B är konstant i din kod kan du förberäkna den reciproka och helt enkelt beräkna med de två sista formlerna. Jag är säker på att detta kommer att gå snabbare än att byta bitar.

Hoppas det här hjälper.

osignerad 128_mod(uint64_tövre,uint64_tlägre,uint64_tdiv) {uint64_tresultat =0;uint64_ta = (~0%div)+1; övre %= div;// den resulterande bitlängden bestämmer antalet cykler som krävs // först räknar vi ut modulär multiplikation av (2^64*övre)%div medan(övre !=0){om(övre&1==1){ resultat += a;om(resultat >= div){resultat -= div;} } a <<=1;om(a >= div){a -= div;} övre >>=1; }// addera de 2 resultaten och returnera modulen om(lower>div){lower -= div;}lämna tillbaka(lägre+resultat)%div;}

Jag vet inte hur man kompilerar assemblerkoderna, all hjälp uppskattas för att kompilera och testa dem.

Jag löste detta problem genom att jämföra med gmplib "mpz_mod()" och summera 1 miljon loopresultat. Det var en lång resa att gå från avmattning (seedup 0,12) till speedup 1,54 — det är anledningen till att jag tror att C-koderna i den här tråden kommer att vara långsamma.

Detaljer inklusive testsele i denna tråd:
https://www.raspberrypi.org/forums/viewtopic.php?f=33&t=311893&p=1873122#p1873122

Detta är "mod_256()" med snabbare användning över att använda gmplib "mpz_mod()", användning av __builtin_clzll() för längre skift var avgörande:

typdef __uint128_t uint256_t[2];#definieramin(x, y) ((xint clz(__uint128_tu){// unsigned long long h = ((osigned long long *)&u)[1]; osignerad lång långh = u >>64;lämna tillbaka(h!=0) ? __builtin_clzll(h) :64+ __builtin_clzll(u);}__uint128_t mod_256(uint256_tx,__uint128_tn){om(x[1] ==0)lämna tillbakax[0] % n;annan{__uint128_tr = x[1] % n;intF = clz(n);intR = clz(r);för(inti=0; i<128; ++i) {om(R>F+1) {inth = min(R-(F+1),128-i); r <<= h; R-=h; i+=(h-1);Fortsätta; } r <<=1;om(r >= n) {r -= n; R=clz(r); } } r += (x[0] % n);om(r >= n) r -= n;lämna tillbakar; }}

Om du har en ny x86-maskin finns det 128-bitarsregister för SSE2+. Jag har aldrig försökt skriva montering för något annat än grundläggande x86, men jag misstänker att det finns några guider där ute.

Jag är 9 år efter striden men här är ett intressant O(1)-kantfall för krafter på 2 som är värt att nämna.

#omfatta // exempel med 32 bitar och 8 bitar.int huvud(){intjag =930;osignerad rödingb = (osignerad röding) i;printf("%d", (int) b);// 162, samma som 930 % 256}

Eftersom det inte finns någon fördefinierad 128-bitars heltalstyp i C, måste bitar av A representeras i en array. Även om B (64-bitars heltal) kan lagras i enosignerad lång lång intvariabel, behövs det att lägga in bitar av B i en annan array för att kunna arbeta på A och B effektivt.

Därefter inkrementeras B som Bx2, Bx3, Bx4, … tills det är det största B mindre än A. Och sedan kan (A-B) beräknas med hjälp av viss subtraktionskunskap för bas 2.

Är det den här typen av lösning du letar efter?

References

Top Articles
Latest Posts
Article information

Author: Aron Pacocha

Last Updated: 11/06/2023

Views: 5321

Rating: 4.8 / 5 (48 voted)

Reviews: 95% of readers found this page helpful

Author information

Name: Aron Pacocha

Birthday: 1999-08-12

Address: 3808 Moen Corner, Gorczanyport, FL 67364-2074

Phone: +393457723392

Job: Retail Consultant

Hobby: Jewelry making, Cooking, Gaming, Reading, Juggling, Cabaret, Origami

Introduction: My name is Aron Pacocha, I am a happy, tasty, innocent, proud, talented, courageous, magnificent person who loves writing and wants to share my knowledge and understanding with you.