How does the cost of EVM memory scale?

Question

How does the cost of EVM memory scales?

I wonder how does the cost of EVM memory (not storage) scale?

Indeed, having read other QA's, the cost of memory is not a linear function, but should increase rapidly as more of it is allocated?

Having written the Floyd-Warshall algorithm in Solidity, which allocates an (n x n) matrix of uint's, I see that the gas consumption increases extremly fast. This is the case even for n = 24, where the gas usage is approximately 19085414.

pragma solidity ^0.4.11;

contract APSPBenchmark is usingOraclize {

    event OraclizeEvent0(string desc);        
    event OraclizeEvent1(string desc, int[] apsp);

    int constant INF = -1;

    function APSPBenchmark() public payable {}    

    /*
    * Local all-pairs shortest path
    */
    function localAPSP(int[] w) public {

        int[] memory apsp  = allPairsShortestPath(w);

        OraclizeEvent0("local");
        //OraclizeEvent1("local", apsp); // Infinity encoded as -1
    }

    /*
    * All-pairs shortest path
    */
    function allPairsShortestPath(int[] w) private constant returns(int[]) {
        uint i; uint j; uint k;
        uint n = babylonian(w.length);
        int[] memory d = new int[](n * n);

        for (i = 0; i < n; ++i) {
            for (j = 0; j < n; ++j) {
                d[i * n +j] = (w[i * n + j] >= 100 ? INF : w[i * n + j]);
            }   
            d[i * n + i] = 0;
        }

        for (k = 0; k < n; ++k) {
            for (i = 0; i < n; ++i) {
                for (j = 0; j < n; ++j) {
                    if (d[i * n + k] == INF || d[k * n + j] == INF) continue;
                    else if (d[i * n + j] == INF) d[i * n + j] = d[i * n + k] + d[k * n + j];
                    else d[i * n + j] = min(d[i * n +j], d[i * n + k] + d[k * n + j]);
                }
            }
        }

        return d;
    }

    /*
    * Minimum of two values
    */
    function min(int a, int b) private constant returns(int)  {
        return a < b ? a : b;
    }    

    /*
    * Babylonian sqrt 
    */
    function babylonian(uint n) private constant returns(uint) {
        uint x = n;
        uint y = 1;
        while (x > y) {
            x = (x + y) / 2;
            y = n / x;
        }
        return x;
    }
}

Answer 1

The gas cost of memory expansion is defined in the Yellow Paper as follows,

G_memory is 3, and a is the number of 256-bit words (ints) allocated.

For a=24^2 this comes to 3 * 576 + 648 = 2,376. So it looks like memory cost is not your main problem.

Answer 2

The formula in this answer is actually not the cost of expansion, it is the memory cost.

The cost of expansion is the additional memory cost for the previously untouched memory words. See Solidity docs and quote below from Appendix H in the yellow paper.

Note also that Cmem is the memory cost function (the expansion function being the difference between the cost before and after). It is a polynomial, with the higher-order coefficient divided and floored, and thus linear up to 724B of memory used, after which it costs substantially more.

(Not enough rep to comment)