C(GCC) 220 216 214 212 bytes

credit to @ceilingcat for 2 bytes

#define L(v)for(int v=2;~v--;)

#define P l/C+r<0|l/C+r>=R|l%C+c<0|l%C+c>=C

f(int R,int C,int*m){for(int*i=m+R*C;~*i;) {int*M,l=*i+++C**i++,a=0,b;L(r)L(c)P?:++a;M=m+l;b=*M/a;b+=!b;*M- =b*a;L(r)L(c)M[r*C+c]+=P?0:b;}}

Run it here

a slightly less golfed version

#define L(v)for(int v=2;~v--;)

#define P l/C+r<0|l/C+r>=R|l%C+c<0|l%C+c>=C

f(int R, int C, int*m) {

    for(int*i=m+R*C;~*i;) {

        int*M,l=*i+++C**i++,a=0,b;

        L(r)

            L(c)

                P?:++a;

        M=m+l;

        b=*M/a;

        b+=!b;

        *M-=b*a;

        L(r)

            L(c)

                M[r*C+c]+=P?0:b;

    }

}

The calling code with an example

int main()

{

  int matrix[] = {3,3,3,3,3,3,3,3,3,0,0,0,1,0,2,-1};

  int rows = 3;

  int columns = 3;

  f(rows,columns,matrix);

  for(int r = 0; r < rows; ++r) {

    for(int c = 0; c < columns; ++c) {

      printf("%03d,",matrix[r*columns + c]);

    }

    printf("n");

  }

}

and the output

001,005,003,

000,006,003,

001,005,003,

JavaScript (ES7),  126 125 123  121 bytes

Saved 2 bytes thanks to @Shaggy

Takes input as (matrix)(list). Outputs by modifying the matrix.

m=>a=>a.map(([Y,X])=>(g=n=>m[m.map((r,y)=>r.map((_,x)=>(x-X)**2+(y-Y)**2<3&&r[n++,x]++)),(m[Y][X]+=~n)<n||g``,Y][X]++)``)

Try it online!

How?

If we were to strictly follow the algorithm described in the challenge, we'd have to perform the following operations for each position pair $(x,y)$:

1. walk through the matrix to compute the number of neighbors $n$
   


2. compute $lfloor m(x,y) / nrfloor$ and deduce the quantity $q$ that needs to be added to each neighbor


3. walk through the matrix again to update each neighbor


4. update $m(x,y)$
   

Instead of that, we use a recursive function which executes a simpler flow of operations, repeated as many times as needed:

1. for each neighbor and for the reference cell itself: increment the cell and increment a counter $n$ (initialized to $0$)


2. subtract $n+1$ from the reference cell


3. if the above result is greater than or equal to $n$, do a recursive call


4. increment the reference cell (all steps of this kind are executed in succession when the last recursive call has completed)

The main benefit is that we only need one loop over the matrix. The second benefit is that we don't have to compute any quotient at all.

Example

$$M=begin{pmatrix}
0&0&0\
0&26&0\
0&0&0
end{pmatrix}text{ and }(x,y)=(1,1)$$

After step 1 of the first iteration, we have:

$$M=begin{pmatrix}
1&1&1\
1&27&1\
1&1&1
end{pmatrix}text{ and }n=9$$

And after step 2 of the first iteration:

$$M=begin{pmatrix}
1&1&1\
1&17&1\
1&1&1
end{pmatrix}$$

The key point here is that the reference cell is updated by taking only the cost ($-9$) into account without the gain ($+1$), so that we know how much we can still draw from it.

Consequently, the sum of the cells is no longer equal to $26$ at this point. But this will be corrected at the end of the process, where all gains on the reference cell are accounted at once.

Step 3 of the first iteration: because $17$ is greater than $9$, we do a recursive call.

After step 1 of the second iteration, we have:

$$M=begin{pmatrix}
2&2&2\
2&18&2\
2&2&2
end{pmatrix}text{ and }n=9$$

And after step 2 of the second iteration:

$$M=begin{pmatrix}
2&2&2\
2&8&2\
2&2&2
end{pmatrix}$$

Step 3 of the second iteration: this time, we have $8<9$, so we stop there.

We now increment the reference cell twice (step 4 of both iterations), leading to the final result:

$$M=begin{pmatrix}
2&2&2\
2&10&2\
2&2&2
end{pmatrix}$$

Commented

m => a =>                     // m[] = input matrix, a[] = list of positions

  a.map(([Y, X]) => (         // for each pair (X, Y) in a[]:

    g = n =>                  //   g = recursive function expecting n = 0

      m[                      //

        m.map((r, y) =>       //     for each row r[] at position y in m[]:

          r.map((_, x) =>     //       for each value at position x in r[]:

            (x - X) ** 2 +    //         if the quadrance between (x, y)

            (y - Y) ** 2 < 3  //         and (X, Y) is less than 3:

            && r[n++, x]++    //           increment n and increment r[x]

          )                   //       end

        ),                    //     end

        (m[Y][X] += ~n)       //     subtract n + 1 from m[Y][X]

        < n                   //     if the result is greater than or equal to n:

        || g``,               //       do a recursive call

        Y                     //     

      ][X]++                  //     increment m[Y][X]

    )``                       //   initial call to g

  )                           // end

Clean, 181 167 bytes

import StdEnv;

foldlm(x,y)={{if(d>2)0b+e-if(d>0)0b*n\e<-:l&v<-[0..],let{b=max m.[y,x]n/n;$a b=2+sign a-(a+1)/size b;n= $x l* $y m;d=(v-x)^2+(u-y)^2}}\l<-:m&u<-[0..]}

Try it online!

In the form of a partially-applied function literal.

Expanded (first version):

f // functinon f on {{Int}} and [(Int,Int)]

    = foldl m (x, y) // fold :: (a -> b -> a) a [b] -> a with first argument  {{Int}} (Int,Int) -> {{Int}} giving  {{Int}} [(Int,Int)] -> {{Int}}

        = {                     // an array of

            {                   // arrays of

                if(d > 2) 0 b   // the amount we give to the neighbors

                + e             // plus the current entry

                - if(d > 0) 0 b // minus the amount taken from the target entry

                * n             // times the number of neighbors, if we're on the target

            \                  // for each

                e <-: l         // element of row l

                & v <- [0..]    // and x-index v

                , let           // local definitions:

                    b           // the amount given to the neighbors

                        = max   // we need at least 1 each, so take the largest of

                            m.[y, x] // the target entry

                            n   // or the number of neighbors

                        / n     // divide it by the number of neighbors

                    n           // the number of neighbors

                        = (     // sum of

                            1   // one

                            + s x // if x is at the left edge = 0 else 1

                            + s ( // if x is at the right edge = 0 else 1

                                size l

                                - x 

                                - 1

                            )

                        ) * (   // times the sum of

                            1   // one

                            + s y // if y is at the top edge = 0 else 1

                            + s ( // if y is at the bottom edge = 0 else 1

                                size m

                                - y

                                - 1

                            )

                        )

                    d           // distance from the target point

                        = (v - x)^2

                        + (u - y)^2

            }

        \                      // for each

            l <-: m             // row l in matrix m

            & u <- [0..]        // and y-index u

        }

Rust - 295 bytes

fn explode(p:(i8,i8),v:&mut Vec<Vec<i8>>){let x=v[p.0 as usize][p.1 as usize];let q=|x,y|x*x+y*y;loop{let mut t=0;for i in 0..v.len(){for j in 0..v[i].len(){if q(i as i8-p.0,j as i8-p.1)<3{v[i][j]+=1;v[p.0 as usize][p.1 as usize]-=1;t+=1;}}}if v[p.0 as usize][p.1 as usize]<=(x/t+x%t){break;}}}

This is pretty long due to Rust requiring unsigned integer indexing of vectors, but requiring signed integers to do subtraction resulting in negatives. However I believe my algorithm is the "shortest algorithm" so far. There is actually no need to deal with detecting edges, bottom, etc.

Notice three things: One, the sum of all cells is always constant. Two, this is a division / remainder situation, so we can apply Bresenham-algorithm style thinking. Three, the question always adds the same number to all cells within a certain distance of the special position cell, before dealing with the "extra" stuff in the special position.

Algorithm:

Store original value of cell at position P into M.

Begin Loop:

Iterate over each cell I in the matrix. If the position of cell I is within 3 Quadrance (squared distance) of the position P, then subtract 1 from cell P and add 1 to the cell I. Count how many times this is done in one iteration through the matrix.

If the value leftover in the cell at position P is less than or equal to M / Count + M modulo Count, then break the loop. Otherwise perform the loop again.

The resulting matrix will be the exploded version. Count is basically a way to count neighbors without dealing with edges. Looping is a way to break down the division/addition stuff into a repeated single addition/subtraction of one. The modulo check ensures we will have appropriate remainder left at position P to deal with 'explosions' that are not evenly divisible amongst neighbors. The do/while loop structure allows P<0 to work properly.

Ungolfed version on the Rust Playground

R, 163 162 161 159 bytes

function(m,l){for(e in l){i=e[1];j=e[2];x=-1:1+i;y=-1:1+j;s=m[x<-x[x<=dim(m)],y<-y[y<=ncol(m)]];z=sum(1|s);d=max(1,m[i,j]%/%z);m[x,y]=s+d;m[i,j]=m[i,j]-d*z};m}

Try it online!

Explanation

function(m,l) {          # Take input as matrix m and 1-indexed list of explosion points l

  for(e in l) {          # Loop over the list of explosion points

    i=e[1]; j=e[2]       # Assign current coordinates to (i,j) for brevity

    x=-1:1+i             # Assign the ranges of neighboring cells: (i-1) to (i+1),

    y=-1:1+j             # and (j-1) to (j+1)

    s=                   # Take the submatrix s=m[x,y]

      m[x<-x[x<=dim(m)]  # But first trim x and y from above to prevent out of bounds errors,

     ,y<-y[y<=ncol(m)]]  # trimming from below isn't necessary, as R tolerates index 0

    z=sum(1|s)           # Count the neighbors

    d=max(1,m[i,j]%/%z)  # Estimate, how much we'll distribute to each neighbor

    m[x,y]=s+d           # Add the distributed amount to each cell of the submatrix

    m[i,j]=m[i,j]-d*z    # Subtract the total amount from the exploded cell

  }

  m                      # Return the modified matrix

}

Common Lisp, 498 bytes

(defmacro s(l c x)`(incf(aref m,l,c),x))

(defmacro w(a &rest f)`(if(,a(or(= l 0)(= l(d 0)))(or(= c 0)(= c(d 1)))),@f))

(defmacro d(n)`(1-(array-dimension m,n)))

(defmacro p(i l m &rest f)`(loop for,i from(1-,l)to(1+,l)when(and(>=,i 0)(<=,i,m))do,@f))

(defmacro g()`(or(p i l(d 0)(p j c(d 1)(s i j 1)))(s l c(- v))))

(defun f(m l c)(let((v(w and 4(w or 6 9))))(if (<(/(s l c 0)v)1)(g)(loop for i to(1-(/(s l c 0)v))do(g)))))

(defun c(m n)(dotimes(i(length n))(f m(nth 0(nth i n))(nth 1(nth i n))))m)

Try it online!

Use this function as (print (c #2A((3 3 3) (3 3 3) (3 3 3)) '((0 0)(0 1)(0 2))))

Better readable version:

(defmacro s (l c x)

  `(incf (aref m ,l ,c) ,x))



(defmacro w (a &rest f)

  `(if (,a (or (= l 0)

           (= l (d 0)))

       (or (= c 0)

           (= c (d 1))))

       ,@f))



(defmacro d (n)

  `(1- (array-dimension m ,n)))



(defmacro p (i l m &rest f)

  `(loop for ,i from (1- ,l) to (1+ ,l)

     when (and (>= ,i 0) (<= ,i ,m))

     do ,@f))



(defmacro g ()

  `(or(p i l (d 0)

     (p j c (d 1)

        (s i j 1)))

      (s l c (- v))))



(defun f (m l c)

  (let ((v (w and 4 (w or 6 9))))

    (if (< (/ (s l c 0) v) 1)

    (g)

      (loop for i to (1- (/ (s l c 0) v))

        do (g)))))



(defun c (m n)

  (dotimes (i (length n))

    (f m (nth 0 (nth i n))

       (nth 1 (nth i n))))

  m)

Output example:

(print (c #2A((3 3 3) (3 3 3) (3 3 3) (3 3 3) (3 3 3) (3 3 3)) '((5 0)(4 1)(0 2))))

;; #2A((3 4 0) (3 4 4) (3 3 3) (4 4 4) (5 -4 4) (1 5 4))



(print (c #2A((3 3 3) (3 3 3) (3 3 3)) '((0 0)(0 1)(0 2))))

; #2A((1 0 1) (5 6 5) (3 3 3))  => #2A((1 0 1) (5 6 5) (3 3 3))



(print (c #2A((9 8 7) (8 9 7) (8 7 9)) '((0 0)(1 1)(2 2))))

;; #2A((4 11 8) (11 5 10) (9 10 4))  => #2A((4 11 8) (11 5 10) (9 10 4))



(print (c #2A((0 0) (0 0)) '((0 0)(0 0)(0 0))))

;; #2A((-9 3) (3 3))  => #2A((-9 3) (3 3))



(print (c #2A((10 20 30)(30 20 10)(40 50 60)) '((0 2)(2 0)(1 1)(1 0))))

;; #2A((21 38 13) (9 12 21) (21 71 64))  => #2A((21 38 13) (9 12 21) (21 71 64))

Java 10, 194 193 bytes

M->C->{for(var q:C){int f=1,n=0,X=q[0],Y=q[1],x,y,c=1,t;for(;f>0;f=(M[X][Y]+=~n)<n?0:c++)for(n=0,x=M.length;x-->0;)for(y=M[x].length;y-->0;)if((t=x-X)*t+(t=y-Y)*t<3)M[x][y]+=++n/n;M[X][Y]+=c;}}

Iterative port of @Arnauld's JavaScript answer.

Modifies the input-matrix instead of returning a new one to save bytes.

Try it online.

Explanation:

M->C->{                     // Method with two integer-matrix parameters and no return-type

  for(var q:C){             //  Loop over the coordinates:

    int f=1,                //   Flag integer, starting at 1 for every coordinate

        n=0,                //   Count integer, starting at 0 for every coordinate

        X=q[0],Y=q[1],      //   The current X,Y coordinate

        x,y,                //   Temp x,y coordinates

        c=1,                //   Counter, starting at 1

        t;                  //   Temp integer

    for(;f>0                //   Loop as long as the flag is not 0:

        ;                   //     After every iteration:

         f=(M[X][Y]+=~n)    //      Decrease the value at X,Y by n+1

            <n?             //      If it's now smaller than the count `n`

               0            //       Change the flag to 0

              :             //      Else:

               c++)         //       Change the flag to counter `c`, and increase `c` by 1

      for(n=0,              //    Reset the count `n` to 0

          x=M.length;x-->0;)//    Loop `x` over the rows:

        for(y=M[x].length;y-->0;)

                            //     Inner loop `y` over the columns:

          if((t=x-X)*t      //      If the difference between `x` and `X` squared

             +(t=y-Y)*t     //      and the difference between `y` and `Y` squared together

              <3){          //      Is smaller than 3:

            M[x][y]+=++n/n; //       Increase count `n` and the value at x,y both by 1

    M[X][Y]+=c;}}           //   Increase the value at X,Y with counter `c`