Continuous EQ defined by a curve
This is a very complicated topic. A good overview on how to do this in MATLAB is here:
https://ccrma.stanford.edu/~jos/vguitar/Fitting_Filters_Matlab.html

An FIR fit can typically done by interpolating up to a full FFT grid and adding either a minimum or a zero phase and then doing an inverse FFT to get the impulse response. You can then try different filter lengths to determine which lengths is required to create an acceptable fit.

IIR fitting is much more tricky. The above example uses MATLAB's invfreqz() which in some cases gives good results and in other only works very poorly. In general the methods are iterative, i.e. there are many design cycles and there are small modifications between design steps that are supposed to "improve" things. These often use "gradient" methods, i.e. you need to evaluate the derivatives of the error function with respect to zeros and poles (or whatever other representation you may choose). They also involve least square error methods. A lot of care has to be taken to correctly formulate the error criteria in terms of frequency spacing, weight functions, relative vs. absolute error etc. The formulation of the error function is a direct representation of what the specific application can tolerate (or not).

EDIT based on OPs question edit:
If you want to do frequency domain filtering, (i.e. FFT->mpy with Filter->inverse FFT) then you need to get all the framing, zero-padding and overlap handling correct. This due to the difference of circular vs. linear convolution and that the FFT "as is" assumes by definition periodic signals. Google "overlap-add" for more info. 

You still need to decide on the phase for your filter. Leaving the phase zero creates a non-causal filter that will create lots of wrap-around problem. At the very minimum you need to shift the impulse response so it's causal.