For the Maxima timings, I did the following:

sage: mp1 = maxima(p1)
sage: mp2 = maxima(p2)
sage: time a = mp1.gcd(mp2)

1 variable over QQ

|| || ezgcd || modular || maxima ||

p - 100

0.02

0.03

0.20

q - 1000

12.10

9.84

7.66

r - 10000

95.03

62.80

30.27

(Univariate)

p - 100

3.07

q - 1000

-

r - 2000

-

2 variables over QQ

ezgcd

modular

p - 20

0.05

0.33

q - 40

0.33

9.47

r - 100

8.71

-

s - 160

89.57

-

3 variables over QQ

ezgcd

modular

p - 20

0.30

-

q - 40

8.40

-

r - 60

92.03

-

4 variables over QQ

ezgcd

modular

maxima

p - 10

0.13

-

1.04

q - 16

0.93

-

8.08

r - 20

2.82

-

s - 30

40.06

-