As far as I can tell, this isn’t done purely within sage, and requires a call to maxima.
If you want to divide the polynomial 
a.maxima_methods().divide(b)
The first element of the output is the quotient and the second is the remainder.
So for example
sage: R=QQ['x']
sage: a=x^210-1
sage: b=R.cyclotomic_polynomial(210)*(x-1)
sage: q,r=a.maxima_methods().divide(b)
sage: q
x^161 + 2*x^160 + 2*x^159 + x^158 - x^156 - x^155 - x^154 - x^153 - x^152 + x^150 + 2*x^149 + 2*x^148 + 2*x^147 + 2*x^146 + 2*x^145 + x^144 - x^142 - x^141 + x^139 + x^138 + x^137 + x^136 + x^135 + x^134 + x^133 + x^132 + x^131 + x^130 + x^129 + x^128 + x^127 + 2*x^126 + 3*x^125 + 3*x^124 + 2*x^123 + x^122 + x^116 + 2*x^115 + 3*x^114 + 3*x^113 + 3*x^112 + 3*x^111 + 3*x^110 + 2*x^109 + x^108 + x^105 + 2*x^104 + 2*x^103 + 2*x^102 + 2*x^101 + 2*x^100 + 2*x^99 + 2*x^98 + 2*x^97 + 2*x^96 + 2*x^95 + 2*x^94 + 2*x^93 + 2*x^92 + 2*x^91 + 2*x^90 + 2*x^89 + 2*x^88 + 2*x^87 + 2*x^86 + 2*x^85 + 2*x^84 + 2*x^83 + 2*x^82 + 2*x^81 + 2*x^80 + 2*x^79 + 2*x^78 + 2*x^77 + 2*x^76 + 2*x^75 + 2*x^74 + 2*x^73 + 2*x^72 + 2*x^71 + 2*x^70 + 2*x^69 + 2*x^68 + 2*x^67 + 2*x^66 + 2*x^65 + 2*x^64 + 2*x^63 + 2*x^62 + 2*x^61 + 2*x^60 + 2*x^59 + 2*x^58 + 2*x^57 + x^56 + x^53 + 2*x^52 + 3*x^51 + 3*x^50 + 3*x^49 + 3*x^48 + 3*x^47 + 2*x^46 + x^45 + x^39 + 2*x^38 + 3*x^37 + 3*x^36 + 2*x^35 + x^34 + x^33 + x^32 + x^31 + x^30 + x^29 + x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 - x^20 - x^19 + x^17 + 2*x^16 + 2*x^15 + 2*x^14 + 2*x^13 + 2*x^12 + x^11 - x^9 - x^8 - x^7 - x^6 - x^5 + x^3 + 2*x^2 + 2*x + 1
sage: r
0
hat tip: this question