| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 1995, 1996, 1997 Robert Gentleman and Ross Ihaka |
| * Copyright (C) 2000-2016 The R Core Team |
| * Copyright (C) 2005 The R Foundation |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| #ifdef HAVE_CONFIG_H |
| #include <config.h> |
| #endif |
| |
| /* Note: gcc -peantic may warn in several places about C99 features |
| as extensions. |
| This was a very-long-standing GCC bug, http://gcc.gnu.org/PR7263 |
| The system <complex.h> header can work around it: some do. |
| It should have been resolved (after a decade) in 2012. |
| */ |
| |
| #if defined(HAVE_CTANH) && !defined(HAVE_WORKING_CTANH) |
| #undef HAVE_CTANH |
| #endif |
| |
| #if 0 |
| /* For testing substitute fns */ |
| #undef HAVE_CARG |
| #undef HAVE_CABS |
| #undef HAVE_CPOW |
| #undef HAVE_CEXP |
| #undef HAVE_CLOG |
| #undef HAVE_CSQRT |
| #undef HAVE_CSIN |
| #undef HAVE_CCOS |
| #undef HAVE_CTAN |
| #undef HAVE_CASIN |
| #undef HAVE_CACOS |
| #undef HAVE_CATAN |
| #undef HAVE_CSINH |
| #undef HAVE_CCOSH |
| #undef HAVE_CTANH |
| #endif |
| |
| #ifdef __SUNPRO_C |
| /* segfaults in Solaris Studio 12.3 */ |
| #undef HAVE_CPOW |
| #endif |
| |
| #include <Defn.h> /* -> ../include/R_ext/Complex.h */ |
| #include <Internal.h> |
| #include <Rmath.h> |
| |
| #include "arithmetic.h" /* complex_* */ |
| #include <complex.h> |
| #include "Rcomplex.h" /* I, SET_C99_COMPLEX, toC99 */ |
| #include <R_ext/Itermacros.h> |
| |
| |
| /* interval at which to check interrupts, a guess */ |
| #define NINTERRUPT 10000000 |
| |
| |
| SEXP attribute_hidden complex_unary(ARITHOP_TYPE code, SEXP s1, SEXP call) |
| { |
| R_xlen_t i, n; |
| SEXP ans; |
| |
| switch(code) { |
| case PLUSOP: |
| return s1; |
| case MINUSOP: |
| ans = NO_REFERENCES(s1) ? s1 : duplicate(s1); |
| Rcomplex *pans = COMPLEX(ans); |
| const Rcomplex *ps1 = COMPLEX_RO(s1); |
| n = XLENGTH(s1); |
| for (i = 0; i < n; i++) { |
| Rcomplex x = ps1[i]; |
| pans[i].r = -x.r; |
| pans[i].i = -x.i; |
| } |
| return ans; |
| default: |
| errorcall(call, _("invalid complex unary operator")); |
| } |
| return R_NilValue; /* -Wall */ |
| } |
| |
| static R_INLINE double complex R_cpow_n(double complex X, int k) |
| { |
| if(k == 0) return (double complex) 1.; |
| else if(k == 1) return X; |
| else if(k < 0) return 1. / R_cpow_n(X, -k); |
| else {/* k > 0 */ |
| double complex z = (double complex) 1.;; |
| while (k > 0) { |
| if (k & 1) z = z * X; |
| if (k == 1) break; |
| k >>= 1; /* efficient division by 2; now have k >= 1 */ |
| X = X * X; |
| } |
| return z; |
| } |
| } |
| |
| #if defined(Win32) |
| # undef HAVE_CPOW |
| #endif |
| /* reason for this: |
| 1) X^n (e.g. for n = +/- 2, 3) is unnecessarily inaccurate in glibc; |
| cut-off 65536 : guided from empirical speed measurements |
| |
| 2) On Mingw (but not Mingw-w64) the system cpow is explicitly linked |
| against the (slow) MSVCRT pow, and gets (0+0i)^Y as 0+0i for all Y. |
| |
| 3) PPC macOS crashed on powers of 0+0i (at least under Rosetta). |
| Really 0i^-1 should by Inf+NaNi, but getting that portably seems too hard. |
| (C1x's CMPLX will eventually be possible.) |
| */ |
| |
| static double complex mycpow (double complex X, double complex Y) |
| { |
| double complex Z; |
| double yr = creal(Y), yi = cimag(Y); |
| int k; |
| if (X == 0.0) { |
| if (yi == 0.0) Z = R_pow(0.0, yr); else Z = R_NaN + R_NaN*I; |
| } else if (yi == 0.0 && yr == (k = (int) yr) && abs(k) <= 65536) |
| Z = R_cpow_n(X, k); |
| else |
| #ifdef HAVE_CPOW |
| Z = cpow(X, Y); |
| #else |
| { |
| /* Used for FreeBSD and MingGW, hence mainly with gcc */ |
| double rho, r, i, theta; |
| r = hypot(creal(X), cimag(X)); |
| i = atan2(cimag(X), creal(X)); |
| theta = i * yr; |
| if (yi == 0.0) |
| rho = pow(r, yr); |
| else { |
| /* rearrangement of cexp(X * clog(Y)) */ |
| r = log(r); |
| theta += r * yi; |
| rho = exp(r * yr - i * yi); |
| } |
| #ifdef __GNUC__ |
| __real__ Z = rho * cos(theta); |
| __imag__ Z = rho * sin(theta); |
| #else |
| Z = rho * cos(theta) + (rho * sin(theta)) * I; |
| #endif |
| } |
| #endif |
| return Z; |
| } |
| |
| |
| |
| SEXP attribute_hidden complex_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2) |
| { |
| R_xlen_t i, i1, i2, n, n1, n2; |
| SEXP ans; |
| |
| /* Note: "s1" and "s2" are protected in the calling code. */ |
| n1 = XLENGTH(s1); |
| n2 = XLENGTH(s2); |
| /* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */ |
| if (n1 == 0 || n2 == 0) return(allocVector(CPLXSXP, 0)); |
| |
| n = (n1 > n2) ? n1 : n2; |
| ans = R_allocOrReuseVector(s1, s2, CPLXSXP, n); |
| PROTECT(ans); |
| |
| Rcomplex *pans = COMPLEX(ans); |
| const Rcomplex *ps1 = COMPLEX_RO(s1); |
| const Rcomplex *ps2 = COMPLEX_RO(s2); |
| |
| switch (code) { |
| case PLUSOP: |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| Rcomplex x1 = ps1[i1]; |
| Rcomplex x2 = ps2[i2]; |
| pans[i].r = x1.r + x2.r; |
| pans[i].i = x1.i + x2.i; |
| }); |
| break; |
| case MINUSOP: |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| Rcomplex x1 = ps1[i1]; |
| Rcomplex x2 = ps2[i2]; |
| pans[i].r = x1.r - x2.r; |
| pans[i].i = x1.i - x2.i; |
| }); |
| break; |
| case TIMESOP: |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| SET_C99_COMPLEX(pans, i, |
| toC99(&ps1[i1]) * toC99(&ps2[i2])); |
| }); |
| break; |
| case DIVOP: |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| SET_C99_COMPLEX(pans, i, |
| toC99(&ps1[i1]) / toC99(&ps2[i2])); |
| }); |
| break; |
| case POWOP: |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| SET_C99_COMPLEX(pans, i, |
| mycpow(toC99(&ps1[i1]), toC99(&ps2[i2]))); |
| }); |
| break; |
| default: |
| error(_("unimplemented complex operation")); |
| } |
| UNPROTECT(1); |
| |
| /* quick return if there are no attributes */ |
| if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue) |
| return ans; |
| |
| /* Copy attributes from longer argument. */ |
| |
| if (ans != s2 && n == n2 && ATTRIB(s2) != R_NilValue) |
| copyMostAttrib(s2, ans); |
| if (ans != s1 && n == n1 && ATTRIB(s1) != R_NilValue) |
| copyMostAttrib(s1, ans); /* Done 2nd so s1's attrs overwrite s2's */ |
| |
| return ans; |
| } |
| |
| SEXP attribute_hidden do_cmathfuns(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP x, y = R_NilValue; /* -Wall*/ |
| R_xlen_t i, n; |
| |
| checkArity(op, args); |
| check1arg(args, call, "z"); |
| if (DispatchGroup("Complex", call, op, args, env, &x)) |
| return x; |
| x = CAR(args); |
| if (isComplex(x)) { |
| n = XLENGTH(x); |
| const Rcomplex *px = COMPLEX_RO(x); |
| switch(PRIMVAL(op)) { |
| case 1: /* Re */ |
| { |
| y = allocVector(REALSXP, n); |
| double *py = REAL(y); |
| for(i = 0 ; i < n ; i++) |
| py[i] = px[i].r; |
| } |
| break; |
| case 2: /* Im */ |
| { |
| y = allocVector(REALSXP, n); |
| double *py = REAL(y); |
| for(i = 0 ; i < n ; i++) |
| py[i] = px[i].i; |
| } |
| break; |
| case 3: /* Mod */ |
| case 6: /* abs */ |
| { |
| y = allocVector(REALSXP, n); |
| double *py = REAL(y); |
| for(i = 0 ; i < n ; i++) |
| #if HAVE_CABS |
| py[i] = cabs(toC99(&px[i])); |
| #else |
| py[i] = hypot(px[i].r, px[i].i); |
| #endif |
| } |
| break; |
| case 4: /* Arg */ |
| { |
| y = allocVector(REALSXP, n); |
| double *py = REAL(y); |
| for(i = 0 ; i < n ; i++) |
| #if HAVE_CARG |
| py[i] = carg(toC99(&px[i])); |
| #else |
| py[i] = atan2(px[i].i, px[i].r); |
| #endif |
| } |
| break; |
| case 5: /* Conj */ |
| { |
| y = NO_REFERENCES(x) ? x : allocVector(CPLXSXP, n); |
| Rcomplex *py = COMPLEX(y); |
| for(i = 0 ; i < n ; i++) { |
| py[i].r = px[i].r; |
| py[i].i = -px[i].i; |
| } |
| } |
| break; |
| } |
| } |
| else if(isNumeric(x)) { /* so no complex numbers involved */ |
| n = XLENGTH(x); |
| if(isReal(x)) PROTECT(x); |
| else PROTECT(x = coerceVector(x, REALSXP)); |
| y = NO_REFERENCES(x) ? x : allocVector(REALSXP, n); |
| double *py = REAL(y); |
| const double *px = REAL_RO(x); |
| |
| switch(PRIMVAL(op)) { |
| case 1: /* Re */ |
| case 5: /* Conj */ |
| for(i = 0 ; i < n ; i++) |
| py[i] = px[i]; |
| break; |
| case 2: /* Im */ |
| for(i = 0 ; i < n ; i++) |
| py[i] = 0.0; |
| break; |
| case 4: /* Arg */ |
| for(i = 0 ; i < n ; i++) |
| if(ISNAN(px[i])) |
| py[i] = px[i]; |
| else if (px[i] >= 0) |
| py[i] = 0; |
| else |
| py[i] = M_PI; |
| break; |
| case 3: /* Mod */ |
| case 6: /* abs */ |
| for(i = 0 ; i < n ; i++) |
| py[i] = fabs(px[i]); |
| break; |
| } |
| UNPROTECT(1); |
| } |
| else errorcall(call, _("non-numeric argument to function")); |
| |
| if (x != y && ATTRIB(x) != R_NilValue) { |
| PROTECT(x); |
| PROTECT(y); |
| SHALLOW_DUPLICATE_ATTRIB(y, x); |
| UNPROTECT(2); |
| } |
| return y; |
| } |
| |
| /* used in format.c and printutils.c */ |
| #define MAX_DIGITS 22 |
| void attribute_hidden z_prec_r(Rcomplex *r, const Rcomplex *x, double digits) |
| { |
| double m = 0.0, m1, m2; |
| int dig, mag; |
| |
| r->r = x->r; r->i = x->i; |
| m1 = fabs(x->r); m2 = fabs(x->i); |
| if(R_FINITE(m1)) m = m1; |
| if(R_FINITE(m2) && m2 > m) m = m2; |
| if (m == 0.0) return; |
| if (!R_FINITE(digits)) { |
| if(digits > 0) return; else {r->r = r->i = 0.0; return ;} |
| } |
| dig = (int)floor(digits+0.5); |
| if (dig > MAX_DIGITS) return; else if (dig < 1) dig = 1; |
| mag = (int)floor(log10(m)); |
| dig = dig - mag - 1; |
| if (dig > 306) { |
| double pow10 = 1.0e4; |
| digits = (double)(dig - 4); |
| r->r = fround(pow10 * x->r, digits)/pow10; |
| r->i = fround(pow10 * x->i, digits)/pow10; |
| } else { |
| digits = (double)(dig); |
| r->r = fround(x->r, digits); |
| r->i = fround(x->i, digits); |
| } |
| } |
| |
| /* These substitute functions are rarely used, and not extensively |
| tested, e.g. over CRAN. Please do not change without very good |
| reason! |
| |
| Currently (Feb 2011) they are used on FreeBSD. |
| */ |
| |
| #ifndef HAVE_CLOG |
| #define clog R_clog |
| /* FIXME: maybe add full IEC60559 support */ |
| static double complex clog(double complex x) |
| { |
| double xr = creal(x), xi = cimag(x); |
| return log(hypot(xr, xi)) + atan2(xi, xr)*I; |
| } |
| #endif |
| |
| #ifndef HAVE_CSQRT |
| #define csqrt R_csqrt |
| /* FreeBSD does have this one */ |
| static double complex csqrt(double complex x) |
| { |
| return mycpow(x, 0.5+0.0*I); |
| } |
| #endif |
| |
| #ifndef HAVE_CEXP |
| #define cexp R_cexp |
| /* FIXME: check/add full IEC60559 support */ |
| static double complex cexp(double complex x) |
| { |
| double expx = exp(creal(x)), y = cimag(x); |
| return expx * cos(y) + (expx * sin(y)) * I; |
| } |
| #endif |
| |
| #ifndef HAVE_CCOS |
| #define ccos R_ccos |
| static double complex ccos(double complex x) |
| { |
| double xr = creal(x), xi = cimag(x); |
| return cos(xr)*cosh(xi) - sin(xr)*sinh(xi)*I; /* A&S 4.3.56 */ |
| } |
| #endif |
| |
| #ifndef HAVE_CSIN |
| #define csin R_csin |
| static double complex csin(double complex x) |
| { |
| double xr = creal(x), xi = cimag(x); |
| return sin(xr)*cosh(xi) + cos(xr)*sinh(xi)*I; /* A&S 4.3.55 */ |
| } |
| #endif |
| |
| #ifndef HAVE_CTAN |
| #define ctan R_ctan |
| static double complex ctan(double complex z) |
| { |
| /* A&S 4.3.57 */ |
| double x2, y2, den, ri; |
| x2 = 2.0 * creal(z); |
| y2 = 2.0 * cimag(z); |
| den = cos(x2) + cosh(y2); |
| /* any threshold between -log(DBL_EPSILON) and log(DBL_XMAX) will do*/ |
| if (ISNAN(y2) || fabs(y2) < 50.0) ri = sinh(y2)/den; |
| else ri = (y2 < 0 ? -1.0 : 1.0); |
| return sin(x2)/den + ri * I; |
| } |
| #endif |
| |
| #ifndef HAVE_CASIN |
| #define casin R_casin |
| static double complex casin(double complex z) |
| { |
| /* A&S 4.4.37 */ |
| double alpha, t1, t2, x = creal(z), y = cimag(z), ri; |
| t1 = 0.5 * hypot(x + 1, y); |
| t2 = 0.5 * hypot(x - 1, y); |
| alpha = t1 + t2; |
| ri = log(alpha + sqrt(alpha*alpha - 1)); |
| /* This comes from |
| 'z_asin() is continuous from below if x >= 1 |
| and continuous from above if x <= -1.' |
| */ |
| if(y < 0 || (y == 0 && x > 1)) ri *= -1; |
| return asin(t1 - t2) + ri*I; |
| } |
| #endif |
| |
| #ifndef HAVE_CACOS |
| #define cacos R_cacos |
| static double complex cacos(double complex z) |
| { |
| return M_PI_2 - casin(z); |
| } |
| #endif |
| |
| #ifndef HAVE_CATAN |
| #define catan R_catan |
| static double complex catan(double complex z) |
| { |
| double x = creal(z), y = cimag(z), rr, ri; |
| rr = 0.5 * atan2(2 * x, (1 - x * x - y * y)); |
| ri = 0.25 * log((x * x + (y + 1) * (y + 1)) / |
| (x * x + (y - 1) * (y - 1))); |
| return rr + ri*I; |
| } |
| #endif |
| |
| #ifndef HAVE_CCOSH |
| #define ccosh R_ccosh |
| static double complex ccosh(double complex z) |
| { |
| return ccos(z * I); /* A&S 4.5.8 */ |
| } |
| #endif |
| |
| #ifndef HAVE_CSINH |
| #define csinh R_csinh |
| static double complex csinh(double complex z) |
| { |
| return -I * csin(z * I); /* A&S 4.5.7 */ |
| } |
| #endif |
| |
| static double complex z_tan(double complex z) |
| { |
| double y = cimag(z); |
| double complex r = ctan(z); |
| if(R_FINITE(y) && fabs(y) > 25.0) { |
| /* at this point the real part is nearly zero, and the |
| imaginary part is one: but some OSes get the imag as NaN */ |
| #if __GNUC__ |
| __imag__ r = y < 0 ? -1.0 : 1.0; |
| #else |
| r = creal(r) + (y < 0 ? -1.0 : 1.0) * I; |
| #endif |
| } |
| return r; |
| } |
| |
| #ifndef HAVE_CTANH |
| #define ctanh R_ctanh |
| static double complex ctanh(double complex z) |
| { |
| return -I * z_tan(z * I); /* A&S 4.5.9 */ |
| } |
| #endif |
| |
| |
| /* Don't rely on the OS at the branch cuts */ |
| |
| static double complex z_asin(double complex z) |
| { |
| if(cimag(z) == 0 && fabs(creal(z)) > 1) { |
| double alpha, t1, t2, x = creal(z), ri; |
| t1 = 0.5 * fabs(x + 1); |
| t2 = 0.5 * fabs(x - 1); |
| alpha = t1 + t2; |
| ri = log(alpha + sqrt(alpha*alpha - 1)); |
| if(x > 1) ri *= -1; |
| return asin(t1 - t2) + ri*I; |
| } |
| return casin(z); |
| } |
| |
| static double complex z_acos(double complex z) |
| { |
| if(cimag(z) == 0 && fabs(creal(z)) > 1) return M_PI_2 - z_asin(z); |
| return cacos(z); |
| } |
| |
| static double complex z_atan(double complex z) |
| { |
| if(creal(z) == 0 && fabs(cimag(z)) > 1) { |
| double y = cimag(z), rr, ri; |
| rr = (y > 0) ? M_PI_2 : -M_PI_2; |
| ri = 0.25 * log(((y + 1) * (y + 1))/((y - 1) * (y - 1))); |
| return rr + ri*I; |
| } |
| return catan(z); |
| } |
| |
| static double complex z_acosh(double complex z) |
| { |
| return z_acos(z) * I; |
| } |
| |
| static double complex z_asinh(double complex z) |
| { |
| return -I * z_asin(z * I); |
| } |
| |
| static double complex z_atanh(double complex z) |
| { |
| return -I * z_atan(z * I); |
| } |
| |
| static Rboolean cmath1(double complex (*f)(double complex), |
| const Rcomplex *x, Rcomplex *y, R_xlen_t n) |
| { |
| R_xlen_t i; |
| Rboolean naflag = FALSE; |
| for (i = 0 ; i < n ; i++) { |
| if (ISNA(x[i].r) || ISNA(x[i].i)) { |
| y[i].r = NA_REAL; y[i].i = NA_REAL; |
| } else { |
| SET_C99_COMPLEX(y, i, f(toC99(x + i))); |
| if ( (ISNAN(y[i].r) || ISNAN(y[i].i)) && |
| !(ISNAN(x[i].r) || ISNAN(x[i].i)) ) naflag = TRUE; |
| } |
| } |
| return naflag; |
| } |
| |
| SEXP attribute_hidden complex_math1(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP x, y; |
| R_xlen_t n; |
| Rboolean naflag = FALSE; |
| |
| PROTECT(x = CAR(args)); |
| n = XLENGTH(x); |
| PROTECT(y = allocVector(CPLXSXP, n)); |
| |
| const Rcomplex *px = COMPLEX_RO(x); |
| Rcomplex *py = COMPLEX(y); |
| |
| switch (PRIMVAL(op)) { |
| case 10003: naflag = cmath1(clog, px, py, n); break; |
| case 3: naflag = cmath1(csqrt, px, py, n); break; |
| case 10: naflag = cmath1(cexp, px, py, n); break; |
| case 20: naflag = cmath1(ccos, px, py, n); break; |
| case 21: naflag = cmath1(csin, px, py, n); break; |
| case 22: naflag = cmath1(z_tan, px, py, n); break; |
| case 23: naflag = cmath1(z_acos, px, py, n); break; |
| case 24: naflag = cmath1(z_asin, px, py, n); break; |
| case 25: naflag = cmath1(z_atan, px, py, n); break; |
| case 30: naflag = cmath1(ccosh, px, py, n); break; |
| case 31: naflag = cmath1(csinh, px, py, n); break; |
| case 32: naflag = cmath1(ctanh, px, py, n); break; |
| case 33: naflag = cmath1(z_acosh, px, py, n); break; |
| case 34: naflag = cmath1(z_asinh, px, py, n); break; |
| case 35: naflag = cmath1(z_atanh, px, py, n); break; |
| |
| default: |
| /* such as sign, gamma */ |
| errorcall(call, _("unimplemented complex function")); |
| } |
| if (naflag) |
| warningcall(call, "NaNs produced in function \"%s\"", PRIMNAME(op)); |
| SHALLOW_DUPLICATE_ATTRIB(y, x); |
| UNPROTECT(2); |
| return y; |
| } |
| |
| static void z_rround(Rcomplex *r, Rcomplex *x, Rcomplex *p) |
| { |
| r->r = fround(x->r, p->r); |
| r->i = fround(x->i, p->r); |
| } |
| |
| static void z_prec(Rcomplex *r, Rcomplex *x, Rcomplex *p) |
| { |
| z_prec_r(r, x, p->r); |
| } |
| |
| static void z_logbase(Rcomplex *r, Rcomplex *z, Rcomplex *base) |
| { |
| double complex dz = toC99(z), dbase = toC99(base); |
| SET_C99_COMPLEX(r, 0, clog(dz)/clog(dbase)); |
| } |
| |
| static void z_atan2(Rcomplex *r, Rcomplex *csn, Rcomplex *ccs) |
| { |
| double complex dr, dcsn = toC99(csn), dccs = toC99(ccs); |
| if (dccs == 0) { |
| if(dcsn == 0) { |
| r->r = NA_REAL; r->i = NA_REAL; /* Why not R_NaN? */ |
| return; |
| } else { |
| double y = creal(dcsn); |
| if (ISNAN(y)) dr = y; |
| else dr = ((y >= 0) ? M_PI_2 : -M_PI_2); |
| } |
| } else { |
| dr = catan(dcsn / dccs); |
| if(creal(dccs) < 0) dr += M_PI; |
| if(creal(dr) > M_PI) dr -= 2 * M_PI; |
| } |
| SET_C99_COMPLEX(r, 0, dr); |
| } |
| |
| |
| /* Complex Functions of Two Arguments */ |
| |
| typedef void (*cm2_fun)(Rcomplex *, Rcomplex *, Rcomplex *); |
| SEXP attribute_hidden complex_math2(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| R_xlen_t i, n, na, nb, ia, ib; |
| Rcomplex ai, bi, *y; |
| const Rcomplex *a, *b; |
| SEXP sa, sb, sy; |
| Rboolean naflag = FALSE; |
| cm2_fun f; |
| |
| switch (PRIMVAL(op)) { |
| case 0: /* atan2 */ |
| f = z_atan2; break; |
| case 10001: /* round */ |
| f = z_rround; break; |
| case 2: /* passed from do_log1arg */ |
| case 10: |
| case 10003: /* passed from do_log */ |
| f = z_logbase; break; |
| case 10004: /* signif */ |
| f = z_prec; break; |
| default: |
| error_return(_("unimplemented complex function")); |
| } |
| |
| PROTECT(sa = coerceVector(CAR(args), CPLXSXP)); |
| PROTECT(sb = coerceVector(CADR(args), CPLXSXP)); |
| na = XLENGTH(sa); nb = XLENGTH(sb); |
| if ((na == 0) || (nb == 0)) { |
| UNPROTECT(2); |
| return(allocVector(CPLXSXP, 0)); |
| } |
| n = (na < nb) ? nb : na; |
| PROTECT(sy = allocVector(CPLXSXP, n)); |
| a = COMPLEX_RO(sa); b = COMPLEX_RO(sb); |
| y = COMPLEX(sy); |
| MOD_ITERATE2(n, na, nb, i, ia, ib, { |
| ai = a[ia]; bi = b[ib]; |
| if(ISNA(ai.r) && ISNA(ai.i) && |
| ISNA(bi.r) && ISNA(bi.i)) { |
| y[i].r = NA_REAL; y[i].i = NA_REAL; |
| } else { |
| f(&y[i], &ai, &bi); |
| if ( (ISNAN(y[i].r) || ISNAN(y[i].i)) && |
| !(ISNAN(ai.r) || ISNAN(ai.i) || ISNAN(bi.r) || ISNAN(bi.i)) ) |
| naflag = TRUE; |
| } |
| }); |
| if (naflag) |
| warning("NaNs produced in function \"%s\"", PRIMNAME(op)); |
| if(n == na) { |
| SHALLOW_DUPLICATE_ATTRIB(sy, sa); |
| } else if(n == nb) { |
| SHALLOW_DUPLICATE_ATTRIB(sy, sb); |
| } |
| UNPROTECT(3); |
| return sy; |
| } |
| |
| SEXP attribute_hidden do_complex(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| /* complex(length, real, imaginary) */ |
| SEXP ans, re, im; |
| R_xlen_t i, na, nr, ni; |
| |
| checkArity(op, args); |
| na = asInteger(CAR(args)); |
| if(na == NA_INTEGER || na < 0) |
| error(_("invalid length")); |
| PROTECT(re = coerceVector(CADR(args), REALSXP)); |
| PROTECT(im = coerceVector(CADDR(args), REALSXP)); |
| nr = XLENGTH(re); |
| ni = XLENGTH(im); |
| /* is always true: if (na >= 0) {*/ |
| na = (nr > na) ? nr : na; |
| na = (ni > na) ? ni : na; |
| /* }*/ |
| ans = allocVector(CPLXSXP, na); |
| Rcomplex *pans = COMPLEX(ans); |
| for(i=0 ; i<na ; i++) { |
| pans[i].r = 0; |
| pans[i].i = 0; |
| } |
| UNPROTECT(2); |
| if(na > 0 && nr > 0) { |
| const double *p_re = REAL_RO(re); |
| for(i=0 ; i<na ; i++) |
| pans[i].r = p_re[i%nr]; |
| } |
| if(na > 0 && ni > 0) { |
| const double *p_im = REAL_RO(im); |
| for(i=0 ; i<na ; i++) |
| pans[i].i = p_im[i%ni]; |
| } |
| return ans; |
| } |
| |
| static void R_cpolyroot(double *opr, double *opi, int *degree, |
| double *zeror, double *zeroi, Rboolean *fail); |
| |
| SEXP attribute_hidden do_polyroot(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| SEXP z, zr, zi, r, rr, ri; |
| Rboolean fail; |
| int degree, i, n; |
| |
| checkArity(op, args); |
| z = CAR(args); |
| switch(TYPEOF(z)) { |
| case CPLXSXP: |
| PROTECT(z); |
| break; |
| case REALSXP: |
| case INTSXP: |
| case LGLSXP: |
| PROTECT(z = coerceVector(z, CPLXSXP)); |
| break; |
| default: |
| UNIMPLEMENTED_TYPE("polyroot", z); |
| } |
| #ifdef LONG_VECTOR_SUPPORT |
| R_xlen_t nn = XLENGTH(z); |
| if (nn > R_SHORT_LEN_MAX) error("long vectors are not supported"); |
| n = (int) nn; |
| #else |
| n = LENGTH(z); |
| #endif |
| const Rcomplex *pz = COMPLEX_RO(z); |
| degree = 0; |
| for(i = 0; i < n; i++) { |
| if(pz[i].r!= 0.0 || pz[i].i != 0.0) degree = i; |
| } |
| n = degree + 1; /* omit trailing zeroes */ |
| if(degree >= 1) { |
| PROTECT(rr = allocVector(REALSXP, n)); |
| PROTECT(ri = allocVector(REALSXP, n)); |
| PROTECT(zr = allocVector(REALSXP, n)); |
| PROTECT(zi = allocVector(REALSXP, n)); |
| |
| double *p_rr = REAL(rr); |
| double *p_ri = REAL(ri); |
| double *p_zr = REAL(zr); |
| double *p_zi = REAL(zi); |
| |
| for(i = 0 ; i < n ; i++) { |
| if(!R_FINITE(pz[i].r) || !R_FINITE(pz[i].i)) |
| error(_("invalid polynomial coefficient")); |
| p_zr[degree-i] = pz[i].r; |
| p_zi[degree-i] = pz[i].i; |
| } |
| R_cpolyroot(p_zr, p_zi, °ree, p_rr, p_ri, &fail); |
| if(fail) error(_("root finding code failed")); |
| UNPROTECT(2); |
| r = allocVector(CPLXSXP, degree); |
| Rcomplex *pr = COMPLEX(r); |
| for(i = 0 ; i < degree ; i++) { |
| pr[i].r = p_rr[i]; |
| pr[i].i = p_ri[i]; |
| } |
| UNPROTECT(3); |
| } |
| else { |
| UNPROTECT(1); |
| r = allocVector(CPLXSXP, 0); |
| } |
| return r; |
| } |
| |
| /* Formerly src/appl/cpoly.c: |
| * |
| * Copyright (C) 1997-1998 Ross Ihaka |
| * Copyright (C) 1999-2001 R Core Team |
| * |
| * cpoly finds the zeros of a complex polynomial. |
| * |
| * On Entry |
| * |
| * opr, opi - double precision vectors of real and |
| * imaginary parts of the coefficients in |
| * order of decreasing powers. |
| * |
| * degree - int degree of polynomial. |
| * |
| * |
| * On Return |
| * |
| * zeror, zeroi - output double precision vectors of |
| * real and imaginary parts of the zeros. |
| * |
| * fail - output int parameter, true only if |
| * leading coefficient is zero or if cpoly |
| * has found fewer than degree zeros. |
| * |
| * The program has been written to reduce the chance of overflow |
| * occurring. If it does occur, there is still a possibility that |
| * the zerofinder will work provided the overflowed quantity is |
| * replaced by a large number. |
| * |
| * This is a C translation of the following. |
| * |
| * TOMS Algorithm 419 |
| * Jenkins and Traub. |
| * Comm. ACM 15 (1972) 97-99. |
| * |
| * Ross Ihaka |
| * February 1997 |
| */ |
| |
| #include <R_ext/Arith.h> /* for declaration of hypot */ |
| #include <R_ext/Memory.h> /* for declaration of R_alloc */ |
| |
| #include <float.h> /* for FLT_RADIX */ |
| |
| #include <Rmath.h> /* for R_pow_di */ |
| |
| static void calct(Rboolean *); |
| static Rboolean fxshft(int, double *, double *); |
| static Rboolean vrshft(int, double *, double *); |
| static void nexth(Rboolean); |
| static void noshft(int); |
| |
| static void polyev(int, double, double, |
| double *, double *, double *, double *, double *, double *); |
| static double errev(int, double *, double *, double, double, double, double); |
| static double cpoly_cauchy(int, double *, double *); |
| static double cpoly_scale(int, double *, double, double, double, double); |
| static void cdivid(double, double, double, double, double *, double *); |
| |
| /* Global Variables (too many!) */ |
| |
| static int nn; |
| static double *pr, *pi, *hr, *hi, *qpr, *qpi, *qhr, *qhi, *shr, *shi; |
| static double sr, si; |
| static double tr, ti; |
| static double pvr, pvi; |
| |
| static const double eta = DBL_EPSILON; |
| static const double are = /* eta = */DBL_EPSILON; |
| static const double mre = 2. * M_SQRT2 * /* eta, i.e. */DBL_EPSILON; |
| static const double infin = DBL_MAX; |
| |
| static void R_cpolyroot(double *opr, double *opi, int *degree, |
| double *zeror, double *zeroi, Rboolean *fail) |
| { |
| static const double smalno = DBL_MIN; |
| static const double base = (double)FLT_RADIX; |
| static int d_n, i, i1, i2; |
| static double zi, zr, xx, yy; |
| static double bnd, xxx; |
| Rboolean conv; |
| int d1; |
| double *tmp; |
| static const double cosr =/* cos 94 */ -0.06975647374412529990; |
| static const double sinr =/* sin 94 */ 0.99756405025982424767; |
| xx = M_SQRT1_2;/* 1/sqrt(2) = 0.707.... */ |
| |
| yy = -xx; |
| *fail = FALSE; |
| |
| nn = *degree; |
| d1 = nn - 1; |
| |
| /* algorithm fails if the leading coefficient is zero. */ |
| |
| if (opr[0] == 0. && opi[0] == 0.) { |
| *fail = TRUE; |
| return; |
| } |
| |
| /* remove the zeros at the origin if any. */ |
| |
| while (opr[nn] == 0. && opi[nn] == 0.) { |
| d_n = d1-nn+1; |
| zeror[d_n] = 0.; |
| zeroi[d_n] = 0.; |
| nn--; |
| } |
| nn++; |
| /*-- Now, global var. nn := #{coefficients} = (relevant degree)+1 */ |
| |
| if (nn == 1) return; |
| |
| /* Use a single allocation as these as small */ |
| const void *vmax = vmaxget(); |
| tmp = (double *) R_alloc((size_t) (10*nn), sizeof(double)); |
| pr = tmp; pi = tmp + nn; hr = tmp + 2*nn; hi = tmp + 3*nn; |
| qpr = tmp + 4*nn; qpi = tmp + 5*nn; qhr = tmp + 6*nn; qhi = tmp + 7*nn; |
| shr = tmp + 8*nn; shi = tmp + 9*nn; |
| |
| /* make a copy of the coefficients and shr[] = | p[] | */ |
| for (i = 0; i < nn; i++) { |
| pr[i] = opr[i]; |
| pi[i] = opi[i]; |
| shr[i] = hypot(pr[i], pi[i]); |
| } |
| |
| /* scale the polynomial with factor 'bnd'. */ |
| bnd = cpoly_scale(nn, shr, eta, infin, smalno, base); |
| if (bnd != 1.) { |
| for (i=0; i < nn; i++) { |
| pr[i] *= bnd; |
| pi[i] *= bnd; |
| } |
| } |
| |
| /* start the algorithm for one zero */ |
| |
| while (nn > 2) { |
| |
| /* calculate bnd, a lower bound on the modulus of the zeros. */ |
| |
| for (i=0 ; i < nn ; i++) |
| shr[i] = hypot(pr[i], pi[i]); |
| bnd = cpoly_cauchy(nn, shr, shi); |
| |
| /* outer loop to control 2 major passes */ |
| /* with different sequences of shifts */ |
| |
| for (i1 = 1; i1 <= 2; i1++) { |
| |
| /* first stage calculation, no shift */ |
| |
| noshft(5); |
| |
| /* inner loop to select a shift */ |
| for (i2 = 1; i2 <= 9; i2++) { |
| |
| /* shift is chosen with modulus bnd */ |
| /* and amplitude rotated by 94 degrees */ |
| /* from the previous shift */ |
| |
| xxx= cosr * xx - sinr * yy; |
| yy = sinr * xx + cosr * yy; |
| xx = xxx; |
| sr = bnd * xx; |
| si = bnd * yy; |
| |
| /* second stage calculation, fixed shift */ |
| |
| conv = fxshft(i2 * 10, &zr, &zi); |
| if (conv) |
| goto L10; |
| } |
| } |
| |
| /* the zerofinder has failed on two major passes */ |
| /* return empty handed */ |
| |
| *fail = TRUE; |
| vmaxset(vmax); |
| return; |
| |
| /* the second stage jumps directly to the third stage iteration. |
| * if successful, the zero is stored and the polynomial deflated. |
| */ |
| L10: |
| d_n = d1+2 - nn; |
| zeror[d_n] = zr; |
| zeroi[d_n] = zi; |
| --nn; |
| for (i=0; i < nn ; i++) { |
| pr[i] = qpr[i]; |
| pi[i] = qpi[i]; |
| } |
| }/*while*/ |
| |
| /* calculate the final zero and return */ |
| cdivid(-pr[1], -pi[1], pr[0], pi[0], &zeror[d1], &zeroi[d1]); |
| |
| vmaxset(vmax); |
| return; |
| } |
| |
| |
| /* Computes the derivative polynomial as the initial |
| * polynomial and computes l1 no-shift h polynomials. */ |
| |
| static void noshft(int l1) |
| { |
| int i, j, jj, n = nn - 1, nm1 = n - 1; |
| |
| double t1, t2, xni; |
| |
| for (i=0; i < n; i++) { |
| xni = (double)(nn - i - 1); |
| hr[i] = xni * pr[i] / n; |
| hi[i] = xni * pi[i] / n; |
| } |
| |
| for (jj = 1; jj <= l1; jj++) { |
| |
| if (hypot(hr[n-1], hi[n-1]) <= |
| eta * 10.0 * hypot(pr[n-1], pi[n-1])) { |
| /* If the constant term is essentially zero, */ |
| /* shift h coefficients. */ |
| |
| for (i = 1; i <= nm1; i++) { |
| j = nn - i; |
| hr[j-1] = hr[j-2]; |
| hi[j-1] = hi[j-2]; |
| } |
| hr[0] = 0.; |
| hi[0] = 0.; |
| } |
| else { |
| cdivid(-pr[nn-1], -pi[nn-1], hr[n-1], hi[n-1], &tr, &ti); |
| for (i = 1; i <= nm1; i++) { |
| j = nn - i; |
| t1 = hr[j-2]; |
| t2 = hi[j-2]; |
| hr[j-1] = tr * t1 - ti * t2 + pr[j-1]; |
| hi[j-1] = tr * t2 + ti * t1 + pi[j-1]; |
| } |
| hr[0] = pr[0]; |
| hi[0] = pi[0]; |
| } |
| } |
| } |
| |
| |
| /* Computes l2 fixed-shift h polynomials and tests for convergence. |
| * initiates a variable-shift iteration and returns with the |
| * approximate zero if successful. |
| */ |
| static Rboolean fxshft(int l2, double *zr, double *zi) |
| { |
| /* l2 - limit of fixed shift steps |
| * zr,zi - approximate zero if convergence (result TRUE) |
| * |
| * Return value indicates convergence of stage 3 iteration |
| * |
| * Uses global (sr,si), nn, pr[], pi[], .. (all args of polyev() !) |
| */ |
| |
| Rboolean pasd, bool, test; |
| static double svsi, svsr; |
| static int i, j, n; |
| static double oti, otr; |
| |
| n = nn - 1; |
| |
| /* evaluate p at s. */ |
| |
| polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi); |
| |
| test = TRUE; |
| pasd = FALSE; |
| |
| /* calculate first t = -p(s)/h(s). */ |
| |
| calct(&bool); |
| |
| /* main loop for one second stage step. */ |
| |
| for (j=1; j<=l2; j++) { |
| |
| otr = tr; |
| oti = ti; |
| |
| /* compute next h polynomial and new t. */ |
| |
| nexth(bool); |
| calct(&bool); |
| *zr = sr + tr; |
| *zi = si + ti; |
| |
| /* test for convergence unless stage 3 has */ |
| /* failed once or this is the last h polynomial. */ |
| |
| if (!bool && test && j != l2) { |
| if (hypot(tr - otr, ti - oti) >= hypot(*zr, *zi) * 0.5) { |
| pasd = FALSE; |
| } |
| else if (! pasd) { |
| pasd = TRUE; |
| } |
| else { |
| |
| /* the weak convergence test has been */ |
| /* passed twice, start the third stage */ |
| /* iteration, after saving the current */ |
| /* h polynomial and shift. */ |
| |
| for (i = 0; i < n; i++) { |
| shr[i] = hr[i]; |
| shi[i] = hi[i]; |
| } |
| svsr = sr; |
| svsi = si; |
| if (vrshft(10, zr, zi)) { |
| return TRUE; |
| } |
| |
| /* the iteration failed to converge. */ |
| /* turn off testing and restore */ |
| /* h, s, pv and t. */ |
| |
| test = FALSE; |
| for (i=1 ; i<=n ; i++) { |
| hr[i-1] = shr[i-1]; |
| hi[i-1] = shi[i-1]; |
| } |
| sr = svsr; |
| si = svsi; |
| polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi); |
| calct(&bool); |
| } |
| } |
| } |
| |
| /* attempt an iteration with final h polynomial */ |
| /* from second stage. */ |
| |
| return(vrshft(10, zr, zi)); |
| } |
| |
| |
| /* carries out the third stage iteration. |
| */ |
| static Rboolean vrshft(int l3, double *zr, double *zi) |
| { |
| /* l3 - limit of steps in stage 3. |
| * zr,zi - on entry contains the initial iterate; |
| * if the iteration converges it contains |
| * the final iterate on exit. |
| * Returns TRUE if iteration converges |
| * |
| * Assign and uses GLOBAL sr, si |
| */ |
| Rboolean bool, b; |
| static int i, j; |
| static double r1, r2, mp, ms, tp, relstp; |
| static double omp; |
| |
| b = FALSE; |
| sr = *zr; |
| si = *zi; |
| |
| /* main loop for stage three */ |
| |
| for (i = 1; i <= l3; i++) { |
| |
| /* evaluate p at s and test for convergence. */ |
| polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi); |
| |
| mp = hypot(pvr, pvi); |
| ms = hypot(sr, si); |
| if (mp <= 20. * errev(nn, qpr, qpi, ms, mp, /*are=*/eta, mre)) { |
| goto L_conv; |
| } |
| |
| /* polynomial value is smaller in value than */ |
| /* a bound on the error in evaluating p, */ |
| /* terminate the iteration. */ |
| |
| if (i != 1) { |
| |
| if (!b && mp >= omp && relstp < .05) { |
| |
| /* iteration has stalled. probably a */ |
| /* cluster of zeros. do 5 fixed shift */ |
| /* steps into the cluster to force */ |
| /* one zero to dominate. */ |
| |
| tp = relstp; |
| b = TRUE; |
| if (relstp < eta) |
| tp = eta; |
| r1 = sqrt(tp); |
| r2 = sr * (r1 + 1.) - si * r1; |
| si = sr * r1 + si * (r1 + 1.); |
| sr = r2; |
| polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi); |
| for (j = 1; j <= 5; ++j) { |
| calct(&bool); |
| nexth(bool); |
| } |
| omp = infin; |
| goto L10; |
| } |
| else { |
| |
| /* exit if polynomial value */ |
| /* increases significantly. */ |
| |
| if (mp * .1 > omp) |
| return FALSE; |
| } |
| } |
| omp = mp; |
| |
| /* calculate next iterate. */ |
| |
| L10: |
| calct(&bool); |
| nexth(bool); |
| calct(&bool); |
| if (!bool) { |
| relstp = hypot(tr, ti) / hypot(sr, si); |
| sr += tr; |
| si += ti; |
| } |
| } |
| return FALSE; |
| |
| L_conv: |
| *zr = sr; |
| *zi = si; |
| return TRUE; |
| } |
| |
| static void calct(Rboolean *bool) |
| { |
| /* computes t = -p(s)/h(s). |
| * bool - logical, set true if h(s) is essentially zero. */ |
| |
| int n = nn - 1; |
| double hvi, hvr; |
| |
| /* evaluate h(s). */ |
| polyev(n, sr, si, hr, hi, |
| qhr, qhi, &hvr, &hvi); |
| |
| *bool = hypot(hvr, hvi) <= are * 10. * hypot(hr[n-1], hi[n-1]); |
| if (!*bool) { |
| cdivid(-pvr, -pvi, hvr, hvi, &tr, &ti); |
| } |
| else { |
| tr = 0.; |
| ti = 0.; |
| } |
| } |
| |
| static void nexth(Rboolean bool) |
| { |
| /* calculates the next shifted h polynomial. |
| * bool : if TRUE h(s) is essentially zero |
| */ |
| int j, n = nn - 1; |
| double t1, t2; |
| |
| if (!bool) { |
| for (j=1; j < n; j++) { |
| t1 = qhr[j - 1]; |
| t2 = qhi[j - 1]; |
| hr[j] = tr * t1 - ti * t2 + qpr[j]; |
| hi[j] = tr * t2 + ti * t1 + qpi[j]; |
| } |
| hr[0] = qpr[0]; |
| hi[0] = qpi[0]; |
| } |
| else { |
| /* if h(s) is zero replace h with qh. */ |
| |
| for (j=1; j < n; j++) { |
| hr[j] = qhr[j-1]; |
| hi[j] = qhi[j-1]; |
| } |
| hr[0] = 0.; |
| hi[0] = 0.; |
| } |
| } |
| |
| /*--------------------- Independent Complex Polynomial Utilities ----------*/ |
| |
| static |
| void polyev(int n, |
| double s_r, double s_i, |
| double *p_r, double *p_i, |
| double *q_r, double *q_i, |
| double *v_r, double *v_i) |
| { |
| /* evaluates a polynomial p at s by the horner recurrence |
| * placing the partial sums in q and the computed value in v_. |
| */ |
| int i; |
| double t; |
| |
| q_r[0] = p_r[0]; |
| q_i[0] = p_i[0]; |
| *v_r = q_r[0]; |
| *v_i = q_i[0]; |
| for (i = 1; i < n; i++) { |
| t = *v_r * s_r - *v_i * s_i + p_r[i]; |
| q_i[i] = *v_i = *v_r * s_i + *v_i * s_r + p_i[i]; |
| q_r[i] = *v_r = t; |
| } |
| } |
| |
| static |
| double errev(int n, double *qr, double *qi, |
| double ms, double mp, double a_re, double m_re) |
| { |
| /* bounds the error in evaluating the polynomial by the horner |
| * recurrence. |
| * |
| * qr,qi - the partial sum vectors |
| * ms - modulus of the point |
| * mp - modulus of polynomial value |
| * a_re,m_re - error bounds on complex addition and multiplication |
| */ |
| double e; |
| int i; |
| |
| e = hypot(qr[0], qi[0]) * m_re / (a_re + m_re); |
| for (i=0; i < n; i++) |
| e = e*ms + hypot(qr[i], qi[i]); |
| |
| return e * (a_re + m_re) - mp * m_re; |
| } |
| |
| |
| static |
| double cpoly_cauchy(int n, double *pot, double *q) |
| { |
| /* Computes a lower bound on the moduli of the zeros of a polynomial |
| * pot[1:nn] is the modulus of the coefficients. |
| */ |
| double f, x, delf, dx, xm; |
| int i, n1 = n - 1; |
| |
| pot[n1] = -pot[n1]; |
| |
| /* compute upper estimate of bound. */ |
| |
| x = exp((log(-pot[n1]) - log(pot[0])) / (double) n1); |
| |
| /* if newton step at the origin is better, use it. */ |
| |
| if (pot[n1-1] != 0.) { |
| xm = -pot[n1] / pot[n1-1]; |
| if (xm < x) |
| x = xm; |
| } |
| |
| /* chop the interval (0,x) unitl f le 0. */ |
| |
| for(;;) { |
| xm = x * 0.1; |
| f = pot[0]; |
| for (i = 1; i < n; i++) |
| f = f * xm + pot[i]; |
| if (f <= 0.0) { |
| break; |
| } |
| x = xm; |
| } |
| |
| dx = x; |
| |
| /* do Newton iteration until x converges to two decimal places. */ |
| |
| while (fabs(dx / x) > 0.005) { |
| q[0] = pot[0]; |
| for(i = 1; i < n; i++) |
| q[i] = q[i-1] * x + pot[i]; |
| f = q[n1]; |
| delf = q[0]; |
| for(i = 1; i < n1; i++) |
| delf = delf * x + q[i]; |
| dx = f / delf; |
| x -= dx; |
| } |
| return x; |
| } |
| |
| static |
| double cpoly_scale(int n, double *pot, |
| double eps, double BIG, double small, double base) |
| { |
| /* Returns a scale factor to multiply the coefficients of the polynomial. |
| * The scaling is done to avoid overflow and to avoid |
| * undetected underflow interfering with the convergence criterion. |
| * The factor is a power of the base. |
| |
| * pot[1:n] : modulus of coefficients of p |
| * eps,BIG, |
| * small,base - constants describing the floating point arithmetic. |
| */ |
| |
| int i, ell; |
| double x, high, sc, lo, min_, max_; |
| |
| /* find largest and smallest moduli of coefficients. */ |
| high = sqrt(BIG); |
| lo = small / eps; |
| max_ = 0.; |
| min_ = BIG; |
| for (i = 0; i < n; i++) { |
| x = pot[i]; |
| if (x > max_) max_ = x; |
| if (x != 0. && x < min_) |
| min_ = x; |
| } |
| |
| /* scale only if there are very large or very small components. */ |
| |
| if (min_ < lo || max_ > high) { |
| x = lo / min_; |
| if (x <= 1.) |
| sc = 1. / (sqrt(max_) * sqrt(min_)); |
| else { |
| sc = x; |
| if (BIG / sc > max_) |
| sc = 1.0; |
| } |
| ell = (int) (log(sc) / log(base) + 0.5); |
| return R_pow_di(base, ell); |
| } |
| else return 1.0; |
| } |
| |
| |
| static |
| void cdivid(double ar, double ai, double br, double bi, |
| double *cr, double *ci) |
| { |
| /* complex division c = a/b, i.e., (cr +i*ci) = (ar +i*ai) / (br +i*bi), |
| avoiding overflow. */ |
| |
| double d, r; |
| |
| if (br == 0. && bi == 0.) { |
| /* division by zero, c = infinity. */ |
| *cr = *ci = R_PosInf; |
| } |
| else if (fabs(br) >= fabs(bi)) { |
| r = bi / br; |
| d = br + r * bi; |
| *cr = (ar + ai * r) / d; |
| *ci = (ai - ar * r) / d; |
| } |
| else { |
| r = br / bi; |
| d = bi + r * br; |
| *cr = (ar * r + ai) / d; |
| *ci = (ai * r - ar) / d; |
| } |
| } |
| |
| /* static double cpoly_cmod(double *r, double *i) |
| * --> replaced by hypot() everywhere |
| */ |