There are two atomic procedures here.
- Convert the projected coordinates from the State Plane zone (in US survey feet) to latitude-longitude on NAD83. (specifically, Lambert conformal conic)
- Convert the NAD83-based latitude-longitude values to WGS84 latitude-longitude values.
As commenters have said, try using existing software like GDAL, PROJ.4. National Geodetic Survey also has programs to do this.
For a very good document on map projections and their equations, see John P. Snyder's Map Projections: A Working Manual. (PDF)
Another document is the IOGP Guidance Note 7.2: Coordinate Conversions and Transformations including Formulas (PDF). It also talks about geographic/datum transformations like what you would need to do for step 2.
You should be aware of how any package or program you use handles geographic/datum transformations. Some hard-code a transformation directly to a geographic coordinate reference system (CRS)/datum object, others have a default, and still others have you specify it during the transformation step.
NAD83 and WGS84 were originally considered equivalent back when they were defined, but they aren't any more. There are also multiple realizations / re-adjustments for both CRS. Today they're about a meter apart in the contiguous US. Your data could really be on the original NAD83, but is more likely to be on a re-adjustment like the current NAD83 (2011) or the earlier NAD83 (HARN), etc. There are file-based methods like the US NGS NADCON and GEOCON that convert between the US realizations. There are more methods that use equations directly rather than calculating offsets from an on-disk file.
Part of deciding which transformation and method to use depends on how accurate the input data is, and how accurate you want the output data to be.
Disclosure: I work for Esri and I am on the subcommittee that maintains the EPSG Geodetic Parameter Registry.
proj.4
folder ongit
OSGeo/proj.4/src/conversions/ which can be found on Wikipedia article Geographic coordinate conversion.