How long will this modular exponentiaton need?

I have an algorithm for prime numbers, coded in python.
There I need to calculate exponentiations modulo large numbers.
To be more precise: For h and k let n = h* 2^k + 1 and I want to calculate 17^((n-1)/2) mod n.

It works, for example I can calculate this with k = 10^4 and h = 123456 in 2 second.
With k = 10^5 and h = 123.456, the calculation needs 20 minutes.

But I want even larger exponents, k = 10^9 for example. But I would like to know whether this is possible in a day or a week.

Can I estimate the time needed?
Or can python even show me the progress?

Here is my code:

import time
start = time.time()
k = 10**4
h = 123456

n = h*2**k+1

print(pow(17, (n-1)//2,n))
print("Time needed:", time.time() - start, "seconds.")

Edit:
I took the exponents k from 10^3 to 10^4 and measured the times it needed for the modular calculation. It’s an exponential growth. How do I estimate now how long it will take for k = 10^9?

Here is the result

How do I estimate now how long it will take for k = 10^9?

You can use numpy.polyfit to produce an equation that estimates the time required for any k.

First, we generate some data, which you’ve already done:

import time
import matplotlib.pyplot as plt
import numpy as np

x = []
y = []

i = 4 # Average each time sample across this many iterations
step = 100 # Increase k by this much between samples
end = 8001 # End of our range for k

h = 123456

# Measure the execution time for lots of k values
for k in range(1, end, step):
    total_time = 0
    # Find the average time across `i` runs to reduce some noise
    for _ in range(i):
        start = time.time()
        n = h*2**k+1
        result = pow(17, (n-1)//2,n)
        stop = time.time()
        total_time += stop-start

    # Collect results
    x.append(k)
    y.append(total_time/i) # Average time across i runs

The execution time look quadratic in your plot, so let’s fit a 2nd degree polynomial to the data using numpy.polyfit and numpy.poly1d

# Find polynomial fit for 
poly_degree = 2 # We'll plot this against the actual data to compare
xarr = np.array(x)
yarr = np.array(y)

# Determine polynomial coefficients
z = np.polyfit(xarr, yarr, poly_degree) 
print(f"Coefficients: {z}")
# Coefficients: [ 3.52708730e-08 -1.24620082e-04  8.82253129e-02]

# Make a calculator for our polynomial
p = np.poly1d(z) 
xp = np.linspace(1,x[-1],step)

# Plot the results
fig, ax = plt.subplots()
ax.plot(x, y) # Empiricle data
ax.plot(xp,p(xp),'--') # Our polynomial curve
ax.set(xlabel="k", ylabel="time (s)")
ax.grid()
plt.show()

How does the quadratic fit against the measured data?

Not bad!

Now we plug in 10**9 for k in the polynomial

k = 10**9
seconds =  z[0]*k**2 + z[1]*k + z[2]

hours = seconds/3600
days = hours/24
years = days/365

print(f"years required when k == 10**9: {years}")

Output: about 1118 years on my computer with these specific parameters. What will really happen is that your computer will fall flat on its face when trying to calculate n for k = 10**9 and the polynomial won’t actually hold.

Edit: Another thing to note is that I sampled k on the order of 10^3 so any small jitter or inaccuracy can produce an error of many thousands of years once you project all the way out to 10^9