It looks like your program is using an approximation based on
$q \approx w = w_s*RH$ with an approximation of Clausius-Clapeyron to find $w_s$. Looking at a few values of RH,T and P, your approximation is pretty close (+/- 5%) to an analytic answer. Based on the output you quoted it looks like you are providing incorrect values of RH. Note in the comments to your routine it says:
@param rh relative humidity (proportion, not %)
This means you need to provide the RH proportion, not the percentage. E.g. divide by 100 -- RH=1 for 100%, RH=0.5 for 50%, etc.
If you adjust your input data you should be able to use your code as-is. If you wish to compare it to something, you can reference the solution below.
If you are given $RH$ (in the range [0,1]), $T$ (K) and $p$ (Pa) you can proceed as follows.
Knowing that $$RH = \dfrac{e}{e_s},$$ $$w = \dfrac{e\ R_d}{R_v(p-e)},$$ and $$q = \dfrac{w}{w+1}$$
Then we can solve for specific humidity $q$.
Rather than combining this into a single formula and solving, it is more straightforward to present this incrementally.
First, find $e_s(T)$ where $$e_s(T) = e_{s0}\exp\left[\left(\dfrac{L_v(T)}{R_v}\right)\left(\dfrac{1}{T_0}-\dfrac{1}{T}\right)\right]$$ and then find $e$ from the first formula ($e = RH*e_s$). Then plug $e$ into the formula for $w$ and then that result into the formula for $q$.
Variables used:
$q$ specific humidity or the mass mixing ratio of water vapor to total air (dimensionless)
$w$ mass mixing ratio of water vapor to dry air (dimensionless)
$e_s(T)$ saturation vapor pressure (Pa)
$e_{s0}$ saturation vapor pressure at $T_0$ (Pa)
$R_d$ specific gas constant for dry air (J kg$^{-1}$ K$^{-1}$)
$R_v$ specific gas constant for water vapor (J kg$^{-1}$ K$^{-1}$)
$p$ pressure (Pa)
$L_v(T)$ specific enthalpy of vaporization (J kg$^{-1}$)
$T$ temperature (K)
$T_0$ reference temperature (typically 273.16 K) (K)