Program to find expected value of given equation for random numbers in Python

Suppose we have a number n. Consider x = rand() mod n, where rand() function generates integers between 0 and 10^100 (both inclusive) uniformly at random. And

$$Y = \sqrt{x+\sqrt{x+\sqrt{x+\sqrt{x+...}}}}$$

We have to find the expected value of Y. The value of n will be in range 1 and 5*10^6.

So, if the input is like n = 5, then the output will be 1.696

Mathematical Background

For a nested radical expression $Y = \sqrt{x+\sqrt{x+\sqrt{x+\sqrt{x+...}}}}$, we can solve this algebraically. If we let $Y = \sqrt{x + Y}$, then squaring both sides gives us $Y^2 = x + Y$, which simplifies to $Y^2 - Y - x = 0$. Using the quadratic formula and taking the positive root: $Y = \frac{1 + \sqrt{1 + 4x}}{2}$.

Algorithm Steps

To solve this, we will follow these steps −

• err := 2235.023971557617
• max_n := 5 * 10^6
• pref := a list initially contains a single 0
• for i in range 1 to 5 * 10^6, do
  • insert (last item of pref + (1 +(4*i + 1)^0.5) * 0.5 at the end of pref
• if n < max_n, then
  • return pref[n - 1] / n
• otherwise,
  • total :=(4 *(n - 1) + 5)^1.5 / 6 - 5^1.5 / 6 - err
  • ans := 0.5 + total /(2 * n)
  • return ans

Implementation

Let us see the following implementation to get better understanding −

def solve(n):
    err = 2235.023971557617
    max_n = 5 * 10**6

    pref = [0]
    for i in range(1, 5 * 10**6):
        pref.append(pref[-1] + (1 + (4 * i + 1)**0.5) * 0.5)

    if n < max_n:
        return pref[n - 1] / n
    else:
        total = (4 * (n - 1) + 5)**1.5 / 6 - 5**1.5 / 6 - err
        ans = 0.5 + total / (2 * n)
        return ans

n = 5
print(f"Expected value for n = {n}: {solve(n):.11f}")

Expected value for n = 5: 1.69647248786

How It Works

The algorithm uses two approaches based on the value of n:

• For small n (< 5×10^6): Pre-computes prefix sums using the formula $(1 + \sqrt{4i + 1}) \times 0.5$ for each possible value of x from 0 to i-1, then returns the average.
• For large n: Uses an analytical approximation formula to avoid memory issues with large arrays.

Testing with Different Values

def solve(n):
    err = 2235.023971557617
    max_n = 5 * 10**6

    pref = [0]
    for i in range(1, min(n + 1, 5 * 10**6)):
        pref.append(pref[-1] + (1 + (4 * i + 1)**0.5) * 0.5)

    if n < max_n:
        return pref[n - 1] / n
    else:
        total = (4 * (n - 1) + 5)**1.5 / 6 - 5**1.5 / 6 - err
        ans = 0.5 + total / (2 * n)
        return ans

# Test with different values
test_cases = [1, 2, 5, 10, 100]
for n in test_cases:
    result = solve(n)
    print(f"n = {n:3d}: Expected value = {result:.6f}")

n =   1: Expected value = 1.000000
n =   2: Expected value = 1.366025
n =   5: Expected value = 1.696472
n =  10: Expected value = 1.949359
n = 100: Expected value = 2.838568

Conclusion

This algorithm efficiently computes the expected value of nested radical expressions by using prefix sums for smaller values and analytical approximations for larger ones. The approach balances accuracy with computational efficiency for the given constraints.