Why does this innocent function cause a segfault?

I'm trying to code metaballs in C++/SFML, my program works just fine in a single thread. I tried to write an MRE to find the problem and here's what I got:

main.cpp

// main

#include <iostream>
#include "threader.h"

float func(float x, float y, float a, float b, float r) {
    return 1.f / sqrt((x - a)*(x - a) + (y - b)*(x - b));
}

int main () {

    sf::RenderWindow window(sf::VideoMode(2800, 1800), "");
    window.setFramerateLimit(20);
    sf::Event event{};
    threader tf(window);
    while (window.isOpen()) {
        while (window.pollEvent(event)) {
            switch (event.type) {
            case sf::Event::Closed: {
                window.close();
            }
            }
        }
        window.clear();

        tf.parallel(func);

        window.display();
    }
}

threader.h

//threader_H

#pragma once

#include <array>
#include <thread>
#include <cmath>
#include <functional>
#include <SFML/Graphics.hpp>

class threader {
public:
    int threadCount = 4;
    std::array<std::thread, 4> threads;
    sf::RenderWindow& window;
public:
    explicit threader(sf::RenderWindow& w) : window(w) {};

    void strip(int start, int end, const std::function<float(float, float, float, float, float)>& func);
    void parallel(const std::function<float(float, float, float, float, float)>& func);

};

threader.cpp

// threader_CPP

#include "threader.h"
#include <iostream>


void threader::strip (int start, int end, const std::function<float (float, float, float, float, float)> &func) {
    for (int X = start; X < end; X += 10) {
        for (int Y = 0; Y < window.getSize().y; Y += 10) {
            auto x = static_cast<float>(X);
            auto y = static_cast<float>(Y);

            x = x / 2800.f * 4 + 0 - 4 / 2.f;
            y = -((1800.f - y) / 1800 * 4 + 0 - 4 / 2.f);

            if (func(x, y, 1, 2, 3) > 1) {
                sf::CircleShape circle(20);
                circle.setPointCount(3);
                circle.setFillColor(sf::Color::Cyan);
                circle.setPosition(X, Y);
                window.draw(circle);
            }
        }
    }
}

void threader::parallel (const std::function<float (float, float, float, float, float)> &func) {
    int start = 0;
    int end = window.getSize().x / threadCount;
    for (auto& t : threads) {
        t = std::thread(&threader::strip, this, start, end, func);
        start = end;
        end += window.getSize().x / threadCount;
    }
    for (auto& t : threads) {
        t.join();
    }

}

Now for the explanation. threader is a class which has two methods: strip that calculates a function for a given strip of the window and parallel that creates threads to separately calculate my function for every strip. This code doesn't work:

But here's the catch: if I adjust the function void func(...) in main to return 1.f / sqrt((x - a) + (y - b)), everything works just fine. What is happening? How a simple calculation can cause a segfault? help please...

EDIT 1: Written in CLion, C++ 20.

EDIT 2: If anything here makes sense, please explain it to me.

Solution

Rather than try to draw to the window on multiple threads, I would instead use std algorithms to filter the points to draw circles, and draw them all on the main thread.

std::vector<std::pair<int, int>> getPoints(const sf::RenderWindow& window) {
    std::vector<std::pair<int, int>> points;
    for (int X = 0; X < window.getSize().x; X += 10) {
        for (int Y = 0; Y < window.getSize().y; Y += 10) {
            points.emplace_back(X, Y);
        }
    }
    return points;
}

template<typename F>
auto filter(F f) {
    return [f](const std::pair<int, int> & point) {
        auto x = static_cast<float>(point.first);
        auto y = static_cast<float>(point.second);

        x = x / 2800.f * 4 + 0 - 4 / 2.f;
        y = -((1800.f - y) / 1800 * 4 + 0 - 4 / 2.f);

        return (func(x, y, 1, 2, 3) > 1);
    }
}

sf::CircleShape toCircle(int X, int Y) {
    sf::CircleShape circle(20);
    circle.setPointCount(3);
    circle.setFillColor(sf::Color::Cyan);
    circle.setPosition(X, Y);
    return circle;
}

template <typename F>
void drawCircles(sf::RenderWindow& window, F f) {
    auto points = getPoints(window);
    auto end = std::remove_if(std::execution::par, points.begin(), points.end(), filter(f));
    points.erase(end, points.end());
    for (auto & [X, Y] : points) {
        window.draw(toCircle(X, Y));
    }
}

float func(float x, float y, float a, float b, float r) {
    return 1.f / sqrt((x - a)*(x - a) + (y - b)*(x - b));
}

int main () {

    sf::RenderWindow window(sf::VideoMode(2800, 1800), "");
    window.setFramerateLimit(20);
    sf::Event event{};
    threader tf(window);
    while (window.isOpen()) {
        while (window.pollEvent(event)) {
            switch (event.type) {
            case sf::Event::Closed: {
                window.close();
            }
            }
        }
        window.clear();

        drawCircles(window, func);

        window.display();
    }
}