CIE XY chromaticities into CCT and back

I'm trying to convert CIE XY chromaticity coordinates 1931 2deg color space into CCT and back. I followed linked formulas.

McCamy function:

func xyToCCT(x: Double, y: Double) -> Int {
    let n = (x - 0.3320) / (y - 0.1858)
    return Int( (-449 * pow(n, 3)) + (3525 * pow(n, 2)) - (6823.3 * n) + 5520.33 )
}

// TEST RESULTS
print(xyToCCT(x: 0.44758, y: 0.40745)) // A   2857.13110 round to 2857 instead 2856
print(xyToCCT(x: 0.34567, y: 0.35850)) // D50 5002.09743 round to 5002 instead 5003
print(xyToCCT(x: 0.31271, y: 0.32902)) // D65 6504.38938 round to 6504 correct 6504
print(xyToCCT(x: 0.31379, y: 0.34531)) // F5  6345.90175 round to 6346 instead 6350

Kang 2002 function:

func cctToXY(k: Int) -> (x: Double, y: Double)? {
    let k = k
    
    guard k > 1667 && k < 25000 else {
        return nil
    }
    
    let cct_3 = pow(Double(k), 3)
    let cct_2 = pow(Double(k), 2)
    
    var x: Double {
        if k <= 4000 {
            return -0.2661239 * pow(10, 9) / cct_3 - 0.2343589 * pow(10, 6) / cct_2 + 0.8776956 * pow(10, 3) / Double(k) + 0.17991
        } else {
            return -3.0258469 * pow(10, 9) / cct_3 + 2.1070379 * pow(10, 6) / cct_2 + 0.2226347 * pow(10, 3) / Double(k) + 0.24039
        }
    }
    
    let x_3 = pow(x, 3)
    let x_2 = pow(x, 2)
    
    var y: Double {
        if k <= 2222 {
            return -1.1063814 * x_3 - 1.34811020 * x_2 + 2.18555832 * x - 0.20219683
        } else if k <= 4000 {
            return -0.9549476 * x_3 - 1.37418593 * x_2 + 2.09137015 * x - 0.16748867
        } else {
            return 3.0817580 * x_3 - 5.8733867 * x_2 + 3.75112997 * x - 0.37001483
        }
    }
    
    return (x: x, y: y)
}

// TEST RESULTS
print(cctToXY(k: 2856)!) //   x: 0.4470706750966438 y: 0.40750879834778053
                         // A x: 0.44758            y: 0.40745

print(cctToXY(k: 5003)!) //     x: 0.3449074537491241 y: 0.35151935324713224
                         // D50 x: 0.34567            y: 0.35850

print(cctToXY(k: 6504)!) //     x: 0.313432036002229 y: 0.323601871509382
                         // D65 x: 0.31271           y: 0.32902

print(cctToXY(k: 6350)!) //    x: 0.3158877226651461 y: 0.3259845718669351
                         // F5 x: 0.31379            y: 0.34531

My source for coordinates and CCT - Wikipedia Standard Illuminant.

Ideas about why results are incorrect? Maybe rounding errors in the official values? I haven't found an official CIE table that can refer for more precise coordinates or the same with more decimals.

Then, exist equivalent formulas for CIE XY chromaticity coordinates 1964 10deg?

As I understand it, there are a couple of different notions of CCT that are commonly used:

• CCT calculated according to the smooth spectral power distribution of an "ideal" blackbody radiator
• CCT calculated according to the uneven spectral power distribution of daylight (illuminant series D)

For blackbody calculations, xy chromaticity can be calculated precisely using an observer's colour matching functions, whereas standard illuminant D provides its own definition.

e.g. CCT 6504 K of...

• blackbody with 1931 2 degree observer = (0.313465, 0.323569)
• blackbody with 1964 10 degree observer = (0.313895, 0.324473)
• daylight according to CIE illuminant D = (0.312714, 0.329119)

Note that:

• "6504" is usually shorthand for "6500 * 1.4388 / 1.4380" due to constant c2 being revised, which gives CIE illuminant D CCT of (0.312720, 0.329125)
• I believe the white point of (0.31272, 0.32903) listed on Wikipedia is calculated according to XYZ -> xyY conversion formula, not CCT -> xy using spectral power distribution equations.

This is all to say: half the battle is knowing which values you are actually trying to match.

With my own library I've found this approach, Ohno (2013), produces accurate roundtrip conversions between xy chromaticity and CCT & Duv (within 1 K CCT and -0.03 to 0.03 Duv for values between 1,000 K - 20,000 K).