1import numpy
2import algopy
3
4def f_fcn(x):
5    A = algopy.zeros((2,2), dtype=x)
6    A[0,0] = x[0]
7    A[1,0] = x[1] * x[0]
8    A[0,1] = x[1]
9    Q,R = algopy.qr(A)
10    return R[0,0]
11
12
13# Method 1: Complex-step derivative approximation (CSDA)
14h = 10**-20
15x0 = numpy.array([3,2],dtype=float)
16x1 = numpy.array([1,0])
17yc = numpy.imag(f_fcn(x0 + 1j * h * x1) - f_fcn(x0 - 1j * h * x1))/(2*h)
18
19# Method 2: univariate Taylor polynomial arithmetic (UTP)
20ax = algopy.UTPM(numpy.zeros((2,1,2)))
21ax.data[0,0] = x0
22ax.data[1,0] = x1
23ay = f_fcn(ax)
24
25# Method 3: finite differences
26h = 10**-8
27yf = (f_fcn(x0 + h * x1) - f_fcn(x0))/h
28
29# Print results
30print('CSDA result =',yc)
31print('UTP result  =',ay.data[1,0])
32print('FD  result  =',yf)
33
34