/*

 *  R : A Computer Language for Statistical Data Analysis

 *  Copyright (C) 2000-2025  The R Core Team

 *  Copyright (C) 2005       The R Foundation

 *  Copyright (C) 1995-1997  Robert Gentleman and Ross Ihaka

 *

 *  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 -pedantic 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>            /* incl "complex math.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

attribute_hidden SEXP 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;

}

attribute_hidden SEXP 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;

}

attribute_hidden SEXP 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;

}

/* Implementing  signif(<complex>)  *and* used in format.c and printutils.c */

#define MAX_DIGITS 22

attribute_hidden void z_prec_r(Rcomplex *r, const Rcomplex *x, double digits)

{

    // Implement    r <- signif(x, digits)

    r->r = x->r; r->i = x->i;

    double m = 0.0,

  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 ;}

    }

    int dig = (int)floor(digits+0.5);

    if (dig > MAX_DIGITS) return; else if (dig < 1) dig = 1;

    int 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 x = creal(z), y = cimag(z),

  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);

}

#ifdef HAVE_CABS

# define R_CABS(Z) cabs(Z)

#else

# define R_CABS(Z) hypot(creal(Z), cimag(Z))

#endif

static bool cmath1(double complex (*f)(double complex),

           const Rcomplex *x, Rcomplex *y, R_xlen_t n)

{

    R_xlen_t i;

    bool 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;

}

attribute_hidden SEXP complex_math1(SEXP call, SEXP op, SEXP args, SEXP env)

{

    SEXP x, y;

    R_xlen_t n;

    bool 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;

  // 1: floor

  // 2: ceil[ing]

    case 3: naflag = cmath1(csqrt, px, py, n); break;

  // 4: sign

    case 10: naflag = cmath1(cexp, px, py, n); break;

  // 11: expm1

  // 12: log1p

    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 *);

attribute_hidden SEXP complex_math2(SEXP call, SEXP op, SEXP args, SEXP env)

{

    cm2_fun f;

    switch (PRIMVAL(op)) {

    case 0: /* atan2 */

  f = z_atan2; break;

    case 10001: /* round */

  f = z_rround; break;

    case 10002: /* passed from do_log1arg */

    case 10010: /*   "     "    " */

    case 10003: /* passed from do_log */

  f = z_logbase; break;

    case 10004: /* signif */

  f = z_prec; break;

    default:

  error_return(_("unimplemented complex function"));

    }

    SEXP

  sa = PROTECT(coerceVector(CAR(args), CPLXSXP)),

  sb = PROTECT(coerceVector(CADR(args), CPLXSXP));

    R_xlen_t na = XLENGTH(sa), nb = XLENGTH(sb), n;

    if ((na == 0) || (nb == 0)) {

  UNPROTECT(2);

  return(allocVector(CPLXSXP, 0));

    }

    n = (na < nb) ? nb : na;

    SEXP sy = PROTECT(allocVector(CPLXSXP, n));

    const Rcomplex

  *a = COMPLEX_RO(sa),

  *b = COMPLEX_RO(sb);

    Rcomplex ai, bi, *y = COMPLEX(sy);

    R_xlen_t i, ia, ib;

    bool naflag = false;

    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;

}

attribute_hidden SEXP 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, bool *fail);

attribute_hidden SEXP do_polyroot(SEXP call, SEXP op, SEXP args, SEXP rho)

{

    checkArity(op, args);

    SEXP 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);

    }

    int i, n;

#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);

    int degree = 0; // := max{i; z[i] != 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 */

    SEXP r;

    if(degree >= 1) {

  SEXP zr, zi, rr, ri;

  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;

  }

  bool fail;

  R_cpolyroot(p_zr, p_zi, &degree, 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.

 *