1 files changed, 91 insertions, 0 deletions
diff --git a/reference/C/CONTRIB/SNIP/fraction.c b/reference/C/CONTRIB/SNIP/fraction.c
new file mode 100755
index 0000000..35a4b4c
--- /dev/null
+++ b/reference/C/CONTRIB/SNIP/fraction.c
@@ -0,0 +1,91 @@
+/*
+**  FRACTION.C - Compute continued fraction series
+**
+**  cfrac() donated to the public domain by the author, Thad Smith
+**  original Fraction.C, public domain by Bob Stout, modified to use cfrac()
+*/
+
+#include <math.h>
+#include <float.h>
+
+#define MAX_LENGTH 100
+
+long double cfrac(long double x, long double *p, long double *q, int bits)
+{
+      double v;                                 /* integer in series    */
+      long double del;                          /* x - v                */
+      long double z;    /* approximated value from truncated series     */
+      long double t;                            /* temp                 */
+      long double p0 = 0.0, q0 = 0.0;           /* last p, q            */
+      long double imax;                         /* max for p, q         */
+      static long double cf[MAX_LENGTH]; /* continued fraction integers */
+      int i, j, ntimes = MAX_LENGTH;;
+
+      if (x < 0)
+            x = -x;
+      imax = floor(pow(2.0, bits)) - 1.0;
+      for (i = 0; i < ntimes; i++)
+      {
+            v = floor((double)x);
+            cf[i] = v;
+            z = cf[i];
+            *p = z; *q = 1;
+            for (j = i; j--; )
+            {
+                  z = cf[j] + 1.0/z;
+                  t = *p;
+                  *p = cf[j] * (*p) + (*q);
+                  *q = t;
+            }
+            del = x-v;
+            if (del < DBL_EPSILON)
+                  break;
+            if ((*p > imax) || (*q > imax))
+            {
+                  *p = p0;
+                  *q = q0;
+                  break;
+            }
+            else
+            {
+                  p0 = *p;
+                  q0 = *q;
+            }
+            x = 1.0 / del;
+      }
+      return (*p)/(*q);
+}
+
+/*
+**  Remove everything below this to use cfrac() as a stand-alone function
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+main (int argc, char *argv[])
+{
+      long double x;      /* value to be approximated */
+      long double r,p,q;  /* approx ratio r = p/q     */
+      int bits;           /* bits of precision        */
+
+      if (argc < 2 || argc > 3)
+      {
+            puts ("Use: FRACTION value [precision]");
+            puts ("where value = floating point value to generate "
+                  "continued fraction");
+            puts ("      precision (optional) = bits in "
+                  "numerator/denominator");
+            return 1;
+      }
+      sscanf (argv[1], "%Lf", &x);
+      if (argc == 3)
+            bits = atoi(argv[2]);
+      else  bits = 32;
+
+      cfrac(x, &p, &q, bits);
+      printf("\n[%.20Lf]\n%.0Lf/%.0Lf = %lXh/%lXh = %.20Lf\n",
+            x, p, q, (long)p, (long)q, r = p/q);
+      printf("Error = %.10Lg, (%.10Lf%%)\n", r - x, 100. * (r - x) / x);
+      return EXIT_SUCCESS;
+}