rational_approximation.cpp

#include <stdio.h>
#include <math.h>
#include <vector>
#include <utility>
#include <algorithm>

typedef std::pair<int,int> intpair_t;

void findRational(std::vector<intpair_t>& result, double v, double e, int qmax) {

    int q = 2;
    int p = (int) (v * (double) q);

    while (q<qmax) {

        int h1, h2;
        double hv1, hv2;

        h1 = p;
        h2 = p+1;
        hv1 = (double) h1/q;
        hv2 = (double) h2/q;
        //printf("hyp: %d/%d %d/%d %f %f\n",p,q,p+1,q,hv1,hv2);

        if (hv2 >= v && hv1 <=v) {

            if (fabs(hv1-v)<fabs(hv2-v)) {
                if (fabs(hv1-v)<e) {
                    result.push_back(intpair_t(h1,q));
                }
            } else {
                if (fabs(hv2-v)<e) {
                    result.push_back(intpair_t(h2,q));
                }
            }

            q++;
            p = (int) (v * (double) q);
        } /*else {
            if (hv2 > v) {
                p--;
            } else {
                p++;
            }
            h1 = p;
            h2 = p+1;
            hv1 = (double) h1/q;
            hv2 = (double) h2/q;
            //printf("hyp: %d/%d %d/%d %f %f\n",p,q,p+1,q,hv1,hv2);
            printf("hey!\n");
        }*/
    }

}

class RationalComparer {
public:
    RationalComparer(double v) {this->v = v;}
    bool operator() (const intpair_t& v1,const intpair_t& v2) {
        double vv1 = (double) v1.first / v1.second;
        double vv2 = (double) v2.first / v2.second;
        return fabs(vv1 - v) < fabs(vv2 - v);
    }
    double v;
};

int main(int argc, char** argv) {


    if (argc < 4) {
        printf("usage: %s value error maximum_denominator\n",argv[0]);
        printf("example: %s 3.1415925 0.000001 1000\n",argv[0]);
        return 0;
    }

    double v = atof(argv[1]);
    double e = atof(argv[2]);
    int qmax = atoi(argv[3]);

    std::vector<intpair_t> result;

    findRational(result,v,e,qmax);

    RationalComparer comparer(v);
    std::sort(result.begin(),result.end(),comparer);

    for (int i=0;i<result.size();i++) {
        intpair_t pq = result.at(i);
        double vi = (double) pq.first / pq.second;
        printf("%d/%d %f (%g)\n",pq.first,pq.second,vi,fabs(vi-v));
    }

    return 0;
}