third-party/mruby/mrbgems/mruby-complex/src/complex.c

   1 #include <mruby.h>
   2 #include <mruby/class.h>
   3 #include <mruby/numeric.h>
   4 #include <math.h>
   5
   6 #ifdef MRB_WITHOUT_FLOAT
   7 # error Complex conflicts 'MRB_WITHOUT_FLOAT' configuration in your 'build_config.rb'
   8 #endif
   9
  10 struct mrb_complex {
  11   mrb_float real;
  12   mrb_float imaginary;
  13 };
  14
  15 #ifdef MRB_USE_FLOAT
  16 #define F(x) x##f
  17 #else
  18 #define F(x) x
  19 #endif
  20
  21 #if defined(MRB_64BIT) || defined(MRB_USE_FLOAT)
  22
  23 #define COMPLEX_USE_ISTRUCT
  24 /* use TT_ISTRUCT */
  25 #include <mruby/istruct.h>
  26
  27 #define complex_ptr(mrb, v) (struct mrb_complex*)mrb_istruct_ptr(v)
  28
  29 static struct RBasic*
  30 complex_alloc(mrb_state *mrb, struct RClass *c, struct mrb_complex **p)
  31 {
  32   struct RIStruct *s;
  33
  34   s = (struct RIStruct*)mrb_obj_alloc(mrb, MRB_TT_ISTRUCT, c);
  35   *p = (struct mrb_complex*)s->inline_data;
  36
  37   return (struct RBasic*)s;
  38 }
  39
  40 #else
  41 /* use TT_DATA */
  42 #include <mruby/data.h>
  43
  44 static const struct mrb_data_type mrb_complex_type = {"Complex", mrb_free};
  45
  46 static struct RBasic*
  47 complex_alloc(mrb_state *mrb, struct RClass *c, struct mrb_complex **p)
  48 {
  49   struct RData *d;
  50
  51   Data_Make_Struct(mrb, c, struct mrb_complex, &mrb_complex_type, *p, d);
  52
  53   return (struct RBasic*)d;
  54 }
  55
  56 static struct mrb_complex*
  57 complex_ptr(mrb_state *mrb, mrb_value v)
  58 {
  59   struct mrb_complex *p;
  60
  61   p = DATA_GET_PTR(mrb, v, &mrb_complex_type, struct mrb_complex);
  62   if (!p) {
  63     mrb_raise(mrb, E_ARGUMENT_ERROR, "uninitialized complex");
  64   }
  65   return p;
  66 }
  67 #endif
  68
  69 static mrb_value
  70 complex_new(mrb_state *mrb, mrb_float real, mrb_float imaginary)
  71 {
  72   struct RClass *c = mrb_class_get(mrb, "Complex");
  73   struct mrb_complex *p;
  74   struct RBasic *comp = complex_alloc(mrb, c, &p);
  75   p->real = real;
  76   p->imaginary = imaginary;
  77   MRB_SET_FROZEN_FLAG(comp);
  78
  79   return mrb_obj_value(comp);
  80 }
  81
  82 static mrb_value
  83 complex_real(mrb_state *mrb, mrb_value self)
  84 {
  85   struct mrb_complex *p = complex_ptr(mrb, self);
  86   return mrb_float_value(mrb, p->real);
  87 }
  88
  89 static mrb_value
  90 complex_imaginary(mrb_state *mrb, mrb_value self)
  91 {
  92   struct mrb_complex *p = complex_ptr(mrb, self);
  93   return mrb_float_value(mrb, p->imaginary);
  94 }
  95
  96 static mrb_value
  97 complex_s_rect(mrb_state *mrb, mrb_value self)
  98 {
  99   mrb_float real, imaginary = 0.0;
 100
 101   mrb_get_args(mrb, "f|f", &real, &imaginary);
 102   return complex_new(mrb, real, imaginary);
 103 }
 104
 105 static mrb_value
 106 complex_to_f(mrb_state *mrb, mrb_value self)
 107 {
 108   struct mrb_complex *p = complex_ptr(mrb, self);
 109
 110   if (p->imaginary != 0) {
 111     mrb_raisef(mrb, E_RANGE_ERROR, "can't convert %v into Float", self);
 112   }
 113
 114   return mrb_float_value(mrb, p->real);
 115 }
 116
 117 static mrb_value
 118 complex_to_i(mrb_state *mrb, mrb_value self)
 119 {
 120   struct mrb_complex *p = complex_ptr(mrb, self);
 121
 122   if (p->imaginary != 0) {
 123     mrb_raisef(mrb, E_RANGE_ERROR, "can't convert %v into Float", self);
 124   }
 125   return mrb_int_value(mrb, p->real);
 126 }
 127
 128 static mrb_value
 129 complex_to_c(mrb_state *mrb, mrb_value self)
 130 {
 131   return self;
 132 }
 133
 134 /* Arithmetic on (significand, exponent) pairs avoids premature overflow in
 135    complex division */
 136 struct float_pair {
 137   mrb_float s;
 138   int x;
 139 };
 140
 141 static void
 142 add_pair(struct float_pair *s, struct float_pair const *a,
 143          struct float_pair const *b)
 144 {
 145   if (b->s == 0.0F) {
 146     *s = *a;
 147   } else if (a->s == 0.0F) {
 148     *s = *b;
 149   } else if (a->x >= b->x) {
 150     s->s = a->s + F(ldexp)(b->s, b->x - a->x);
 151     s->x = a->x;
 152   } else {
 153     s->s = F(ldexp)(a->s, a->x - b->x) + b->s;
 154     s->x = b->x;
 155   }
 156 }
 157
 158 static void
 159 mul_pair(struct float_pair *p, struct float_pair const *a,
 160          struct float_pair const *b)
 161 {
 162   p->s = a->s * b->s;
 163   p->x = a->x + b->x;
 164 }
 165
 166 static void
 167 div_pair(struct float_pair *q, struct float_pair const *a,
 168          struct float_pair const *b)
 169 {
 170   q->s = a->s / b->s;
 171   q->x = a->x - b->x;
 172 }
 173
 174 static mrb_value
 175 complex_div(mrb_state *mrb, mrb_value self)
 176 {
 177   mrb_value rhs = mrb_get_arg1(mrb);
 178   struct mrb_complex *a, *b;
 179   struct float_pair ar, ai, br, bi;
 180   struct float_pair br2, bi2;
 181   struct float_pair div;
 182   struct float_pair ar_br, ai_bi;
 183   struct float_pair ai_br, ar_bi;
 184   struct float_pair zr, zi;
 185
 186   a = complex_ptr(mrb, self);
 187   b = complex_ptr(mrb, rhs);
 188
 189   /* Split floating point components into significand and exponent */
 190   ar.s = F(frexp)(a->real, &ar.x);
 191   ai.s = F(frexp)(a->imaginary, &ai.x);
 192   br.s = F(frexp)(b->real, &br.x);
 193   bi.s = F(frexp)(b->imaginary, &bi.x);
 194
 195   /* Perform arithmetic on (significand, exponent) pairs to produce
 196      the result: */
 197
 198   /* the divisor */
 199   mul_pair(&br2, &br, &br);
 200   mul_pair(&bi2, &bi, &bi);
 201   add_pair(&div, &br2, &bi2);
 202
 203   /* real component */
 204   mul_pair(&ar_br, &ar, &br);
 205   mul_pair(&ai_bi, &ai, &bi);
 206   add_pair(&zr, &ar_br, &ai_bi);
 207   div_pair(&zr, &zr, &div);
 208
 209   /* imaginary component */
 210   mul_pair(&ai_br, &ai, &br);
 211   mul_pair(&ar_bi, &ar, &bi);
 212   ar_bi.s = -ar_bi.s;
 213   add_pair(&zi, &ai_br, &ar_bi);
 214   div_pair(&zi, &zi, &div);
 215
 216   /* assemble the result */
 217   return complex_new(mrb, F(ldexp)(zr.s, zr.x), F(ldexp)(zi.s, zi.x));
 218 }
 219
 220 void mrb_mruby_complex_gem_init(mrb_state *mrb)
 221 {
 222   struct RClass *comp;
 223
 224 #ifdef COMPLEX_USE_ISTRUCT
 225   mrb_assert(sizeof(struct mrb_complex) < ISTRUCT_DATA_SIZE);
 226 #endif
 227   comp = mrb_define_class(mrb, "Complex", mrb_class_get(mrb, "Numeric"));
 228 #ifdef COMPLEX_USE_ISTRUCT
 229   MRB_SET_INSTANCE_TT(comp, MRB_TT_ISTRUCT);
 230 #else
 231   MRB_SET_INSTANCE_TT(comp, MRB_TT_DATA);
 232 #endif
 233   mrb_undef_class_method(mrb, comp, "new");
 234   mrb_define_class_method(mrb, comp, "rectangular", complex_s_rect, MRB_ARGS_REQ(1)|MRB_ARGS_OPT(1));
 235   mrb_define_class_method(mrb, comp, "rect", complex_s_rect, MRB_ARGS_REQ(1)|MRB_ARGS_OPT(1));
 236   mrb_define_method(mrb, mrb->kernel_module, "Complex", complex_s_rect, MRB_ARGS_REQ(1)|MRB_ARGS_OPT(1));
 237   mrb_define_method(mrb, comp, "real", complex_real, MRB_ARGS_NONE());
 238   mrb_define_method(mrb, comp, "imaginary", complex_imaginary, MRB_ARGS_NONE());
 239   mrb_define_method(mrb, comp, "to_f", complex_to_f, MRB_ARGS_NONE());
 240   mrb_define_method(mrb, comp, "to_i", complex_to_i, MRB_ARGS_NONE());
 241   mrb_define_method(mrb, comp, "to_c", complex_to_c, MRB_ARGS_NONE());
 242   mrb_define_method(mrb, comp, "__div__", complex_div, MRB_ARGS_REQ(1));
 243 }
 244
 245 void
 246 mrb_mruby_complex_gem_final(mrb_state* mrb)
 247 {
 248 }