Back to index

plt-scheme  4.2.1
complex.c
Go to the documentation of this file.
00001 /*
00002   MzScheme
00003   Copyright (c) 2004-2009 PLT Scheme Inc.
00004   Copyright (c) 1995-2001 Matthew Flatt
00005 
00006     This library is free software; you can redistribute it and/or
00007     modify it under the terms of the GNU Library General Public
00008     License as published by the Free Software Foundation; either
00009     version 2 of the License, or (at your option) any later version.
00010 
00011     This library is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014     Library General Public License for more details.
00015 
00016     You should have received a copy of the GNU Library General Public
00017     License along with this library; if not, write to the Free
00018     Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
00019     Boston, MA 02110-1301 USA.
00020 
00021   libscheme
00022   Copyright (c) 1994 Brent Benson
00023   All rights reserved.
00024 */
00025 
00026 #include "schpriv.h"
00027 #include <ctype.h>
00028 #include <math.h>
00029 
00030 #define zero scheme_exact_zero
00031 
00032 static Scheme_Object *make_complex(const Scheme_Object *r, const Scheme_Object *i,
00033                                int normalize)
00034 {
00035   Scheme_Complex *c;
00036 
00037   c = (Scheme_Complex *)scheme_malloc_small_dirty_tagged(sizeof(Scheme_Complex));
00038   CLEAR_KEY_FIELD(&c->so);
00039   c->so.type = scheme_complex_type;
00040   c->r = (Scheme_Object *)r;
00041   c->i = (Scheme_Object *)i;
00042 
00043   if (normalize)
00044     return scheme_complex_normalize((Scheme_Object *)c);
00045   else
00046     return (Scheme_Object *)c;
00047 }
00048 
00049 Scheme_Object *scheme_make_complex(const Scheme_Object *r, const Scheme_Object *i)
00050 {
00051   return make_complex(r, i, 1);
00052 }
00053 
00054 Scheme_Object *scheme_real_to_complex(const Scheme_Object *n)
00055 {
00056   return make_complex(n, zero, 0);
00057 }
00058 
00059 #ifdef MZ_XFORM
00060 START_XFORM_SKIP;
00061 #endif
00062 
00063 Scheme_Object *scheme_make_small_complex(const Scheme_Object *n, Small_Complex *s)
00064 {
00065   s->so.type = scheme_complex_type;
00066   s->r = (Scheme_Object *)n;
00067   s->i = zero;
00068 
00069   return (Scheme_Object *)s;
00070 }
00071 
00072 #ifdef MZ_XFORM
00073 END_XFORM_SKIP;
00074 #endif
00075 
00076 int scheme_is_complex_exact(const Scheme_Object *o)
00077 {
00078   Scheme_Complex *c = (Scheme_Complex *)o;
00079 
00080   return !SCHEME_FLOATP(c->r) && !SCHEME_FLOATP(c->i);
00081 }
00082 
00083 Scheme_Object *scheme_complex_normalize(const Scheme_Object *o)
00084 {
00085   Scheme_Complex *c = (Scheme_Complex *)o;
00086 
00087   if (c->i == zero)
00088     return c->r;
00089   if (c->r == zero) {
00090     /* No coercions */
00091     return (Scheme_Object *)c; 
00092   }
00093 
00094   if (SCHEME_DBLP(c->i)) {
00095     if (!SCHEME_DBLP(c->r)) {
00096       Scheme_Object *r;
00097       r = scheme_make_double(scheme_get_val_as_double(c->r));
00098       c->r = r;
00099     }
00100   } else if (SCHEME_DBLP(c->r)) {
00101     Scheme_Object *i;
00102     i = scheme_make_double(scheme_get_val_as_double(c->i));
00103     c->i = i;
00104   }
00105 
00106 #ifdef MZ_USE_SINGLE_FLOATS
00107   if (SCHEME_FLTP(c->i)) {
00108     if (!SCHEME_FLTP(c->r)) {
00109       if (SCHEME_DBLP(c->r)) {
00110        c->i = scheme_make_double(SCHEME_FLT_VAL(c->i));
00111       } else {
00112        c->r = scheme_make_float(scheme_get_val_as_float(c->r));
00113       }
00114     }
00115   }
00116 #endif
00117 
00118   return (Scheme_Object *)c;
00119 }
00120 
00121 Scheme_Object *scheme_complex_real_part(const Scheme_Object *n)
00122 {
00123   return ((Scheme_Complex *)n)->r;
00124 }
00125 
00126 Scheme_Object *scheme_complex_imaginary_part(const Scheme_Object *n)
00127 {
00128   return ((Scheme_Complex *)n)->i;
00129 }
00130 
00131 int scheme_complex_eq(const Scheme_Object *a, const Scheme_Object *b)
00132 {
00133   Scheme_Complex *ca = (Scheme_Complex *)a;
00134   Scheme_Complex *cb = (Scheme_Complex *)b;
00135   return scheme_bin_eq(ca->r, cb->r) && scheme_bin_eq(ca->i, cb->i);
00136 }
00137 
00138 Scheme_Object *scheme_complex_negate(const Scheme_Object *o)
00139 {
00140   Scheme_Complex *c = (Scheme_Complex *)o;
00141 
00142   return make_complex(scheme_bin_minus(scheme_make_integer(0),
00143                                    c->r), 
00144                     scheme_bin_minus(scheme_make_integer(0),
00145                                    c->i),
00146                     0);
00147 }
00148 
00149 Scheme_Object *scheme_complex_add(const Scheme_Object *a, const Scheme_Object *b)
00150 {
00151   Scheme_Complex *ca = (Scheme_Complex *)a;
00152   Scheme_Complex *cb = (Scheme_Complex *)b;
00153 
00154   return scheme_make_complex(scheme_bin_plus(ca->r, cb->r),
00155                           scheme_bin_plus(ca->i, cb->i));
00156 }
00157 
00158 Scheme_Object *scheme_complex_subtract(const Scheme_Object *a, const Scheme_Object *b)
00159 {
00160   Scheme_Complex *ca = (Scheme_Complex *)a;
00161   Scheme_Complex *cb = (Scheme_Complex *)b;
00162 
00163   return scheme_make_complex(scheme_bin_minus(ca->r, cb->r),
00164                           scheme_bin_minus(ca->i, cb->i));
00165 }
00166 
00167 Scheme_Object *scheme_complex_add1(const Scheme_Object *n)
00168 {
00169   Small_Complex s;
00170 
00171   return scheme_complex_add(scheme_make_small_complex(scheme_make_integer(1), &s), 
00172                          n);
00173 }
00174 
00175 Scheme_Object *scheme_complex_sub1(const Scheme_Object *n)
00176 {
00177   Small_Complex s;
00178 
00179   return scheme_complex_add(n, scheme_make_small_complex(scheme_make_integer(-1), 
00180                                                   &s));
00181 }
00182 
00183 Scheme_Object *scheme_complex_multiply(const Scheme_Object *a, const Scheme_Object *b)
00184 {
00185   Scheme_Complex *ca = (Scheme_Complex *)a;
00186   Scheme_Complex *cb = (Scheme_Complex *)b;
00187 
00188   return scheme_make_complex(scheme_bin_minus(scheme_bin_mult(ca->r, cb->r),
00189                                          scheme_bin_mult(ca->i, cb->i)),
00190                           scheme_bin_plus(scheme_bin_mult(ca->r, cb->i),
00191                                         scheme_bin_mult(ca->i, cb->r)));
00192   
00193 }
00194 
00195 Scheme_Object *scheme_complex_divide(const Scheme_Object *_n, const Scheme_Object *_d)
00196 { 
00197   Scheme_Complex *cn = (Scheme_Complex *)_n;
00198   Scheme_Complex *cd = (Scheme_Complex *)_d;
00199   Scheme_Object *den, *r, *i, *a, *b, *c, *d, *cm, *dm, *aa[1];
00200   int swap;
00201   
00202   if ((cn->r == zero) && (cn->i == zero))
00203     return zero;
00204 
00205   a = cn->r;
00206   b = cn->i;
00207   c = cd->r;
00208   d = cd->i;
00209 
00210   /* Check for exact-zero simplifications in d: */
00211   if (c == zero) {
00212     i = scheme_bin_minus(zero, scheme_bin_div(a, d));
00213     r = scheme_bin_div(b, d);
00214     return scheme_make_complex(r, i);
00215   } else if (d == zero) {
00216     r = scheme_bin_div(a, c);
00217     i = scheme_bin_div(b, c);
00218     return scheme_make_complex(r, i);
00219   }
00220 
00221   if (!SCHEME_FLOATP(c) && !SCHEME_FLOATP(d)) {
00222     /* The simple way: */
00223     cm = scheme_bin_plus(scheme_bin_mult(c, c), 
00224                          scheme_bin_mult(d, d));
00225     
00226     r = scheme_bin_div(scheme_bin_plus(scheme_bin_mult(c, a),
00227                                        scheme_bin_mult(d, b)),
00228                        cm);
00229     i = scheme_bin_div(scheme_bin_minus(scheme_bin_mult(c, b),
00230                                         scheme_bin_mult(d, a)),
00231                        cm);
00232     
00233     return scheme_make_complex(r, i);
00234   }
00235 
00236   if (scheme_is_zero(d)) {
00237     /* This is like dividing by a real number, except that
00238        the inexact 0 imaginary part can interact with +inf.0 and +nan.0 */
00239     r = scheme_bin_plus(scheme_bin_div(a, c),
00240                      /* Either 0.0 or +nan.0: */
00241                      scheme_bin_mult(d, b));
00242     i = scheme_bin_minus(scheme_bin_div(b, c),
00243                       /* Either 0.0 or +nan.0: */
00244                       scheme_bin_mult(d, a));
00245     
00246     return scheme_make_complex(r, i);
00247   }
00248   if (scheme_is_zero(c)) {
00249     r = scheme_bin_plus(scheme_bin_div(b, d),
00250                      /* Either 0.0 or +nan.0: */
00251                      scheme_bin_mult(c, a));
00252     i = scheme_bin_minus(scheme_bin_mult(c, b),  /* either 0.0 or +nan.0 */
00253                       scheme_bin_div(a, d));
00254 
00255     return scheme_make_complex(r, i);
00256   }
00257 
00258   aa[0] = c;
00259   cm = scheme_abs(1, aa);
00260   aa[0] = d;
00261   dm = scheme_abs(1, aa);
00262 
00263   if (scheme_bin_lt(cm, dm)) {
00264     cm = a;
00265     a = b;
00266     b = cm;
00267     cm = c;
00268     c = d;
00269     d = cm;
00270     swap = 1;
00271   } else
00272     swap = 0;
00273 
00274   r = scheme_bin_div(c, d);
00275 
00276   den = scheme_bin_plus(d, scheme_bin_mult(c, r));
00277 
00278   if (swap)
00279     i = scheme_bin_div(scheme_bin_minus(a, scheme_bin_mult(b, r)), den);
00280   else
00281     i = scheme_bin_div(scheme_bin_minus(scheme_bin_mult(b, r), a), den);
00282 
00283   r = scheme_bin_div(scheme_bin_plus(b, scheme_bin_mult(a, r)), den);
00284 
00285   return scheme_make_complex(r, i);
00286 }
00287 
00288 Scheme_Object *scheme_complex_power(const Scheme_Object *base, const Scheme_Object *exponent)
00289 {
00290   Scheme_Complex *cb = (Scheme_Complex *)base;
00291   Scheme_Complex *ce = (Scheme_Complex *)exponent;
00292   double a, b, c, d, bm, ba, nm, na, r1, r2;
00293 
00294   if ((ce->i == zero) && !SCHEME_FLOATP(ce->r)) {
00295     if (SCHEME_INTP(ce->r) || SCHEME_BIGNUMP(ce->r))
00296       return scheme_generic_integer_power(base, ce->r);
00297   }
00298 
00299   a = scheme_get_val_as_double(cb->r);
00300   b = scheme_get_val_as_double(cb->i);
00301   c = scheme_get_val_as_double(ce->r);
00302   d = scheme_get_val_as_double(ce->i);
00303 
00304   bm = sqrt(a * a + b * b);
00305   ba = atan2(b, a);
00306 
00307   /* New mag & angle */
00308   nm = pow(bm, c) * exp(-(ba * d));
00309   na = log(bm) * d + ba * c;
00310 
00311   r1 = nm * cos(na);
00312   r2 = nm * sin(na);
00313 
00314 #ifdef MZ_USE_SINGLE_FLOATS
00315   /* Coerce to double or float? */
00316 #ifdef USE_SINGLE_FLOATS_AS_DEFAULT
00317   if (!SCHEME_DBLP(cb->r) && !SCHEME_DBLP(cb->i)
00318       && !SCHEME_DBLP(ce->r) && !SCHEME_DBLP(ce->i))
00319 #else
00320   if (SCHEME_FLTP(cb->r) && SCHEME_FLTP(cb->i)
00321       && SCHEME_FLTP(ce->r) && SCHEME_FLTP(ce->i))
00322 #endif
00323     return scheme_make_complex(scheme_make_float((float)r1), 
00324                             scheme_make_float((float)r2));
00325 #endif
00326 
00327   return scheme_make_complex(scheme_make_double(r1), 
00328                           scheme_make_double(r2));
00329 }
00330 
00331 Scheme_Object *scheme_complex_sqrt(const Scheme_Object *o)
00332 {
00333   Scheme_Complex *c = (Scheme_Complex *)o;
00334   Scheme_Object *r, *i, *ssq, *srssq, *nrsq, *prsq, *nr, *ni;
00335 
00336   r = c->r;
00337   i = c->i;
00338 
00339   if (scheme_is_zero(i)) {
00340     /* Special case for x+0.0i: */
00341     r = scheme_sqrt(1, &r);
00342     if (!SCHEME_COMPLEXP(r))
00343       return scheme_make_complex(r, i);
00344     else
00345       return r;
00346   }
00347 
00348   ssq = scheme_bin_plus(scheme_bin_mult(r, r),
00349                      scheme_bin_mult(i, i));
00350 
00351   srssq = scheme_sqrt(1, &ssq);
00352 
00353   if (SCHEME_FLOATP(srssq)) {
00354     /* We may have lost too much precision, if i << r.  The result is
00355        going to be inexact, anyway, so switch to using expt. */
00356     Scheme_Object *a[2];
00357     a[0] = (Scheme_Object *)o;
00358     a[1] = scheme_make_double(0.5);
00359     return scheme_expt(2, a);
00360   }
00361 
00362   nrsq = scheme_bin_div(scheme_bin_minus(srssq, r),
00363                      scheme_make_integer(2));
00364 
00365   nr = scheme_sqrt(1, &nrsq);
00366   if (scheme_is_negative(i))
00367     nr = scheme_bin_minus(zero, nr);
00368     
00369   prsq = scheme_bin_div(scheme_bin_plus(srssq, r),
00370                      scheme_make_integer(2));
00371 
00372   ni = scheme_sqrt(1, &prsq);
00373 
00374   return scheme_make_complex(ni, nr);
00375 }