[Question]
  [
[sandbox (deleted)](https://codegolf.meta.stackexchange.com/a/17288/78039) 
Lets define a matrix of 9s as:
$$ N = \begin{bmatrix} 9&9&9\\9&9&9\\9&9&9 \end{bmatrix} $$
Lets define an exploding number as a number at position \$(x,y)\$ that can be decomposed into equal integers between all its adjacent neighbors (including itself) and the absolute value of each portion is greater than 0.
From the previous matrix, lets explode the number at position \$(1,1)\$ (0 indexed)
$$ N = \begin{bmatrix} 9&9&9\\9&\color{red}9&9\\9&9&9 \end{bmatrix} $$
$$ N = \begin{bmatrix} 9+\color{red}1&9+\color{red}1&9+\color{red}1\\9+\color{red}1&\color{blue}0+\color{red}1&9+\color{red}1\\9+\color{red}1&9+\color{red}1&9+\color{red}1 \end{bmatrix} $$
$$ N = \begin{bmatrix} 10&10&10\\10&\color{red}1&10\\10&10&10 \end{bmatrix} $$
Sometimes, decomposing result into a rational number greater than 1. This is something we need to avoid when exploding numbers. In this cases the remainder will be assigned to the exploded number.
To demonstrate it, lets continue working with our previous matrix. This time we will explode the number at position \$(0,0)\$
$$ N = \begin{bmatrix} \color{red}{10}&10&10\\10&1&10\\10&10&10 \end{bmatrix} $$
Here we have 3 neightbors and the number itself. Here the equation is something like \$10/4\$ which give us **2** for each and **2** as remainder.
$$ N = \begin{bmatrix} \color{blue}2+\color{red}{2}&\color{red}{10+2}&10\\\color{red}{10+2}&\color{red}{1+2}&10\\10&10&10 \end{bmatrix} $$
$$ N = \begin{bmatrix} \color{red}{4}&12&10\\12&3&10\\10&10&10 \end{bmatrix} $$
As well, sometimes a number wont be big enough to be decomposed in equal parts (where the abs is greater than 0) between his neighbors (|rational number| < 1). In this cases we need to *"borrow"* from the exploded number in order to maintain the *"greater than 0"* condition. Lets continue with our previous example and explode the number at position \$(1,1)\$.
$$ N = \begin{bmatrix} 4&12&10\\12&\color{red}3&10\\10&10&10 \end{bmatrix} $$
$$ N = \begin{bmatrix} 4+\color{red}1&12+\color{red}1&10+\color{red}1\\12+\color{red}1&\color{blue}0+\color{red}{1}-\color{green}6&10+\color{red}1\\10+\color{red}1&10+\color{red}1&10+\color{red}1 \end{bmatrix} $$
$$ N = \begin{bmatrix} 5&13&11\\13&\color{red}{-5}&11\\11&11&11 \end{bmatrix} $$
---
**The challenge** is, given a list of \$(x,y)\$ positions and an finite non-empty array of natural numbers, return the exploded form after each number from the positions list has been exploded.
---
**Test cases**
Input: `initial matrix: [[3, 3, 3], [3, 3, 3], [3, 3, 3]], numbers: [[0,0],[0,1],[0,2]]`
Output: `[[1, 0, 1], [5, 6, 5], [3, 3, 3]]`
---
Input: `Initial matrix: [[9, 8, 7], [8, 9, 7], [8, 7, 9]], numbers: [[0,0],[1,1],[2,2]]`
Output: `[[4, 11, 8],[11, 5, 10],[9, 10, 4]]`
---
Input: `Initial matrix: [[0, 0], [0, 0]], numbers: [[0,0],[0,0],[0,0]]`
Output: `[[-9, 3],[3, 3]]`
---
Input: `Initial Matrix: [[10, 20, 30],[30, 20, 10],[40, 50, 60]], numbers: [[0,2],[2,0],[1,1],[1,0]]`
Output: `[[21, 38, 13], [9, 12, 21], [21, 71, 64]]`
---
Input: `Initial Matrix: [[1]], numbers: [[0,0]]`
Output: `[[1]]`
---
Input: `Initial Matrix: [[1, 2, 3]], numbers: [[0,0], [0, 1]]`
Output: `[[1, 1, 4]]`
---
**Notes**
* [Input/Output rules](https://codegolf.meta.stackexchange.com/questions/2447/default-for-code-golf-input-output-methods) apply
* You can assume input matrix will never be empty
* You can assume coordinates are always going to be valid
* Input coord in test cases is given as (row, column). If you need it to be (x, y) you can swap the values. If so, please state that in your answer
  
[Answer]
# 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](https://tio.run/##XZDBbsIwDIbvfQqPCSmJU1HojRA49EolxJVxaAOdIrVlyhibVMqrd0loYZqixPKf77edqPBdqa57PRwLXR9hTS60OBmi6zNc5EzcLmEoaDBcb6CcJGgW0dXHpdxey3GCygkuLmUSFN685e5M3Mkq2vQ1mZYVblkibkwL2jgl5aVkGhET5gLPZMRzsSaGromim9UcMROptZUilyydZCJH@ZILloYyZ9mDTHeG2Qn2KDeraJ6Ltu3cBFWma0KDJgC4p2ejf3Z7kNDE/P@K/JraPePhtBW9yZy@P60hHnJ1Kr@q@iEVxAG8V/m9BfU3/U8ai0bChoWvJQDRUHAzPRl1Z5Rl@koOUwMG8GEsVpDROIoPfNS3sY8epkFQe98VoA3@Gt7qkdfboO1@AQ)
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,
```
[Answer]
# 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!](https://tio.run/##lZJvb4IwEMbf@ynulWnl0Lb4B6Llc2iaZhCnZouAQbNo5vbVWUuHOF1cBiV9uB73u6flNX1L98vyZXfw8@J5Va1llck4NaOfpTtC1ALnmsqYbGQu40xlLlziyQRL9/KER5tx9Oe01xMeOfkLK2ZBt1uq3PPwqD2PUiSZWmg11578zOksP583SYIu4tEkodVhtT@ABJIhpBRkDO@wJhklKZ3Cssj3xXbV3xYbE5rCR6dj04lSAZpbI9wKqxRDZmeG3E1Cawq312BgMjkCQ6jTRghjhJErBXZo3dAiDHFiV0KMGjHB6AeNO5p4QBsalCGGdboRBsnrbyMrEIYt8WKB3Vhqw/eQe6IfWR/fltrqBiYYBnWtoNaujyHDEcPxhSmco2t/3KHr6sJ4CELTe9CYEAjCbYNZmphnfOWJt1b@av6RJ35VEcUvZ/6f4s1fwOvdr74A "JavaScript (Node.js) – 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) / n\rfloor\$ 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
```
[Answer]
# [R](https://www.r-project.org/), ~~163~~ ~~162~~ ~~161~~ ~~159~~ ~~155~~ 146 bytes
```
function(m,l){for(e in l){v=m[i<-e[1],j<-e[2]];s=m[x<--1:(i<dim(m))+i,y<--1:(j<ncol(m))+j];z=sum(1|s);d=max(1,v%/%z);m[x,y]=s+d;m[i,j]=v+d-d*z};m}
```
[Try it online!](https://tio.run/##bZDRboMgFIbv9xTeNIF5zETWdQ55i901XDRYE4xootZYtz67O9KuqWwJcA7n/Pwf0M5FIOfiVOveNDWxUNGvomnJMTB1gPkg7d5k0XHPFJRLTJQSHRbHLIrYBzFZbiyxlIYGztdSmdW6qVytVGKS3ckS9t1RkUt7GAmDYfOymahADzgr2YU5pgZKJYcwj/Ln6SLsZS6IPfStGUkKHDiFynQ90SSBhIJGE7aEZUfp013LV9qbiP0e4SutRud32OFM3bqDlC7H4dM3uDG5ZxCDa/iwe1jDGMpj4G5ggtvXGLYxvMV/qfyKe4QnniEClv4D3MPh5dCG/fce5n5t/gE "R – Try It Online")
### Explanation
(Corresponds to a previous version of the code)
```
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
}
```
[Answer]
# [Clean](https://github.com/Ourous/curated-clean-linux), ~~181~~ 167 bytes
```
import StdEnv;
```
# 
```
foldl\m(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!](https://tio.run/##ZY69boMwGEXn5Ck@CRTZwSbA1jrOlA6VumUEKhkwkSXbVOFHUMSr13WrbtWV7j3bubWWwjrTNaOWYISyzCnz0T0GuA3Ni52Yuw3iMex3AbTAoe10owuDZrJgvq6qRc0lw0kVSfrLieejLQp5ps/6MJ1pnsRxSbQc1oobMYOJ84XMpT1ZFgqoeBb16m5BUCSiFJ969SmhYpZDOIM@QriAYQ1HE53xexahkS5@t60otDeYw/hn2PY7fw7WYA2eiM9G/sEGOUpJignKSOb7h0vmvupWi3vv6Oubuy5WGFX33w)
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
        }
```
[Answer]
# 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](https://play.rust-lang.org/?version=stable&mode=debug&edition=2018&gist=a9b034495766443fdafea121d4fcd367)
[Answer]
# Java 10, ~~194~~ ~~193~~ ~~191~~ ~~190~~ ~~184~~ ~~182~~ 171 bytes
```
M->C->{for(var q:C){int n,X=q[0],Y=q[1],x,y,c=0;do{c++;x=n=0;for(var r:M){y=0;for(int $:r)r[y]+=Math.hypot(x-X,y++-Y)<2?++n/n:0;x++;}}while((M[X][Y]+=~n)>=n);M[X][Y]+=c;}}
```
Iterative port of [*@Arnauld*'s JavaScript answer](https://codegolf.stackexchange.com/a/180205/52210).  
-17 bytes thanks to *@Arnauld*.
Modifies the input-matrix instead of returning a new one to save bytes.
[Try it online.](https://tio.run/##nVTRjppAFH33KyakDxCuFHDb3cJi05g06YNpk33REB6miBWLAzsMroSwv25HGAUsbkmDUe@de8859wwzW7zH4@3q99GPcJqiOQ5JMUIoZZiFPtryVS1jYaStM@KzMCbaV/HnMSTM9VwP@opmMUmzXUDPRdMpWiPnOB9PZ@NpsY6pvMcUPVszpeAViMDCeXZ1D5b8x/DgADn4jm6v4sJXVfvgEB6cu6g1V4pcJE7d7yyqUDf3VGeO2Ubb5EnM5MN4AbmqjpfKo/lZVcl7Yun2gYOV5csmjAJZnrsLz13yrleiTB2i2JeEz4uO9ojbkGQ/I26DcGMfhyu04w7JT4yG5JfrIayc3EKIBSmTSfCCxMB1FqHiE/CnhDfDElC7tTDAKEvFfhvY0KH6NGA8RMZ1pq65ZtBB/zcDugPD7OKZgCbDGAbN8L/mQDUADCOZAH8a1N6wz58TSU1lDpnkAe4b1AcuvRPe90zSzACFOYSkamgHt2Sfv4e8RKYOkxbspMq0t/hOhw86fOxjMyvlzRTGhbMcNZdIdWwqCaKVnyF@fg5whkJ@HNNVSDAvEufpKU9ZsNPijGkJP2osIrL0jSQZs5AkhkpowFj@47Qq14BiQdwLHVxkdUjEbNcsslRIalPHbyRVgk7G4JnSvmjokSlWaiUaTpIoP8vTsO8HCZPbSm7iSN8zNmzeWyL@2oM2QncrhOuNc5gbJpZueMVZmov/C6U4TzUW1xejjBWhoDz@AQ)
**Explanation:**
```
M->C->{                      // Method with two integer-matrix parameters and no return-type
  for(var q:C){              //  Loop over the coordinates:
    int n,                   //   Count integer
        X=q[0],Y=q[1],       //   The current X,Y coordinate
        x,y,                 //   Temp x,y coordinates
        c=0;                 //   Counter, starting at 0
    do{                      //   Do-while:
      c++;                   //    Increase the counter `c` by 1
      x=n=0;                 //    (Re)set both `x` and the count `n` to 0
      for(var r:M)           //    Loop over the rows `r`:
        y=0;                 //     (Re)set `y` to 0
        for(int $:r)         //     Loop over the cells of the current row:
          r[y]+=             //      Increase the value at x,y by:
            Math.hypot(      //       If the hypot (builtin for `sqrt(a*a, b*b)`) of:
              x-X,           //        the difference between `x` and `X`,
                  y++-Y)     //        and difference between `y` and `Y`
                             //        (and increase `y` by 1 afterwards with `y++`)
              <2?            //       Is smaller than 2:
                 ++n/n       //        Increase count `n` and the value at x,y both by 1
                :            //       Else:
                 0;          //        Leave the value at x,y the same by increasing by 0
       x++;}}                //     Increase `x` by 1
    while((M[X][Y]+=~n)      //    Decrease the value at X,Y by n+1
          >=n);              //    Continue the do-while if this new value is still larger
                             //    than or equal to count `n`
    M[X][Y]+=c;}}            //   Increase the value at X,Y with counter `c`
```
[Answer]
# [Common Lisp](http://www.clisp.org/), 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!](https://tio.run/##XVHRToNAEHz3KzYx0dkIEaoxMWmN/onk4Npr4I4c2Navx7m22mJI2GEYZuYW07qhnybUje0qE4MMaMXIQT/hvLGoYmOly9rMaHZQvbkI96jkLjbDKDaJLbIKIWIlrRSaBmoC1RNpTqQhWZLU7N3O3Gp4upQ5A2P1ndeua/zggme215myh2NEd5XdhtCLDTFzYmPo6JK1OgaUD5z7TeNR@RpvK75niSVn1qnW4X@HNWjGtr24c3vC7bkzBrLbY3lC7gi57PS32pcXiy7RirYZgR32wlR55gxRXuQ1abkmwRKPJwcG7JTOa/07AjNS8XwmYVVqrrMMs7yiDiP3NMAx06/HDTkFfxc8cXG8uyOXUHl5Vu10Qh@dH2HkdvEBPAkvnQ@9B4q0hSKdv5BF@nL6AQ "Common Lisp – 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))
```
[Answer]
# [Python 2](https://docs.python.org/2/), 171 bytes
```
M,p=input()
l=len
for i,j in p:
 y=[(i+y,j+x)for x in-1,0,1for y in-1,0,1if l(M)>j+x>-1<i+y<l(M[0])];x=max(M[i][j]/l(y),1);M[i][j]-=x*l(y)
 for i,j in y:M[i][j]+=x
print M
```
[Try it online!](https://tio.run/##TY1BDoMgEEX3nGKWTB1T6FLFG3gCwqKLmmIoEqMJnJ5Caho3M//9Pz8T0v5e/SPniYKyPhw7R@aUe3k2rxtYWsB6CB2DpDS3TaKliVijWIJWkiBZKf3JzuD4hGO5G1s5lMpQWAuDpo/q84wFrNGLuTuekCT2J7cq3qrF4PI5dWfaqMjCZv0OU85ca0EgDMFvV1FUMS7T4Bc "Python 2 – Try It Online")
]