Death By Shock Probe: That's a lot of dice

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 - 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";

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