matrix(SCALE/MAX_SCALE,0,(WIDTH*SCALE-WIDTH*MAX_SCALE)/2,0,SCALE/MAX_SCALE,(HEIGHT*SCALE-HEIGHT*MAX_SCALE)/2)

var/sfac = SCALE/MAX_SCALE
var/ofac = SCALE-MAX_SCALE
var/HW = WIDTH/2, var/HH = HEIGHT/2
matrix(sfac,0,HW*ofac,0,sfac,HH*ofac)

var/matrix/m = matrix()
var/sfac = SCALE/MAX_SCALE
m.Scale(sfac,sfac)
var/HW = WIDTH/2
var/HH = HEIGHT/2
var/ofac = SCALE-MAX_SCALE
m.Translate(HW*ofac,HH*ofac)

Square()
    // Returns the square of the src object, that is src**2. There are a few slight simplifications
    // that can be made in this method that are not possible in the multiply method, due to different
    // assumptions. In most situations, the computational difference between these two should be quite
    // slim and largely irrelevant, but if you're doing a lot of computations (especially exponentiations)
    // at once, they should be helpful.

    var
        // The final result.
        pif_LongInt/UnsignedDouble/Square = new(0)

        // Bytes of src.

        src0 =  src.block_1 & 0x00FF
        src1 = (src.block_1 & 0xFF00) >> 8
        src2 =  src.block_2 & 0x00FF
        src3 = (src.block_2 & 0xFF00) >> 8

        // Bytes of the final result's blocks. As in the multiply method, the first byte of each
        // of these is the actual value, while the second byte is a buffer until that data is
        // pushed (i.e., added) to the next block.

        byte_1 = 0
        byte_2 = 0
        byte_3 = 0
        byte_4 = 0

        // Buffer to temporarily store multiplications.

        buffer

    /*
     *
     * Performing the computation.
     *

    The idea behind the computation basically follows as above. That is, we let

        A = r**3 a_3 + r**2 a_2 + r a_1 + a_0

    where 0x00 <= a_n <= 0xFF. Then,

        A*A = r**6 (a_3*a_3) + r**5 (2 a_2*a_3) + r**4 (2 a_1*a_3 + a_2*a_2) + r**3 (2 a_0*a_3 +
              2 a_1*a_2) + r**2 (2 a_0*a_2 + a_1*a_1) + r (2 a_0*a_1) + a_0*a_0.

    and we can disregard the terms with coefficients of r**6, r**5, or r**4. Now, when performing
    these computations we have several fewer multiplications, as well as some multiplications by
    2. The multiplications by two can simply be treated as a bitshift by one to the left.

    Upper bounds at each step (i.e., after summing the terms of each given coefficient) are exactly
    the same as in multiplication. This is because when doing multiplication, we tacitly assume (for
    the purposes of upper bounds) that we're computing 0xFFFFFFFF*0xFFFFFFFF = (0xFFFFFFFF)**2.

    */

    // Terms with coefficient r**0 == 1.

    byte_1 = src0*src0
    byte_2 = (byte_1 & 0xFF00) >> 8
    FLUSH(byte_1)

    // Terms with coefficient r**1 == r.

    buffer = src0*src1 // *2

    byte_2 += (buffer & 0x00FF) << 1 // By bitshifting by one, we multiply by two.
    byte_3 =  ((byte_2 & 0xFF00) >> 8) + ((buffer & 0xFF00) >> 7) // Likewise, except we're shifting to
                                                                  // to the left by seven instead of eight.
    byte_4 = (byte_3 & 0xFF00) >> 8 // We may have pushed into byte_3's buffer, so shift accordingly.

    FLUSH(byte_2)
    FLUSH(byte_3)

    // Terms with coefficient r**2.

    buffer = src1*src1
    byte_3 += (buffer & 0x00FF)
    byte_4 += ((buffer & 0xFF00) >> 8) + ((byte_3 & 0xFF00) >> 8)

    FLUSH(byte_3)
    FLUSH(byte_4)

    buffer = src0*src2 // *2.
    byte_3 += (buffer & 0x00FF) << 1.
    byte_4 += ((byte_3 & 0xFF00) >> 8) + ((buffer & 0xFF00) >> 7)

    FLUSH(byte_3)
    FLUSH(byte_4)

    // Terms with coefficient r**3.
    byte_4 += (src0*src3) << 1 // Multiply by 2.
    FLUSH(byte_4)

    byte_4 += (src1*src2) << 1 // Multiply by 2.
    FLUSH(byte_4)

    /*
     * Glue the pieces back together, and deliver the result.
     */

    Square._SetBlock(1, byte_1 | (byte_2 << 8))
    Square._SetBlock(2, byte_3 | (byte_4 << 8))

    return Square

Square()
    return src. Multiply (src)