First of all, thanks for all answers and comments. Unfortunetaly, the existing tools were not fully compatible with the latest versions of QGIS and ArcGIS. Therefore I made my own solution using the tool indicated by @polygeo, the QGIS plugin from @Alexandre and the name of the algorithm (four color map) from @Jens.
Here is my code for those interested (for ArcGIS but the second part could be used in QGIS as well).
arcpy.MakeFeatureLayer_management(fc, fc[:-4]+ "_lyr" )
arcpy.AddField_management(fc[:-4] + "_lyr", "color", "SHORT")
print "field alread exists"
arcpy.CalculateField_management(fc[:-4] + "_lyr", "color", "10" , "PYTHON")
arcpy.PolygonNeighbors_analysis(fc[:-4] + "_lyr", fc[:-4] + "_tb.dbf" )
graph = 
cursor=arcpy.da.SearchCursor( fc[:-4] + "_tb.dbf" , ("src_FID","nbr_FID") )
for row in cursor:
pols = arcpy.da.UpdateCursor(fc[:-4] + "_lyr", ("OID@","color"))
colored = 
for pol in pols:
nbrs = [ second for first, second in graph if first == pol]
usedcolors = 
for nbr in nbrs:
usedcolors += [second for first, second in colored if first == nbr]
pol=[color for color in range(10) if color not in usedcolors]
Note that the algorithm does not guarantee that only 4 colors are used: though it has been proven that the solution exists, the "brute force" is necessary to achieve it. In my case, I got 7 colors which is small enough. The script could have an additional loop until the solution is found, but I need to do it for hundreds of maps and 7 colors is OK.