Different Pillow versions yielding different grayscale values for the same image

While converting PNG (RGB) to grayscale using two different versions of Pillow and Python, different grayscale values were observed.

I used Pillow version 5.4.1 on Python 2.7.18, and Pillow version 11.2.1 on Python 3.12.0

On checking the code on Github it mentions in the comment that for grayscale conversion following formula is used:

Link: https://github.com/python-pillow/Pillow/blob/11.2.x/src/PIL/Image.py

I need the same grayscale values on Python 3.12.0 as that in the Python 2.7.18. Pillow version 5.4.1 is not supported in Python 3.12.0.

Pillow uses L = R * 299/1000 + G * 587/1000 + B * 114/1000

But it does not yield same grayscale pixel value for two different version of Pillow (5.4.1 and 11.2.1)

I used the formula itself instead of Pillow for grayscale conversion L = floor(R * 0.299 + G * 0.587 + B * 0.114) in Python 3.12.0

But I encountered the following problem: For pixel value 4,4,4 Pillow 5.4.1 converted it to grayscale value of 4 But the formula yielded the grayscale value of 3

How can be the result of Pillow 5.4.1 replicated in Python 3.12.0 for grayscale conversion?

Note: Since I'm using a PNG, when image is first read, the pixel values are same in case of Pillow 5.4.1 and Pillow 11.2.1

Here is the relevant part of convert.c version 5.4.x

#define L24(rgb)\
    ((rgb)[0]*19595 + (rgb)[1]*38470 + (rgb)[2]*7471)


static void
rgb2l(UINT8* out, const UINT8* in, int xsize)
{
    int x;
    for (x = 0; x < xsize; x++, in += 4)
        /* ITU-R Recommendation 601-2 (assuming nonlinear RGB) */
        *out++ = L24(in) >> 16;
}

It's similar to, but not exactly the same as, the expression given in the documentation. Using a divisor of 65536 (a power of 2) instead of 1000 lets you replace the slow integer divide with a really quick right-shift, which adds up when it's performed per-pixel.

Here's the equivalent in Python:

from PIL import Image
import numpy as np

def rgb2l(img):
    array = np.array(img.convert('RGB'))
    red = array[:,:,0]
    green = array[:,:,1]
    blue = array[:,:,2]
    gray = red.astype(np.uint32)*19595 + green.astype(np.uint32)*38470 + blue.astype(np.uint32)*7471
    gray8 = (gray >> 16).astype(np.uint8)
    return Image.fromarray(gray8)

I was curious about why the two versions differed, so I compared the sources. The function rgb2l is identical, but the macro L24 is different. Here's the one from 11.2.x:

#define L24(rgb) ((rgb)[0] * 19595 + (rgb)[1] * 38470 + (rgb)[2] * 7471 + 0x8000)

Can you spot the difference? It's that added + 0x8000 on the end, which is used for rounding - an RGB value of (0,1,0) will result in a gray of 1 instead of 0. I went through all 16777216 RGB combinations, and 8388586 were different - almost exactly 50%. In no case was the difference greater than 1.