Uiua, 10 bytes & 16 bytes

/ט-1ⁿ⇡₁99

Try it

Explanation of the simple approach:

      ⇡₁99 # generate array of 1-99
     ⁿ     # take the input to the power of each element in array
  ˜-1      # 1- elements
/×         # multiply iterations of the product

An iterative approach (which is slightly longer), but more elegant in my opinion (16 bytes):

⍥(⊙⊸××⊙⤚-⊙⊙1)∞1.

Try it

                  . # duplicate p
                 1  # initialize the product with 1
⍥(⊙⊸××⊙⤚-⊙⊙1)∞   # fixed point iteration on top stack element
           ⊙⊙1     # push a 1 to below current value and p^n
       ⊙⤚-         # subtract 1 from p^n but keep p^n below on the stack
      ×             # multiply the current value with 1-p^n
  ⊙⊸×              # compute p^(n+1) by multiplying p^n * p

Haskell (48 34 bytes)

Thanks to @Wheat Wizard for saving 14 bytes.

f p n|n>0=(1-(p^n))*f p(n-1)|1>0=1

The function f takes the arguments p (probability) and n (the number of iterations).

$$f(p, n) = \prod_{i=1}^{n}(1 - p^i)$$

MMIX, 36 bytes (9 instrs)

double inf_pochhammer(double p)

Hexdump (jxd)

00000000: c1020000 e0ff3ff0 e0013ff0 0603ff02  Ḋ£¡¡ṭ”?ṅṭ¢?ṅ©¤”£
00000010: 10020200 10010103 e403c010 5b03fffc  Ñ££¡Ñ¢¢¤ỵ¤ĊÑ[¤”‘
00000020: f8020000                             ẏ£¡¡

Disassembled:

inf_pochammer IS @
    SET  $2,$0          // m = p
    SETH $255,#3FF0     // one = 1.
    SETH $1,#3FF0       // prod = 1.
0H  FSUB $3,$255,$2     // do {x = one -. m
    FMUL $2,$2,$0       // m *.= p
    FMUL $1,$1,$3       // prod *.= x
    INCH $3,#C010
    PBNZ $3,0B          // } while x != 1.
    POP  2,0            // return prod

The reason for the termination condition being x == 1. rather than m == 0 is that it's entirely possible for the rounding mode to prevent m from ever rounding to 0, but it must eventually get small enough that x rounds to 1.0.

CASIO BASIC (CASIO fx-9750GIII), 44 bytes

?→P
P-Not P→P
Σ((-1)^N×P^(N(3N-1)÷2),N,-100,500

accurate for the most part

uses this cool algorithm i think: \$ \sum_{k=-100}^{500}(-1^N\times (P-not(P)) ^{\frac{N(3N-1)}{2}} ) \$

uses the pentagonal number theorem, N ranges from -100 to 500 and the outputs are SLIGHTLY different (due to rounding errors and stuff) so I will list them out.

input output
0.0 1
0.1 0.8900101
0.2 0.7603327959
0.3 0.6126481542
0.4 0.4518605508
0.5 0.2887880951
0.6 0.1431214821
0.7 0.04231589738
0.8 3.368005852e-03
0.9 1.286067435e-06

SAKO, 72 bytes

CZYTAJ:P
W=1
*1)W=W×(P*I-1)
POWTORZ:I=1(1)99
DRUKUJ(2,7):-W
STOP1
KONIEC

Iterates only 99 levels down, but it looks alright to me.

BQN, 12 bytes

×´1-×⍟(↕99)˜

BQN online REPL

MATL, 7 6 bytes

1 byte saved thanks to @Neil!

9W:^qp

Try it online! Or verify all test cases.

Explanation

9    % Push 9
W    % 2 raised to that
:    % Range (1-based)
^    % Implicit input. Power, elementwise
q    % Subtract 1, elementwise
p    % Product. Implicit input

05AB1E, 6 5 bytes

₄Lm<P

-1 byte thanks to @Neil.

Could be made more precise at the cost of 1 byte by changing ₄ (=1000) to a larger even integer (e.g. žm=9,876,543,210 or T°=\$10^{10}\$=10,000,000,000). (If it's an odd number, the result will be negative due to the <.)

Try it online or verify all test cases.

Explanation:

₄      # Push 1000
 L     # Pop and push a list in the range [1,1000]
  m    # Take the (implicit) input to the power of each of these integers
   <   # Subtract each by 1
    P  # Take the product of the list
       # (after which the result is output implicitly)

Charcoal, 9 8 bytes

ＩΠ⊖ＸＮ⊕…φ

Try it online! Link is to verbose version of code. Works for p up to about 0.99 with the given margin of error (up to about 0.97 with calculation error less than floating-point error). Explanation:

    Ｎ       Input `p` as a number
   Ｘ        Vectorised raised to power
      …     Range up to
       φ    Predefined variable `1000`
     ⊕      Vectorised incremented
  ⊖         Vectorised decremented
 Π          Take the product
Ｉ           Cast to string
            Implicitly print

Fig, \$8\log_{256}(96)\approx\$ 6.585 bytes

r-1^xa@3

Try it online!

Explanation: it computes, this time, a different kind of function:

$$\prod_{n=1}^{1000}(p^n-1)$$

Note that \$p<1\implies p^n-1<0\$ but if it works, it works. Since the following is true:

$$y>0,\space y\in\mathbb{N},\space\forall x(n_x\in\mathbb{R}\space\land\space n_x\nless0),\space\prod_{x=1}^{10^y}(-n_x)=\prod_{x=1}^{10^y}(n_x)$$

Explanation in code:

r-1^xa@3
r        Product of:
     a@3 range [1...1000]
   ^x    each, raised to the power of input
 -1      then for each element, do x-1

Python 3, 45, 38 bytes

f=lambda x,y=1:y>99or(1-x**y)*f(x,y+1)

Try it online!

Explanation: simply computes this function like:

$$f(x)=\prod_{y=1}^{99}1-x^y$$

-7 thanks to Neil!

PARI/GP, 21 bytes

p->prodinf(n=1,1-p^n)

Attempt This Online!

Google Sheets, 46 bytes

=reduce(1,sequence(99),lambda(p,n,p*(1-A1^n)))

Put \$p\$ in cell A1 and the formula in B1. The formula only iterates 99 levels down but that seems plenty given the updated requirements.

APL(NARS), 44 chars

{p←⍵⋄+/1,{(1+p*⍵)×(¯1*⍵)×p*2÷⍨⍵ׯ1+⍵×3}¨⍳99}

One has to calculate the sum from the positive side, from the negative side, take care that in the 0 there is only a value 1. It is possible to simplify.

  m←{p←⍵⋄+/1,{(1+p*⍵)×(¯1*⍵)×p*2÷⍨⍵ׯ1+⍵×3}¨⍳99}
  ⍪{(1⍕⍵),m⍵}¨10÷⍨¯1+⍳10
0.0 1                 
0.1 0.8900101         
0.2 0.7603327959      
0.3 0.6126481542      
0.4 0.4518605508      
0.5 0.2887880951      
0.6 0.1431214821      
0.7 0.04231589738     
0.8 0.003368005852    
0.9 0.000001286067434 

JavaScript (V8), 32 31 30 bytes

f=(p,q=p)=>p*p?-f(p*q,q)*--p:1

Try it online!

-1 byte thanks to l4m2
-1 more byte thanks to Arnauld

This code is a simple solution using the given infinite product. The only problem with this code is that it needs to stop when p is too low, and sadly, JS is not the greatest when it comes to calculate small values (e.g. 1e-323 * 0.8 == 1e-323), so you have to manipulate the value before checking it.

Jelly, 5 bytes

*Ɱ⁹CP

A monadic Link that accepts a number from \$[0,0.9]\$ and yields an estimate of the probability.

Try it online!

How?

Takes advantage of the loose requirements by taking the product of the first 256 terms of the Euler function.

*Ɱ⁹CP - Link: initial probability, P, (from [0,0.9])
  ⁹   - 256
 Ɱ    - map across {[1,256]} with:
*     -   {P} exponentiate {that}
   C  - complement {those}
    P - produce of {that}

Wolfram Language (Mathematica), 23 bytes

Product[1-#^i,{i,∞}]&

Try it online!

It doesn't get much more straightforward than this.

There is, of course, a built-in:

Wolfram Language (Mathematica), 11 bytes

QPochhammer

Try it online!

APL+WIN, 10 bytes

Prompts for probability.

×/1-⎕*⍳1E4

Try it online! Thanks to Dyalog Classic