Build bipartite graph structure – connect a[i] with all its divisors from b[].
Then find maximum matching and check whether it is perfect matching (number of edges in matching is equal to the number of pairs (if graph is directed) or to doubled number).
Arbitrary chosen Kuhn algorithm implementation here.
Upd:
@Eric Duminil made great concise Python implementation here
This approach has polynomial complexity from O(n^2) to O(n^3) depending on chosen matching algorithm and number of edges (division pairs) against factorial complexity for brute-force algorithm.