Created
August 2, 2020 14:03
-
-
Save jhoenicke/fca701e9179b3d8d4ab45de6afa2b42f to your computer and use it in GitHub Desktop.
Remez Algorithm to approximate arbitrary functions by low-level polynomials on a fixed interval
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <math.h> | |
#include <gmp.h> | |
#include <stdio.h> | |
#include <stdlib.h> | |
#define DEBUG | |
/*#define MOREDEBUG*/ | |
#define BITS 1000 | |
#define DIGITS 30 //(BITS/3 - 5) | |
#define ROUND_COEFFS 48 | |
#define RELATIVE_ERROR 0 | |
#define ERROR_AT_LAST_IS_ZERO 1 | |
char number[DIGITS+3]; | |
mp_exp_t expt; | |
mpf_t PI_2; | |
mpf_t LN2; | |
#define FIRST_POLY 1 | |
#define STEP_POLY 1 | |
#define NUM_POLY 3 | |
void newton(mpf_t m, void (*f) (mpf_t, mpf_t), void (*df) (mpf_t, mpf_t)) | |
{ | |
mpf_t tmp1, tmp2; | |
int retry; | |
mpf_init(tmp1); | |
mpf_init(tmp2); | |
for (retry = 0; retry < 4; retry++) | |
{ | |
f(tmp1, m); | |
df(tmp2, m); | |
mpf_div(tmp1, tmp1, tmp2); | |
mpf_sub(m, m, tmp1); | |
} | |
} | |
void findzero(mpf_t m, mpf_t a, mpf_t b, void (*f) (mpf_t, mpf_t)) | |
{ | |
int retry; | |
int sign; | |
mpf_t val; | |
mpf_init(val); | |
f(val, b); | |
sign = (mpf_cmp_si(val, 0) > 0); | |
for (retry = 0; retry < BITS/10; retry++) | |
{ | |
mpf_add(m, a, b); | |
mpf_div_ui(m, m, 2); | |
f(val, m); | |
if ((mpf_cmp_si(val, 0) > 0) ^ sign) | |
a = m; | |
else | |
b = m; | |
} | |
} | |
void find_all_zeros(mpf_t *roots, int n, | |
void (*f) (mpf_t, mpf_t), void (*df) (mpf_t, mpf_t), | |
mpf_t a, mpf_t b, mpf_t eps) | |
{ | |
mpf_t tmp1, tmp2, tmp3; | |
int retry, i, j, done; | |
mpf_init(tmp1); | |
mpf_init(tmp2); | |
mpf_init(tmp3); | |
done = 0; | |
for (retry = 0; retry < 50 && !done; retry++) | |
{ | |
done = 1; | |
for (i = 0; i < n; i++) | |
{ | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, roots[i]); | |
printf("%d: root %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
f(tmp1, roots[i]); | |
df(tmp2, roots[i]); | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, tmp1); | |
printf("%d: f %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
mpf_get_str(number, &expt, 10, DIGITS, tmp2); | |
printf("%d: df %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
mpf_div(tmp2, tmp2, tmp1); | |
#if 0 | |
mpf_ui_div(tmp3, 1, tmp2); | |
mpf_get_str(number, &expt, 10, DIGITS, tmp3); | |
printf("%d: delta %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
for (j = 0; j < n; j++) | |
{ | |
if (i != j) | |
{ | |
mpf_sub(tmp3, roots[i], roots[j]); | |
if (mpf_cmp_si(tmp3, 0) == 0) | |
{ | |
mpf_set(tmp2, eps); | |
break; | |
} | |
mpf_ui_div(tmp3, 1, tmp3); | |
mpf_sub(tmp2, tmp2, tmp3); | |
#if STEP_POLY == 2 | |
mpf_add(tmp3, roots[i], roots[j]); | |
mpf_ui_div(tmp3, 1, tmp3); | |
mpf_sub(tmp2, tmp2, tmp3); | |
#endif | |
} | |
} | |
#if FIRST_POLY > 1 | |
mpf_ui_div(tmp3, FIRST_POLY-1, roots[i]); | |
mpf_sub(tmp2, tmp2, tmp3); | |
#endif | |
mpf_ui_div(tmp2, 1, tmp2); | |
mpf_abs(tmp3, tmp2); | |
if (mpf_cmp(tmp3, eps) > 0) | |
done = 0; | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, roots[i]); | |
printf("%d: root %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
mpf_sub(roots[i], roots[i], tmp2); | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, roots[i]); | |
printf("%d: root %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
if (mpf_cmp(roots[i], a) < 0) | |
mpf_set(roots[i], a); | |
if (mpf_cmp(roots[i], b) > 0) | |
mpf_set(roots[i], b); | |
} | |
} | |
if (done) | |
printf("Newton succeeds after %d iterations!\n", retry); | |
else | |
printf("Newton didn't succeed!\n"); | |
{ | |
int compare(const mpf_t *a, const mpf_t *b) | |
{ | |
return mpf_cmp(*a, *b); | |
} | |
qsort(roots, n, sizeof(mpf_t), | |
(int (*) (const void*, const void*)) compare); | |
} | |
} | |
void approx(int n, | |
void (*f) (mpf_t, mpf_t), | |
void (*T) (mpf_t, int, mpf_t), | |
mpf_t *xcof, | |
mpf_t *a, mpf_t error) | |
{ | |
int i, j, k; | |
int swapR[n+1]; | |
mpf_t M[n+1][n+1]; | |
mpf_t b[n+1]; | |
mpf_t fac, tmp, tmp2; | |
#ifdef MOREDEBUG | |
char number[DIGITS+3]; | |
mp_exp_t expt; | |
#endif | |
mpf_init(fac); | |
mpf_init(tmp); | |
mpf_init(tmp2); | |
/* build an equation system M*a = b */ | |
/* Here a_0 is error, a_1...,a_{n+1} are the coefficients for poly */ | |
/* The equation system M_ij is SUM a_j*T_j(xi) +/- a_{n+1} = f(xi) */ | |
/* initialize b to the f(xi) */ | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_init(b[i]); | |
f(b[i], xcof[i]); | |
} | |
for (j = 0; j < n; j++) | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_init(M[i][n-j]); | |
T(M[i][n-j], j, xcof[i]); | |
} | |
j = 1; | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_t factor; | |
mpf_init(factor); | |
f(factor, xcof[i]); | |
mpf_add_ui(factor, factor, 1); | |
mpf_init_set_si(M[i][0], j); | |
mpf_mul(M[i][0], M[i][0], factor); | |
j = -j; | |
} | |
mpf_t t; | |
mpf_set_d(t, 1/65536.); | |
mpf_add(b[1], b[1], t); | |
//mpf_set_d(t, 1/65536.); | |
//mpf_add(b[3], b[3], t); | |
#if ERROR_AT_LAST_IS_ZERO | |
mpf_init_set_si(M[n][0], 0); | |
#endif | |
//mpf_init_set_si(M[1][0], 0); | |
for (i = 0; i <= n; i++) | |
swapR[i] = i; | |
#ifdef MOREDEBUG | |
printf("INITIAL EQUATION: \n"); | |
for (k = 0; k <=n; k++) | |
{ | |
for (j = 0; j <= n; j++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, M[swapR[k]][j]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
mpf_get_str(number, &expt, 10, DIGITS, b[swapR[k]]); | |
printf("= %*sE%2d\n", DIGITS+2, number, (int)expt); | |
} | |
#endif | |
/* Produce upper triangle matrix */ | |
/* after this step M[swapR[i]][j] is an upper triangle */ | |
for (i = 0; i < n; i++) | |
{ | |
int best = i; | |
int t; | |
mpf_abs(tmp, M[swapR[i]][i]); | |
for (k = i+1; k <= n; k++) | |
{ | |
mpf_abs(tmp2, M[swapR[k]][i]); | |
if (mpf_cmp(tmp2, tmp) > 0) | |
{ | |
mpf_set(tmp, tmp2); | |
best = k; | |
} | |
} | |
#ifdef MOREDEBUG | |
printf("SELECTING %d (%d)\n", swapR[best], best); | |
#endif | |
t = swapR[i]; | |
swapR[i] = swapR[best]; | |
swapR[best] = t; | |
for (k = i+1; k <= n; k++) | |
{ | |
mpf_div(fac, M[swapR[k]][i], M[swapR[i]][i]); | |
for (j = i; j <= n; j++) | |
{ | |
mpf_mul(tmp, fac, M[swapR[i]][j]); | |
mpf_sub(M[swapR[k]][j], M[swapR[k]][j], tmp); | |
} | |
mpf_mul(tmp, fac, b[swapR[i]]); | |
mpf_sub(b[swapR[k]], b[swapR[k]], tmp); | |
} | |
#ifdef MOREDEBUG | |
printf("EQUATION step %d: \n", i); | |
for (k = 0; k <=n; k++) | |
{ | |
for (j = 0; j <= n; j++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, M[swapR[k]][j]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
mpf_get_str(number, &expt, 10, DIGITS, b[swapR[k]]); | |
printf("= %*sE%2d\n", DIGITS+2, number, (int)expt); | |
} | |
#endif | |
} | |
for (i = n; i >= 0; i--) | |
{ | |
for (j = n; j > i; j--) | |
{ | |
mpf_mul(tmp, j==0 ? error : a[n-j], M[swapR[i]][j]); | |
mpf_sub(b[swapR[i]], b[swapR[i]], tmp); | |
} | |
mpf_div(tmp, b[swapR[i]], M[swapR[i]][i]); | |
if (i == 0) | |
mpf_set(error, tmp); | |
else | |
{ | |
#if ROUND_COEFFS > 0 | |
/* round to FixedPoint representable number */ | |
mpf_mul_2exp(tmp, tmp, (ROUND_COEFFS - 16*(n-i)) + 1); | |
mpf_add_ui(tmp, tmp, 1); | |
mpf_div_ui(tmp, tmp, 2); | |
mpf_floor(tmp, tmp); | |
mpf_div_2exp(tmp, tmp, (ROUND_COEFFS - 16*(n-i))); | |
#endif | |
mpf_set(a[n-i], tmp); | |
} | |
} | |
mpf_abs(error, error); | |
} | |
void newX(mpf_t eps, int n, mpf_t a, mpf_t b, | |
void (*f) (mpf_t, mpf_t), | |
void (*df) (mpf_t, mpf_t), | |
void (*ddf) (mpf_t, mpf_t), | |
void (*T) (mpf_t, int, mpf_t), | |
void (*dT) (mpf_t, int, mpf_t), | |
void (*ddT) (mpf_t, int, mpf_t), | |
mpf_t *xcof, mpf_t *acof, mpf_t error) | |
{ | |
int i; | |
mpf_t nxcof[n]; | |
mpf_t tmp, tmp1; | |
void fun(mpf_t r, mpf_t x) | |
{ | |
int k; | |
mpf_t val; | |
mpf_init(val); | |
f(r, x); | |
for (k = 0; k < n; k++) | |
{ | |
T(val, k, x); | |
mpf_mul(val, acof[k], val); | |
mpf_sub(r, r, val); | |
} | |
} | |
void dev(mpf_t r, mpf_t x) | |
{ | |
int k; | |
mpf_t val; | |
mpf_init(val); | |
df(r, x); | |
for (k = 0; k < n; k++) | |
{ | |
dT(val, k, x); | |
mpf_mul(val, acof[k], val); | |
mpf_sub(r, r, val); | |
} | |
} | |
void ddev(mpf_t r, mpf_t x) | |
{ | |
int k; | |
mpf_t val; | |
mpf_init(val); | |
ddf(r, x); | |
for (k = 0; k < n; k++) | |
{ | |
ddT(val, k, x); | |
mpf_mul(val, acof[k], val); | |
mpf_sub(r, r, val); | |
} | |
} | |
mpf_init(tmp); | |
mpf_init(tmp1); | |
#if 1 | |
find_all_zeros(xcof, n, dev, ddev, a, b, eps); | |
#else | |
if (useNewton) | |
{ | |
for (i=0; i<n; i++) | |
{ | |
mpf_init(nxcof[i]); | |
newton(xcof[i], dev, ddev); | |
} | |
} | |
else | |
{ | |
mpf_t del; | |
mpf_t last; | |
mpf_t tmp2; | |
int sign, p; | |
mpf_init(del); | |
mpf_init(tmp2); | |
mpf_init(last); | |
mpf_sub(del, b, a); | |
mpf_div_ui(del, del, 1000); | |
mpf_mul_ui(tmp, del, 10); | |
dev(tmp1, tmp); | |
sign = (mpf_cmp_si(tmp1, 0) > 0); | |
p = 0; | |
/* mpf_set(xcof[p], a); */ | |
for (i = 10; i < 1000; i++) | |
{ | |
/* char number[DIGITS+3]; */ | |
/* mp_exp_t expt; */ | |
mpf_set(last, tmp); | |
mpf_add(tmp, tmp, del); | |
/* mpf_get_str(number, &expt, 10, DIGITS, tmp); */ | |
/* printf("d'(%*sE%2d) = ", DIGITS+2, number, (int)expt); */ | |
dev(tmp1, tmp); | |
/* mpf_get_str(number, &expt, 10, DIGITS, tmp1); */ | |
/* printf("%*sE%2d\n", DIGITS+2, number, (int)expt); */ | |
if ((mpf_cmp_si(tmp1, 0) > 0) ^ sign) { | |
/* int signtmp; */ | |
findzero(tmp1, last, tmp, dev); | |
/* fun(tmp2, tmp1); */ | |
/* signtmp = mpf_cmp_si(tmp2, 0) > 0; */ | |
/* fun(tmp2, xcof[p]); */ | |
/* if (signtmp != (mpf_cmp_si(tmp2, 0) > 0)) */ | |
/* p++; */ | |
mpf_set(xcof[p++], tmp1); | |
sign ^=1; | |
} | |
} | |
while (p < n) | |
mpf_set(xcof[p++], b); | |
} | |
#endif | |
mpf_set_si(error, 0); | |
for (i=0; i< n; i++) | |
{ | |
fun(tmp, xcof[i]); | |
mpf_abs(tmp, tmp); | |
if (mpf_cmp(tmp, error) > 0) | |
mpf_set(error, tmp); | |
} | |
} | |
void remez(int n, mpf_t a, mpf_t b, mpf_t eps, | |
void (*f) (mpf_t, mpf_t), | |
void (*df) (mpf_t, mpf_t), | |
void (*ddf) (mpf_t, mpf_t), | |
void (*T) (mpf_t, int, mpf_t), | |
void (*dT) (mpf_t, int, mpf_t), | |
void (*ddT) (mpf_t, int, mpf_t), | |
mpf_t *xcof, mpf_t *acof) | |
{ | |
int i, j; | |
mpf_t tmp; | |
mpf_t error, maxerror; | |
mpf_init(tmp); | |
mpf_init(error); | |
mpf_init(maxerror); | |
#ifdef DEBUG | |
printf(" XCOF:"); | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, xcof[i]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
printf("\n"); | |
#endif | |
for (j = 0; j < 100; j++) | |
{ | |
#ifdef DEBUG | |
printf("REMEZ iter %d: \n", j); | |
#endif | |
approx(n, f, T, xcof, acof, error); | |
#ifdef DEBUG | |
printf(" ACOF:"); | |
for (i = 0; i < n; i++) | |
{ | |
#if ROUND_COEFFS | |
mpf_mul_2exp(tmp, acof[i], ROUND_COEFFS - 16*i); | |
printf(" %ld", mpf_get_si(tmp)); | |
#else | |
mpf_get_str(number, &expt, 10, DIGITS, acof[i]); | |
printf(" %*sE%02d", DIGITS+2, number, (int)expt); | |
#endif | |
} | |
printf("\n"); | |
mpf_get_str(number, &expt, 10, DIGITS, error); | |
printf(" ERR: %*sE%02d\n", DIGITS+2, number, (int)expt); | |
#endif | |
newX(eps, n, a, b, f, df, ddf, T,dT,ddT, xcof, acof, maxerror); | |
#ifdef DEBUG | |
printf(" XCOF:"); | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, xcof[i]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
printf("\n"); | |
mpf_get_str(number, &expt, 10, DIGITS, maxerror); | |
printf("MAXERR: %*sE%2d\n", DIGITS+2, number, (int) expt); | |
#endif | |
mpf_sub(tmp, maxerror, error); | |
mpf_div(tmp, tmp, maxerror); | |
if (mpf_cmp(tmp, eps) < 0) | |
break; | |
} | |
} | |
void poly(mpf_t r, int i, mpf_t x) { | |
i = STEP_POLY*i + FIRST_POLY; | |
mpf_set_ui(r, 1); | |
while (i-- > 0) | |
mpf_mul(r, r, x); | |
} | |
void dpoly(mpf_t r, int i, mpf_t x) { | |
i = STEP_POLY*i + FIRST_POLY; | |
mpf_set_ui(r, i); | |
i--; | |
while (i-- > 0) | |
mpf_mul(r, r, x); | |
} | |
void ddpoly(mpf_t r, int i, mpf_t x) { | |
i = STEP_POLY*i + FIRST_POLY; | |
mpf_set_ui(r, i * (i-1)); | |
i -= 2; | |
while (i-- > 0) | |
mpf_mul(r, r, x); | |
} | |
void mpf_sin(mpf_t r, mpf_t x) { | |
mpf_t w, p; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_mul(w, x, x); | |
mpf_set(p, x); | |
mpf_set(r, p); | |
for (i=3; i< 40; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_div_ui(p, p,i*(i-1)); | |
mpf_neg(p, p); | |
mpf_add(r, r, p); | |
} | |
} | |
void loghelp(mpf_t r, mpf_t x) { | |
/* compute ln (1+x) - ln (1-x) - 2x */ | |
mpf_t w, p,q; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_init(q); | |
mpf_mul(w, x, x); | |
mpf_mul_ui(p, x,2); | |
mpf_set_ui(r, 0); | |
mpf_add(r, r, p); | |
for (i=3; i< 800; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_div_ui(q, p, i); | |
mpf_add(r, r, q); | |
} | |
} | |
void mpf_log(mpf_t r, mpf_t x) { | |
mpf_t p, q; | |
int i; | |
mpf_init(p); | |
mpf_init(q); | |
mpf_add_ui(p, x, 1); | |
mpf_sub_ui(q, x, 1); | |
mpf_div(p, p, q); | |
mpf_get_str(number, &expt, 10, DIGITS, p); | |
printf("p: %*sE%2d\n", DIGITS+2, number, (int)expt); | |
// p = (x+1)/(x-1) | |
loghelp(r, p); | |
mpf_get_str(number, &expt, 10, DIGITS, r); | |
printf("r: %*sE%2d\n", DIGITS+2, number, (int)expt); | |
// r = ln((1+p) / (1-p)) - 2p | |
// r = ln(x) - 2p | |
mpf_mul_ui(p, p, 2); | |
mpf_add(r, r, p); | |
} | |
void dloghelp(mpf_t r, mpf_t x) { | |
mpf_mul(r, x, x); | |
mpf_ui_sub(r, 1, r); | |
mpf_ui_div(r, 2, r); | |
//mpf_sub_ui(r, r, 2); | |
} | |
void ddloghelp(mpf_t r, mpf_t x) { | |
mpf_mul(r, x, x); | |
mpf_ui_sub(r, 1, r); | |
mpf_mul(r, r, r); | |
mpf_div(r, x, r); | |
mpf_mul_ui(r, r, 4); | |
} | |
void mpf_cos(mpf_t r, mpf_t x) { | |
mpf_t w, p; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_mul(w, x, x); | |
mpf_set_si(p, 1); | |
mpf_set(r, p); | |
for (i=2; i< 40; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_div_ui(p, p, i*(i-1)); | |
mpf_neg(p,p); | |
mpf_add(r, r, p); | |
} | |
} | |
void mpf_exp(mpf_t r, mpf_t x) { | |
mpf_t p; | |
int i; | |
mpf_init(p); | |
mpf_set_si(p, 1); | |
mpf_set(r, p); | |
for (i=1; i< 1000; i++) | |
{ | |
mpf_mul(p, p, x); | |
mpf_div_ui(p, p, i); | |
mpf_add(r, r, p); | |
} | |
} | |
void mpf_atan(mpf_t r, mpf_t x) { | |
mpf_t w, p, tmp; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_init(tmp); | |
mpf_mul(w, x, x); | |
mpf_set(p, x); | |
mpf_set(r, p); | |
for (i=3; i< 80; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_neg(p,p); | |
mpf_div_ui(tmp, p, i); | |
mpf_add(r, r, tmp); | |
} | |
} | |
void c(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_div_ui(w, w, 1048576); | |
mpf_mul_ui(w, w, 0); | |
/* mpf_div_ui(w, w, 1024); */ | |
mpf_mul(r, x, PI_2); | |
mpf_cos(r, r); | |
mpf_sub_ui(r, r, 1); | |
mpf_sub(r, r, w); | |
} | |
void dc(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, x); | |
mpf_div_ui(w, w, 1048576); | |
mpf_mul_ui(w, w, (0)*8); | |
mpf_mul(r, x, PI_2); | |
mpf_sin(r, r); | |
mpf_neg(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void ddc(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, x); | |
mpf_mul(w, w, w); | |
mpf_div_ui(w, w, 1048576); | |
mpf_mul_ui(w, w, (0)*8*7); | |
mpf_mul(r, x, PI_2); | |
mpf_cos(r, r); | |
mpf_neg(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void s(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, x); | |
mpf_div_ui(w, w, 8192*4); | |
mpf_mul_ui(w, w, 5); | |
mpf_mul(r, x, PI_2); | |
mpf_sin(r, r); | |
mpf_sub(r, r, w); | |
mpf_set_d(w, 4.6566128730773e-10/2); | |
mpf_sub(r, r, w); | |
} | |
void ds(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_div_ui(w, w, 8192*4); | |
mpf_mul_ui(w, w, 9*5); | |
mpf_mul(r, x, PI_2); | |
mpf_cos(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void dds(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_div(w, w, x); | |
mpf_div_ui(w, w, 8192*4); | |
mpf_mul_ui(w, w, 9*8*5); | |
mpf_mul(r, x, PI_2); | |
mpf_sin(r, r); | |
mpf_neg(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void at(mpf_t r, mpf_t x) { | |
mpf_atan(r, x); | |
mpf_div(r, r, PI_2); | |
} | |
void dat(mpf_t r, mpf_t x) { | |
mpf_mul(r, x, x); | |
mpf_add_ui(r, r, 1); | |
mpf_mul(r, r, PI_2); | |
mpf_ui_div(r, 1, r); | |
} | |
void ddat(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_add_ui(w, w, 1); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, PI_2); | |
mpf_div(r, x, w); | |
mpf_mul_ui(r, r, 2); | |
mpf_neg(r, r); | |
} | |
void sn(mpf_t r, mpf_t x) { | |
mpf_sin(r, x); | |
} | |
void dsn(mpf_t r, mpf_t x) { | |
mpf_cos(r, x); | |
} | |
void ddsn(mpf_t r, mpf_t x) { | |
mpf_sin(r, x); | |
mpf_neg(r, r); | |
} | |
void exphelp(mpf_t r, mpf_t x) { | |
/* x*(e+1)/(e-1) - 2 */ | |
mpf_t e, t, xln2; | |
mpf_init(e); | |
mpf_init(t); | |
mpf_exp(e, x); | |
mpf_sub_ui(t, e, 1); | |
mpf_add_ui(e, e, 1); | |
mpf_mul(e, e, x); | |
mpf_div(t, e, t); | |
mpf_sub_ui(r, t, 2); | |
} | |
void dexphelp(mpf_t r, mpf_t x) { | |
/* (ee - 2xe - 1)/(e-1)^2 */ | |
if (mpf_sgn(x) == 0) { | |
mpf_set_ui(r, 0); | |
return; | |
} | |
mpf_t e, n, xln2; | |
mpf_init(e); | |
mpf_init(n); | |
mpf_exp(e, x); | |
mpf_sub_ui(n, e, 1); | |
mpf_mul(n, n, n); | |
mpf_mul_ui(r, x, 2); | |
mpf_sub(r, e, r); | |
mpf_mul(r, e, r); | |
mpf_sub_ui(r, r, 1); | |
mpf_div(r, r, n); | |
} | |
void ddexphelp(mpf_t r, mpf_t x) { | |
/* 2e (-2e + xe + 2 + x) / (e-1)^3 */ | |
if (mpf_sgn(x) == 0) { | |
mpf_set_ui(r, 2); | |
mpf_div_ui(r, r, 3); | |
return; | |
} | |
mpf_t e, n; | |
mpf_t xln2; | |
mpf_init(e); | |
mpf_init(n); | |
mpf_exp(e, x); | |
mpf_sub_ui(n, e, 1); | |
mpf_mul(r, n, n); | |
mpf_mul(n, r, n); | |
// n = (e-1)^3 | |
mpf_mul(r, e, x); | |
mpf_add_ui(r, r, 2); | |
// r = xe + 2 | |
mpf_add(r, r, x); | |
// r = xe + 2 + x | |
mpf_add(e, e, e); | |
// e = 2e | |
mpf_sub(r, r, e); | |
mpf_mul(r, r, e); | |
mpf_div(r, r, n); | |
} | |
void exphelp2(mpf_t r, mpf_t x) { | |
mpf_t xln2; | |
mpf_init(xln2); | |
mpf_mul(xln2, x, LN2); | |
exphelp(r, xln2); | |
mpf_div(r, r, LN2); | |
} | |
void dexphelp2(mpf_t r, mpf_t x) { | |
mpf_t xln2; | |
mpf_init(xln2); | |
mpf_mul(xln2, x, LN2); | |
dexphelp(r, xln2); | |
} | |
void ddexphelp2(mpf_t r, mpf_t x) { | |
mpf_t xln2; | |
mpf_init(xln2); | |
mpf_mul(xln2, x, LN2); | |
ddexphelp(r, xln2); | |
mpf_mul(r, r, LN2); | |
} | |
int main() { | |
int i, n = NUM_POLY; | |
mpf_t xcof[n+1]; | |
mpf_t acof[n]; | |
mpf_t a, b, eps, tmp; | |
mpf_set_default_prec(BITS); | |
mpf_init(PI_2); | |
mpf_set_str(PI_2, "1.57079632679489661923132169163975144209858469968755291048747229615390820314310", 10); | |
mpf_init(LN2); | |
{ | |
mpf_t a; | |
mpf_init(a); | |
mpf_set_ui(LN2, 0); | |
mpf_set_ui(a, 1); | |
for (int n = 1; n < BITS; n++) { | |
mpf_div_ui(a, a, 2); | |
mpf_add(LN2, LN2, a); | |
mpf_mul_ui(a, a, n); | |
mpf_div_ui(a, a, n+1); | |
} | |
} | |
mpf_get_str(number, &expt, 10, DIGITS, LN2); | |
printf("LN2: %*sE%2d\n", DIGITS+2, number, (int)expt); | |
for (i = 0; i <= n; i++) | |
mpf_init(xcof[i]); | |
for (i = 0; i < n; i++) | |
mpf_init(acof[i]); | |
mpf_init(eps); | |
mpf_init(a); | |
mpf_init(b); | |
mpf_init(tmp); | |
#if 0 | |
mpf_set_d(a, -.01); | |
mpf_set_d(b, .34657359); | |
#endif | |
mpf_set_d(a, 0.0001); | |
mpf_set_d(b, 1); | |
mpf_set_d(eps, 1E-16); | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_mul_ui(xcof[i], b, i+1); | |
mpf_mul_ui(tmp, a, n-i); | |
mpf_add(xcof[i], xcof[i], tmp); | |
mpf_div_ui(xcof[i], xcof[i], n+1); | |
} | |
/* mpf_set_d(xcof[0], .10); */ | |
/* mpf_set_d(xcof[1], .2); */ | |
/* mpf_set_d(xcof[2], .3); */ | |
/* mpf_set(xcof[3], b); */ | |
{ | |
/* | |
void f(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
loghelp(r,x); | |
mpf_mul(tmp, x, x); | |
mpf_mul(tmp, tmp, tmp); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul_ui(tmp, tmp, 5); | |
mpf_div_ui(tmp, tmp, 16); | |
mpf_sub(r, r, tmp); | |
} | |
void df(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
dloghelp(r,x); | |
mpf_mul(tmp, x, x); | |
mpf_mul(tmp, tmp, tmp); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul_ui(tmp, tmp, 7*5); | |
mpf_div_ui(tmp, tmp, 16); | |
mpf_sub(r, r, tmp); | |
} | |
void ddf(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
ddloghelp(r,x); | |
mpf_mul(tmp, x, x); | |
mpf_mul(tmp, tmp, tmp); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul_ui(tmp, tmp, 7*6*5); | |
mpf_div_ui(tmp, tmp, 16); | |
mpf_sub(r, r,tmp); | |
} | |
*/ | |
void f(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
mpf_mul(tmp, x, LN2); | |
mpf_exp(r, tmp); | |
mpf_sub_ui(r, r, 1); | |
} | |
void df(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
mpf_mul(tmp, x, LN2); | |
mpf_exp(r, tmp); | |
mpf_mul(r, r, LN2); | |
} | |
void ddf(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
mpf_mul(tmp, x, LN2); | |
mpf_exp(r, tmp); | |
mpf_mul(r, r, LN2); | |
mpf_mul(r, r, LN2); | |
} | |
remez(n, a, b, eps, | |
f, df, ddf, poly, dpoly, ddpoly, xcof, acof); | |
} | |
return 0; | |
} | |
#include <math.h> | |
#include <gmp.h> | |
#include <stdio.h> | |
#include <stdlib.h> | |
#define DEBUG | |
/*#define MOREDEBUG*/ | |
#define BITS 1000 | |
#define DIGITS 30 //(BITS/3 - 5) | |
#define ROUND_COEFFS 48 | |
#define RELATIVE_ERROR 0 | |
#define ERROR_AT_LAST_IS_ZERO 1 | |
char number[DIGITS+3]; | |
mp_exp_t expt; | |
mpf_t PI_2; | |
mpf_t LN2; | |
#define FIRST_POLY 1 | |
#define STEP_POLY 1 | |
#define NUM_POLY 3 | |
void newton(mpf_t m, void (*f) (mpf_t, mpf_t), void (*df) (mpf_t, mpf_t)) | |
{ | |
mpf_t tmp1, tmp2; | |
int retry; | |
mpf_init(tmp1); | |
mpf_init(tmp2); | |
for (retry = 0; retry < 4; retry++) | |
{ | |
f(tmp1, m); | |
df(tmp2, m); | |
mpf_div(tmp1, tmp1, tmp2); | |
mpf_sub(m, m, tmp1); | |
} | |
} | |
void findzero(mpf_t m, mpf_t a, mpf_t b, void (*f) (mpf_t, mpf_t)) | |
{ | |
int retry; | |
int sign; | |
mpf_t val; | |
mpf_init(val); | |
f(val, b); | |
sign = (mpf_cmp_si(val, 0) > 0); | |
for (retry = 0; retry < BITS/10; retry++) | |
{ | |
mpf_add(m, a, b); | |
mpf_div_ui(m, m, 2); | |
f(val, m); | |
if ((mpf_cmp_si(val, 0) > 0) ^ sign) | |
a = m; | |
else | |
b = m; | |
} | |
} | |
void find_all_zeros(mpf_t *roots, int n, | |
void (*f) (mpf_t, mpf_t), void (*df) (mpf_t, mpf_t), | |
mpf_t a, mpf_t b, mpf_t eps) | |
{ | |
mpf_t tmp1, tmp2, tmp3; | |
int retry, i, j, done; | |
mpf_init(tmp1); | |
mpf_init(tmp2); | |
mpf_init(tmp3); | |
done = 0; | |
for (retry = 0; retry < 50 && !done; retry++) | |
{ | |
done = 1; | |
for (i = 0; i < n; i++) | |
{ | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, roots[i]); | |
printf("%d: root %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
f(tmp1, roots[i]); | |
df(tmp2, roots[i]); | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, tmp1); | |
printf("%d: f %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
mpf_get_str(number, &expt, 10, DIGITS, tmp2); | |
printf("%d: df %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
mpf_div(tmp2, tmp2, tmp1); | |
#if 0 | |
mpf_ui_div(tmp3, 1, tmp2); | |
mpf_get_str(number, &expt, 10, DIGITS, tmp3); | |
printf("%d: delta %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
for (j = 0; j < n; j++) | |
{ | |
if (i != j) | |
{ | |
mpf_sub(tmp3, roots[i], roots[j]); | |
if (mpf_cmp_si(tmp3, 0) == 0) | |
{ | |
mpf_set(tmp2, eps); | |
break; | |
} | |
mpf_ui_div(tmp3, 1, tmp3); | |
mpf_sub(tmp2, tmp2, tmp3); | |
#if STEP_POLY == 2 | |
mpf_add(tmp3, roots[i], roots[j]); | |
mpf_ui_div(tmp3, 1, tmp3); | |
mpf_sub(tmp2, tmp2, tmp3); | |
#endif | |
} | |
} | |
#if FIRST_POLY > 1 | |
mpf_ui_div(tmp3, FIRST_POLY-1, roots[i]); | |
mpf_sub(tmp2, tmp2, tmp3); | |
#endif | |
mpf_ui_div(tmp2, 1, tmp2); | |
mpf_abs(tmp3, tmp2); | |
if (mpf_cmp(tmp3, eps) > 0) | |
done = 0; | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, roots[i]); | |
printf("%d: root %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
mpf_sub(roots[i], roots[i], tmp2); | |
#if 0 | |
mpf_get_str(number, &expt, 10, DIGITS, roots[i]); | |
printf("%d: root %*sE%2d\n", i, DIGITS+2, number, (int)expt); | |
#endif | |
if (mpf_cmp(roots[i], a) < 0) | |
mpf_set(roots[i], a); | |
if (mpf_cmp(roots[i], b) > 0) | |
mpf_set(roots[i], b); | |
} | |
} | |
if (done) | |
printf("Newton succeeds after %d iterations!\n", retry); | |
else | |
printf("Newton didn't succeed!\n"); | |
{ | |
int compare(const mpf_t *a, const mpf_t *b) | |
{ | |
return mpf_cmp(*a, *b); | |
} | |
qsort(roots, n, sizeof(mpf_t), | |
(int (*) (const void*, const void*)) compare); | |
} | |
} | |
void approx(int n, | |
void (*f) (mpf_t, mpf_t), | |
void (*T) (mpf_t, int, mpf_t), | |
mpf_t *xcof, | |
mpf_t *a, mpf_t error) | |
{ | |
int i, j, k; | |
int swapR[n+1]; | |
mpf_t M[n+1][n+1]; | |
mpf_t b[n+1]; | |
mpf_t fac, tmp, tmp2; | |
#ifdef MOREDEBUG | |
char number[DIGITS+3]; | |
mp_exp_t expt; | |
#endif | |
mpf_init(fac); | |
mpf_init(tmp); | |
mpf_init(tmp2); | |
/* build an equation system M*a = b */ | |
/* Here a_0 is error, a_1...,a_{n+1} are the coefficients for poly */ | |
/* The equation system M_ij is SUM a_j*T_j(xi) +/- a_{n+1} = f(xi) */ | |
/* initialize b to the f(xi) */ | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_init(b[i]); | |
f(b[i], xcof[i]); | |
} | |
for (j = 0; j < n; j++) | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_init(M[i][n-j]); | |
T(M[i][n-j], j, xcof[i]); | |
} | |
j = 1; | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_t factor; | |
mpf_init(factor); | |
f(factor, xcof[i]); | |
mpf_add_ui(factor, factor, 1); | |
mpf_init_set_si(M[i][0], j); | |
mpf_mul(M[i][0], M[i][0], factor); | |
j = -j; | |
} | |
mpf_t t; | |
mpf_set_d(t, 1/65536.); | |
mpf_add(b[1], b[1], t); | |
//mpf_set_d(t, 1/65536.); | |
//mpf_add(b[3], b[3], t); | |
#if ERROR_AT_LAST_IS_ZERO | |
mpf_init_set_si(M[n][0], 0); | |
#endif | |
//mpf_init_set_si(M[1][0], 0); | |
for (i = 0; i <= n; i++) | |
swapR[i] = i; | |
#ifdef MOREDEBUG | |
printf("INITIAL EQUATION: \n"); | |
for (k = 0; k <=n; k++) | |
{ | |
for (j = 0; j <= n; j++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, M[swapR[k]][j]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
mpf_get_str(number, &expt, 10, DIGITS, b[swapR[k]]); | |
printf("= %*sE%2d\n", DIGITS+2, number, (int)expt); | |
} | |
#endif | |
/* Produce upper triangle matrix */ | |
/* after this step M[swapR[i]][j] is an upper triangle */ | |
for (i = 0; i < n; i++) | |
{ | |
int best = i; | |
int t; | |
mpf_abs(tmp, M[swapR[i]][i]); | |
for (k = i+1; k <= n; k++) | |
{ | |
mpf_abs(tmp2, M[swapR[k]][i]); | |
if (mpf_cmp(tmp2, tmp) > 0) | |
{ | |
mpf_set(tmp, tmp2); | |
best = k; | |
} | |
} | |
#ifdef MOREDEBUG | |
printf("SELECTING %d (%d)\n", swapR[best], best); | |
#endif | |
t = swapR[i]; | |
swapR[i] = swapR[best]; | |
swapR[best] = t; | |
for (k = i+1; k <= n; k++) | |
{ | |
mpf_div(fac, M[swapR[k]][i], M[swapR[i]][i]); | |
for (j = i; j <= n; j++) | |
{ | |
mpf_mul(tmp, fac, M[swapR[i]][j]); | |
mpf_sub(M[swapR[k]][j], M[swapR[k]][j], tmp); | |
} | |
mpf_mul(tmp, fac, b[swapR[i]]); | |
mpf_sub(b[swapR[k]], b[swapR[k]], tmp); | |
} | |
#ifdef MOREDEBUG | |
printf("EQUATION step %d: \n", i); | |
for (k = 0; k <=n; k++) | |
{ | |
for (j = 0; j <= n; j++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, M[swapR[k]][j]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
mpf_get_str(number, &expt, 10, DIGITS, b[swapR[k]]); | |
printf("= %*sE%2d\n", DIGITS+2, number, (int)expt); | |
} | |
#endif | |
} | |
for (i = n; i >= 0; i--) | |
{ | |
for (j = n; j > i; j--) | |
{ | |
mpf_mul(tmp, j==0 ? error : a[n-j], M[swapR[i]][j]); | |
mpf_sub(b[swapR[i]], b[swapR[i]], tmp); | |
} | |
mpf_div(tmp, b[swapR[i]], M[swapR[i]][i]); | |
if (i == 0) | |
mpf_set(error, tmp); | |
else | |
{ | |
#if ROUND_COEFFS > 0 | |
/* round to FixedPoint representable number */ | |
mpf_mul_2exp(tmp, tmp, (ROUND_COEFFS - 16*(n-i)) + 1); | |
mpf_add_ui(tmp, tmp, 1); | |
mpf_div_ui(tmp, tmp, 2); | |
mpf_floor(tmp, tmp); | |
mpf_div_2exp(tmp, tmp, (ROUND_COEFFS - 16*(n-i))); | |
#endif | |
mpf_set(a[n-i], tmp); | |
} | |
} | |
mpf_abs(error, error); | |
} | |
void newX(mpf_t eps, int n, mpf_t a, mpf_t b, | |
void (*f) (mpf_t, mpf_t), | |
void (*df) (mpf_t, mpf_t), | |
void (*ddf) (mpf_t, mpf_t), | |
void (*T) (mpf_t, int, mpf_t), | |
void (*dT) (mpf_t, int, mpf_t), | |
void (*ddT) (mpf_t, int, mpf_t), | |
mpf_t *xcof, mpf_t *acof, mpf_t error) | |
{ | |
int i; | |
mpf_t nxcof[n]; | |
mpf_t tmp, tmp1; | |
void fun(mpf_t r, mpf_t x) | |
{ | |
int k; | |
mpf_t val; | |
mpf_init(val); | |
f(r, x); | |
for (k = 0; k < n; k++) | |
{ | |
T(val, k, x); | |
mpf_mul(val, acof[k], val); | |
mpf_sub(r, r, val); | |
} | |
} | |
void dev(mpf_t r, mpf_t x) | |
{ | |
int k; | |
mpf_t val; | |
mpf_init(val); | |
df(r, x); | |
for (k = 0; k < n; k++) | |
{ | |
dT(val, k, x); | |
mpf_mul(val, acof[k], val); | |
mpf_sub(r, r, val); | |
} | |
} | |
void ddev(mpf_t r, mpf_t x) | |
{ | |
int k; | |
mpf_t val; | |
mpf_init(val); | |
ddf(r, x); | |
for (k = 0; k < n; k++) | |
{ | |
ddT(val, k, x); | |
mpf_mul(val, acof[k], val); | |
mpf_sub(r, r, val); | |
} | |
} | |
mpf_init(tmp); | |
mpf_init(tmp1); | |
#if 1 | |
find_all_zeros(xcof, n, dev, ddev, a, b, eps); | |
#else | |
if (useNewton) | |
{ | |
for (i=0; i<n; i++) | |
{ | |
mpf_init(nxcof[i]); | |
newton(xcof[i], dev, ddev); | |
} | |
} | |
else | |
{ | |
mpf_t del; | |
mpf_t last; | |
mpf_t tmp2; | |
int sign, p; | |
mpf_init(del); | |
mpf_init(tmp2); | |
mpf_init(last); | |
mpf_sub(del, b, a); | |
mpf_div_ui(del, del, 1000); | |
mpf_mul_ui(tmp, del, 10); | |
dev(tmp1, tmp); | |
sign = (mpf_cmp_si(tmp1, 0) > 0); | |
p = 0; | |
/* mpf_set(xcof[p], a); */ | |
for (i = 10; i < 1000; i++) | |
{ | |
/* char number[DIGITS+3]; */ | |
/* mp_exp_t expt; */ | |
mpf_set(last, tmp); | |
mpf_add(tmp, tmp, del); | |
/* mpf_get_str(number, &expt, 10, DIGITS, tmp); */ | |
/* printf("d'(%*sE%2d) = ", DIGITS+2, number, (int)expt); */ | |
dev(tmp1, tmp); | |
/* mpf_get_str(number, &expt, 10, DIGITS, tmp1); */ | |
/* printf("%*sE%2d\n", DIGITS+2, number, (int)expt); */ | |
if ((mpf_cmp_si(tmp1, 0) > 0) ^ sign) { | |
/* int signtmp; */ | |
findzero(tmp1, last, tmp, dev); | |
/* fun(tmp2, tmp1); */ | |
/* signtmp = mpf_cmp_si(tmp2, 0) > 0; */ | |
/* fun(tmp2, xcof[p]); */ | |
/* if (signtmp != (mpf_cmp_si(tmp2, 0) > 0)) */ | |
/* p++; */ | |
mpf_set(xcof[p++], tmp1); | |
sign ^=1; | |
} | |
} | |
while (p < n) | |
mpf_set(xcof[p++], b); | |
} | |
#endif | |
mpf_set_si(error, 0); | |
for (i=0; i< n; i++) | |
{ | |
fun(tmp, xcof[i]); | |
mpf_abs(tmp, tmp); | |
if (mpf_cmp(tmp, error) > 0) | |
mpf_set(error, tmp); | |
} | |
} | |
void remez(int n, mpf_t a, mpf_t b, mpf_t eps, | |
void (*f) (mpf_t, mpf_t), | |
void (*df) (mpf_t, mpf_t), | |
void (*ddf) (mpf_t, mpf_t), | |
void (*T) (mpf_t, int, mpf_t), | |
void (*dT) (mpf_t, int, mpf_t), | |
void (*ddT) (mpf_t, int, mpf_t), | |
mpf_t *xcof, mpf_t *acof) | |
{ | |
int i, j; | |
mpf_t tmp; | |
mpf_t error, maxerror; | |
mpf_init(tmp); | |
mpf_init(error); | |
mpf_init(maxerror); | |
#ifdef DEBUG | |
printf(" XCOF:"); | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, xcof[i]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
printf("\n"); | |
#endif | |
for (j = 0; j < 100; j++) | |
{ | |
#ifdef DEBUG | |
printf("REMEZ iter %d: \n", j); | |
#endif | |
approx(n, f, T, xcof, acof, error); | |
#ifdef DEBUG | |
printf(" ACOF:"); | |
for (i = 0; i < n; i++) | |
{ | |
#if ROUND_COEFFS | |
mpf_mul_2exp(tmp, acof[i], ROUND_COEFFS - 16*i); | |
printf(" %ld", mpf_get_si(tmp)); | |
#else | |
mpf_get_str(number, &expt, 10, DIGITS, acof[i]); | |
printf(" %*sE%02d", DIGITS+2, number, (int)expt); | |
#endif | |
} | |
printf("\n"); | |
mpf_get_str(number, &expt, 10, DIGITS, error); | |
printf(" ERR: %*sE%02d\n", DIGITS+2, number, (int)expt); | |
#endif | |
newX(eps, n, a, b, f, df, ddf, T,dT,ddT, xcof, acof, maxerror); | |
#ifdef DEBUG | |
printf(" XCOF:"); | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_get_str(number, &expt, 10, DIGITS, xcof[i]); | |
printf(" %*sE%2d", DIGITS+2, number, (int)expt); | |
} | |
printf("\n"); | |
mpf_get_str(number, &expt, 10, DIGITS, maxerror); | |
printf("MAXERR: %*sE%2d\n", DIGITS+2, number, (int) expt); | |
#endif | |
mpf_sub(tmp, maxerror, error); | |
mpf_div(tmp, tmp, maxerror); | |
if (mpf_cmp(tmp, eps) < 0) | |
break; | |
} | |
} | |
void poly(mpf_t r, int i, mpf_t x) { | |
i = STEP_POLY*i + FIRST_POLY; | |
mpf_set_ui(r, 1); | |
while (i-- > 0) | |
mpf_mul(r, r, x); | |
} | |
void dpoly(mpf_t r, int i, mpf_t x) { | |
i = STEP_POLY*i + FIRST_POLY; | |
mpf_set_ui(r, i); | |
i--; | |
while (i-- > 0) | |
mpf_mul(r, r, x); | |
} | |
void ddpoly(mpf_t r, int i, mpf_t x) { | |
i = STEP_POLY*i + FIRST_POLY; | |
mpf_set_ui(r, i * (i-1)); | |
i -= 2; | |
while (i-- > 0) | |
mpf_mul(r, r, x); | |
} | |
void mpf_sin(mpf_t r, mpf_t x) { | |
mpf_t w, p; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_mul(w, x, x); | |
mpf_set(p, x); | |
mpf_set(r, p); | |
for (i=3; i< 40; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_div_ui(p, p,i*(i-1)); | |
mpf_neg(p, p); | |
mpf_add(r, r, p); | |
} | |
} | |
void loghelp(mpf_t r, mpf_t x) { | |
/* compute ln (1+x) - ln (1-x) - 2x */ | |
mpf_t w, p,q; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_init(q); | |
mpf_mul(w, x, x); | |
mpf_mul_ui(p, x,2); | |
mpf_set_ui(r, 0); | |
mpf_add(r, r, p); | |
for (i=3; i< 800; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_div_ui(q, p, i); | |
mpf_add(r, r, q); | |
} | |
} | |
void mpf_log(mpf_t r, mpf_t x) { | |
mpf_t p, q; | |
int i; | |
mpf_init(p); | |
mpf_init(q); | |
mpf_add_ui(p, x, 1); | |
mpf_sub_ui(q, x, 1); | |
mpf_div(p, p, q); | |
mpf_get_str(number, &expt, 10, DIGITS, p); | |
printf("p: %*sE%2d\n", DIGITS+2, number, (int)expt); | |
// p = (x+1)/(x-1) | |
loghelp(r, p); | |
mpf_get_str(number, &expt, 10, DIGITS, r); | |
printf("r: %*sE%2d\n", DIGITS+2, number, (int)expt); | |
// r = ln((1+p) / (1-p)) - 2p | |
// r = ln(x) - 2p | |
mpf_mul_ui(p, p, 2); | |
mpf_add(r, r, p); | |
} | |
void dloghelp(mpf_t r, mpf_t x) { | |
mpf_mul(r, x, x); | |
mpf_ui_sub(r, 1, r); | |
mpf_ui_div(r, 2, r); | |
//mpf_sub_ui(r, r, 2); | |
} | |
void ddloghelp(mpf_t r, mpf_t x) { | |
mpf_mul(r, x, x); | |
mpf_ui_sub(r, 1, r); | |
mpf_mul(r, r, r); | |
mpf_div(r, x, r); | |
mpf_mul_ui(r, r, 4); | |
} | |
void mpf_cos(mpf_t r, mpf_t x) { | |
mpf_t w, p; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_mul(w, x, x); | |
mpf_set_si(p, 1); | |
mpf_set(r, p); | |
for (i=2; i< 40; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_div_ui(p, p, i*(i-1)); | |
mpf_neg(p,p); | |
mpf_add(r, r, p); | |
} | |
} | |
void mpf_exp(mpf_t r, mpf_t x) { | |
mpf_t p; | |
int i; | |
mpf_init(p); | |
mpf_set_si(p, 1); | |
mpf_set(r, p); | |
for (i=1; i< 1000; i++) | |
{ | |
mpf_mul(p, p, x); | |
mpf_div_ui(p, p, i); | |
mpf_add(r, r, p); | |
} | |
} | |
void mpf_atan(mpf_t r, mpf_t x) { | |
mpf_t w, p, tmp; | |
int i; | |
mpf_init(w); | |
mpf_init(p); | |
mpf_init(tmp); | |
mpf_mul(w, x, x); | |
mpf_set(p, x); | |
mpf_set(r, p); | |
for (i=3; i< 80; i+=2) | |
{ | |
mpf_mul(p, p, w); | |
mpf_neg(p,p); | |
mpf_div_ui(tmp, p, i); | |
mpf_add(r, r, tmp); | |
} | |
} | |
void c(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_div_ui(w, w, 1048576); | |
mpf_mul_ui(w, w, 0); | |
/* mpf_div_ui(w, w, 1024); */ | |
mpf_mul(r, x, PI_2); | |
mpf_cos(r, r); | |
mpf_sub_ui(r, r, 1); | |
mpf_sub(r, r, w); | |
} | |
void dc(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, x); | |
mpf_div_ui(w, w, 1048576); | |
mpf_mul_ui(w, w, (0)*8); | |
mpf_mul(r, x, PI_2); | |
mpf_sin(r, r); | |
mpf_neg(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void ddc(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, x); | |
mpf_mul(w, w, w); | |
mpf_div_ui(w, w, 1048576); | |
mpf_mul_ui(w, w, (0)*8*7); | |
mpf_mul(r, x, PI_2); | |
mpf_cos(r, r); | |
mpf_neg(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void s(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, x); | |
mpf_div_ui(w, w, 8192*4); | |
mpf_mul_ui(w, w, 5); | |
mpf_mul(r, x, PI_2); | |
mpf_sin(r, r); | |
mpf_sub(r, r, w); | |
mpf_set_d(w, 4.6566128730773e-10/2); | |
mpf_sub(r, r, w); | |
} | |
void ds(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_div_ui(w, w, 8192*4); | |
mpf_mul_ui(w, w, 9*5); | |
mpf_mul(r, x, PI_2); | |
mpf_cos(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void dds(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, w); | |
mpf_div(w, w, x); | |
mpf_div_ui(w, w, 8192*4); | |
mpf_mul_ui(w, w, 9*8*5); | |
mpf_mul(r, x, PI_2); | |
mpf_sin(r, r); | |
mpf_neg(r, r); | |
mpf_mul(r, r, PI_2); | |
mpf_mul(r, r, PI_2); | |
mpf_sub(r, r, w); | |
} | |
void at(mpf_t r, mpf_t x) { | |
mpf_atan(r, x); | |
mpf_div(r, r, PI_2); | |
} | |
void dat(mpf_t r, mpf_t x) { | |
mpf_mul(r, x, x); | |
mpf_add_ui(r, r, 1); | |
mpf_mul(r, r, PI_2); | |
mpf_ui_div(r, 1, r); | |
} | |
void ddat(mpf_t r, mpf_t x) { | |
mpf_t w; | |
mpf_init(w); | |
mpf_mul(w, x, x); | |
mpf_add_ui(w, w, 1); | |
mpf_mul(w, w, w); | |
mpf_mul(w, w, PI_2); | |
mpf_div(r, x, w); | |
mpf_mul_ui(r, r, 2); | |
mpf_neg(r, r); | |
} | |
void sn(mpf_t r, mpf_t x) { | |
mpf_sin(r, x); | |
} | |
void dsn(mpf_t r, mpf_t x) { | |
mpf_cos(r, x); | |
} | |
void ddsn(mpf_t r, mpf_t x) { | |
mpf_sin(r, x); | |
mpf_neg(r, r); | |
} | |
void exphelp(mpf_t r, mpf_t x) { | |
/* x*(e+1)/(e-1) - 2 */ | |
mpf_t e, t, xln2; | |
mpf_init(e); | |
mpf_init(t); | |
mpf_exp(e, x); | |
mpf_sub_ui(t, e, 1); | |
mpf_add_ui(e, e, 1); | |
mpf_mul(e, e, x); | |
mpf_div(t, e, t); | |
mpf_sub_ui(r, t, 2); | |
} | |
void dexphelp(mpf_t r, mpf_t x) { | |
/* (ee - 2xe - 1)/(e-1)^2 */ | |
if (mpf_sgn(x) == 0) { | |
mpf_set_ui(r, 0); | |
return; | |
} | |
mpf_t e, n, xln2; | |
mpf_init(e); | |
mpf_init(n); | |
mpf_exp(e, x); | |
mpf_sub_ui(n, e, 1); | |
mpf_mul(n, n, n); | |
mpf_mul_ui(r, x, 2); | |
mpf_sub(r, e, r); | |
mpf_mul(r, e, r); | |
mpf_sub_ui(r, r, 1); | |
mpf_div(r, r, n); | |
} | |
void ddexphelp(mpf_t r, mpf_t x) { | |
/* 2e (-2e + xe + 2 + x) / (e-1)^3 */ | |
if (mpf_sgn(x) == 0) { | |
mpf_set_ui(r, 2); | |
mpf_div_ui(r, r, 3); | |
return; | |
} | |
mpf_t e, n; | |
mpf_t xln2; | |
mpf_init(e); | |
mpf_init(n); | |
mpf_exp(e, x); | |
mpf_sub_ui(n, e, 1); | |
mpf_mul(r, n, n); | |
mpf_mul(n, r, n); | |
// n = (e-1)^3 | |
mpf_mul(r, e, x); | |
mpf_add_ui(r, r, 2); | |
// r = xe + 2 | |
mpf_add(r, r, x); | |
// r = xe + 2 + x | |
mpf_add(e, e, e); | |
// e = 2e | |
mpf_sub(r, r, e); | |
mpf_mul(r, r, e); | |
mpf_div(r, r, n); | |
} | |
void exphelp2(mpf_t r, mpf_t x) { | |
mpf_t xln2; | |
mpf_init(xln2); | |
mpf_mul(xln2, x, LN2); | |
exphelp(r, xln2); | |
mpf_div(r, r, LN2); | |
} | |
void dexphelp2(mpf_t r, mpf_t x) { | |
mpf_t xln2; | |
mpf_init(xln2); | |
mpf_mul(xln2, x, LN2); | |
dexphelp(r, xln2); | |
} | |
void ddexphelp2(mpf_t r, mpf_t x) { | |
mpf_t xln2; | |
mpf_init(xln2); | |
mpf_mul(xln2, x, LN2); | |
ddexphelp(r, xln2); | |
mpf_mul(r, r, LN2); | |
} | |
int main() { | |
int i, n = NUM_POLY; | |
mpf_t xcof[n+1]; | |
mpf_t acof[n]; | |
mpf_t a, b, eps, tmp; | |
mpf_set_default_prec(BITS); | |
mpf_init(PI_2); | |
mpf_set_str(PI_2, "1.57079632679489661923132169163975144209858469968755291048747229615390820314310", 10); | |
mpf_init(LN2); | |
{ | |
mpf_t a; | |
mpf_init(a); | |
mpf_set_ui(LN2, 0); | |
mpf_set_ui(a, 1); | |
for (int n = 1; n < BITS; n++) { | |
mpf_div_ui(a, a, 2); | |
mpf_add(LN2, LN2, a); | |
mpf_mul_ui(a, a, n); | |
mpf_div_ui(a, a, n+1); | |
} | |
} | |
mpf_get_str(number, &expt, 10, DIGITS, LN2); | |
printf("LN2: %*sE%2d\n", DIGITS+2, number, (int)expt); | |
for (i = 0; i <= n; i++) | |
mpf_init(xcof[i]); | |
for (i = 0; i < n; i++) | |
mpf_init(acof[i]); | |
mpf_init(eps); | |
mpf_init(a); | |
mpf_init(b); | |
mpf_init(tmp); | |
#if 0 | |
mpf_set_d(a, -.01); | |
mpf_set_d(b, .34657359); | |
#endif | |
mpf_set_d(a, 0.0001); | |
mpf_set_d(b, 1); | |
mpf_set_d(eps, 1E-16); | |
for (i = 0; i <= n; i++) | |
{ | |
mpf_mul_ui(xcof[i], b, i+1); | |
mpf_mul_ui(tmp, a, n-i); | |
mpf_add(xcof[i], xcof[i], tmp); | |
mpf_div_ui(xcof[i], xcof[i], n+1); | |
} | |
/* mpf_set_d(xcof[0], .10); */ | |
/* mpf_set_d(xcof[1], .2); */ | |
/* mpf_set_d(xcof[2], .3); */ | |
/* mpf_set(xcof[3], b); */ | |
{ | |
/* | |
void f(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
loghelp(r,x); | |
mpf_mul(tmp, x, x); | |
mpf_mul(tmp, tmp, tmp); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul_ui(tmp, tmp, 5); | |
mpf_div_ui(tmp, tmp, 16); | |
mpf_sub(r, r, tmp); | |
} | |
void df(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
dloghelp(r,x); | |
mpf_mul(tmp, x, x); | |
mpf_mul(tmp, tmp, tmp); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul_ui(tmp, tmp, 7*5); | |
mpf_div_ui(tmp, tmp, 16); | |
mpf_sub(r, r, tmp); | |
} | |
void ddf(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
ddloghelp(r,x); | |
mpf_mul(tmp, x, x); | |
mpf_mul(tmp, tmp, tmp); | |
mpf_mul(tmp, tmp, x); | |
mpf_mul_ui(tmp, tmp, 7*6*5); | |
mpf_div_ui(tmp, tmp, 16); | |
mpf_sub(r, r,tmp); | |
} | |
*/ | |
void f(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
mpf_mul(tmp, x, LN2); | |
mpf_exp(r, tmp); | |
mpf_sub_ui(r, r, 1); | |
} | |
void df(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
mpf_mul(tmp, x, LN2); | |
mpf_exp(r, tmp); | |
mpf_mul(r, r, LN2); | |
} | |
void ddf(mpf_t r, mpf_t x) { | |
mpf_t tmp; | |
mpf_init(tmp); | |
mpf_mul(tmp, x, LN2); | |
mpf_exp(r, tmp); | |
mpf_mul(r, r, LN2); | |
mpf_mul(r, r, LN2); | |
} | |
remez(n, a, b, eps, | |
f, df, ddf, poly, dpoly, ddpoly, xcof, acof); | |
} | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Compile with
gcc remez.c -lgmp
. Warning: this code is an undocumented mess and will break if you touch it too much.