diff --git a/activitysim/abm/tables/skims.py b/activitysim/abm/tables/skims.py
index 3a2245278c..877f9dd8bd 100644
--- a/activitysim/abm/tables/skims.py
+++ b/activitysim/abm/tables/skims.py
@@ -15,6 +15,8 @@
 import multiprocessing
 
 from collections import OrderedDict
+from functools import reduce
+from operator import mul
 
 import numpy as np
 import openmatrix as omx
@@ -90,7 +92,7 @@ def get_skim_info(omx_file_path, tags_to_load=None):
 
     if MAX_BLOCK_BYTES:
         max_block_items = MAX_BLOCK_BYTES // np.dtype(skim_dtype).itemsize
-        max_skims_per_block = max_block_items // np.prod(omx_shape)
+        max_skims_per_block = max_block_items // multiply_large_numbers(omx_shape)
     else:
         max_skims_per_block = num_skims
 
@@ -149,14 +151,14 @@ def block_name(block):
 def buffers_for_skims(skim_info, shared=False):
 
     skim_dtype = skim_info['dtype']
-    omx_shape = [np.float64(x) for x in skim_info['omx_shape']]
+    omx_shape = skim_info['omx_shape']
     blocks = skim_info['blocks']
 
     skim_buffers = {}
     for block_name, block_size in iteritems(blocks):
 
         # buffer_size must be int (or p2.7 long), not np.int64
-        buffer_size = int(np.prod(omx_shape) * block_size)
+        buffer_size = int(multiply_large_numbers(omx_shape) * block_size)
 
         csz = buffer_size * np.dtype(skim_dtype).itemsize
         logger.info("allocating shared buffer %s for %s (%s) matrices (%s)" %
@@ -191,7 +193,7 @@ def skim_data_from_buffers(skim_buffers, skim_info):
     for block_name, block_size in iteritems(blocks):
         skims_shape = omx_shape + (block_size,)
         block_buffer = skim_buffers[block_name]
-        assert len(block_buffer) == int(np.prod(skims_shape))
+        assert len(block_buffer) == int(multiply_large_numbers(skims_shape))
         block_data = np.frombuffer(block_buffer, dtype=skim_dtype).reshape(skims_shape)
         skim_data.append(block_data)
 
@@ -261,6 +263,10 @@ def skim_dict(data_dir, settings):
     return skim_dict
 
 
+def multiply_large_numbers(list_of_numbers):
+    return reduce(mul, list_of_numbers)
+
+
 @inject.injectable(cache=True)
 def skim_stack(skim_dict):
 
diff --git a/activitysim/abm/test/test_skims.py b/activitysim/abm/test/test_skims.py
new file mode 100644
index 0000000000..bdc786b5d8
--- /dev/null
+++ b/activitysim/abm/test/test_skims.py
@@ -0,0 +1,67 @@
+from collections import OrderedDict
+from future.utils import iteritems
+
+import numpy as np
+import pytest
+
+from activitysim.abm.tables import skims
+
+
+@pytest.fixture(scope="session")
+def matrix_dimension():
+    return 5922
+
+
+@pytest.fixture(scope="session")
+def num_of_matrices():
+    return 845
+
+
+@pytest.fixture(scope="session")
+def skim_info(num_of_matrices, matrix_dimension):
+    time_periods = ['EA', 'AM', 'MD', 'PM', 'NT']
+
+    omx_keys = OrderedDict()
+    omx_key1_block_offsets = OrderedDict()
+    omx_block_offsets = OrderedDict()
+    omx_blocks = OrderedDict()
+    omx_blocks['skim_arc_skims_0'] = num_of_matrices
+
+    for i in range(0, num_of_matrices + 1):
+        key1_name = 'm{}'.format(i // len(time_periods) + 1)
+        time_period = time_periods[i % len(time_periods)]
+
+        omx_keys[(key1_name, time_period)] = '{}__{}'.format(key1_name, time_period)
+        omx_block_offsets[(key1_name, time_period)] = (0, i)
+
+        if 0 == i % len(time_periods):
+            omx_key1_block_offsets[key1_name] = (0, i)
+
+    skim_info = {
+        'omx_name': 'arc_skims',
+        'omx_shape': (matrix_dimension, matrix_dimension),
+        'num_skims': num_of_matrices,
+        'dtype': np.float32,
+        'omx_keys': omx_keys,
+        'key1_block_offsets': omx_key1_block_offsets,
+        'block_offsets': omx_block_offsets,
+        'blocks': omx_blocks
+    }
+
+    return skim_info
+
+
+def test_multiply_large_numbers(skim_info, num_of_matrices, matrix_dimension):
+    omx_shape = skim_info['omx_shape']
+    blocks = skim_info['blocks']
+
+    for block_name, block_size in iteritems(blocks):
+        # If overflow, this number will go negative
+        assert int(skims.multiply_large_numbers(omx_shape) * block_size) == \
+               num_of_matrices * matrix_dimension ** 2
+
+
+def test_multiple_large_floats():
+    calculated_value = skims.multiply_large_numbers([6205.1, 5423.2, 932.4, 15.4])
+    actual_value = 483200518316.9472
+    assert abs(calculated_value - actual_value) < 0.0001