Perl 6, 26 22 bytes

{sum (1..4).pick xx 1+<20}

{sum roll 1..4: 1+<20}

It's not fast...takes about nine seconds on my system.

Try it online!

Six years ago, I apparently didn't know about the roll method.

Vyxal s, 38 bits^v2, 4.75 bytes

k₂ƛ4℅

Try it Online!

Explained

k₂ƛ4℅­⁡​‎‎⁡⁠⁡‏⁠‎⁡⁠⁢‏⁠‎⁡⁠⁣‏‏​⁡⁠⁡‌⁢​‎‎⁡⁠⁤‏⁠‎⁡⁠⁢⁡‏‏​⁡⁠⁡‌⁣​‎‏​⁢⁠⁡‌­
k₂ƛ    # ‎⁡For each number in the range [1, 2**20]:
   4℅  # ‎⁢  Choose a random number between 1 and 4
# ‎⁣s flag sums
💎

Created with the help of Luminespire.

Thunno 2 S, 6 bytes

T4@ı4ɼ

Try it online!

Explanation

T4@ı4ɼ  # Implicit input
T4@     # Push 4 ** 10
   ı    # Map over the range:
    4ɼ  #  Random number from [1..4]
        # Implicit output of sum

R, 32 24 23 21 bytes

Edit: Got rid of the as.integer and used integer division %/%. Speed it up slightly.

Thanks to Alex A for the sample tip ... and Giuseppe for removing the r=

sum(sample(4,2^20,T))

Tested with

i = s = 0
repeat {
i = i + 1
print(sum(sample(4,2^20,r=T)))
s = s + system.time(sum(sample(4,2^20,r=T)))[3]
if (i == 10) break
}
print (s/10)

Outputs

[1] 2621936
[1] 2620047
[1] 2621004
[1] 2621783
[1] 2621149
[1] 2619777
[1] 2620428
[1] 2621840
[1] 2621458
[1] 2620680
elapsed 
   0.029 

For pure speed the following completes in microseconds. However I'm not sure I've got my logic correct for it. The results appear consistent with the random method. Shame it's a longer length.

sum(rmultinom(1,2^20,rep(1,4))*1:4)

Here's a timing run I did on my machine

system.time(for(i in 1:1000000)sum(rmultinom(1,2^20,rep(1,4))*1:4))
                   user                  system                 elapsed 
7.330000000000040927262 0.000000000000000000000 7.370000000000345607987 

Ti-84 Basic, 17 bytes

Total footprint - Size of program header = 17 bytes

Run Time: Unknown, estimated at 5-6 hours based on performance for smaller numbers of rolls (so, basically, not very good)

Σ(randInt(1,4),A,1,2^20

VBA, 40 Bytes

for i=1 to 2^20:q=q+int(rnd*4)+1:next:?q

Javascript, 55 53 50 47 41 bytes

for(a=i=1<<20;i--;)a+=(Math.random()*4)|0

I didn't realize that non-random numbers were a known irritant, so I figure that I ought to post a real solution. Meant no disrespect.

Commentary: as noted by others above you can skip the +1 to each roll by starting off with the number of rolls in your answer, and by not having to write a=0,i=1<<20 you save two bytes, and another 2 because you don't add +1 to each roll. The parseInt function does the same thing as Math.floor but is 2 characters shorter.

Python 2, 69 67 Bytes

import random
print sum([random.randint(1,4)for _ in range(4**10)])

11 9 bytes more that the existing Python 2 answer :( I'm sure I can find some places to shave off a few bytes though.

C++, 154 bytes

Golfed:

#include<iostream>
#include <stdlib.h>
#include <time.h>
int main(){int c=0;srand(time(NULL));for(int x=0;x<=1048576;x++){c=c+(rand()%4+1);}std::cout<<c;}

Ungolfed:

#include<iostream>
#include <stdlib.h>
#include <time.h>

int main(){
    int c = 0;
    srand (time(NULL));
    for(int x = 0; x <= 1048576; x++){
        c = c + (rand() % 4 + 1);
    }

    std::cout<<c;
}

PARI/GP, 25 bytes

Really, no need for golfing here -- this is the straightforward way of doing the calculation in GP. It runs in 90 milliseconds on my machine. Hoisting the +1 saves about 20 milliseconds.

sum(i=1,2^20,random(4)+1)

Just for fun: if one were optimizing for performance in PARI,

inline long sum32d4(void) {
  long n = rand64();
  // Note: __builtin_popcountll could replace hamming_word if using gcc
  return hamming_word(n) + hamming_word(n & 0xAAAAAAAAAAAAAAAALL);
}

long sum1048576d4(void) {
  long total = 0;
  int i;
  for(i=0; i<32768; i++) total += sum32d4();
  return total;
}

has a very small total operation count -- if xorgens needs ~27 cycles per 64-bit word (can anyone verify this?), then a processor with POPCNT should take only about 0.5 cycle/bit, or a few hundred microseconds for the final number.

This should have close-to-optimal worst-case performance among methods using random numbers of similar or higher quality. It should be possible to greatly increase average speed by combining cases -- maybe a million rolls at a time -- and selecting with (essentially) arithmetic coding.

C# - 91

I went the function route for this

Func<int,int>G=s=>{int e=s=0;var r=new Random();for(;s<1<<20;s++)e+=r.Next(1,5);return e;};

Java (Java 8) - 52

int f(){return new Random().ints(1<<20,1,5).sum();}

Swift, 64 bytes

Nothing clever, golfing in Swift is hard...

func r()->Int{var x=0;for _ in 0..<(2<<19) {x+=Int(arc4random()%4)+1;};return x;}

Version 2 (too late)

var x=0;for _ in 0..<(2<<19){x+=Int(arc4random()%4)+1;};print(x)

Swift, 69 bytes

var b={return[UInt32](0..<8<<17).reduce(0){$1*0+$0+arc4random()%4+1}}

$1*0 is required for the compiler to infer the type of the closure.

Ruby, 51 47 chars

x=[];(2**20).times{x<<rand(4)+1};p x.inject(:+)

I looked at all of the answers before I did this, and the sum(2**20 times {randInt(4)}) strategy really stuck out, so I used that.

><>, 76 chars

012a*&>2*&1v
|.!33&^?&:-<
3.v < >-:v >
   vxv1v^<;3
  1234    n+
  >>>> >?!^^

I'm not sure if this one works, because my browser crashed when I tried to test it, but here's the online interpreter.

Scratch, 102 bytes

The consensus seems to be to convert it to the plaintext form used by the scratchblocks code on the scratch forums, so here we go.

when green flag clicked
set [a v] to [0]
repeat (1048576)
change [a v] by (pick random (1) to (4))
end

Result will be stored in the variable a.

Microscript, 7 bytes

First, a disclaimer: This is not actually a competing entry for this challenge, as this language is significantly too new.

Now, here's the program: 20ec1r4

K, 13 bytes

+/1+(_2^20)?4

Take the sum (+/) of one plus a vector of 2^20 (_2^20) numbers from [0,4) (?4).

The floor _ is necessary because ^ returns a float, which is slightly inconvenient for this problem.

Tested with Kona. On my machine, this runs in roughly 24 milliseconds.

edit:

To test timing, you can prefix the program with \t and get the runtime in milliseconds.

  \t +/1+(_2^20)?4
25

PowerShell, 41 37 bytes

1..1mb|%{(get-random)%4+1}|measure -s

Took my machine 2 minutes 40 seconds

Ruby, 32 bytes

(1..2**20).inject{|x|x-~rand(4)}

In a more readable form:

(1..2**20).inject(0) do |x|
  x + rand(4) + 1
end

It creates a range from 1 to 1048576 and then iterates over the block that many times. Each time the block is executed the value from the previous iteration is passed in as x (initially 0, the default for inject). Each iteration it calculates a random number between 0 and 3 (inclusive), adds one so it simulates rolling a d4 and adds that to the total.

On my machine it's pretty fast to run (0.25 real, 0.22 user, 0.02 sys).

If you've got Ruby installed you can run it with ruby -e 'p (1..2**20).inject{|x|x+rand(4)+1}' (the p is necessary to see the output when run in this manner, omit it if you don't care for that or just run it inside IRB where the result is printed to the screen for you). I've tested it on Ruby 2.1.6.

Thanks to histocrat for the bit twiddling hack that replaces x + rand(4) + 1 with x-~rand(4).

Mathematica, 30 27 bytes

Tr[RandomInteger[3,2^20]+1]

Mathematica has quite long function names...

JavaScript (ES6), 54 bytes

Average time < 100 msec. Run snippet to test (in Firefox)

// This is the answer
f=t=>(i=>{for(t=i;i--;)t+=Math.random()*4|0})(1<<20)|t

// This is the test
test();

function test(){
  var time = ~new Date;
  var tot = f();
  time -= ~new Date;
  
  Out.innerHTML = "Tot: " + tot + " in msec: " + time + "\n" + Out.innerHTML;
}

<button onclick="test()">Repeat test</button><br>
<pre id=Out></pre>

Explanation

With no statistical package built-in, in Javascript the shortest way to obtain the sum of 1 million random number is to call random() for a million times. So simply

f=()=>{
   var t = 0, r, i
   for (i=1<<20; i--; ) 
   {
      r = Math.random()*4 // random number between 0 and 3.9999999
      r = r + 1 // range 1 ... 4.999999
      r = r | 0 // truncate to int, so range 1 ... 4
      t = t+r
   }
   return t
}

Now, adding 1 for a million times is exactly the same than adding 1 million, or even better, start the sum with 1 million and then add the rest:

f=()=>{
   var t, r, i
   for (t = i = 1<<20; i--; ) 
   {
      r = Math.random()*4 // random number between 0 and 3.9999999
      r = r | 0 // truncate to int, so range 0 ... 3
      t = t+r
   }
   return t
}

Then golf, drop the temp variable r and drop the declaration of local variables. t is a parameter, as one is needed to shorten the declaration of f. i is global (bad thing)

f=t=>{
   for(t=i=1<<20;i--;) 
      t+=Math.random()*4|0
   return t
}

Then find a way to avoid 'return' using a nameless inner function. As a side effect, we gain another parameter so no globals used

f=t=>(
  (i=>{ // start inner function body
     for(t=i;i--;)t=t+Math.random()*4|0 // assign t without returning it
   })(1<<20) // value assigned to parameter i
  | t // the inner function returns 'undefined', binary ored with t gives t again
) // and these open/close bracket can be removed too

Commodore Basic, 37 bytes

1F┌I=1TO2↑20:C=C+INT(R/(1)*4+1):N─:?C

PETSCII substitutions: ─ = SHIFT+E, / = SHIFT+N, ┌ = SHIFT+O

Estimated runtime based on runs with lower dice counts: 4.25 hours.

It's tempting to try to golf off two bytes by making C an integer, getting implicit rounding of the random numbers. However, the range on integers in Commodore Basic is -32678 to 32767 -- not enough, when the median answer is 2621440.

Javascript + ActionScript2, 67 bytes

According to the ActionScript3 documentation, it is only required to use the keyword var when setting a type.

Other than that, it pretty much looks like Javascript.
Since I'm a fan of polyglots, I did this one.

The code:

function r(){i=a=0;while(i++<1<<20)a+=1+Math.random()*4&7;return a}

The &7 has a function to convert the number to an integer value. I could use other methods, but the gain in size would be null.

Test it here:

function r(){i=a=0;while(i++<1<<20)a+=1+Math.random()*4&7;return a}

document.write(r());

This answer is based on Brian's answer in ActionScript3, based in a comment I left with a solution for him.

PostgreSQL, 59 bytes

select sum(ceil(random()*4)) from generate_series(1,1<<20);

I'll admit to the slight problem that random() could, in theory, produce exactly zero, in which case the die roll would be zero.

PHP, 38 37 bytes

This uses a very simple idea: sum them all!

Also, I've noticed that 1048576 is 10000000000000000000 in binary, equivalent to 1<<20.

Here's the code:

while($i++<1<<20)$v+=rand(1,4);echo$v

Test in your browser (with VERY LITTLE changes):

$i=$v=0;while($i++<1<<20)$v+=rand(1,4);printf($v);

//'patch' to make a TRUE polyglot, to work in JS

if('\0'=="\0")//this will be false in PHP, true in JS
{
  //rand function, takes 2 parameters
  function rand($m,$n){
    return ((Math.random()*($n-$m+1))+$m)>>0;
    
    /*
     *returns an integer number between $m and $n
     *example run, with rand(4,9):
     *(Math.random()*(9-4+1))+4
     *
     *if you run Math.random()*(9-4+1),
     *it will returns numbers between 0 and 5 (9-4+1=6, but Math.random() never returns 1)
     *but the minimum is 4, so, we add it in the end
     *this results in numbers between 0+4 and 4+5
     *
     */
    
  }
  
  //for this purpose, this is enough
  function printf($s){
    document.write($s);
  }
}


/*
 *Changes:
 *
 *- instead of echo, use printf
 *    printf outputs a formatted string in php
 *    if I used echo instead, PHP would complain
 *- set an initial value on $i and $v
 *    this avoids errors in Javascript
 *
 */


$i=$v=0;while($i++<1<<20)$v+=rand(1,4);printf($v);

All the changes in the code are explained in comments.

C#: 105 bytes

using System.Linq;class C{int D(){var a=new System.Random();return new int[1<<20].Sum(i=>a.Next(1,5));}}

Lua 5.1, 51 bytes :)

a=0 for _=1,2^20 do a=a+math.random(4) end print(a)

math.random(n) already returns a number between 1 and n, so I didn't use @edc65's optimization (clever though it is). Another 26 bytes is required to properly seed the RNG, however....

Perl - 48 44 37 39 34 bytes

$-+=rand(4)+1for(1..2**20);print$-

Prints the sum without a trailing newline.
Saved 4 bytes by substituting for 2**20 (thanks Maltysen) and removing quotes around print.
Saved another 7 bytes by rearranging the code (thanks Thaylon!)
Lost 2 bytes because my old code generated 0-4 (it should be 1-4).
Once again, saved 5 bytes thanks to Caek and nutki.

Ungolfed, properly written code:

my $s = 0
$s += int( rand(4) + 1 ) for (1 .. 2**20);
print "$s";

C#, 148 (console app)

using R=System.Random;using C=System.Console;class P{static void Main(){int o=0;R r=new R();for(int i=0;i<1048576;i++){o+=r.Next(1,4);}C.Write(o);}}

Actionscript 3, 81 77 75 bytes:

function r(){var a=0,i=1<<20;while(i--){a+=int(Math.random()*4+1)}return a}

CJam, 12 11 bytes

YK#_{4mr+}*

This is pretty straight foward:

YK                  e# Y is 2, K is 20
  #                 e# 2 to the power 20
   _                e# Copy this 2 to the power 20. The first one acts as a base value
    {    }*         e# Run this code block 2 to the power 20 times
     4mr            e# Get a random int from 0 to 3. 0 to 3 works because we already have
                    e# 2 to the power 20 as base value for summation.
        +           e# Add it to the current sum (initially 2 to the power 20)

But the beauty of this is that its really fast too! On my machine (and using the Java compiler) it takes on an average of 70 milliseconds.

The online version takes around 1.7 seconds on my machine.

Update: 1 byte saved thanks to DocMax

J (12 bytes, about 9.8 milliseconds)

+/>:?4$~2^20

I suspect this is mostly memory bandwith-limited: I can't even get it to max out a single core...

You can test this with the following code:

   timeit =: 13 : '(1000 * >./ ($/x) 6!:2"0 1 y)'
   4 20 timeit '+/>:?4$~2^20'
9.90059

This runs it in 4 groups of 20 trails, and returns the number of milliseconds of the avarage time in the quickest group. An interpreter can be found here.

C, 57 bytes

main(a,b){for(b=a=1<<20;a--;b+=rand()%4);printf("%d",b);}

This code works... once. If you ever need to roll those dice again, you'll need to put srand(time(0)) in there, adding 14 bytes.

Haskell, 73 bytes

import System.Random
f=fmap sum.(sequence.replicate(2^20))$randomRIO(1,4)

Usage:

$ ghci sourcefile.hs
ghci> f
2622130

Pyth - 9 8 bytes

Uses obvious simple method of summation of randint. Took me minute to realize 1048576 was 2^20, now I feel really stupid. Thanks to @Jakube for saving me a byte by pointing out 2^20 = 4^10.

smhO4^4T

The runtime is horrible, it has yet to finish on my computer, so there is no point running it online so here is the 2^10 one: Try it online here.

s        Summation
 m       Map
  h      Incr (accounts for 0-indexed randint)
   O4    Randint 4
  ^4T    Four raised to ten

Perl, 29

Generates a table of the required length.

print~~map{0..rand 4}1..2**20

APL, 11 10 bytes

+/?4⍴⍨2*20

This just takes the sum of an array of 2²⁰ = 1048576 random integers between 1 and 4.

+/           ⍝ Reduce by summing a
  ?          ⍝ random integer
   4⍴⍨       ⍝ array with values between 1 and 4
      2*20   ⍝ of length 2^20

You can benchmark this on TryAPL by printing the timestamp before and after. It takes about 0.02 seconds.

Saved a byte thanks to marinus and FUZxxl!

Matlab, 24

First submission ever!

sum(randi([1,4],1,2^20))

I had hoped to make use of randi([1,4],1024), which gives a matrix of 1048576 elements, but then I needed a double sum, which takes more characters than this.

Regarding the running speed mentioned in the question, timeit tells me the runtime is about 0.031 seconds. So, pretty much instant.

Java, 65

Since we have scores listed by language, why not throw Java into the mix? There's not much to golf here, just a simple loop, but I was able to squeeze a couple out of my initial attempt:

int f(){int i=1<<20,s=i;while(i-->0)s+=Math.random()*4;return s;}

Go, 87 bytes

Naive solution

import"math/rand";func r(){o,n:=0,2<<19;for i:=0;i<n;i++{o+=rand.Intn(4)};println(o+n)}

Run online here: http://play.golang.org/p/gwP5Os7_Sq

Due to the way the Go playground works you have to manually change the seed (time is always the same)

Go, 78 bytes

Golfed

import."math/rand";func r()(o int){for i:=2<<19;i>=0;i--{o+=Intn(4)+1};return}

Still working on it

Run online here http://play.golang.org/p/pCliUpu9Eq

Pyth, 10 bytes

u+GhO4^4TZ

This has slightly more bytes than @Maltysen's Pyth solution. But it runs in 8.5 seconds on my laptop, while @Maltysen's solution produced no solution in 20 minutes running time.

But still a little bit too slow for the online compiler.

Explanation

u     ^4TZ   start with G = 0, for H in 0, ... 4^10-1:
                G = 
 +GhO4              G + (rand_int(4) + 1)
             result is printed implicitly 

Clip 10, 12 bytes

r+`m[)r4}#WT

         #4T    .- 4^10 = 1048576             -.
   m[   }       .- that many...               -.
     )r4        .-          ...random numbers -.
r+`             .- sum                        -.

It takes approximately 0.6 seconds to run on my machine.

Python 2, 58 bytes

We get 1048576 random characters from the operating system, take 2 bits of each, and add them up. Using the os library seems to save a few characters over using the random library.

import os
print sum(1+ord(c)%4 for c in os.urandom(1<<20))

This takes about 0.2 seconds on my PC.