[C] Compilatore e divisione

[apatriarca]

Ho comunque provato la formula più precisa invece della approssimazione e viene che supporta tutti i numeri fino a \(2^{32}\) con soli \(34\) bit invece di averne bisogno di \(36\) come nella formula approssimata.

[utente__medio]

Ho scritto il seguente piccolo script in Julia:

Codice:

δ(d, n) = ceil(1/d, digits=n, base=2)
ϵ(d, n) = δ(d, n) - 1/d
log₂ϵ(d, n) = log2(ϵ(d, n))
log₂M(d, n) = - log₂ϵ(d, n) - log2(d)
M(d, n) = exp2(log₂M(d, n))
M(7, 11)


che mi restituisce \(682.6666666666407\) per cui sembra che la formula \(\log_2(M) + \log_2(\epsilon) < \log_2(\delta)\) sia abbastanza accurata (la differenza di due unità potrebbe essere dovuto a qualche errore di calcolo immagino).

Ma se la condizione originale da rispettare è \( k\,d\epsilon < \delta \), penso si possa fare anche a meno dei logaritmi, infatti con semplici passaggi si ottiene che \( k< \delta/(d \delta -1) \), e quindi a rigore l'unica inferenza che possiamo fare (posto \( A = [ \delta/(d \delta -1)] \)) è che $ A*d <= D_{MAX} < (A+1)*d $ . O sbaglio?

Nel nostro caso facendo i conti viene $ 679 <= D_{MAX} < 686 $.

[apatriarca]

La formula \(k\,d\,\epsilon < \delta\) considera solo i multipli di \(d\) ed è la ragione per cui hai ottenuto solo l'intervallo di valori all'interno del quale trovi il massimo. Bella trovata comunque, non ci avevo pensato ed è effettivamente più semplice. Nota che il passaggio a \(d\,\epsilon = d\delta - 1\) non è necessario perché sia \(\epsilon\) che \(d\) sono delle costanti per cui è solo un modo diverso di ottenere lo stesso valore.

Si può a questo punto ricavare anche \(M = A\,d + a\) osservando che devi avere \(A\,d\,\epsilon + a\,\delta < 1,\) da cui puoi semplicemente ricavare
\[ A\,(d\,\delta - 1) + a\,\delta = (A\,d + a)\,\delta - A = M\delta - A < 1 \implies M < (1 + A)/\delta \]

Ho modificato il codice con questa nuova logica:

Codice:

function max_dividend(d, n)
    δ = ceil(1/d, digits=n, base=2)
    A = floor(δ/(d*δ - 1))
    return floor(Int64, (1 + A)/δ)
end


Restituisce 684.

[utente__medio]

Ottimo, passare da un intervallo di ampiezza $d$ al valore esatto è un passo in avanti non da poco!

Però non capisco da dove viene la condizione \( A\,d\,\epsilon + a\,\delta < 1 \).

[apatriarca]

Viene da
\[ A = \lfloor(A\,d + a)\,\delta\rfloor = \lfloor A\,d\,\delta + a\,\delta\rfloor = \lfloor A\,d\,(1/d + \epsilon) + a\,\delta\rfloor = A + \lfloor A\,d\,\epsilon + a\,\delta\rfloor\]
A questo punto rimuovendo A da entrambi i lati ti rimane che quello che hai nel floor è positivo e deve essere minore di uno perché il suo floor sia zero.

[utente__medio]

Chiaro, grazie!

Stavo cercando di ricapitolare il tutto, ma mi è sorto un dubbio. Il nostro scopo è dimostrare che \(D \div d = \lfloor D \delta \rfloor \) per determinati valori di $D$. Nell'uguaglianza appena scritta ho utilizzato il floor, ma l'operazione di troncamento \( [\cdot] \) sarebbe stata in questo caso del tutto equivalente, dal momento che stiamo considerando valori positivi.
Il discorso cambia invece quando consideriamo la condizione

\( \lfloor kd \epsilon - \delta\rfloor = -1\Rightarrow \begin{cases} kd \epsilon - \delta < 0 \\ kd \epsilon - \delta \geq -1 \end{cases} \)

o la condizione

\( [ kd \epsilon - \delta] = -1\Rightarrow \begin{cases} kd \epsilon - \delta \leq -1 \\ kd \epsilon - \delta > -2 \end{cases} \)

Su quali basi dovremmo decidere a priori se utilizzare il floor o il troncamento?

[apatriarca]

Con floor possiamo estrarre valori interi dal floor senza che il risultato cambi, ma questo non è possibile nel caso del troncamento se, come in questo caso, passiamo da valori positivi a negativi.

[utente__medio]

Ho provato a ricapitolare il tutto riformulando alcune cose:

Testo nascosto, fai click qui per vederlo

La divisione intera \( D \backslash d = \lfloor D / d \rfloor = \lfloor D \cdot \frac{1}{d} \rfloor \) (con \( D > 0 \) variabile intera e \( d > 1 \) costante intera) è una funzione costante a tratti con gradini unitari in corrispondenza dei multipli del divisore, ossia per quei valori di \( D \) esprimibili nella forma \( k \ d \) (con \( k \in N^+ \)).
Detto \( a \) un intero appartenente all'intervallo \( [ 1 \ ; \ d ) \), sarà quidi:

\( f(D) = \lfloor D \cdot \frac{1}{d} \rfloor = \begin{cases} k , \ \ \ \ \ \ \ \ \ \ \ D = k \ d \\ k - 1 , \ \ \ \ D = k \ d - a \end{cases} \)

Sia \( q = \lceil \frac{1}{d} \rceil _{n, \ b} \) un arrotondamento per eccesso di \( \frac{1}{d} \) all'\( n \)-esima cifra frazionaria nella base numerica \( b \) (ovviamente sarà \( \frac{1}{d} \leq q < 1 \)), affinché \( \lfloor D \ q \rfloor \) restituisca gli stessi risultati della divisione intera dovranno quindi essere rispettate le seguenti condizioni:

\( \begin{cases} \lfloor ( k \ d ) \ q \rfloor = k \\ \lfloor ( k \ d - a ) \ q \rfloor = k - 1 \end{cases} \Rightarrow \begin{cases} k \ d \ q \geq k \\ k \ d \ q < k + 1 \\ ( k \ d - a ) \ q \geq k - 1 \\ ( k \ d - a ) \ q < k \end{cases} \)

Introdotta la differenza \( \epsilon = q - \frac{1}{d} \geq 0 \), si evince che la prima e la terza disequazione risultano sempre verificate:

1) \( \ \ k \ d \ q \geq k \ \Rightarrow \ k \ d \ ( \epsilon + \frac{1}{d} ) - k \geq 0 \ \Rightarrow \ k \ d \ \epsilon + k - k \geq 0 \ \Rightarrow \ k \ d \ \epsilon \geq 0 \)

3) \( \ \ ( k \ d - a ) \ q \geq k - 1 \ \Rightarrow \ (k \ d - a ) ( \epsilon + \frac{1}{d} ) - k + 1 \geq 0 \ \Rightarrow \ k \ d \ \epsilon + k - a \ \epsilon - \frac{a}{d} - k + 1 \geq 0 \ \Rightarrow \ \epsilon \ ( k \ d - a ) + ( 1 - \frac{a}{d} ) \geq 0 \)

Restano quindi la seconda e la quarta disequazione:

\( \begin{cases} k \ d \ q < k + 1 \ \ \ \ \Rightarrow \ k \ d \ q - k < 1 \\ ( k \ d - a ) \ q < k \ \Rightarrow \ k \ d \ q - k < a \ q \end{cases} \)

La condizione più restrittiva è data dal minore tra i secondi membri delle due disequazioni, ossia \( q \), che si ottiene da \( a \ q \) per \( a = a_{min} = 1 \), e che risulta minore dell'unità.
Quindi in definitiva:

\( \begin{cases} \lfloor ( k \ d ) \ q \rfloor = k \\ \lfloor ( k \ d - a ) \ q \rfloor = k - 1 \end{cases} \ \Rightarrow \ k \ d \ q - k < q \ \Rightarrow \ k < \frac{q}{d \ q - 1} \ \Rightarrow \ k_{max} = \lfloor \frac{q}{d \ q - 1} \rfloor \)

A questo punto è possibile restringere il campo passando dall'intervallo \( [ k_{max} \ d \ ; \ ( k_{max} + 1 ) \ d ) \) ad un singolo valore esprimibile come \( D_{max} = k_{max} \ d + a \). Per farlo notiamo che dovrà essere:

\( \lfloor D_{max} \ q \rfloor = k_{max} \ \Rightarrow \begin{cases} D_{max} \ q \geq k_{max} \ \Rightarrow \ ( k_{max} \ d + a ) ( \epsilon + \frac{1}{d} ) - k_{max} \geq 0 \ \Rightarrow \ k_{max} \ d \ \epsilon + a \ \epsilon + \frac{a}{d} \geq 0 \\ D_{max} \ q < k_{max} + 1 \ \Rightarrow \ D_{max} < \frac{k_{max} + 1}{q} \end{cases} \)

La prima disequazione è sempre verificata e quindi dalla seconda risulta che \( D_{max} = \lfloor \frac{k_{max} + 1}{q} \rfloor \).

In questo modo viene anche meno la questione di distinguere tra floor e troncamento, in quanto non avendo a che fare con numeri negativi le due funzioni risultano del tutto equivalenti.

Passiamo ora dalla parte più "matematica" a quella più "informatica":

Testo nascosto, fai click qui per vederlo

Volendo operare con gli interi, ed essendo in questo caso \( q = \lceil \frac{1}{d} \rceil _{n, \ 2} \), si ricava che

\( \lfloor D \ q \rfloor = \lfloor \frac {D \ ( q \ \cdot \ 2^n )}{2^n} \rfloor = \lfloor \frac {D \ \cdot \ Q}{2^n} \rfloor = D \cdot Q \gg n \)

dove è bene ricordare che l'operatore di moltiplicazione ha la precedenza su quello di right-shift.
Per farla semplice utilizziamo per \( D \) e \( Q \) delle variabili intere senza segno a \(32 \) bit (uint32_t), in modo che il loro prodotto possa essere contenuto in un uint64_t. Sarà quindi \( n = n_B + n_Z \), dove \( n_B = 32 \) è il numero di bit utilizzati per \( Q \) e \( n_Z = \lfloor \log_2(d) \rfloor \) è il numero di zeri frazionari iniziali di \( q \).
Ovviamente per l'applicabilità del metodo bisognerà sempre controllare che \( D \leq D_{max} \), con

\( k_{max} = \lfloor \frac{q}{d \ q - 1} \rfloor = \lfloor \frac{q \ \cdot \ 2^n}{( d \ q - 1 ) \ 2^n} \rfloor = \lfloor \frac{Q}{d \ Q - 2^n} \rfloor \),

\( D_{max} = \lfloor \frac{k_{max} + 1}{q} \rfloor = \lfloor \frac{ ( k_{max} + 1 ) \ 2^n}{q \ \cdot \ 2^n} \rfloor = \lfloor \frac{ ( k_{max} + 1 ) \ 2^n}{Q} \rfloor \).

Quello che ho scritto è corretto o ho commesso qualche errore?

[utente__medio]

Una delle due questioni che ponevo nel post iniziale può oramai ritenersi abbondantemente risolta (o almeno credo), passo quindi all'altra.

Relativamente alla divisione nativa, nei casi in cui il divisore è una costante, i coefficienti (che io chiamo $Q$ e $n$) per quel determinato valore di $d$ possono essere tranquillamente calcolati in fase di compilazione; se invece per esempio d = rand() + 1 il discorso cambia e diventa necessario calcolare i coefficienti per RAND_MAX + 1 valori (che io vado a salvare negli array v_Q[RAND_MAX + 1] e v_n[RAND_MAX + 1]). La domanda quindi sorge spontanea: come si comporta il compilatore in questo caso?
Nella mia ignoranza, dando un'occhiata ai tempi quasi sovrapponibili

Codice:

50000000x [ 941787252 / (rand() + 1) ] IN 0.732000s (NATIVA)
50000000x [ 941787252 / (rand() + 1) ] IN 0.766000s (PERSONALE)

mi verrebbe da pensare che anche il compilatore (ammesso che utilizzi lo stesso procedimento di eseguire la divisione mediante una moltiplicazione ed uno shift) prelevi i coefficienti da qualche parte, e che questo prelievo sia leggermente più veloce rispetto a quello che io faccio dagli array globali v_Q e v_n.
Concordate sul fatto che anche il compilatore deve avere già a disposizione da qualche parte questi valori? E in tal caso dove?

[apatriarca]

Molto probabilmente fa semplicemente una divisione intera invece che usare il prodotto e uno shift. In effetti qui puoi vedere che fa quello che dico.