Skip to content

Commit

Permalink
Accuracy improvements for C.LN1+X, and complex LN, ATAN, and ATANH.
Browse files Browse the repository at this point in the history
  • Loading branch information
@thomasokken
thomasokken committed Mar 26, 2024
1 parent e494995 commit 6c95e99
Show file tree
Hide file tree
Showing 7 changed files with 153 additions and 134 deletions.
2 changes: 1 addition & 1 deletion common/core_commands2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@ int docmd_comb(arg_struct *arg) {
#ifdef BCD_MATH
s = x == 0 ? 1 : pow(10, 1 + floor(log10(x)));
#elif ANDROID
s = x == 0 ? 1 : pow(2, 1 + floor(log(x) / log(2)));
s = x == 0 ? 1 : pow(2, 1 + floor(log(x) / log(phloat(2))));
#else
s = x == 0 ? 1 : pow(2, 1 + floor(log2(x)));
#endif
Expand Down
24 changes: 7 additions & 17 deletions common/core_commands4.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -393,23 +393,13 @@ static int mappable_ln_1_x_c(phloat xre, phloat xim, phloat *yre, phloat *yim) {
return ERR_NONE;
} else {
phloat x1re = xre + 1;
phloat h = hypot(x1re, xim);
phloat bre;
phloat s;
if (p_isinf(h)) {
s = 10000;
h = hypot(x1re / s, xim / s);
bre = log(h) + log(s);
} else {
s = 1;
bre = log(h);
}
phloat bim = atan2(xim, x1re);
phloat cre = x1re - 1 - xre;
phloat dre = cre * x1re / h / h * s;
phloat dim = -cre * xim / h / h * s;
*yre = bre - dre;
*yim = bim - dim;
phloat are, aim;
math_ln(x1re, xim, &are, &aim);
phloat bre = x1re - 1 - xre;
phloat cre, cim;
math_inv(x1re, xim, &cre, &cim);
*yre = are - bre * cre;
*yim = aim - bre * cim;
return ERR_NONE;
}
}
Expand Down
107 changes: 12 additions & 95 deletions common/core_commands6.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,27 +392,23 @@ static int mappable_log_c(phloat xre, phloat xim, phloat *yre, phloat *yim) {
*yim = 0;
} else {
*yre = log10(-xre);
*yim = PI / log(10.0);
*yim = PI / log(phloat(10.0));
}
return ERR_NONE;
} else if (xre == 0) {
if (xim > 0) {
*yre = log10(xim);
*yim = PI / log(100.0);
*yim = PI / log(phloat(100.0));
} else {
*yre = log10(-xim);
*yim = -PI / log(100.0);
*yim = -PI / log(phloat(100.0));
}
return ERR_NONE;
} else {
phloat h = hypot(xre, xim);
if (p_isinf(h)) {
const phloat s = 10000;
h = hypot(xre / s, xim / s);
*yre = log10(h) + 4;
} else
*yre = log10(h);
*yim = atan2(xim, xre) / log(10.0);
math_ln(xre, xim, yre, yim);
phloat s = log(phloat(10.0));
*yre /= s;
*yim /= s;
return ERR_NONE;
}
}
Expand All @@ -426,7 +422,7 @@ int docmd_log(arg_struct *arg) {
if (flags.f.real_result_only)
return ERR_INVALID_DATA;
else {
vartype *r = new_complex(log10(-x->x), PI / log(10.0));
vartype *r = new_complex(log10(-x->x), PI / log(phloat(10.0)));
if (r == NULL)
return ERR_INSUFFICIENT_MEMORY;
else {
Expand Down Expand Up @@ -465,7 +461,7 @@ static int mappable_10_pow_x_r(phloat x, phloat *y) {
static int mappable_10_pow_x_c(phloat xre, phloat xim, phloat *yre, phloat *yim){
int inf;
phloat h;
xim *= log(10.0);
xim *= log(phloat(10.0));
if ((inf = p_isinf(xim)) != 0)
xim = inf < 0 ? NEG_HUGE_PHLOAT : POS_HUGE_PHLOAT;
h = pow(10, xre);
Expand Down Expand Up @@ -511,40 +507,6 @@ static int mappable_ln_r(phloat x, phloat *y) {
}
}

static int mappable_ln_c(phloat xre, phloat xim, phloat *yre, phloat *yim) {
if (xim == 0) {
if (xre == 0)
return ERR_INVALID_DATA;
if (xre > 0) {
*yre = log(xre);
*yim = 0;
} else {
*yre = log(-xre);
*yim = PI;
}
return ERR_NONE;
} else if (xre == 0) {
if (xim > 0) {
*yre = log(xim);
*yim = PI / 2;
} else {
*yre = log(-xim);
*yim = -PI / 2;
}
return ERR_NONE;
} else {
phloat h = hypot(xre, xim);
if (p_isinf(h)) {
const phloat s = 10000;
h = hypot(xre / s, xim / s);
*yre = log(h) + log(s);
} else
*yre = log(h);
*yim = atan2(xim, xre);
return ERR_NONE;
}
}

int docmd_ln(arg_struct *arg) {
if (stack[sp]->type == TYPE_REAL) {
vartype_real *x = (vartype_real *) stack[sp];
Expand Down Expand Up @@ -573,7 +535,7 @@ int docmd_ln(arg_struct *arg) {
}
} else {
vartype *v;
int err = map_unary(stack[sp], &v, mappable_ln_r, mappable_ln_c);
int err = map_unary(stack[sp], &v, mappable_ln_r, math_ln);
if (err == ERR_NONE)
unary_result(v);
return err;
Expand Down Expand Up @@ -725,51 +687,6 @@ static int mappable_inv_r(phloat x, phloat *y) {
return ERR_NONE;
}

static int mappable_inv_c(phloat xre, phloat xim, phloat *yre, phloat *yim) {
phloat r, t, rre, rim;
int inf;
if (xre == 0 && xim == 0)
return ERR_DIVIDE_BY_0;
if (fabs(xim) <= fabs(xre)) {
r = xim / xre;
t = 1 / (xre + xim * r);
if (r == 0) {
rre = t;
rim = -xim * (1 / xre) * t;
} else {
rre = t;
rim = -r * t;
}
} else {
r = xre / xim;
t = 1 / (xre * r + xim);
if (r == 0) {
rre = xre * (1 / xim) * t;
rim = -t;
} else {
rre = r * t;
rim = -t;
}
}
inf = p_isinf(rre);
if (inf != 0) {
if (flags.f.range_error_ignore)
rre = inf == 1 ? POS_HUGE_PHLOAT : NEG_HUGE_PHLOAT;
else
return ERR_OUT_OF_RANGE;
}
inf = p_isinf(rim);
if (inf != 0) {
if (flags.f.range_error_ignore)
rim = inf == 1 ? POS_HUGE_PHLOAT : NEG_HUGE_PHLOAT;
else
return ERR_OUT_OF_RANGE;
}
*yre = rre;
*yim = rim;
return ERR_NONE;
}

int docmd_inv(arg_struct *arg) {
vartype *x = stack[sp];
vartype *v;
Expand All @@ -781,7 +698,7 @@ int docmd_inv(arg_struct *arg) {
err = unit_div(x, one, &v);
free_vartype(one);
} else
err = map_unary(stack[sp], &v, mappable_inv_r, mappable_inv_c);
err = map_unary(stack[sp], &v, mappable_inv_r, math_inv);
if (err == ERR_NONE)
unary_result(v);
return err;
Expand Down Expand Up @@ -982,7 +899,7 @@ int docmd_y_pow_x(arg_struct *arg) {
res = new_complex(0, 0);
goto done;
}
err = mappable_ln_c(yre, yim, &lre, &lim);
err = math_ln(yre, yim, &lre, &lim);
if (err != ERR_NONE)
return err;
tmp = lre * xre - lim * xim;
Expand Down
144 changes: 123 additions & 21 deletions common/core_math2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,34 +243,91 @@ int math_atanh(phloat xre, phloat xim, phloat *yre, phloat *yim) {
return ERR_NONE;
}

phloat are, aim, bre, bim, cre, cim, h;
phloat x = fabs(xim);
phloat y = fabs(xre);

/* TODO: review, and deal with overflows in intermediate results */
#ifdef BCD_MATH
const Phloat BIG(1000000000000000000LL);
#else
const double BIG = 0x8000000; // 2^27
#endif

/* a = 1 + x */
are = 1 + xre;
aim = xim;

/* b = 1 - x */
bre = 1 - xre;
bim = - xim;

/* c = a / b */
h = hypot(bre, bim);
cre = (are * bre + aim * bim) / h / h;
cim = (aim * bre - are * bim) / h / h;
if (x >= BIG || y >= BIG) { // atan(1/z) ~= 1/z
math_inv(x, y, &x, &y); // atan(x+y*I) ~= pi/2 + t
x += PI / 2; // |re(t^3/3)/re(t)| <= 0x1p-54
// |im(t^3/3)/im(t)| <= 0x1p-54
} else {
phloat x2 = x * x;
phloat ym = 1 - y;
x = atan2(2 * x, (1 + y) * ym - x2) / 2;
y = log1p(4 * y / (ym * ym + x2)) / 4;
}

/* y = log(c) / 2 */
*yre = log(hypot(cre, cim)) / 2;
*yim = atan2(cim, cre) / 2;
x = fabs(x);
if (xim < 0)
x = -x;
y = fabs(y);
if (xre < 0)
y = -y;

/* Theoretically, you could go out of range, but in practice,
* you can't get close enough to the critical values to cause
* trouble.
*/
*yre = y;
*yim = x;
return ERR_NONE;
}

int math_ln(phloat xre, phloat xim, phloat *yre, phloat *yim) {
if (xim == 0) {
if (xre == 0)
return ERR_INVALID_DATA;
if (xre > 0) {
*yre = log(xre);
*yim = 0;
} else {
*yre = log(-xre);
*yim = PI;
}
return ERR_NONE;
} else if (xre == 0) {
if (xim > 0) {
*yre = log(xim);
*yim = PI / 2;
} else {
*yre = log(-xim);
*yim = -PI / 2;
}
return ERR_NONE;
} else {
phloat h = xre * xre + xim * xim;
phloat a = atan2(xim, xre);
if (h > 0.5 && h < 3.0) {
if ((h = fabs(xre)) < (xim = fabs(xim))) {
h = xim;
xim = xre;
}
h -= 1; /* x^2-1 = 2*(x-1) + (x-1)^2 */
h = log1p(2 * h + h * h + xim * xim) / 2;
} else if (p_isnormal(h)) {
h = log(h) / 2;
} else {
#ifdef BCD_MATH
phloat m = scalbn(phloat(1), 38);
phloat b = -log(phloat(10)) * 38;
#else
phloat m = 0x2000000000000000; // 2^61
phloat b = -log(phloat(2)) * 61;
#endif
if (p_isinf(h)) {
m = 1 / m;
b = -b;
}
h = log(hypot(m * xre, m * xim)) + b;
}
*yre = h;
*yim = a;
return ERR_NONE;
}
}

int math_sqrt(phloat xre, phloat xim, phloat *yre, phloat *yim) {
if (xre == 0) {
if (xim == 0) {
Expand Down Expand Up @@ -326,3 +383,48 @@ int math_sqrt(phloat xre, phloat xim, phloat *yre, phloat *yim) {
}
return ERR_NONE;
}

int math_inv(phloat xre, phloat xim, phloat *yre, phloat *yim) {
phloat r, t, rre, rim;
int inf;
if (xre == 0 && xim == 0)
return ERR_DIVIDE_BY_0;
if (fabs(xim) <= fabs(xre)) {
r = xim / xre;
t = 1 / (xre + xim * r);
if (r == 0) {
rre = t;
rim = -xim * (1 / xre) * t;
} else {
rre = t;
rim = -r * t;
}
} else {
r = xre / xim;
t = 1 / (xre * r + xim);
if (r == 0) {
rre = xre * (1 / xim) * t;
rim = -t;
} else {
rre = r * t;
rim = -t;
}
}
inf = p_isinf(rre);
if (inf != 0) {
if (flags.f.range_error_ignore)
rre = inf == 1 ? POS_HUGE_PHLOAT : NEG_HUGE_PHLOAT;
else
return ERR_OUT_OF_RANGE;
}
inf = p_isinf(rim);
if (inf != 0) {
if (flags.f.range_error_ignore)
rim = inf == 1 ? POS_HUGE_PHLOAT : NEG_HUGE_PHLOAT;
else
return ERR_OUT_OF_RANGE;
}
*yre = rre;
*yim = rim;
return ERR_NONE;
}
2 changes: 2 additions & 0 deletions common/core_math2.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ int math_tan(phloat x, phloat *y, bool rad);
int math_asinh(phloat xre, phloat xim, phloat *yre, phloat *yim);
int math_acosh(phloat xre, phloat xim, phloat *yre, phloat *yim);
int math_atanh(phloat xre, phloat xim, phloat *yre, phloat *yim);
int math_ln(phloat xre, phloat xim, phloat *yre, phloat *yim);
int math_sqrt(phloat xre, phloat xim, phloat *yre, phloat *yim);
int math_inv(phloat xre, phloat xim, phloat *yre, phloat *yim);

#endif
6 changes: 6 additions & 0 deletions common/core_phloat.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -405,6 +405,12 @@ int p_isnan(Phloat p) {
return r;
}

int p_isnormal(Phloat p) {
int r;
bid128_isNormal(&r, &p.val);
return r;
}

int to_digit(Phloat p) {
BID_UINT128 ten, res;
int d10 = 10;
Expand Down

0 comments on commit 6c95e99

Please sign in to comment.