I presume you're using ArcGIS Desktop which uses python 2x, not ArcGIS Pro which uses python 3x. Side note: this is one example of why it's important to specify the software you are using.
In python 2x, when both operands are integers, division returns an integer result rounded towards minus infinity. This was changed in Python 3x so division returns a float unless you force integer division using the //
operator.
Here's a worked example to demonstrate:
import arcpy
from arcpy.sa import Raster, Int, Float, CreateConstantRaster
arcpy.CheckOutExtension('spatial')
red = Int(CreateConstantRaster(500))
nir = Int(CreateConstantRaster(1000))
#Or load from file
#Note load files into Raster objects just once,
#instead of each time you use them as you were
#in your original code
#red = Raster(r'path\to\red.tif')
#nir = Raster(r'path\to\nir.tif')
#Calculate NDVI using the python division operator instead
#of the ArcGIS "Divide" function so code is easier to read
#Integer division
ndvi = (nir - red) / (nir + red)
print ndvi.pixelType, ndvi.maximum, ndvi.minimum
#prints S32 (signed 32bit int) 0.0 0.0
#One operand floating point
ndvi = Float((nir - red)) / (nir + red)
print ndvi.pixelType, ndvi.maximum, ndvi.minimum
#prints F32 (32bit float) 0.333333343267 0.333333343267
Alternatively, you can force floating point division by using from __future__ import division
at the very top of your script:
from __future__ import division
#imports etc...
red = Raster(r'path\to\red.tif')
nir = Raster(r'path\to\nir.tif')
#Integer division
ndvi = (nir - red) / (nir + red)
print ndvi.pixelType, ndvi.maximum, ndvi.minimum
#prints F32 (32bit float) 0.333333343267 0.333333343267