Как вариант, можно использовать fractions (документация) - Ваш коэффициент по сути дробь, а эта библиотека (встроенная) обеспечивает поддержку арифметики рациональных чисел. С ним можно провести расчеты в формате дробей и вывести итог только в конце, чтобы минимизировать проблемы с округлением:

from itertools import count
from fractions import Fraction
from tqdm.auto import tqdm

def postponed_sieve():
    yield 2; yield 3; yield 5; yield 7;
    sieve = {}
    ps = postponed_sieve()
    p = next(ps) and next(ps)
    q = p*p
    for c in count(9, 2):
        if c in sieve:
            s = sieve.pop(c)
        elif c < q:
            yield c
            continue
        else:
            s = count(q+2*p, 2*p)
            p = next(ps)
            q = p*p
        for m in s:
            if m not in sieve:
                break
        sieve[m] = s

N = int(input('Enter limit (N): '))
prime_gen = postponed_sieve()

result = Fraction(1)
pbar = tqdm(range(1, N+1), ncols=80)
for n in pbar:
    pbar.set_description(f'Process: {float(result):.17f}')
    p = next(prime_gen)
    result *= Fraction((p-1), p)

print(f'\nFinal result: {float(result)}')

UPD: Добавил вывод процесса расчета - можно наглядно увидеть как результат медленно уменьшается к нулю.

Должен вас разочаровать, предел равен нулю. Другими словами произведение, которое вы хотите вычислить, уменьшается, хотя и очень медленно, до нуля.

Рассмотрим ln∏^(p_i-1)/_pi = ln∏1 - ¹/_pi = ∑ ln(1 - ¹/_pi).

С ростом p_i имеем ln(1 - ¹/_pi) ≈ ^-1/_pi. Это один из "замечательных" пределов.

В статье Ряд обратных простых чисел находим факт: сумма обратных простых расходится. Следовательно и сумма логарифмов стремится к -∞. А раз логарифм числа неограниченно уменьшается, то само число стремиться к нулю.

Зайдём ещё немного дальше. Из той же статьи:

∑_pi≤n ¹/_pi ≥ ln ln(n + 1) - ln ^π2/₆.

Правую часть обозначим за f(n).

Заметим, что ln 1 = 0, ln' 1 = 1, ln'' < 0. Следовательно ln(1 + z) ≤ z, если z > -1.

Тогда ln(1 - ¹/_pi) ≤ - ¹/_pi.

Перепишем логарифм произведения ещё раз:

ln∏_pi≤n ^(p_i-1)/_pi = ln∏_pi≤n (1 - ¹/_pi) = ∑_pi≤n ln(1 - ¹/_pi) ≤ ∑_pi≤n - ¹/_pi = -∑_pi≤n ¹/_pi ≤ -f(n)

Если убрать логарифм, то

∏_pi≤n ^(p_i-1)/_pi ≤ e^-f(n) = e^{ln π2/₆ - ln ln(n + 1)} = e^{ln(π2/_{6·ln(n + 1)})} = ^π2/_{6·ln(n + 1)}

P.S. Заодно я посчитал произведение для пятидесяти миллионов простых: 0.0271163310942829.... В этом числе все знаки верные. Могу рассказать как, если вам ещё интересно.

Предел произведения равен нулю. Это доказывается элементарно, если вы знаете, что такое дзета-функция Римана

Частичное определение этой функции даёт тождество Эйлера

Произведение, которое вам нужно вычислить, легко свести к тождеству Эйлера:

По определению дзета-функция в единице равна сумме гармонического ряда

Сумма этого ряда равна плюс бесконечности. Следовательно, 1/ζ(1) сходится к нулю.

Вычисление произведения

Если же на самом деле хочется помучить компьютер и повычислять произведение, то можно воспользоваться библиотекой для вычислений произвольной точности GMP. Для неё есть обёртка gmpy2:

pip install gmpy2

import gmpy2


class Calculator:
    """
    Вычисляет $\\prod_{p \\text{ prime}} (1 - \\frac{1}{p})$ для первых N простых чисел.

    Теоретически это произведение стремится к нулю, но очень медленно.

    Пример использования: вычисление произведения для первых 1,000,000 простых чисел
    с точностью 200 битов (примерно 60 десятичных цифр).
    ```python
        import gmpy2
        gmpy2.get_context().precision = 200
        calc = Calculator()
        result = calc.compute(1_000_000)  # Вычисляет произведение для первых 1,000,000 простых чисел.
        print(result) # 0.033913656528660925320555243684117535433335178185542590787103493
    ```
    """

    def __init__(self):
        self.p: gmpy2.mpz = gmpy2.mpz(2)
        self.count: int = 1
        self.x: gmpy2.mpfr = gmpy2.mpfr(0.5)

    def step(self) -> None:
        """Выполняет один шаг вычисления, находя следующее простое число и обновляя произведение."""
        self.p = gmpy2.next_prime(self.p)
        self.x *= gmpy2.mpfr(self.p - 1) / gmpy2.mpfr(self.p)
        self.count += 1

    def compute(self, target_count: int) -> gmpy2.mpfr:
        """Вычисляет произведение для первых target_count простых чисел.

        Параметры:
            target_count (int): Количество простых чисел для включения в произведение.
        """
        if target_count < self.count:
            raise ValueError(
                f"target_count {target_count} не может быть меньше текущего count {self.count}"
            )
        while self.count < target_count:
            self.step()
        return self.x


ctx = gmpy2.get_context()
ctx.precision = 200

calc = Calculator()
for count in (1, 2, 3, 100, 1000, 10_000, 100_000, 1_000_000, 10_000_000, 100_000_000):
    print(f"Вычисление приближения 1/Zeta(1) с использованием {count:_} простых чисел...")

    result = calc.compute(count)
    print(f"Последнее использованное простое число: {int(calc.p):_}")
    print(result)

Результат для первых 100 миллионов простых чисел (7 минут на 12th Gen Intel(R) Core(TM) i7-1260P)

Вычисление приближения 1/Zeta(1) с использованием 1 простых чисел...
Последнее использованное простое число: 2
0.5
Вычисление приближения 1/Zeta(1) с использованием 2 простых чисел...
Последнее использованное простое число: 3
0.33333333333333333333333333333333333333333333333333333333333344
Вычисление приближения 1/Zeta(1) с использованием 3 простых чисел...
Последнее использованное простое число: 5
0.26666666666666666666666666666666666666666666666666666666666681
Вычисление приближения 1/Zeta(1) с использованием 100 простых чисел...
Последнее использованное простое число: 541
0.088749883775329826000183072638107984947697792075234128337628271
Вычисление приближения 1/Zeta(1) с использованием 1_000 простых чисел...
Последнее использованное простое число: 7_919
0.062466592946666309983931797897919132441722728138905212403029476
Вычисление приближения 1/Zeta(1) с использованием 10_000 простых чисел...
Последнее использованное простое число: 104_729
0.048558971171659960911776653489251760039769957383772203098709548
Вычисление приближения 1/Zeta(1) с использованием 100_000 простых чисел...
Последнее использованное простое число: 1_299_709
0.039880656936185428860214645822332742439137456266629879389169669
Вычисление приближения 1/Zeta(1) с использованием 1_000_000 простых чисел...
Последнее использованное простое число: 15_485_863
0.033913656528660925320555243684117535433335178185542590787103493
Вычисление приближения 1/Zeta(1) с использованием 10_000_000 простых чисел...
Последнее использованное простое число: 179_424_673
0.029542146596289014034012454730218671076098574743083600206253598
Вычисление приближения 1/Zeta(1) с использованием 100_000_000 простых чисел...
Последнее использованное простое число: 2_038_074_743
0.026193226431561846176900186868986133139585593442336001422582174

Как показал @StanislavVolodarskiy, эта последовательность убывает как O(1/ln(P)), где P - последнее простое число в произведении. Другими словами, если вам нужно получить n точных десятичных цифр, то нужно считать до простого числа порядка e^(10^n).