Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
branch: master
mjd
236 lines (199 sloc) 4.2 kb
#include <stdlib.h>
#include <math.h>
#include "cf_arith.h"
typedef struct st_cf_arith {
BASE_CF;
int a, b, c, d;
int e, f, g, h;
CF x, y;
char x_done, y_done;
} *CFA;
#define cfa ((CFA)cf)
static int next_arith(CF);
static float bound(int, int);
static void egest(CFA, int);
static void ingest_x(CFA);
static void ingest_y(CFA);
#define MAX(i, j) ((i) > (j) ? (i) : (j))
#define IDIFF(p, q) ((p) == C_INF ? C_INF : \
(q) == C_INF ? C_INF : \
fabs((p)-(q)))
CF
new_arith(int a, int b, int c, int d,
int e, int f, int g, int h,
CF x, CF y)
{
CFA cf = (CFA) new_cf(sizeof(struct st_cf_arith));
if (cf == 0) return 0;
cf->a = a;
cf->b = b;
cf->c = c;
cf->d = d;
cf->e = e;
cf->f = f;
cf->g = g;
cf->h = h;
cf->x = x;
cf->x_done = 0;
cf->y = y;
cf->y_done = 0;
cf->next = next_arith;
return (CF)cf;
}
static int
next_arith(CF cf)
{
int oops = 30;
while (1) {
#ifdef DEBUG_COEFFS
fprintf(stderr, "I'm now %d %d %d %d / %d %d %d %d\n",
cfa->a, cfa->b, cfa->c, cfa->d,
cfa->e, cfa->f, cfa->g, cfa->h
);
#endif
if (cfa->e == 0 && cfa->f == 0 && cfa->g == 0 && cfa->h == 0)
return C_INF;
else {
double
b11 = bound(cfa->a, cfa->e),
b01 = bound(cfa->c, cfa->g),
b10 = bound(cfa->b, cfa->f),
b00 = bound(cfa->d, cfa->h);
int i11 = b11, i10 = b10, i01 = b01, i00 = b00;
if (i11 == i10 && i10 == i01 && i01 == i00) {
egest(cfa, i11);
return i11;
} else {
double
xw = MAX(IDIFF(b11, b01), IDIFF(b10, b00)),
yw = MAX(IDIFF(b11, b10), IDIFF(b01, b00));
#ifdef DEBUG_BOUNDS
fprintf(stderr, "bounds: %f %f %f %f\n", b00, b01, b10, b11);
fprintf(stderr, "x diffs: %f %f\n", IDIFF(b11, b01), IDIFF(b10, b00));
fprintf(stderr, "y diffs: %f %f\n", IDIFF(b11, b10), IDIFF(b01, b00));
fprintf(stderr, "x bound = %f, y bound = %f\n", xw, yw);
#endif
if (xw > yw) {
ingest_x(cfa);
} else {
ingest_y(cfa);
}
}
}
if (--oops < 1) {
fprintf(stderr, "next_arith loop ran too many times\n");
exit(1);
}
}
}
float
bound(int n, int d)
{
return d == 0 ? C_INF : (float)n/(float)d;
}
void
egest(CFA cf, int q)
{
int a = cf->a, b = cf->b, c = cf->c, d = cf->d;
cf->a = cf->e;
cf->b = cf->f;
cf->c = cf->g;
cf->d = cf->h;
cf->e = a - q*cf->e;
cf->f = b - q*cf->f;
cf->g = c - q*cf->g;
cf->h = d - q*cf->h;
#ifdef DEBUG_IO
fprintf(stderr, "Egested %d\n", q);
#endif
}
void
ingest_x(CFA cf)
{
int p;
p = next(cf->x);
if (p == C_INF) {
#ifdef DEBUG_IO
fprintf(stderr, "Ingested oo from x\n");
#endif
cf->x_done = 1;
cf->c = cf->a;
cf->d = cf->b;
cf->g = cf->e;
cf->h = cf->f;
} else {
int a = cf->a, b = cf->b, e = cf->e, f = cf->f;
#ifdef DEBUG_IO
fprintf(stderr, "Ingested %d from x\n", p);
#endif
cf->a = a*p + cf->c;
cf->b = b*p + cf->d;
cf->c = a;
cf->d = b;
cf->e = e*p + cf->g;
cf->f = f*p + cf->h;
cf->g = e;
cf->h = f;
}
}
void
ingest_y(CFA cf)
{
int p;
p = next(cf->y);
if (p == C_INF) {
#ifdef DEBUG_IO
fprintf(stderr, "Ingested oo from y\n");
#endif
cf->y_done = 1;
cf->b = cf->a;
cf->d = cf->c;
cf->f = cf->e;
cf->h = cf->g;
} else {
int a = cf->a, c = cf->c, e = cf->e, g = cf->g;
#ifdef DEBUG_IO
fprintf(stderr, "Ingested %d from y\n", p);
#endif
cf->a = a*p + cf->b;
cf->b = a;
cf->c = c*p + cf->d;
cf->d = c;
cf->e = e*p + cf->f;
cf->f = e;
cf->g = g*p + cf->h;
cf->h = g;
}
}
CF
cadd(CF x, CF y)
{
return
new_arith(0, 1, 1, 0,
0, 0, 0, 1,
x, y);
}
CF
csub(CF x, CF y)
{
return
new_arith(0, 1, -1, 0,
0, 0, 0, 1,
x, y);
}
CF
cmul(CF x, CF y)
{
return
new_arith(1, 0, 0, 0,
0, 0, 0, 1,
x, y);
}
CF
cdiv(CF x, CF y)
{
return
new_arith(0, 1, 0, 0,
0, 0, 1, 0,
x, y);
}
Jump to Line
Something went wrong with that request. Please try again.