The R Project SVN R

Rev

Rev 49001 | Rev 49066 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 49001 Rev 49062
Line 80... Line 80...
80 
    }
 80 
    }
81 
    return R_NilValue; /* -Wall */
 81 
    return R_NilValue; /* -Wall */
82 
}
 82 
}
83 
 
 83 
 
84 
#ifndef HAVE_C99_COMPLEX
 84 
#ifndef HAVE_C99_COMPLEX
- 
 
 85 
/* c := a / b --- where c may overwrite 'b' but not 'a' */
85 
static void complex_div(Rcomplex *c, Rcomplex *a, Rcomplex *b)
 86 
static void complex_div(Rcomplex *c, Rcomplex *a, Rcomplex *b)
86 
{
 87 
{
87 
    double ratio, den;
 88 
    double ratio, den;
88 
    double abr, abi;
 89 
    double abr, abi;
89 
 
 90 
 
Line 106... Line 107...
106 
}
 107 
}
107 
#endif
 108 
#endif
108 
 
 109 
 
109 
#ifndef HAVE_C99_COMPLEX
 110 
#ifndef HAVE_C99_COMPLEX
110 
 
 111 
 
- 
 
 112 
/* "export" this (to be API) : */
- 
 
 113 
void R_cpow_n(Rcomplex *r, Rcomplex *x, int k) {
- 
 
 114 
    if(k == 0) {
- 
 
 115 
	r->r = x->r;
- 
 
 116 
	r->i = x->i;
- 
 
 117 
    } else if(k < 0) {
- 
 
 118 
	Rcomplex h;
- 
 
 119 
	R_cpow_n(r, x, -k);
- 
 
 120 
	/* r := 1/r : */
- 
 
 121 
	h.r = 1.; h.i = 0;
- 
 
 122 
	complex_div(r, &h, r);
- 
 
 123 
    }
- 
 
 124 
    else {/* k > 0 */
- 
 
 125 
	Rcomplex X;
- 
 
 126 
	r->r = X.r = x->r;
- 
 
 127 
	r->i = X.i = x->i;
- 
 
 128 
	k--;
- 
 
 129 
	while (k > 0) {
- 
 
 130 
	    double rr;
- 
 
 131 
	    if (k & 1) { /* r := r * x */
- 
 
 132 
		rr   = r->r * x->r - r->i * x->i;
- 
 
 133 
		r->i = r->r * x->i + r->i * x->r;
- 
 
 134 
		r->r = rr;
- 
 
 135 
	    }
- 
 
 136 
	    if(k == 1)
- 
 
 137 
		break;
- 
 
 138 
	    k >>= 1; /* efficient division by 2; now have k >= 1 */
- 
 
 139 
 
- 
 
 140 
	    /* X := X * X : */
- 
 
 141 
	    rr	= (X.r - X.i) * (X.r + X.i); /* = X.r * X.r - X.i * X.i */
- 
 
 142 
	    X.i = 2 * X.r * X.i;	     /* = X.r * X.i + X.i * X.r */
- 
 
 143 
	    X.r = rr;
- 
 
 144 
	}
- 
 
 145 
	return;
- 
 
 146 
    }
- 
 
 147 
}
- 
 
 148 
 
111 
static void complex_pow(Rcomplex *r, Rcomplex *a, Rcomplex *b)
 149 
static void complex_pow(Rcomplex *r, Rcomplex *a, Rcomplex *b)
112 
{
 150 
{
113 
/* r := a^b */
 151 
/* r := a^b */
114 
    double logr, logi, x, y;
 152 
    double logr, logi, x, y;
115 
    int ib;
 153 
    int ib;
Line 119... Line 157...
119 
	    r->r = a->r; r->i = a->i; return;
 157 
	    r->r = a->r; r->i = a->i; return;
120 
	}
 158 
	}
121 
	if(a->i == 0. && a->r >= 0.) {
 159 
	if(a->i == 0. && a->r >= 0.) {
122 
	    r->r = R_pow(a->r, b->r); r->i = 0.; return;
 160 
	    r->r = R_pow(a->r, b->r); r->i = 0.; return;
123 
	}
 161 
	}
124 
	if(a->r == 0. && b->r == (ib = (int)b->r)) {/* (|a|*i)^b */
 162 
	if(b->r == (ib = (int)b->r)) { /* z ^ k  (k integer) */
- 
 
 163 
	    if(a->r == 0.) {/* (|a|*i)^b */
125 
	    x = R_pow_di(a->i, ib);
 164 
		x = R_pow_di(a->i, ib);
126 
	    if(ib % 2) {	/* ib is odd ==> imaginary r */
 165 
		if(ib % 2) {	/* ib is odd ==> imaginary r */
127 
		r->r = 0.;
 166 
		    r->r = 0.;
128 
		r->i = ((ib>0 && ib %4 == 3)||(ib<0 && (-ib)%4 == 1))? -x : x;
 167 
		    r->i = ((ib>0 && ib %4 == 3)||(ib<0 && (-ib)%4 == 1))? -x : x;
129 
	    } else {		/* even exponent b : real r */
 168 
		} else {	/* even exponent b : real r */
130 
		r->r = (ib %4)? -x : x; r->i = 0.;
 169 
		    r->r = (ib %4)? -x : x; r->i = 0.;
- 
 
 170 
		}
- 
 
 171 
	    } else {
- 
 
 172 
		R_cpow_n(r, a, ib);
131 
	    }
 173 
	    }
132 
	    return;
 174 
	    return;
133 
	}
 175 
	}
134 
    }
 176 
    }
135 
    logr = log(hypot(a->r, a->i) );
 177 
    logr = log(hypot(a->r, a->i) );
Line 140... Line 182...
140 
    r->i = x * sin(y);
 182 
    r->i = x * sin(y);
141 
}
 183 
}
142 
 
 184 
 
143 
#else /* HAVE_C99_COMPLEX */
 185 
#else /* HAVE_C99_COMPLEX */
144 
 
 186 
 
- 
 
 187 
 
- 
 
 188 
/* "export" this (to be API) : */
- 
 
 189 
double complex R_cpow_n(double complex X, int k) {
- 
 
 190 
    if(k == 0)
- 
 
 191 
	return X;
- 
 
 192 
    else if(k < 0)
- 
 
 193 
	return 1. / R_cpow_n(X, -k);
- 
 
 194 
    else {/* k > 0 */
- 
 
 195 
	double complex z = X;
- 
 
 196 
	k--;
- 
 
 197 
	while (k > 0) {
- 
 
 198 
	    if (k & 1)
- 
 
 199 
		z = z * X;
- 
 
 200 
	    if(k == 1)
- 
 
 201 
		break;
- 
 
 202 
	    k >>= 1; /* efficient division by 2; now have k >= 1 */
- 
 
 203 
	    /* x := x * x */
- 
 
 204 
	    X = X * X;
- 
 
 205 
	}
- 
 
 206 
	return z;
- 
 
 207 
    }
- 
 
 208 
}
- 
 
 209 
 
145 
#ifdef Win32
 210 
#ifdef Win32
146 
/* Need this because the system one is explicitly linked
 211 
/* Need this because the system one is explicitly linked
147 
   against MSVCRT's pow, and gets (0+0i)^Y as 0+0i for all Y */
 212 
   against MSVCRT's pow, and gets (0+0i)^Y as 0+0i for all Y */
148 
static double complex mycpow (double complex X, double complex Y)
 213 
static double complex mycpow (double complex X, double complex Y)
149 
{
 214 
{
150 
  double complex Res;
 215 
    double complex Res; int k;
151 
  if (X == 0.0) {
 216 
    if (X == 0.0) {
152 
      __real__ Res = R_pow(0.0, __real__ Y);
 217 
	__real__ Res = R_pow(0.0, __real__ Y);
153 
      __imag__ Res = 0.0;
 218 
	__imag__ Res = 0.0;
- 
 
 219 
    } else if (__imag__ Y == 0.0 && __real__ Y == (k = (int)__real__ Y)
- 
 
 220 
	       && abs(k) <= 65536) {
- 
 
 221 
	return R_cpow_n(X, k);
154 
  } else {
 222 
    } else {
155 
      double rho, r,i, theta;
 223 
	double rho, r,i, theta;
156 
      r = hypot (__real__ X, __imag__ X);
 224 
	r = hypot (__real__ X, __imag__ X);
157 
      i = carg (X);
 225 
	i = carg (X);
158 
      theta = i * __real__ Y;
 226 
	theta = i * __real__ Y;
159 
 
 227 
 
160 
      if (__imag__ Y == 0.0)
 228 
	if (__imag__ Y == 0.0)
161 
	  rho = pow (r, __real__ Y);
 229 
	    rho = pow (r, __real__ Y);
162 
      else {
 230 
	else {
163 
	  r = log (r);
 231 
	    r = log (r);
164 
	  /* rearrangement of cexp(X * clog(Y)) */
 232 
	    /* rearrangement of cexp(X * clog(Y)) */
165 
	  theta += r * __imag__ Y;
 233 
	    theta += r * __imag__ Y;
166 
	  rho = exp (r * __real__ Y - i * __imag__ Y);
 234 
	    rho = exp (r * __real__ Y - i * __imag__ Y);
167 
      }
 235 
	}
168 
      __real__ Res = rho * cos (theta);
 236 
	__real__ Res = rho * cos (theta);
169 
      __imag__ Res = rho * sin (theta);
 237 
	__imag__ Res = rho * sin (theta);
170 
  }
 238 
    }
171 
  return  Res;
 239 
    return  Res;
172 
}
 240 
}
173 
#else /* not Win32 */
 241 
#else /* not Win32 */
- 
 
 242 
/* reason for this:
174 
/* reason for this: glibc gets (0+0i)^y = Inf+NaNi for y < 0
 243 
 * 1) glibc gets (0+0i)^y = Inf+NaNi for y < 0
- 
 
 244 
 * 2)   X ^ n  (e.g. for n = 2, 3)  is unnecessarily inaccurate;
- 
 
 245 
 *	cut-off 65536 : guided from empirical speed measurements
175 
*/
 246 
*/
- 
 
 247 
 
176 
static double complex mycpow (double complex X, double complex Y)
 248 
static double complex mycpow (double complex X, double complex Y)
177 
{
 249 
{
178 
    double tmp = cimag(Y);
 250 
    double iY = cimag(Y);
179 
    if (X == 0.0 && tmp == 0) {
 251 
    if (iY == 0) {/* X ^ <real> */
- 
 
 252 
	double complex Z; int k;
- 
 
 253 
	if (X == 0.) {
180 
	double complex Z = R_pow(0.0, creal(Y));
 254 
	    Z = R_pow(0., creal(Y));
- 
 
 255 
	} else if(creal(Y) == (k = (int)creal(Y)) && abs(k) <= 65536) {
- 
 
 256 
	    Z = R_cpow_n(X, k);
- 
 
 257 
	} else {
- 
 
 258 
	    Z = cpow(X, Y);
- 
 
 259 
	}
181 
	return Z;
 260 
	return Z;
182 
    } else
 261 
    } else
183 
	return cpow(X, Y);
 262 
	return cpow(X, Y);
184 
}
 263 
}
185 
#endif
 264 
#endif