Fix mYPowers/mZPowers typo in ApprPolynomialSpecial4::Evaluate

[CodeReclaimers]

Line 199 of ApprPolynomialSpecial4.h: Evaluate() reads z-powers from mYPowers instead of mZPowers, so the result is correct only when y == z.

Fix: mYPowers[mZDegrees[i]] → mZPowers[mZDegrees[i]]

Test (CMakeLists.txt)

cmake_minimum_required(VERSION 3.14)
project(test_issue_2_1 LANGUAGES CXX)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(GTE_ROOT "${CMAKE_CURRENT_SOURCE_DIR}/../../GTE")
add_executable(test_issue_2_1 test_issue_2_1.cpp)
target_include_directories(test_issue_2_1 PRIVATE "${GTE_ROOT}")

Test (test_issue_2_1.cpp)

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <array>
#include <vector>
#include <Mathematics/ApprPolynomialSpecial4.h>

static bool test_z_vs_y_independence()
{
    using Real = double;
    std::vector<int32_t> xDegs = { 0, 1 };
    std::vector<int32_t> yDegs = { 0, 1 };
    std::vector<int32_t> zDegs = { 0, 1 };
    gte::ApprPolynomialSpecial4<Real> fitter(xDegs, yDegs, zDegs);

    std::vector<std::array<Real, 4>> obs;
    for (Real x = 0.0; x <= 1.0; x += 0.25)
        for (Real y = 0.0; y <= 1.0; y += 0.25)
            for (Real z = 0.0; z <= 1.0; z += 0.25)
                obs.push_back({ x, y, z, 1.0 + 2.0 * x * y * z });

    fitter.Fit(obs.size(), obs.data());

    Real x = 0.5, y = 0.3, z = 0.8;
    Real expected = 1.0 + 2.0 * x * y * z;
    Real actual = fitter.Evaluate(x, y, z);
    return std::fabs(actual - expected) < 0.02;
}

static bool test_z_y_symmetric_bug()
{
    using Real = double;
    std::vector<int32_t> xDegs = { 0, 1 };
    std::vector<int32_t> yDegs = { 0, 1 };
    std::vector<int32_t> zDegs = { 0, 1 };
    gte::ApprPolynomialSpecial4<Real> fitter(xDegs, yDegs, zDegs);

    std::vector<std::array<Real, 4>> obs;
    for (Real x = 0.1; x <= 1.0; x += 0.1)
        for (Real y = 0.1; y <= 1.0; y += 0.1)
            for (Real z = 0.1; z <= 1.0; z += 0.1)
                obs.push_back({ x, y, z, x * y * z });

    fitter.Fit(obs.size(), obs.data());

    Real x = 0.5, y = 0.3, z = 0.9;
    Real val_xyz = fitter.Evaluate(x, y, z);
    Real val_xyy = fitter.Evaluate(x, y, y);
    // With the bug, these are equal because z-powers use y-power table
    return std::fabs(val_xyz - val_xyy) > 1e-10;
}

int main()
{
    int bugs = (!test_z_vs_y_independence() ? 1 : 0)
             + (!test_z_y_symmetric_bug() ? 1 : 0);
    return bugs > 0 ? EXIT_FAILURE : EXIT_SUCCESS;
}