Unfortunately, there isn't really a guaranteed method to figure out what algorithm is being used. By far, the most reliable way to figure it out is by looking through the source code. If you know a particular Sage function, then you can do introspection on that. If you're lucky, then the choice of algorithm or reference will be mentioned in the docstring or the source code for that function. If you don't know exactly where an algorithm might be implemented, then you may have to do a global search on the code base. Since Sage contains other components, you may also have to delve further down to see what they do.

I don't know a good solution for this type of problem. If anyone does, I'd be more than interested in hearing it. The hardest problem is making sure that that documentation actually stays in sync with the code.

As for your example, if we look at the documentation for the factor command on a univariate polynomial over the integers:

```
sage: R.<x> = ZZ[]
sage: f = (x+1)^2
sage: f.factor?
```

we see that a choice is made between either NTL or PARI for factoring. The NTL implementation notes state that it uses either the Zassenhaus method or van Hoeij's method. PARI's documentation also states that it uses van Hoeij's method.

For multivariate factoring, Sage does not currently support it over the integers. However, there is code at #2179 which could be used to add this.