What type of orthogonal polynomials does R use?

poly uses QR factorization, as described in some detail in this answer.

I think that what you really seem to be looking for is how to replicate the output of R's poly using python.

Here I have written a function to do that based on R's implementation. I have also added some comments so that you can see the what the equivalent statements in R look like:

import numpy as npdef poly(x, degree):    xbar = np.mean(x)    x = x - xbar    # R: outer(x, 0L:degree, "^")    X = x[:, None] ** np.arange(0, degree+1)    #R: qr(X)$qr    q, r = np.linalg.qr(X)    #R: r * (row(r) == col(r))    z = np.diag((np.diagonal(r)))      # R: Z = qr.qy(QR, z)    Zq, Zr = np.linalg.qr(q)    Z = np.matmul(Zq, z)    # R: colSums(Z^2)    norm1 = (Z**2).sum(0)    #R: (colSums(x * Z^2)/norm2 + xbar)[1L:degree]    alpha = ((x[:, None] * (Z**2)).sum(0) / norm1 +xbar)[0:degree]    # R: c(1, norm2)    norm2 = np.append(1, norm1)    # R: Z/rep(sqrt(norm1), each = length(x))    Z = Z / np.reshape(np.repeat(norm1**(1/2.0), repeats = x.size), (-1, x.size), order='F')    #R: Z[, -1]    Z = np.delete(Z, 0, axis=1)    return [Z, alpha, norm2];

Checking that this works:

x = np.arange(10) + 1degree = 9poly(x, degree)

The first row of the returned matrix is

[-0.49543369,  0.52223297, -0.45342519,  0.33658092, -0.21483446,  0.11677484, -0.05269379,  0.01869894, -0.00453516],

compared to the same operation in R

poly(1:10, 9)# [1] -0.495433694  0.522232968 -0.453425193  0.336580916 -0.214834462# [6]  0.116774842 -0.052693786  0.018698940 -0.004535159