Your Gauss-Krueger projection uses +datum=potsdam. Up to 2012, this was hard coded in proj4 to a very unprecise value using a 3-parameter-transformation.
You find more exact values for 7-parameter transformations in this topic:
http://forum.openstreetmap.org/viewtopic.php?id=12723
There is an even better ntv2-grid transformation available here (take the binary), that has to be in the same folder as your application and data, unlsss you specify full pathnames.
To compare the different possible values, I made this test batch file:
echo epsg31467-epsg4326 >out.txt
cs2cs +init=epsg:31467 +to +init=epsg:4326 31467.txt >>out.txt
echo proj-Definition epsg >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs +to +init=epsg:4326 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:4326 31467.txt >>out.txt
echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:4326 31467.txt >>out.txt
echo epsg31467-epsg3785 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3785 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3785 31467.txt >>out.txt
echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:3785 31467.txt >>out.txt
echo epsg31467-epsg3857 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3857 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3857 31467.txt >>out.txt
echo epsg31467-epsg900913 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:900913 31467.txt >>out.txt
echo epsg31467-proj900913 >>out.txt
cs2cs +init=epsg:31467 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt
echo proj31467-proj900913 >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt
with any sample Gauss-Krüger lon-lat coordinate pair in 31467.txt
With GDAL 2.2-dev I get results which are closer to those that you consider correct. The Proj.4 version in that build is proj-4.8.
gdaltransform -s_srs epsg:31466 -t_srs epsg:4326
2944378 5705784
12.3756178459526 51.3128099490148 44.7304746098816
I suggest you to check which version your library is and update if it is older as the first step.
You should also find out what is your reference implementation so that others can evaluate if that is something that is known to be good for your area (Germany?).
Best Answer
Did you define the towgs84 parameters, if yes which set you are using? Whiteout parameters the difference can grow up to 230m.
In relation to transformation with gridshift (beta2007), i detect the following possible differences
The transformation descriptions are here, the boundless parameter here, proj4js just use by default a 3 parameter set.
If you don't use the gridshift method, you should use the de_dhdn_2001 Parameter set for whole Germany.