msvcp: fixed complex division
Andrey Zhezherun
zhezherun at yandex.ru
Sun Nov 3 08:58:29 CST 2013
On 1 November 2013 10:40, Piotr Caban <piotr.caban at gmail.com> wrote:
> How about making it more numerically safe?
Makes sense, although I checked what libstdc++ does and it has the
vanilla version, similar to what I had before. Not sure what the
native msvcp does though. A revised patch is attached.
Regards,
Andrey
-------------- next part --------------
--- wine-1.7.5.orig/dlls/msvcp71/math.c 2013-10-25 18:45:30.000000000 +0100
+++ wine-1.7.5/dlls/msvcp71/math.c 2013-11-01 11:34:42.000000000 +0000
@@ -1087,6 +1087,7 @@
/* ??$?KM at std@@YA?AV?$complex at M@0 at AEBV10@0 at Z */
complex_float* __cdecl complex_float_div(complex_float *ret, const complex_float *l, const complex_float *r)
{
+ float tmp, den;
if((!r->real && !r->imag) || _isnan(l->real) || _isnan(l->imag)
|| _isnan(r->real) || _isnan(r->imag)) {
ret->real = std_numeric_limits_float_quiet_NaN();
@@ -1094,15 +1095,16 @@
return ret;
}
- ret->real = 0;
- ret->imag = 0;
- if(r->real) {
- ret->real = l->real / r->real;
- ret->imag = l->imag / r->real;
- }
- if(r->imag) {
- ret->real += l->imag / r->imag;
- ret->imag -= l->real / r->imag;
+ if (fabsf(r->real) >= fabsf(r->imag)) {
+ tmp = r->imag / r->real;
+ den = r->real + r->imag*tmp;
+ ret->real = (l->real + l->imag*tmp) / den;
+ ret->imag = (l->imag - l->real*tmp) / den;
+ } else {
+ tmp = r->real / r->imag;
+ den = r->real*tmp + r->imag;
+ ret->real = (l->real*tmp + l->imag) / den;
+ ret->imag = (l->imag*tmp - l->real) / den;
}
return ret;
}
@@ -1705,6 +1707,7 @@
/* ??$?KO at std@@YA?AV?$complex at O@0 at AEBV10@0 at Z */
complex_double* __cdecl complex_double_div(complex_double *ret, const complex_double *l, const complex_double *r)
{
+ double tmp, den;
if((!r->real && !r->imag) || _isnan(l->real) || _isnan(l->imag)
|| _isnan(r->real) || _isnan(r->imag)) {
ret->real = std_numeric_limits_double_quiet_NaN();
@@ -1712,15 +1715,16 @@
return ret;
}
- ret->real = 0;
- ret->imag = 0;
- if(r->real) {
- ret->real = l->real / r->real;
- ret->imag = l->imag / r->real;
- }
- if(r->imag) {
- ret->real += l->imag / r->imag;
- ret->imag -= l->real / r->imag;
+ if (fabs(r->real) >= fabs(r->imag)) {
+ tmp = r->imag / r->real;
+ den = r->real + r->imag*tmp;
+ ret->real = (l->real + l->imag*tmp) / den;
+ ret->imag = (l->imag - l->real*tmp) / den;
+ } else {
+ tmp = r->real / r->imag;
+ den = r->real*tmp + r->imag;
+ ret->real = (l->real*tmp + l->imag) / den;
+ ret->imag = (l->imag*tmp - l->real) / den;
}
return ret;
}
--- wine-1.7.5.orig/dlls/msvcp80/math.c 2013-10-25 18:45:30.000000000 +0100
+++ wine-1.7.5/dlls/msvcp80/math.c 2013-11-01 11:34:22.000000000 +0000
@@ -1087,6 +1087,7 @@
/* ??$?KM at std@@YA?AV?$complex at M@0 at AEBV10@0 at Z */
complex_float* __cdecl complex_float_div(complex_float *ret, const complex_float *l, const complex_float *r)
{
+ float tmp, den;
if((!r->real && !r->imag) || _isnan(l->real) || _isnan(l->imag)
|| _isnan(r->real) || _isnan(r->imag)) {
ret->real = std_numeric_limits_float_quiet_NaN();
@@ -1094,15 +1095,16 @@
return ret;
}
- ret->real = 0;
- ret->imag = 0;
- if(r->real) {
- ret->real = l->real / r->real;
- ret->imag = l->imag / r->real;
- }
- if(r->imag) {
- ret->real += l->imag / r->imag;
- ret->imag -= l->real / r->imag;
+ if (fabsf(r->real) >= fabsf(r->imag)) {
+ tmp = r->imag / r->real;
+ den = r->real + r->imag*tmp;
+ ret->real = (l->real + l->imag*tmp) / den;
+ ret->imag = (l->imag - l->real*tmp) / den;
+ } else {
+ tmp = r->real / r->imag;
+ den = r->real*tmp + r->imag;
+ ret->real = (l->real*tmp + l->imag) / den;
+ ret->imag = (l->imag*tmp - l->real) / den;
}
return ret;
}
@@ -1705,6 +1707,7 @@
/* ??$?KO at std@@YA?AV?$complex at O@0 at AEBV10@0 at Z */
complex_double* __cdecl complex_double_div(complex_double *ret, const complex_double *l, const complex_double *r)
{
+ double tmp, den;
if((!r->real && !r->imag) || _isnan(l->real) || _isnan(l->imag)
|| _isnan(r->real) || _isnan(r->imag)) {
ret->real = std_numeric_limits_double_quiet_NaN();
@@ -1712,15 +1715,16 @@
return ret;
}
- ret->real = 0;
- ret->imag = 0;
- if(r->real) {
- ret->real = l->real / r->real;
- ret->imag = l->imag / r->real;
- }
- if(r->imag) {
- ret->real += l->imag / r->imag;
- ret->imag -= l->real / r->imag;
+ if (fabs(r->real) >= fabs(r->imag)) {
+ tmp = r->imag / r->real;
+ den = r->real + r->imag*tmp;
+ ret->real = (l->real + l->imag*tmp) / den;
+ ret->imag = (l->imag - l->real*tmp) / den;
+ } else {
+ tmp = r->real / r->imag;
+ den = r->real*tmp + r->imag;
+ ret->real = (l->real*tmp + l->imag) / den;
+ ret->imag = (l->imag*tmp - l->real) / den;
}
return ret;
}
--- wine-1.7.5.orig/dlls/msvcp90/math.c 2013-10-25 18:45:30.000000000 +0100
+++ wine-1.7.5/dlls/msvcp90/math.c 2013-11-01 11:34:06.000000000 +0000
@@ -1087,6 +1087,7 @@
/* ??$?KM at std@@YA?AV?$complex at M@0 at AEBV10@0 at Z */
complex_float* __cdecl complex_float_div(complex_float *ret, const complex_float *l, const complex_float *r)
{
+ float tmp, den;
if((!r->real && !r->imag) || _isnan(l->real) || _isnan(l->imag)
|| _isnan(r->real) || _isnan(r->imag)) {
ret->real = std_numeric_limits_float_quiet_NaN();
@@ -1094,15 +1095,16 @@
return ret;
}
- ret->real = 0;
- ret->imag = 0;
- if(r->real) {
- ret->real = l->real / r->real;
- ret->imag = l->imag / r->real;
- }
- if(r->imag) {
- ret->real += l->imag / r->imag;
- ret->imag -= l->real / r->imag;
+ if (fabsf(r->real) >= fabsf(r->imag)) {
+ tmp = r->imag / r->real;
+ den = r->real + r->imag*tmp;
+ ret->real = (l->real + l->imag*tmp) / den;
+ ret->imag = (l->imag - l->real*tmp) / den;
+ } else {
+ tmp = r->real / r->imag;
+ den = r->real*tmp + r->imag;
+ ret->real = (l->real*tmp + l->imag) / den;
+ ret->imag = (l->imag*tmp - l->real) / den;
}
return ret;
}
@@ -1705,6 +1707,7 @@
/* ??$?KO at std@@YA?AV?$complex at O@0 at AEBV10@0 at Z */
complex_double* __cdecl complex_double_div(complex_double *ret, const complex_double *l, const complex_double *r)
{
+ double tmp, den;
if((!r->real && !r->imag) || _isnan(l->real) || _isnan(l->imag)
|| _isnan(r->real) || _isnan(r->imag)) {
ret->real = std_numeric_limits_double_quiet_NaN();
@@ -1712,15 +1715,16 @@
return ret;
}
- ret->real = 0;
- ret->imag = 0;
- if(r->real) {
- ret->real = l->real / r->real;
- ret->imag = l->imag / r->real;
- }
- if(r->imag) {
- ret->real += l->imag / r->imag;
- ret->imag -= l->real / r->imag;
+ if (fabs(r->real) >= fabs(r->imag)) {
+ tmp = r->imag / r->real;
+ den = r->real + r->imag*tmp;
+ ret->real = (l->real + l->imag*tmp) / den;
+ ret->imag = (l->imag - l->real*tmp) / den;
+ } else {
+ tmp = r->real / r->imag;
+ den = r->real*tmp + r->imag;
+ ret->real = (l->real*tmp + l->imag) / den;
+ ret->imag = (l->imag*tmp - l->real) / den;
}
return ret;
}
More information about the wine-devel
mailing list