How to extract the factor of polynomial in mathematica?

Questions
,
1 
f = 49667625-20256600 a1x+673365 a1x^2 + 61159350 a1y-
1335600 a1x a1y - 3444735 a1y^2 - 2105040 b0x + 566678 a1x b0x +13552 a1x^2 b0x - 3382344 a1y b0x - 61376 a1x a1y b0x +    336000 a1y^2 b0x + 23569 b0x^2 + 560 a1x b0x^2+44296 a1y b0x^2+1120 a1x a1y b0x^2 - 5684 a1y^2 b0x^2 - 9961770 b0y +4904984 a1x b0y - 201344 a1x^2 b0y - 9043902 a1y b0y +288592 a1x a1y b0y + 350400 a1y^2 b0y + 200480 b0x b0y -103352 a1x b0x b0y - 2464 a1x^2 b0x b0y + 268256 a1y b0x b0y +10520 a1x a1y b0x b0y - 25504 a1y^2 b0x b0y + 521809 b0y^2 -228976 a1x b0y^2 + 14348 a1x^2 b0y^2 + 301504 a1y b0y^2-21920 a1x a1y b0y^2 - 1856 a1y^2 b0y^2;

f /. {a1x -> 20, a1y -> 5, b0x -> 0, b0y -> 0}

As shown above, ‘a1x= 20, a1y= 5, b0x=0, b0y=0’ is a set of roots of polynomial f. My question is how to extract the factor ‘(a1x-20)(a1y-5)*b0x*b0y’ from the polynomial f and make polynomial f. be expressed as ‘f=(a1x-20)(a1y-5)*b0x*b0y*F’. In the representation, F is a polynomial related to f.

2

In principle, could use Factor or Cancel or PolynomialReduce, but in this example such a factorization does not exist -- how would one get a constant term of 49667625 if everything is multiplied by b0x*b0y?

3

What

f /. {a1x -> 20, a1y -> 5, b0x -> 0, b0y -> 0}

shows is that

f == p1 * (alx - 20) + p2 * (aly - 5) + p3 * b0x + p4 * b0y

for some polynomials p1, p2, p3, p4. The polynomials are given by

PolynomialReduce[f, {(a1x - 20), (a1y - 5), b0x, b0y}, Variables[f]]