/*
 *  basicnumarray.c -
 *
 *  Author: gotoken <gotoken@notwork.org>
 *  $Id: basicnumarray.c,v 1.3 2000/04/14 04:50:28 gotoken Exp $
 *
 */

#include <math.h>
#include "ruby.h"

VALUE cBasicNumArray;
VALUE cBasicNumArrayDouble;
static ID icoerce;
static ID iadd;
static ID isub;
static ID imul;
static ID idiv;
static ID ipow;
static ID imod;

enum eType {
    tDOUBLE
};

enum eOpt {
    oADD, 
    oSUB, 
    oMUL, 
    oDIV,
    oPOW,
    oMOD,

    oACOS,
    oACOSH,
    oASIN,
    oASINH,
    oATAN,
    oATANH,
    oCBRT,
    oCEIL,
    oCOPYSIGN,
    oCOS,
    oCOSH,
    oERF,
    oERFC,
    oEXP,
    oEXPM1,
    oFABS,
    oFLOOR,
    oILOGB,
    oJ0,
    oJ1,
    oLGAMMA,
    oLOG,
    oLOG10,
    oLOG1P,
    oRINT,
    oSIN,
    oSINH,
    oSQRT,
    oTAN,
    oTANH,
    oY0,
    oY1
};

typedef enum eType etype;
typedef enum eOpt  eopt;

struct BasicNumArray {
    etype  type;
    size_t len;
    void   *ptr;
};

#define BNA struct BasicNumArray

BNA * alloc_bna(etype type, size_t len);
void free_bna(BNA *ary);
VALUE bna_s_new(VALUE mod, VALUE klass, VALUE size);
VALUE na_trim(VALUE obj);
VALUE na_to_basic(VALUE obj);
VALUE bna_setall(VALUE obj, VALUE x);
static VALUE bna_scaler(VALUE x, VALUE y);
VALUE bna_len(VALUE obj);
VALUE bna_aref0(VALUE obj, VALUE idx);
VALUE bna_aref(VALUE obj, VALUE idx);
VALUE bna_aset0(VALUE obj, VALUE idx, VALUE x);
VALUE bna_aset(VALUE obj, VALUE idx, VALUE x);
static etype bna_etype(VALUE obj);
static int bna_is_a(VALUE x, VALUE c);
VALUE bna_coerce(VALUE x, VALUE y);
static void do_coerce(VALUE *x ,VALUE *y);
VALUE bna_binop(VALUE x, VALUE y, ID op, eopt eop);
VALUE bna_add(VALUE x, VALUE y);
VALUE bna_sub(VALUE x, VALUE y);
VALUE bna_mul(VALUE x, VALUE y);
VALUE bna_div(VALUE x, VALUE y);
VALUE bna_pow(VALUE x, VALUE y);
VALUE bna_mod(VALUE x, VALUE y);
VALUE unary_op(VALUE x, eopt op);
VALUE bna_acos(VALUE x);
VALUE bna_acosh(VALUE x);
VALUE bna_asin(VALUE x);
VALUE bna_asinh(VALUE x);
VALUE bna_atan(VALUE x);
VALUE bna_atanh(VALUE x);
VALUE bna_cbrt(VALUE x);
VALUE bna_ceil(VALUE x);
VALUE bna_copysign(VALUE x);
VALUE bna_cos(VALUE x);
VALUE bna_cosh(VALUE x);
VALUE bna_erf(VALUE x);
VALUE bna_erfc(VALUE x);
VALUE bna_exp(VALUE x);
VALUE bna_expm1(VALUE x);
VALUE bna_fabs(VALUE x);
VALUE bna_floor(VALUE x);
VALUE bna_j0(VALUE x);
VALUE bna_j1(VALUE x);
VALUE bna_lgamma(VALUE x);
VALUE bna_log(VALUE x);
VALUE bna_log10(VALUE x);
VALUE bna_1p(VALUE x);
VALUE bna_rint(VALUE x);
VALUE bna_sin(VALUE x);
VALUE bna_sinh(VALUE x);
VALUE bna_sqrt(VALUE x);
VALUE bna_tan(VALUE x);
VALUE bna_tanh(VALUE x);
VALUE bna_y0(VALUE x);
VALUE bna_y1(VALUE x);
VALUE m_cos(VALUE obj, VALUE x); 
VALUE m_sin(VALUE obj, VALUE x);
VALUE m_tan(VALUE obj, VALUE x);
VALUE m_exp(VALUE obj, VALUE x);
VALUE m_log(VALUE obj, VALUE x);
VALUE m_log10(VALUE obj, VALUE x);
VALUE m_sqrt(VALUE obj, VALUE x);
static VALUE bna_inspect_ary(VALUE obj);
static VALUE bna_ary_inspect(VALUE obj);
void Init_basicnumarray(void);

BNA *
alloc_bna(etype type, size_t len)
{
    BNA *ary;
    ary = xmalloc(sizeof(BNA));

    ary->type = type;
    ary->len = len;

    switch(type) {
      case tDOUBLE:
	ary->ptr = (double *)xmalloc(len*sizeof(double));
    }

    return ary;
}

void
free_bna(BNA *ary)
{
    free(ary->ptr);
    free(ary);
}

VALUE
bna_s_new(VALUE mod, VALUE klass, VALUE size)
{
    BNA    *ary;
    VALUE   obj;
    size_t  len;

    len = NUM2ULONG(size);

    if (klass == rb_cFloat) {
	ary = alloc_bna(tDOUBLE, len);
	return Data_Wrap_Struct(cBasicNumArrayDouble, 0, free_bna, ary);
    } else {
	rb_raise(rb_eNotImpError, "non supported class");
    }
}

VALUE
bna_setall(VALUE obj, VALUE x)
{
    BNA *ary;
    long i;
    double d, *dp;

    Data_Get_Struct(obj, BNA, ary);

    switch (ary->type) {
      case tDOUBLE:
	d = RFLOAT(rb_Float(x))->value;
	dp = ary->ptr;
	for (i = 0; i < ary->len; i++) {
	    dp[i] = d;
	}
	return obj;
    }
}

static VALUE
bna_scaler(VALUE x, VALUE y)
{
    VALUE obj;
    BNA *ary0, *ary1;

    Data_Get_Struct(x, BNA, ary0);
    ary1 = alloc_bna(ary0->type, ary0->len);
    switch (ary0->type) {
      case tDOUBLE:
	obj = Data_Wrap_Struct(cBasicNumArrayDouble, 0, free_bna, ary1);
	return bna_setall(obj, y);
    }
}

VALUE
bna_len(VALUE obj)
{
    BNA *ary;

    Data_Get_Struct(obj, BNA, ary);
    return INT2NUM(ary->len);
}


VALUE
bna_aref0(VALUE obj, VALUE idx)
{
    long pos;
    BNA *ary;

    double *dp;

    pos = NUM2INT(idx);
    Data_Get_Struct(obj, BNA, ary);

    switch(ary->type) {
      case tDOUBLE:
	dp = ary->ptr;
	return rb_float_new(dp[pos]);
    }
}

VALUE
bna_aref(VALUE obj, VALUE idx)
{
    if (TYPE(idx) == T_FIXNUM || TYPE(idx) == T_BIGNUM) {
	return bna_aref0(obj, idx);
    } else {
	rb_raise(rb_eArgError, "unsupported kind of index");
    }
}

VALUE
bna_aset0(VALUE obj, VALUE idx, VALUE x)
{
    long pos;
    BNA *ary;
    double *dp;

    pos = NUM2INT(idx);
    Data_Get_Struct(obj, BNA, ary);

    switch(ary->type) {
      case tDOUBLE:
	dp = ary->ptr;
	dp[pos] = RFLOAT(rb_Float(x))->value;
	return x;
    }
}

VALUE
bna_aset(VALUE obj, VALUE idx, VALUE x)
{
    VALUE val;

    if (TYPE(idx) == T_FIXNUM || TYPE(idx) == T_BIGNUM) {
	return bna_aset0(obj, idx, x);
    } else {
	rb_raise(rb_eArgError, "unsupported kind of index");
    }
}

static etype
bna_etype(VALUE obj)
{
    BNA *ary;
    Data_Get_Struct(obj, BNA, ary);
    return ary->type;
}

static int
bna_is_a(VALUE x, VALUE c)
{
    VALUE cl = CLASS_OF(x);

    while (cl) {
        if (cl == c || RCLASS(cl)->m_tbl == RCLASS(c)->m_tbl)
            return 1;
        cl = RCLASS(cl)->super;
    }
    return 0;
}

VALUE
bna_coerce(VALUE x, VALUE y)
{
    if (CLASS_OF(x) == CLASS_OF(y))
	return rb_assoc_new(x, y);
    if (bna_is_a(y, rb_cNumeric)) {
	return rb_assoc_new(bna_scaler(x, y), x);
    }
    rb_raise(rb_eTypeError, "cannot coerce");
}

static void
do_coerce(VALUE *x, VALUE *y)
{
    VALUE ary;

    ary = rb_funcall(*y, icoerce, 1, *x);
    if (TYPE(ary) != T_ARRAY || RARRAY(ary)->len != 2) {
	rb_raise(rb_eTypeError, "coerce must return [x, y]");
    }
    *x = RARRAY(ary)->ptr[0];
    *y = RARRAY(ary)->ptr[1];
}

VALUE
bna_binop(VALUE x, VALUE y, ID op, eopt eop)
{
    VALUE obj, res, tmp;
    BNA *ary0, *ary1, *ary;
    int i;
    double *dp0, *dp1, *dp;

    if (bna_is_a(y, rb_cNumeric)) {
	y = bna_scaler(x, y);
    } else if (!bna_is_a(y, cBasicNumArray)) {
	do_coerce(&x, &y);
	return rb_funcall(x, op, 1, y);
    }

    Data_Get_Struct(x, BNA, ary0);
    Data_Get_Struct(y, BNA, ary1);
    if (ary0->len != ary1->len)
	rb_raise(rb_eArgError, "size mismatch");
    
    ary = alloc_bna(ary0->type, ary0->len);
    dp0 = ary0->ptr;
    dp1 = ary1->ptr;
    dp = ary->ptr;
    switch (eop) {
      case oADD:
	for (i = 0; i < ary0->len; i++) dp[i] = dp0[i] + dp1[i];
	break;
      case oSUB:
	for (i = 0; i < ary0->len; i++) dp[i] = dp0[i] - dp1[i];
	break;
      case oMUL:
	for (i = 0; i < ary0->len; i++) dp[i] = dp0[i] * dp1[i];
	break;
      case oDIV:
	for (i = 0; i < ary0->len; i++) dp[i] = dp0[i] / dp1[i];
	break;
      case oPOW:
	for (i = 0; i < ary0->len; i++) dp[i] = pow(dp0[i], dp1[i]);
	break;
      case oMOD:
	for (i = 0; i < ary0->len; i++) dp[i] = fmod(dp0[i], dp1[i]);
	break;
      defalut:
	rb_raise(rb_eNameError, "[BNA BUG] what's happes?? (eop = %d)", eop);
    }
    return Data_Wrap_Struct(cBasicNumArrayDouble, 0, free_bna, ary);
}

VALUE
bna_add(VALUE x, VALUE y)
{
    return bna_binop(x, y, iadd, oADD);
}

VALUE
bna_sub(VALUE x, VALUE y)
{
    return bna_binop(x, y, isub, oSUB);
}

VALUE
bna_mul(VALUE x, VALUE y)
{
    return bna_binop(x, y, imul, oMUL);
}

VALUE
bna_div(VALUE x, VALUE y)
{
    return bna_binop(x, y, idiv, oDIV);
}

VALUE
bna_pow(VALUE x, VALUE y)
{
    return bna_binop(x, y, ipow, oPOW);
}

VALUE
bna_mod(VALUE x, VALUE y)
{
    return bna_binop(x, y, imod, oMOD);
}

VALUE
unary_op(VALUE x, eopt op)
{
    BNA *ary0, *ary;
    double *dp0, *dp;
    int i;

    Data_Get_Struct(x, BNA, ary0);

    if (ary0->type == tDOUBLE) {
	ary = alloc_bna(ary0->type, ary0->len);
	dp0 = ary0->ptr;
	dp = ary->ptr;
	
	switch (op) {
	  case oACOS:
	    for (i = 0; i < ary0->len; i++) dp[i] = acos(dp0[i]);
	    break;
	  case oACOSH:
	    for (i = 0; i < ary0->len; i++) dp[i] = acosh(dp0[i]);
	    break;
	  case oASIN:
	    for (i = 0; i < ary0->len; i++) dp[i] = asin(dp0[i]);
	    break;
	  case oASINH:
	    for (i = 0; i < ary0->len; i++) dp[i] = asinh(dp0[i]);
	    break;
	  case oATAN:
	    for (i = 0; i < ary0->len; i++) dp[i] = atan(dp0[i]);
	    break;
	  case oATANH:
	    for (i = 0; i < ary0->len; i++) dp[i] = atanh(dp0[i]);
	    break;
	  case oCBRT:
	    for (i = 0; i < ary0->len; i++) dp[i] = cbrt(dp0[i]);
	    break;
	  case oCEIL:
	    for (i = 0; i < ary0->len; i++) dp[i] = ceil(dp0[i]);
	    break;
	  case oCOS:
	    for (i = 0; i < ary0->len; i++) dp[i] = cos(dp0[i]);
	    break;
	  case oCOSH:
	    for (i = 0; i < ary0->len; i++) dp[i] = cosh(dp0[i]);
	    break;
	  case oERF:
	    for (i = 0; i < ary0->len; i++) dp[i] = erf(dp0[i]);
	    break;
	  case oERFC:
	    for (i = 0; i < ary0->len; i++) dp[i] = erfc(dp0[i]);
	    break;
	  case oEXP:
	    for (i = 0; i < ary0->len; i++) dp[i] = exp(dp0[i]);
	    break;
	  case oEXPM1:
	    for (i = 0; i < ary0->len; i++) dp[i] = expm1(dp0[i]);
	    break;
	  case oFABS:
	    for (i = 0; i < ary0->len; i++) dp[i] = fabs(dp0[i]);
	    break;
	  case oFLOOR:
	    for (i = 0; i < ary0->len; i++) dp[i] = floor(dp0[i]);
	    break;
	  case oJ0:
	    for (i = 0; i < ary0->len; i++) dp[i] = j0(dp0[i]);
	    break;
	  case oJ1:
	    for (i = 0; i < ary0->len; i++) dp[i] = j1(dp0[i]);
	    break;
	  case oLGAMMA:
	    for (i = 0; i < ary0->len; i++) dp[i] = lgamma(dp0[i]);
	    break;
	  case oLOG:
	    for (i = 0; i < ary0->len; i++) dp[i] = log(dp0[i]);
	    break;
	  case oLOG10:
	    for (i = 0; i < ary0->len; i++) dp[i] = log10(dp0[i]);
	    break;
	  case oLOG1P:
	    for (i = 0; i < ary0->len; i++) dp[i] = log1p(dp0[i]);
	    break;
	  case oRINT:
	    for (i = 0; i < ary0->len; i++) dp[i] = rint(dp0[i]);
	    break;
	  case oSIN:
	    for (i = 0; i < ary0->len; i++) dp[i] = sin(dp0[i]);
	    break;
	  case oSINH:
	    for (i = 0; i < ary0->len; i++) dp[i] = sinh(dp0[i]);
	    break;
	  case oSQRT:
	    for (i = 0; i < ary0->len; i++) dp[i] = sqrt(dp0[i]);
	    break;
	  case oTAN:
	    for (i = 0; i < ary0->len; i++) dp[i] = tan(dp0[i]);
	    break;
	  case oTANH:
	    for (i = 0; i < ary0->len; i++) dp[i] = tanh(dp0[i]);
	    break;
	  case oY0:
	    for (i = 0; i < ary0->len; i++) dp[i] = y0(dp0[i]);
	    break;
	  case oY1:
	    for (i = 0; i < ary0->len; i++) dp[i] = y1(dp0[i]);
	    break;
	  defalut:
	    rb_raise(rb_eNameError, "[BNA BUG] what's happes? (op = %d)", op);
	}
	return Data_Wrap_Struct(cBasicNumArrayDouble, 0, free_bna, ary);
    }
}

VALUE
bna_acos(VALUE x)
{
    return unary_op(x, oACOS);
}

VALUE
bna_acosh(VALUE x)
{
    return unary_op(x, oACOSH);
}

VALUE
bna_asin(VALUE x)
{
    return unary_op(x, oASIN);
}

VALUE
bna_asinh(VALUE x)
{
    return unary_op(x, oASINH);
}

VALUE
bna_atan(VALUE x)
{
    return unary_op(x, oATAN);
}

VALUE
bna_atanh(VALUE x)
{
    return unary_op(x, oATANH);
}

VALUE
bna_cbrt(VALUE x)
{
    return unary_op(x, oCBRT);
}

VALUE
bna_ceil(VALUE x)
{
    return unary_op(x, oCEIL);
}

VALUE
bna_copysign(VALUE x)
{
    return unary_op(x, oCOPYSIGN);
}

VALUE
bna_cos(VALUE x)
{
    return unary_op(x, oCOS);
}

VALUE
bna_cosh(VALUE x)
{
    return unary_op(x, oCOSH);
}

VALUE
bna_erf(VALUE x)
{
    return unary_op(x, oERF);
}

VALUE
bna_erfc(VALUE x)
{
    return unary_op(x, oERFC);
}

VALUE
bna_exp(VALUE x)
{
    return unary_op(x, oEXP);
}

VALUE
bna_expm1(VALUE x)
{
    return unary_op(x, oEXPM1);
}

VALUE
bna_fabs(VALUE x)
{
    return unary_op(x, oFABS);
}

VALUE
bna_floor(VALUE x)
{
    return unary_op(x, oFLOOR);
}

VALUE
bna_j0(VALUE x)
{
    return unary_op(x, oJ0);
}

VALUE
bna_j1(VALUE x)
{
    return unary_op(x, oJ1);
}

VALUE
bna_lgamma(VALUE x)
{
    return unary_op(x, oLGAMMA);
}

VALUE
bna_log(VALUE x)
{
    return unary_op(x, oLOG);
}

VALUE
bna_log10(VALUE x)
{
    return unary_op(x, oLOG10);
}

VALUE
bna_log1p(VALUE x)
{
    return unary_op(x, oLOG1P);
}

VALUE
bna_rint(VALUE x)
{
    return unary_op(x, oRINT);
}

VALUE
bna_sin(VALUE x)
{
    return unary_op(x, oSIN);
}

VALUE
bna_sinh(VALUE x)
{
    return unary_op(x, oSINH);
}

VALUE
bna_sqrt(VALUE x)
{
    return unary_op(x, oSQRT);
}

VALUE
bna_tan(VALUE x)
{
    return unary_op(x, oTAN);
}

VALUE
bna_tanh(VALUE x)
{
    return unary_op(x, oTANH);
}

VALUE
bna_y0(VALUE x)
{
    return unary_op(x, oY0);
}

VALUE
bna_y1(VALUE x)
{
    return unary_op(x, oY1);
}

VALUE
m_cos(VALUE obj, VALUE x)
{
    if (bna_is_a(x, cBasicNumArray)) {
	return bna_cos(x);
    } else {
	x = rb_Float(x);
	return rb_float_new(cos(RFLOAT(x)->value));
    }
}

VALUE
m_sin(VALUE obj, VALUE x)
{
    if (bna_is_a(x, cBasicNumArray)) {
	return bna_sin(x);
    } else {
	x = rb_Float(x);
	return rb_float_new(sin(RFLOAT(x)->value));
    }
}

VALUE
m_tan(VALUE obj, VALUE x)
{
    if (bna_is_a(x, cBasicNumArray)) {
	return bna_tan(x);
    } else {
	x = rb_Float(x);
	return rb_float_new(tan(RFLOAT(x)->value));
    }
}

VALUE
m_exp(VALUE obj, VALUE x)
{
    if (bna_is_a(x, cBasicNumArray)) {
	return bna_exp(x);
    } else {
	x = rb_Float(x);
	return rb_float_new(exp(RFLOAT(x)->value));
    }
}

VALUE
m_log(VALUE obj, VALUE x)
{
    if (bna_is_a(x, cBasicNumArray)) {
	return bna_log(x);
    } else {
	x = rb_Float(x);
	return rb_float_new(log(RFLOAT(x)->value));
    }
}

VALUE
m_log10(VALUE obj, VALUE x)
{
    if (bna_is_a(x, cBasicNumArray)) {
	return bna_log10(x);
    } else {
	x = rb_Float(x);
	return rb_float_new(log10(RFLOAT(x)->value));
    }
}

VALUE
m_sqrt(VALUE obj, VALUE x)
{
    if (bna_is_a(x, cBasicNumArray)) {
	return bna_sqrt(x);
    } else {
	x = rb_Float(x);
	return rb_float_new(sqrt(RFLOAT(x)->value));
    }
}

static VALUE
bna_inspect_ary(VALUE obj)
{
    long i = 0;
    VALUE s, str;
    double *dp;
    BNA *ary;

    Data_Get_Struct(obj, BNA, ary);
    dp=ary->ptr;
    
    str = rb_str_new2("[");

    for (i=0; i<ary->len; i++) {
	s = rb_inspect(rb_float_new(dp[i]));
	if (i > 0) rb_str_cat(str, ", ", 2);
	rb_str_cat(str, RSTRING(s)->ptr, RSTRING(s)->len);
    }
    rb_str_cat(str, "]", 1);

    return str;
}

static VALUE
bna_ary_inspect(VALUE obj)
{
    BNA *ary;
    
    Data_Get_Struct(obj, BNA, ary);
    
    if (ary->len == 0) return rb_str_new2("[]");
    return bna_inspect_ary(obj);
}

void
Init_basicnumarray()
{
    cBasicNumArray = rb_define_class("BasicNumArray", rb_cObject);
    cBasicNumArrayDouble = rb_define_class("BasicNumArrayDouble", cBasicNumArray);
    rb_define_singleton_method(cBasicNumArray, "new", bna_s_new, 2);
    rb_define_method(cBasicNumArray, "len", bna_len, 0);
    rb_define_method(cBasicNumArray, "size", bna_len, 0);
    rb_define_method(cBasicNumArray, "[]", bna_aref, 1);
    rb_define_method(cBasicNumArray, "[]=", bna_aset, 2);
    rb_define_method(cBasicNumArray, "setall", bna_setall, 1);
    rb_define_method(cBasicNumArray, "coerce", bna_coerce, 1);
    rb_define_method(cBasicNumArray, "+", bna_add, 1);
    rb_define_method(cBasicNumArray, "-", bna_sub, 1);
    rb_define_method(cBasicNumArray, "*", bna_mul, 1);
    rb_define_method(cBasicNumArray, "/", bna_div, 1);
    rb_define_method(cBasicNumArray, "**", bna_pow, 1);
    rb_define_method(cBasicNumArray, "%", bna_mod, 1);
    rb_define_method(cBasicNumArray, "acos", bna_acos, 0);
    rb_define_method(cBasicNumArray, "acosh", bna_acosh, 0);
    rb_define_method(cBasicNumArray, "asin", bna_asin, 0);
    rb_define_method(cBasicNumArray, "asinh", bna_asinh, 0);
    rb_define_method(cBasicNumArray, "atan", bna_atan, 0);
    rb_define_method(cBasicNumArray, "atanh", bna_atanh, 0);
    rb_define_method(cBasicNumArray, "cbrt", bna_cbrt, 0);
    rb_define_method(cBasicNumArray, "ceil", bna_ceil, 0);
    rb_define_method(cBasicNumArray, "cos", bna_cos, 0);
    rb_define_method(cBasicNumArray, "cosh", bna_cosh, 0);
    rb_define_method(cBasicNumArray, "erf", bna_erf, 0);
    rb_define_method(cBasicNumArray, "erfc", bna_erfc, 0);
    rb_define_method(cBasicNumArray, "exp", bna_exp, 0);
    rb_define_method(cBasicNumArray, "expm1", bna_expm1, 0);
    rb_define_method(cBasicNumArray, "fabs", bna_fabs, 0);
    rb_define_method(cBasicNumArray, "floor", bna_floor, 0);
    rb_define_method(cBasicNumArray, "j0", bna_j0, 0);
    rb_define_method(cBasicNumArray, "j1", bna_j1, 0);
    rb_define_method(cBasicNumArray, "lgamma", bna_lgamma, 0);
    rb_define_method(cBasicNumArray, "log", bna_log, 0);
    rb_define_method(cBasicNumArray, "log10", bna_log10, 0);
    rb_define_method(cBasicNumArray, "log1p", bna_log1p, 0);
    rb_define_method(cBasicNumArray, "rint", bna_rint, 0);
    rb_define_method(cBasicNumArray, "sin", bna_sin, 0);
    rb_define_method(cBasicNumArray, "sinh", bna_sinh, 0);
    rb_define_method(cBasicNumArray, "sqrt", bna_sqrt, 0);
    rb_define_method(cBasicNumArray, "tan", bna_tan, 0);
    rb_define_method(cBasicNumArray, "tanh", bna_tanh, 0);
    rb_define_method(cBasicNumArray, "y0", bna_y0, 0);
    rb_define_method(cBasicNumArray, "y1", bna_y1, 0);
    rb_define_method(cBasicNumArray, "inspect", bna_ary_inspect, 0);
    rb_define_module_function(rb_mMath, "cos", m_cos, 1);
    rb_define_module_function(rb_mMath, "sin", m_sin, 1);
    rb_define_module_function(rb_mMath, "tan", m_tan, 1);
    rb_define_module_function(rb_mMath, "exp", m_exp, 1);
    rb_define_module_function(rb_mMath, "log", m_log, 1);
    rb_define_module_function(rb_mMath, "log10", m_log10, 1);
    rb_define_module_function(rb_mMath, "sqrt", m_sqrt, 1);

    iadd = rb_intern("+");
    isub = rb_intern("-");
    imul = rb_intern("*");
    idiv = rb_intern("/");
    ipow = rb_intern("**");
    imod = rb_intern("%");
    icoerce = rb_intern("coerce");
}