Python coding with numpy sympy

Question

New to python, and I am working on a problem with the ultimate goal of being able to evaluate

$$\frac{\sqrt{x}\bigg[\prod_{k=0}^{n-1}(\sqrt{x}+2^k)- \prod_{k=0}^{n-1}(-\sqrt{x}+2^k)\bigg]}{\prod_{k=0}^{n-1}(\sqrt{x}+2^k)+ \prod_{k=0}^{n-1}(-\sqrt{x}+2^k)}$$

at whichever value of $x, n$ I choose.

First, I probably need to expand the function: how should I go about doing so? Perhaps sympy? Then, how would I evaluate the sympy function?

Answer 1

import numpy as np
def prod(n,a):
    p=1
    for k in range(n):
       p*=a+2**k
    return p
def f(n,x):
    return np.sqrt(x)*(prod(np.sqrt(x))+prod(-np.sqrt(x)))/(prod(np.sqrt(x))-prod(-np.sqrt(x)))

The first function prod computes $\Pi_{k=0}^{n-1}(a+2^k)$ for given $n$ and $a$. The function f is the one you specified, which uses prod 4 times, each time with $\pm\sqrt{x}$ as $a$ in the argument.

Answer 2

As the infinite product $\prod_{k=0}^\infty(1 + 2^{-k} a)$ converges for all $a$, it seems better numerically to compute the expression as $$\frac{\sqrt{x}\bigg[\prod_{k=0}^{n-1}(1 + 2^{-k}\sqrt{x})- \prod_{k=0}^{n-1}(1 - 2^{-k}\sqrt{x})\bigg]}{\prod_{k=0}^{n-1}(1 + 2^{-k}\sqrt{x})+ \prod_{k=0}^{n-1}(1 - 2^{-k}\sqrt{x})}$$ The following code does that

import numpy as np

def prod(a, n):
    res = 1.0
    y = a
    for k in range(n):
        res *= (1.0 + y)
        y /= 2.0
    return res

def expr(x, n):
    r = np.sqrt(x)
    u, v = prod(r, n), prod(-r, n)
    return r * (u - v)/(u + v)

if __name__ == '__main__':
    for x in range(10):
        print(expr(x, 30))

""" output ->
0.0
1.0
1.433872292244485
1.7430889640419365
2.0
2.2274334143759584
2.4351764585875912
2.6283079353312218
2.809840294960271
2.9817458083387267
"""

Answer 3

With sympy it's even easier (and more precise):

from sympy.abc import x, k, n
from sympy import Product, sqrt, simplify
nm = sqrt(x) * (Product(sqrt(x) + 2**k, (k, 0, n-1)) - Product(-sqrt(x) + 2**k, (k, 0, n-1)))
dn = (Product(sqrt(x) + 2**k, (k, 0, n-1)) + Product(-sqrt(x) + 2**k, (k, 0, n-1)))
r = nm / dn
# assign values, e.g., x=2, n=30 and simplify the expression r
r = r.subs({x:2, n:30})
simplify(r.doit(), full=True)

The above code will produce the output as the following rational number:

$\displaystyle \frac{15027793389823468799645258066452193394332497033326437886769757423624751185248020668383259645787227100364629147062572788303205899}{10480566136263082489017135219769116375568335802982506381129986143601518100707743247130303530922557209156825461255161211247319516}$

We may want to approximate the output rational number obtained with a float, using the following code:

Float(simplify(r.doit(), full=True), 20) # with precision of 20 digits
# 1.4338722922444847313