Skip to content

Instantly share code, notes, and snippets.

@jhoenicke

jhoenicke/remez.c

Created August 2, 2020 14:03
Show Gist options
  • Save jhoenicke/fca701e9179b3d8d4ab45de6afa2b42f to your computer and use it in GitHub Desktop.
Save jhoenicke/fca701e9179b3d8d4ab45de6afa2b42f to your computer and use it in GitHub Desktop.
Download ZIP
Remez Algorithm to approximate arbitrary functions by low-level polynomials on a fixed interval
Raw
remez.c
#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;
}
@jhoenicke
Copy link
Author

jhoenicke commented Aug 2, 2020

Compile with gcc remez.c -lgmp. Warning: this code is an undocumented mess and will break if you touch it too much.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment