3 files changed, 298 insertions, 0 deletions
diff --git a/mrbgems/mruby-cmath/mrbgem.rake b/mrbgems/mruby-cmath/mrbgem.rake
new file mode 100644
index 000000000..e00725fef
--- /dev/null
+++ b/mrbgems/mruby-cmath/mrbgem.rake
@@ -0,0 +1,8 @@
+# This `mruby-cmath` gem uses C99 _Complex features
+# You need C compler that support C99+
+MRuby::Gem::Specification.new('mruby-cmath') do |spec|
+  spec.license = 'MIT'
+  spec.author  = 'mruby developers'
+  spec.summary = 'standard Math module with complex'
+  spec.add_dependency 'mruby-complex'
+end
diff --git a/mrbgems/mruby-cmath/src/cmath.c b/mrbgems/mruby-cmath/src/cmath.c
new file mode 100644
index 000000000..14e8ba670
--- /dev/null
+++ b/mrbgems/mruby-cmath/src/cmath.c
@@ -0,0 +1,249 @@
+/*
+** cmath.c - Math module with complex numbers
+**
+** See Copyright Notice in mruby.h
+*/
+
+/*
+** This `mruby-cmath` gem uses C99 _Complex features
+** You need C compler that support C99+
+*/
+
+#include <mruby.h>
+
+#ifdef MRB_NO_FLOAT
+# error CMath conflicts with 'MRB_NO_FLOAT' configuration
+#endif
+
+#include <math.h>
+#include <complex.h>
+
+mrb_value mrb_complex_new(mrb_state *mrb, mrb_float real, mrb_float imag);
+void mrb_complex_get(mrb_state *mrb, mrb_value cpx, mrb_float*, mrb_float*);
+
+static mrb_bool
+cmath_get_complex(mrb_state *mrb, mrb_value c, mrb_float *r, mrb_float *i)
+{
+  if (mrb_integer_p(c)) {
+    *r = (mrb_float)mrb_integer(c);
+    *i = 0;
+    return FALSE;
+  }
+  else if (mrb_float_p(c)) {
+    *r = mrb_float(c);
+    *i = 0;
+    return FALSE;
+  }
+  else if (mrb_obj_is_kind_of(mrb, c, mrb_class_get(mrb, "Complex"))) {
+    mrb_complex_get(mrb, c, r, i);
+    return TRUE;
+  }
+  else {
+    mrb_raise(mrb, E_TYPE_ERROR, "Numeric required");
+    return FALSE;
+  }
+}
+
+#ifdef MRB_USE_FLOAT32
+#define F(x) x##f
+#else
+#endif
+#define F(x) x
+
+#ifdef _WIN32
+
+#ifdef MRB_USE_FLOAT32
+typedef _Fcomplex mrb_complex;
+#define CX(r,i) _FCbuild(r,i)
+#else
+typedef _Dcomplex mrb_complex;
+#define CX(r,i) _Cbuild(r,i)
+#endif
+
+static mrb_complex
+CXDIVf(mrb_complex x, mrb_float y)
+{
+  return CX(creal(x)/y, cimag(x)/y);
+}
+
+static mrb_complex
+CXDIVc(mrb_complex a, mrb_complex b)
+{
+  mrb_float ratio, den;
+  mrb_float abr, abi, cr, ci;
+
+  if ((abr = creal(b)) < 0)
+    abr = - abr;
+  if ((abi = cimag(b)) < 0)
+    abi = - abi;
+  if (abr <= abi) {
+    ratio = creal(b) / cimag(b) ;
+    den = cimag(a) * (1 + ratio*ratio);
+    cr = (creal(a)*ratio + cimag(a)) / den;
+    ci = (cimag(a)*ratio - creal(a)) / den;
+  }
+  else {
+    ratio = cimag(b) / creal(b) ;
+    den = creal(a) * (1 + ratio*ratio);
+    cr = (creal(a) + cimag(a)*ratio) / den;
+    ci = (cimag(a) - creal(a)*ratio) / den;
+  }
+  return CX(cr, ci);
+}
+
+#else
+
+#ifdef MRB_USE_FLOAT32
+typedef float _Complex mrb_complex;
+#else
+typedef double _Complex mrb_complex;
+#endif
+    
+#define CX(r,i) (r+i*I)
+#define CXDIVf(x,y) (x)/(y)
+#define CXDIVc(x,y) (x)/(y)
+  
+#endif
+
+#define DEF_CMATH_METHOD(name) \
+static mrb_value \
+cmath_ ## name(mrb_state *mrb, mrb_value self)\
+{\
+  mrb_value z = mrb_get_arg1(mrb);\
+  mrb_float real, imag;\
+  if (cmath_get_complex(mrb, z, &real, &imag)) {\
+    mrb_complex c = CX(real,imag);\
+    c = F(c ## name)(c);\
+    return mrb_complex_new(mrb, creal(c), cimag(c));\
+  }\
+  return mrb_float_value(mrb, F(name)(real));\
+}
+
+/* exp(z): return the exponential of z */
+DEF_CMATH_METHOD(exp)
+
+/* log(z): return the natural logarithm of z, with branch cut along the negative real axis */
+static mrb_value
+cmath_log(mrb_state *mrb, mrb_value self) {
+  mrb_value z;
+  mrb_float base;
+  mrb_float real, imag;
+
+  mrb_int n = mrb_get_args(mrb, "o|f", &z, &base);
+
+#ifndef M_E
+#define M_E F(exp)(1.0)
+#endif
+
+  if (n == 1) base = M_E;
+  if (cmath_get_complex(mrb, z, &real, &imag) || real < 0.0) {
+    mrb_complex c = CX(real,imag);
+    c = F(clog)(c);
+    if (n == 2) c = CXDIVc(c, F(clog)(base));
+    return mrb_complex_new(mrb, creal(c), cimag(c));
+  }
+  if (n == 1) return mrb_float_value(mrb, F(log)(real));
+  return mrb_float_value(mrb, F(log)(real)/F(log)(base));
+}
+
+/* log10(z): return the base-10 logarithm of z, with branch cut along the negative real axis */
+static mrb_value
+cmath_log10(mrb_state *mrb, mrb_value self) {
+  mrb_value z = mrb_get_arg1(mrb);
+  mrb_float real, imag;
+  if (cmath_get_complex(mrb, z, &real, &imag) || real < 0.0) {
+    mrb_complex c = CX(real,imag);
+    c = CXDIVf(F(clog)(c),log(10));
+    return mrb_complex_new(mrb, creal(c), cimag(c));
+  }
+  return mrb_float_value(mrb, F(log10)(real));
+}
+
+/* log2(z): return the base-2 logarithm of z, with branch cut along the negative real axis */
+static mrb_value
+cmath_log2(mrb_state *mrb, mrb_value self) {
+  mrb_value z = mrb_get_arg1(mrb);
+  mrb_float real, imag;
+  if (cmath_get_complex(mrb, z, &real, &imag) || real < 0.0) {
+    mrb_complex c = CX(real,imag);
+    c = CXDIVf(F(clog)(c),log(2));
+    return mrb_complex_new(mrb, creal(c), cimag(c));
+  }
+  return mrb_float_value(mrb, F(log2)(real));
+}
+
+/* sqrt(z): return square root of z */
+static mrb_value
+cmath_sqrt(mrb_state *mrb, mrb_value self) {
+  mrb_value z = mrb_get_arg1(mrb);
+  mrb_float real, imag;
+  if (cmath_get_complex(mrb, z, &real, &imag) || real < 0.0) {
+    mrb_complex c = CX(real,imag);
+    c = F(csqrt)(c);
+    return mrb_complex_new(mrb, creal(c), cimag(c));
+  }
+  return mrb_float_value(mrb, F(sqrt)(real));
+}
+
+/* sin(z): sine function */
+DEF_CMATH_METHOD(sin)
+/* cos(z): cosine function */
+DEF_CMATH_METHOD(cos)
+/* tan(z): tangent function */
+DEF_CMATH_METHOD(tan)
+/* asin(z): arc sine function */
+DEF_CMATH_METHOD(asin)
+/* acos(z): arc cosine function */
+DEF_CMATH_METHOD(acos)
+/* atan(z): arg tangent function */
+DEF_CMATH_METHOD(atan)
+/* sinh(z): hyperbolic sine function */
+DEF_CMATH_METHOD(sinh)
+/* cosh(z): hyperbolic cosine function */
+DEF_CMATH_METHOD(cosh)
+/* tanh(z): hyperbolic tangent function */
+DEF_CMATH_METHOD(tanh)
+/* asinh(z): inverse hyperbolic sine function */
+DEF_CMATH_METHOD(asinh)
+/* acosh(z): inverse hyperbolic cosine function */
+DEF_CMATH_METHOD(acosh)
+/* atanh(z): inverse hyperbolic tangent function */
+DEF_CMATH_METHOD(atanh)
+
+/* ------------------------------------------------------------------------*/
+
+void
+mrb_mruby_cmath_gem_init(mrb_state* mrb)
+{
+  struct RClass *cmath;
+  cmath = mrb_define_module(mrb, "CMath");
+
+  mrb_include_module(mrb, cmath, mrb_module_get(mrb, "Math"));
+
+  mrb_define_module_function(mrb, cmath, "sin", cmath_sin, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "cos", cmath_cos, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "tan", cmath_tan, MRB_ARGS_REQ(1));
+
+  mrb_define_module_function(mrb, cmath, "asin", cmath_asin, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "acos", cmath_acos, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "atan", cmath_atan, MRB_ARGS_REQ(1));
+
+  mrb_define_module_function(mrb, cmath, "sinh", cmath_sinh, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "cosh", cmath_cosh, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "tanh", cmath_tanh, MRB_ARGS_REQ(1));
+
+  mrb_define_module_function(mrb, cmath, "asinh", cmath_asinh, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "acosh", cmath_acosh, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "atanh", cmath_atanh, MRB_ARGS_REQ(1));
+
+  mrb_define_module_function(mrb, cmath, "exp", cmath_exp, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "log", cmath_log, MRB_ARGS_REQ(1)|MRB_ARGS_OPT(1));
+  mrb_define_module_function(mrb, cmath, "log2", cmath_log2, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "log10", cmath_log10, MRB_ARGS_REQ(1));
+  mrb_define_module_function(mrb, cmath, "sqrt", cmath_sqrt, MRB_ARGS_REQ(1));
+}
+
+void
+mrb_mruby_cmath_gem_final(mrb_state* mrb)
+{
+}
diff --git a/mrbgems/mruby-cmath/test/cmath.rb b/mrbgems/mruby-cmath/test/cmath.rb
new file mode 100644
index 000000000..ddd0e5106
--- /dev/null
+++ b/mrbgems/mruby-cmath/test/cmath.rb
@@ -0,0 +1,41 @@
+##
+# CMath Test
+
+def assert_complex(c1, c2)
+  assert('assert_complex') do
+    assert_float(c1.real, c2.real)
+    assert_float(c1.imaginary, c2.imaginary)
+  end
+end
+
+assert('CMath.exp') do
+  assert_float(1.0, CMath.exp(0))
+  assert_complex(-1+0i, CMath.exp(Math::PI.i))
+  assert_complex((-1.1312043837568135+2.4717266720048188i), CMath.exp(1+2i))
+end
+
+assert('CMath.log') do
+  assert_float(0, CMath.log(1))
+  assert_float(3.0, CMath.log(8,2))
+  assert_complex((1.092840647090816-0.42078724841586035i), CMath.log(-8,-2))
+end
+
+assert('CMath.sqrt') do
+  assert_complex(Complex(0,2), CMath.sqrt(-4.0))
+  assert_complex(Complex(0,3), CMath.sqrt(-9.0))
+end
+
+assert('CMath trigonometric_functions') do
+  assert_complex(Math.sinh(2).i, CMath.sin(2i))
+  assert_complex(Math.cosh(2)+0i,   CMath.cos(2i))
+  assert_complex(Math.tanh(2).i, CMath.tan(2i))
+  assert_complex(Math.sin(2).i, CMath.sinh(2i))
+  assert_complex(Math.cos(2)+0i, CMath.cosh(2i))
+  assert_complex(Math.tan(2).i, CMath.tanh(2i))
+  assert_complex(1+1i, CMath.sin(CMath.asin(1+1i)))
+  assert_complex(1+1i, CMath.cos(CMath.acos(1+1i)))
+  assert_complex(1+1i, CMath.tan(CMath.atan(1+1i)))
+  assert_complex(1+1i, CMath.sinh(CMath.asinh(1+1i)))
+  assert_complex(1+1i, CMath.cosh(CMath.acosh(1+1i)))
+  assert_complex(1+1i, CMath.tanh(CMath.atanh(1+1i)))
+end