Missmatch in TensorFlow and my manual calculations in a regressor CNN: rounding error?

I need to build a function that can manually compute the output of a regressor CNN given an input tensor (I cannot use TensorFlow for doing so).

To ensure my manual calculations are accurate, I am comparing them against the predictions provided by TensorFlow. Although I get identical predictions for CNNs with one and two convolutions layers, the results are not identical for CNNs with three or more convolutional layers! To further examine the situation, I am obtaining the intermediate outputs from TensorFlow and comparing it with my manual calculation in each layer (to find where the mismatch is starting from).

To further explain this situation, consider a regressor CNN with 2*30*6 input tensors, where M=2, P=30, and C=6 refers to the number of rows, columns and channels in the input tensor, respectively. This CNN includes three convolutional layers, each followed by an average pooling layer, and four fully connected layers. In the convolution layer, the stride size is 1, the padding is valid, and the kernel size is 5*5 (ks = 5).

I obtain the intermediate output of the first convolutional layer as follows:

import keras
from keras.models import model_from_json
popx = input_processor(pop, inp, result)
json_file = open('CNN_Model_Insti1_j20_k2_o1_m2_inst0.hdf5', 'r')
loaded_model_json = json_file.read()
json_file.close()
model = model_from_json(loaded_model_json)
model.load_weights(filepath='CNN_BestWeight_Insti1_j20_k2_o1_m2_inst0.hdf5')

intermediate_layer_model0 = keras.Model(inputs=model.input,
                                       outputs=model.get_layer('conv2d').output)
intermediate_output0 = intermediate_layer_model0(np.expand_dims(popx, axis=0))

The outcome computed by TensorFlow for the first row, the first column and the first neuron/channel (m=0, p=0, and c=0, respectively) is equal to 0.25817838.

My manual calculation for the first convolutional layer, and for the first row, the first column and the first neuron/channel (m=0, p=0, and c=0, respectively), is also made as follows:

m, p, c = 0, 0, 0
val = 0
for mp, pp, cp in itertools.product(range(ks), range(ks), range(C)):
    if (m-int((ks-1)/2)+mp >= 0) and (p-int((ks-1)/2)+pp >= 0) and \
            (m-int((ks-1)/2)+mp < M) and (p-int((ks-1)/2)+pp < P):
        val += weight[0][mp, pp, cp, c] * popx[m-int((ks-1)/2)+mp, p-int((ks-1)/2)+pp, cp]
val += weight[1][0]

However, my manual prediction is equal to 0.23464933590083903. Is this difference due to rounding? If so, what can I do about it?