If you only need the coefficients, you can try this:
library(data.table)
setDT(df)
dafr <- df[, as.list(lm.fit(cbind(1, b), a)$coef), by=list(c, d)]
setnames(dafr, c("c", "d", "intercept", "slope"))
# c d intercept slope
#1: 1 5 1.869449e-13 0.5
#2: 2 6 5.176935e-13 0.5
#3: 3 7 5.000000e+02 0.5
#4: 4 8 5.000000e+02 0.5