What is the conversion formula GIMP uses to convert RGB to LCh and back?

Don’t know if this is the right place to ask. Anyway I am a programming enthusiast and this is my Stack Overflow profile, and my Code Review profile.

I am extremely intelligent and very experienced in Python and currently learning C++. I want to convert from RGB to LCh and back, as a programming challenge. I want to do this without using any color conversion libraries, and I want to convert BGR float arrays to LCh arrays. The values in my input arrays are scaled down to 0 to 1, and I want the output values to be the same.

I have written Python implementations to convert BGR to HSV/HSL and back using information from Wikipedia:

import numba as nb
import numpy as np
from typing import Callable, Tuple

@nb.njit(cache=True, fastmath=True)
def extrema(a: float, b: float, c: float) -> Tuple[float]:
    i = 2
    if b > c:
        b, c = c, b
        i = 1

    if a > b:
        a, b = b, a

    if b > c:
        b, c = c, b
        i = 0

    return i, a, c

@nb.njit(cache=True, fastmath=True)
def hue(b: float, g: float, r: float, d: float, i: float) -> float:
    if i == 2:
        h = (g - b) / (6 * d)
    elif i:
        h = 1 / 3 + (b - r) / (6 * d)
        h = 2 / 3 + (r - g) / (6 * d)

    return h % 1

@nb.njit(cache=True, fastmath=True)
def HSL_pixel(
    b: float, g: float, r: float, i: float, x: float, z: float
) -> Tuple[float]:
    s = x + z
    d = z - x
    avg = s / 2

    return (hue(b, g, r, d, i), d / (1 - abs(s - 1)), avg) if d else (0, 0, avg)

@nb.njit(cache=True, fastmath=True)
def HSV_pixel(
    b: float, g: float, r: float, i: float, x: float, z: float
) -> Tuple[float]:
    d = z - x
    return (hue(b, g, r, d, i), d / z, z) if d else (0, 0, z)

@nb.njit(cache=True, parallel=True)
def from_BGR(img: np.ndarray, mode: Callable) -> np.ndarray:
    height, width = img.shape[:2]
    out = np.empty_like(img)
    for y in nb.prange(height):
        for x in nb.prange(width):
            b, g, r = img[y, x]
            i, a, c = extrema(b, g, r)
            out[y, x] = mode(b, g, r, i, a, c)

    return out

def BGR_to_HSL(img: np.ndarray) -> np.ndarray:
    return from_BGR(img, HSL_pixel)

def BGR_to_HSV(img: np.ndarray) -> np.ndarray:
    return from_BGR(img, HSV_pixel)

@nb.njit(cache=True, parallel=True)
def to_BGR(val: np.ndarray, mode: Callable) -> np.ndarray:
    height, width = val.shape[:2]
    img = np.zeros_like(val)
    for y in nb.prange(height):
        for x in nb.prange(width):
            a, b, c = val[y, x]
            img[y, x] = mode(a, b, c)

    return img

@nb.njit(cache=True, fastmath=True)
def HSL_helper(h: float, n: float) -> float:
    k = (n + 12 * h) % 12
    return max(-1, min(k - 3, 9 - k, 1))

@nb.njit(cache=True, fastmath=True)
def BGR_from_HSL_pixel(h: float, s: float, l: float):
    a = s * min(l, 1 - l)
    return l - a * HSL_helper(h, 4), l - a * HSL_helper(h, 8), l - a * HSL_helper(h, 0)

@nb.njit(cache=True, fastmath=True)
def HSV_helper(h: float, n: float) -> float:
    k = (n + 6 * h) % 6
    return max(0, min(k, 4 - k, 1))

@nb.njit(cache=True, fastmath=True)
def BGR_from_HSV_pixel(h: float, s: float, v: float):
    a = v * s
    return v - a * HSV_helper(h, 1), v - a * HSV_helper(h, 3), v - a * HSV_helper(h, 5)

def HSL_short(h: float, s: float, l: float) -> Tuple[float]:
    return BGR_from_HSL_pixel(h, s, l) if s else (l, l, l)

def HSV_short(h: float, s: float, v: float) -> Tuple[float]:
    return BGR_from_HSV_pixel(h, s, v) if s else (v, v, v)

def HSL_to_BGR(hsl: np.ndarray) -> np.ndarray:
    return to_BGR(hsl, HSL_short)

def HSV_to_BGR(hsv: np.ndarray) -> np.ndarray:
    return to_BGR(hsv, HSV_short)

But I have failed to find the formula GIMP uses for the conversion. I learned that there are many specifications of LCh and I want the output of my code to agree with GIMP, because I use GIMP.

I have even downloaded the source code of GIMP and tried to find the code for conversion by searching the files, but there are way too many files and I failed to find the code for conversion.

Here is what I have managed so far:

XYZ_MAT = np.array(
    [[0.4124, 0.3576, 0.1805], [0.2126, 0.7152, 0.0722], [0.0193, 0.1192, 0.9505]],

@nb.njit(cache=True, fastmath=True)
def gamma_expand(c: float) -> float:
    return  c / 12.92 if c <= 0.04045 else ((c + 0.055) / 1.055) ** 2.4

@nb.njit(cache=True, fastmath=True)
def BGR_to_XYZ(bgr: np.ndarray) -> np.ndarray:
    xyz = np.zeros(3)
    rgb = xyz.copy()
    s = 0
    for i in (2, 1, 0):
        c = gamma_expand(bgr[i])
        j = 2 - i
        s += XYZ_MAT[0, j] * c
        rgb[j] = c
    xyz[0] = s
    xyz[1] = (XYZ_MAT[1] * rgb).sum()
    xyz[2] = (XYZ_MAT[2] * rgb).sum()
    return xyz

I used D65 matrix to convert BGR to XYZ, without rescaling XYZ values. I read that GIMP uses D50 in this post, so I may be using the wrong matrix. What formula exactly does GIMP use for the conversion?

The color conversions in GIMP are generally done using babl. A quick look there found this.