MolLib/Ruby_bind/ruby_types.c

   1 /*
   2  *  ruby_types.c
   3  *  Molby
   4  *
   5  *  Created by Toshi Nagata on 09/01/24.
   6  *  Copyright 2009 Toshi Nagata. All rights reserved.
   7  *
   8  This program is free software; you can redistribute it and/or modify
   9  it under the terms of the GNU General Public License as published by
  10  the Free Software Foundation version 2 of the License.
  11
  12  This program is distributed in the hope that it will be useful,
  13  but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  GNU General Public License for more details.
  16  */
  17
  18 #include "Molby.h"
  19
  20 #include <errno.h>
  21 #include <string.h>
  22 #include <ctype.h>
  23 #include <limits.h>
  24
  25 #pragma mark ====== Global Values ======
  26
  27 VALUE rb_cVector3D, rb_cTransform, rb_cIntGroup;
  28
  29 #pragma mark ====== Utility functions (Vector/Matrix) ======
  30
  31 /*
  32 static VALUE
  33 s_VectorClass(void)
  34 {
  35         if (rb_cVector == 0)
  36                 rb_cVector = rb_const_get(rb_cObject, rb_intern("Vector"));
  37         return rb_cVector;
  38 }
  39 */
  40
  41 /*
  42 VALUE
  43 ValueFromVector(const Vector *vp)
  44 {
  45         static ID mname = 0;
  46         if (mname == 0)
  47                 mname = rb_intern("[]");
  48         return rb_funcall(rb_cVector3D, mname, 3, rb_float_new(vp->x), rb_float_new(vp->y), rb_float_new(vp->z));
  49 }
  50 */
  51
  52 #pragma mark ====== Vector3D Class ======
  53
  54 void
  55 VectorFromValue(VALUE val, Vector *vp)
  56 {
  57         Vector *vp1;
  58         if (rb_obj_is_kind_of(val, rb_cVector3D)) {
  59                 Data_Get_Struct(val, Vector, vp1);
  60                 *vp = *vp1;
  61         } else {
  62                 static ID mname = 0;
  63                 if (mname == 0)
  64                         mname = rb_intern("[]");
  65                 vp->x = NUM2DBL(rb_Float(rb_funcall(val, mname, 1, INT2FIX(0))));
  66                 vp->y = NUM2DBL(rb_Float(rb_funcall(val, mname, 1, INT2FIX(1))));
  67                 vp->z = NUM2DBL(rb_Float(rb_funcall(val, mname, 1, INT2FIX(2))));
  68         }
  69 }
  70
  71 static VALUE
  72 s_Vector3D_Alloc(VALUE klass)
  73 {
  74         Vector *vp = ALLOC(Vector);
  75         vp->x = vp->y = vp->z = 0.0;
  76         return Data_Wrap_Struct(klass, 0, -1, vp);
  77 }
  78
  79 VALUE
  80 ValueFromVector(const Vector *vp)
  81 {
  82         Vector *vp1;
  83         VALUE retval = s_Vector3D_Alloc(rb_cVector3D);
  84         Data_Get_Struct(retval, Vector, vp1);
  85         *vp1 = *vp;
  86         return retval;
  87 }
  88
  89 /*
  90  *  call-seq:
  91  *     new
  92  *     new(Vector3D)
  93  *     new(Array)
  94  *
  95  *  Returns a new Vector3D object. In the first form, a zero vector
  96  *  is returned. In the second form, the given vector3d is duplicated.
  97  *  In the third form, a vector [ary[0], ary[1], ary[2]] is returned.
  98  *  The argument ary can be anything that responds to '[]' method.
  99  */
 100 static VALUE
 101 s_Vector3D_Initialize(int argc, VALUE *argv, VALUE self)
 102 {
 103         VALUE val;
 104         Vector *vp;
 105         Data_Get_Struct(self, Vector, vp);
 106         rb_scan_args(argc, argv, "01", &val);
 107         if (!NIL_P(val))
 108                 VectorFromValue(val, vp);
 109         return Qnil;
 110 }
 111
 112 /*
 113  *  call-seq:
 114  *     size -> Integer
 115  *
 116  *  Returns 3. This method is present only to be consistent with classes like
 117  *  Array or Vector.
 118  */
 119 static VALUE
 120 s_Vector3D_Size(VALUE self)
 121 {
 122         return INT2FIX(3);
 123 }
 124
 125 /*
 126  *  call-seq:
 127  *     self[index]        -> Float
 128  *
 129  *  Element Reference---Returns the element at the given index. If the index is
 130  *  less than 0 or more than 2, an exception is thrown.
 131  */
 132 static VALUE
 133 s_Vector3D_ElementAtIndex(VALUE self, VALUE val)
 134 {
 135         Vector *vp;
 136         Double w;
 137         int n = NUM2INT(val);
 138         Data_Get_Struct(self, Vector, vp);
 139         if (n < 0 || n >= 3)
 140                 rb_raise(rb_eMolbyError, "index to Vector3D out of range");
 141         w = (n == 0 ? vp->x : (n == 1 ? vp->y : vp->z));
 142         return rb_float_new(w);
 143 }
 144
 145 /*
 146  *  call-seq:
 147  *     self[index] = val
 148  *
 149  *  Element Assignment---Set the element at _index_. If _index_ is
 150  *  less than 0 or more than 2, an exception is thrown.
 151  */
 152 static VALUE
 153 s_Vector3D_SetElementAtIndex(VALUE self, VALUE idx, VALUE val)
 154 {
 155         Vector *vp;
 156         Double w = NUM2DBL(rb_Float(val));
 157         int n = NUM2INT(idx);
 158         Data_Get_Struct(self, Vector, vp);
 159         if (n < 0 || n >= 3)
 160                 rb_raise(rb_eMolbyError, "index to Vector3D out of range");
 161         if (n == 0)
 162                 vp->x = w;
 163         else if (n == 1)
 164                 vp->y = w;
 165         else
 166                 vp->z = w;
 167         return rb_float_new(w);
 168 }
 169
 170 /*
 171  *  call-seq:
 172  *     self == val   ->   bool
 173  *
 174  *  Equality---Two vector3ds are equal if their elements are all equal.
 175  *  Usual caution about comparison between floating point numbers should be
 176  *  paid. Also consider using something like <code>vector3d.length < 1e-10</code>.
 177  */
 178 static VALUE
 179 s_Vector3D_IsEqual(VALUE self, VALUE val)
 180 {
 181         Vector *vp1, v2;
 182         Data_Get_Struct(self, Vector, vp1);
 183         VectorFromValue(val, &v2);
 184         if (vp1->x == v2.x && vp1->y == v2.y && vp1->z == v2.z)
 185                 return Qtrue;
 186         else return Qfalse;
 187 }
 188
 189 /*
 190  *  call-seq:
 191  *     self + val   ->   (new) Vector3D
 192  *
 193  *  Add two vectors element by element. Val is converted to Vector3D.
 194  */
 195 static VALUE
 196 s_Vector3D_Add(VALUE self, VALUE val)
 197 {
 198         Vector *vp1, v2;
 199         VALUE retval;
 200         Data_Get_Struct(self, Vector, vp1);
 201         VectorFromValue(val, &v2);
 202         retval = s_Vector3D_Alloc(rb_cVector3D);
 203         v2.x += vp1->x;
 204         v2.y += vp1->y;
 205         v2.z += vp1->z;
 206         return ValueFromVector(&v2);
 207 }
 208
 209 /*
 210  *  call-seq:
 211  *     self - val   ->   (new) Vector3D
 212  *
 213  *  Subtract two vectors element by element. Val is converted to Vector3D.
 214  */
 215 static VALUE
 216 s_Vector3D_Subtract(VALUE self, VALUE val)
 217 {
 218         Vector *vp1, v2;
 219         VALUE retval;
 220         Data_Get_Struct(self, Vector, vp1);
 221         VectorFromValue(val, &v2);
 222         retval = s_Vector3D_Alloc(rb_cVector3D);
 223         v2.x = vp1->x - v2.x;
 224         v2.y = vp1->y - v2.y;
 225         v2.z = vp1->z - v2.z;
 226         return ValueFromVector(&v2);
 227 }
 228
 229 /*
 230  *  call-seq:
 231  *     self.dot(val) ->  Float
 232  *
 233  *  Calculate the dot (inner) product of the two vectors. Val is converted to Vector3D.
 234  *
 235  *  <b>See Also:</b> Vector3D.*
 236  */
 237 static VALUE
 238 s_Vector3D_Dot(VALUE self, VALUE val)
 239 {
 240         Vector *vp1, v2;
 241         Data_Get_Struct(self, Vector, vp1);
 242         VectorFromValue(val, &v2);
 243         return rb_float_new(vp1->x * v2.x + vp1->y * v2.y + vp1->z * v2.z);
 244 }
 245
 246 /*
 247  *  call-seq:
 248  *     self * numeric           ->  (new) Vector3D
 249  *     self * val    ->  Float
 250  *
 251  *  In the first form, the vector is scaled by the numeric. In the second
 252  *  form, the dot (inner) product of the two vectors are returned, which is equivalent to
 253  *  self.dot(val).
 254  */
 255 static VALUE
 256 s_Vector3D_Multiply(VALUE self, VALUE val)
 257 {
 258         Vector *vp1, v2;
 259         Data_Get_Struct(self, Vector, vp1);
 260         if (rb_obj_is_kind_of(val, rb_cNumeric)) {
 261                 double w = NUM2DBL(rb_Float(val));
 262                 v2.x = vp1->x * w;
 263                 v2.y = vp1->y * w;
 264                 v2.z = vp1->z * w;
 265                 return ValueFromVector(&v2);
 266         } else return s_Vector3D_Dot(self, val);
 267 }
 268
 269 /*
 270  *  call-seq:
 271  *     self / numeric           ->  (new) Vector3D
 272  *
 273  *  The vector is scaled by the inverse of the given numeric.
 274  */
 275 static VALUE
 276 s_Vector3D_Divide(VALUE self, VALUE val)
 277 {
 278         Vector *vp1, v2;
 279         double w = NUM2DBL(rb_Float(val));
 280         Data_Get_Struct(self, Vector, vp1);
 281         v2.x = vp1->x / w;
 282         v2.y = vp1->y / w;
 283         v2.z = vp1->z / w;
 284         return ValueFromVector(&v2);
 285 }
 286
 287 /*
 288  *  call-seq:
 289  *     self.cross(val) ->  (new) Vector3D
 290  *
 291  *  Calculate the cross (outer) product of the two vectors. Val is converted to Vector3D.
 292  */
 293 static VALUE
 294 s_Vector3D_Cross(VALUE self, VALUE val)
 295 {
 296         Vector *vp1, v2, v3;
 297         Data_Get_Struct(self, Vector, vp1);
 298         VectorFromValue(val, &v2);
 299         v3.x = vp1->y * v2.z - vp1->z * v2.y;
 300         v3.y = vp1->z * v2.x - vp1->x * v2.z;
 301         v3.z = vp1->x * v2.y - vp1->y * v2.x;
 302         return ValueFromVector(&v3);
 303 }
 304
 305 /*
 306  *  call-seq:
 307  *     -self                ->  (new) Vector3D
 308  *
 309  *  Calculate the opposite vector.
 310  */
 311 static VALUE
 312 s_Vector3D_UnaryMinus(VALUE self)
 313 {
 314         Vector *vp, v1;
 315         Data_Get_Struct(self, Vector, vp);
 316         v1.x = -vp->x;
 317         v1.y = -vp->y;
 318         v1.z = -vp->z;
 319         return ValueFromVector(&v1);
 320 }
 321
 322 /*
 323  *  call-seq:
 324  *     length                ->  Float
 325  *
 326  *  Calculate the Pythagorean length of the vector.
 327  *  Note that this method is <em>not</em> an alias of Vector3D#size, which returns 3.
 328  */
 329 static VALUE
 330 s_Vector3D_Length(VALUE self)
 331 {
 332         Vector *vp;
 333         Data_Get_Struct(self, Vector, vp);
 334         return rb_float_new(sqrt(vp->x * vp->x + vp->y * vp->y + vp->z * vp->z));
 335 }
 336
 337 /*
 338  *  call-seq:
 339  *     length2                ->  Float
 340  *
 341  *  Calculate the square of the Pythagorean length of the vector.
 342  */
 343 static VALUE
 344 s_Vector3D_Length2(VALUE self)
 345 {
 346         Vector *vp;
 347         Data_Get_Struct(self, Vector, vp);
 348         return rb_float_new(vp->x * vp->x + vp->y * vp->y + vp->z * vp->z);
 349 }
 350
 351 /*
 352  *  call-seq:
 353  *     normalize              ->  (new) Vector3D
 354  *
 355  *  Returns a unit vector with the same direction. Raises an exception when the
 356  *  vector is a zero vector.
 357  */
 358 static VALUE
 359 s_Vector3D_Normalize(VALUE self)
 360 {
 361         Vector *vp, v;
 362         double w;
 363         Data_Get_Struct(self, Vector, vp);
 364         w = 1.0 / sqrt(vp->x * vp->x + vp->y * vp->y + vp->z * vp->z);
 365         if (!isfinite(w))
 366                 rb_raise(rb_eMolbyError, "trying to normalize a (nearly) zero vector");
 367         v.x = vp->x * w;
 368         v.y = vp->y * w;
 369         v.z = vp->z * w;
 370         return ValueFromVector(&v);
 371 }
 372
 373 /*
 374  *  call-seq:
 375  *     to_a                    ->  Array
 376  *
 377  *  Returns [self.x, self.y, self.z].
 378  */
 379 static VALUE
 380 s_Vector3D_ToArray(VALUE self)
 381 {
 382         Vector *vp;
 383         Data_Get_Struct(self, Vector, vp);
 384         return rb_ary_new3(3, rb_float_new(vp->x), rb_float_new(vp->y), rb_float_new(vp->z));
 385 }
 386
 387 /*
 388  *  call-seq:
 389  *     each {|item| ...}
 390  *
 391  *  Calls block for x, y, z elements, passing that element as a parameter.
 392  */
 393 static VALUE
 394 s_Vector3D_Each(VALUE self)
 395 {
 396         Vector *vp;
 397         Data_Get_Struct(self, Vector, vp);
 398         rb_yield(rb_float_new(vp->x));
 399         rb_yield(rb_float_new(vp->y));
 400         rb_yield(rb_float_new(vp->z));
 401         return self;
 402 }
 403
 404 /*
 405  *  call-seq:
 406  *     x       -> Float
 407  *
 408  *  Get the x element of the vector.
 409  */
 410 static VALUE
 411 s_Vector3D_GetX(VALUE self)
 412 {
 413         Vector *vp;
 414         Data_Get_Struct(self, Vector, vp);
 415         return rb_float_new(vp->x);
 416 }
 417
 418 /*
 419  *  call-seq:
 420  *     y       -> Float
 421  *
 422  *  Get the y element of the vector.
 423  */
 424 static VALUE
 425 s_Vector3D_GetY(VALUE self)
 426 {
 427         Vector *vp;
 428         Data_Get_Struct(self, Vector, vp);
 429         return rb_float_new(vp->y);
 430 }
 431
 432 /*
 433  *  call-seq:
 434  *     z       -> Float
 435  *
 436  *  Get the z element of the vector.
 437  */
 438 static VALUE
 439 s_Vector3D_GetZ(VALUE self)
 440 {
 441         Vector *vp;
 442         Data_Get_Struct(self, Vector, vp);
 443         return rb_float_new(vp->z);
 444 }
 445
 446 /*
 447  *  call-seq:
 448  *     x = val
 449  *
 450  *  Set the x element of the vector.
 451  */
 452 static VALUE
 453 s_Vector3D_SetX(VALUE self, VALUE val)
 454 {
 455         Vector *vp;
 456         Data_Get_Struct(self, Vector, vp);
 457         vp->x = NUM2DBL(rb_Float(val));
 458         return rb_float_new(vp->x);
 459 }
 460
 461 /*
 462  *  call-seq:
 463  *     y = val
 464  *
 465  *  Set the y element of the vector.
 466  */
 467 static VALUE
 468 s_Vector3D_SetY(VALUE self, VALUE val)
 469 {
 470         Vector *vp;
 471         Data_Get_Struct(self, Vector, vp);
 472         vp->y = NUM2DBL(rb_Float(val));
 473         return rb_float_new(vp->y);
 474 }
 475
 476 /*
 477  *  call-seq:
 478  *     z = val
 479  *
 480  *  Set the z element of the vector.
 481  */
 482 static VALUE
 483 s_Vector3D_SetZ(VALUE self, VALUE val)
 484 {
 485         Vector *vp;
 486         Data_Get_Struct(self, Vector, vp);
 487         vp->z = NUM2DBL(rb_Float(val));
 488         return rb_float_new(vp->z);
 489 }
 490
 491 /*
 492  *  call-seq:
 493  *     Vector3d[fx, fy, fz]   -> (new) Vector3D
 494  *
 495  *  Create a new vector3d object. Equivalent to Vector3D#new([fx, fy, fz]).
 496  */
 497 static VALUE
 498 s_Vector3D_Create(VALUE klass, VALUE args)
 499 {
 500         VALUE val = s_Vector3D_Alloc(klass);
 501         s_Vector3D_Initialize(1, &args, val);
 502         return val;
 503 }
 504
 505 /*
 506  *  call-seq:
 507  *     inspect        -> String
 508  *
 509  *  Create a readable string like "Vector3D[fx, fy, fz]".
 510  */
 511 static VALUE
 512 s_Vector3D_Inspect(VALUE self)
 513 {
 514         /*  self.class.name << self.to_a.inspect  */
 515 /*      VALUE klass = CLASS_OF(self); */
 516 /*      VALUE val = rb_funcall(klass, rb_intern("name"), 0); */
 517         VALUE val = rb_str_new2("Vector3D");
 518         self = s_Vector3D_ToArray(self);
 519         return rb_funcall(val, rb_intern("<<"), 1, rb_funcall(self, rb_intern("inspect"), 0));
 520 }
 521
 522 #pragma mark ====== Transform Class ======
 523
 524 void
 525 TransformFromValue(VALUE val, Transform *tp)
 526 {
 527         int i, j;
 528         Transform *tp1;
 529         if (rb_obj_is_kind_of(val, rb_cTransform)) {
 530                 Data_Get_Struct(val, Transform, tp1);
 531                 memmove(tp, tp1, sizeof(Transform));
 532         } else if (TYPE(val) == T_ARRAY) {
 533                 /*  Array must be
 534                  (1) [[a11 a21 a31] [a12 a22 a32] [a13 a23 a33] [a14 a24 a34]] or
 535                  (2) [a11 a21 a31 a12 a22 a32 a13 a23 a33 a14 a24 a34] */
 536                 int len = RARRAY_LEN(val);
 537                 VALUE *valp = RARRAY_PTR(val);
 538                 if (len == 12) {
 539                         for (i = 0; i < 12; i++)
 540                                 (*tp)[i] = NUM2DBL(rb_Float(valp[i]));
 541                         return;
 542                 } else if (len == 4) {
 543                         VALUE val2, *valp2;
 544                         for (i = 0; i < 4; i++) {
 545                                 val2 = valp[i];
 546                                 if (TYPE(val2) != T_ARRAY)
 547                                         val2 = rb_funcall(val2, rb_intern("to_a"), 0);
 548                                 if (TYPE(val2) != T_ARRAY || RARRAY_LEN(val2) != 3)
 549                                         goto array_format_error;
 550                                 valp2 = RARRAY_PTR(val2);
 551                                 for (j = 0; j < 3; j++)
 552                                         (*tp)[3 * i + j] = NUM2DBL(rb_Float(valp2[j]));
 553                         }
 554                         return;
 555                 }
 556         array_format_error:
 557                 rb_raise(rb_eMolbyError, "wrong array format; must be an array of either (1) four 3D column vectors or (2) 12 numerics");
 558         } else {
 559                 static ID index_mid = 0, row_size_mid, column_size_mid;
 560                 if (index_mid == 0) {
 561                         index_mid = rb_intern("[]");
 562                         row_size_mid = rb_intern("row_size");
 563                         column_size_mid = rb_intern("column_size");
 564                 }
 565                 if (rb_respond_to(val, row_size_mid) && rb_respond_to(val, column_size_mid)) {
 566                         /*  Matrix-type object  */
 567                         for (i = 0; i < 4; i++) {
 568                                 for (j = 0; j < 3; j++)
 569                                         (*tp)[i * 3 + j] = NUM2DBL(rb_Float(rb_funcall(val, index_mid, 2, INT2FIX(i), INT2FIX(j))));
 570                         }
 571                 } else {
 572                         /*  Other "array-like" object  */
 573                         for (i = 0; i < 12; i++)
 574                                 (*tp)[i] = NUM2DBL(rb_Float(rb_funcall(val, index_mid, 1, INT2FIX(i))));
 575                 }
 576         }
 577 }
 578
 579 static VALUE
 580 s_Transform_Alloc(VALUE klass)
 581 {
 582         Transform *tp = ALLOC(Transform);
 583         memset(tp, 0, sizeof(Transform));
 584         (*tp)[0] = (*tp)[4] = (*tp)[8] = 1.0;
 585         return Data_Wrap_Struct(klass, 0, -1, tp);
 586 }
 587
 588 VALUE
 589 ValueFromTransform(Transform *tp)
 590 {
 591         Transform *tp1;
 592         VALUE val = s_Transform_Alloc(rb_cTransform);
 593         Data_Get_Struct(val, Transform, tp1);
 594         memmove(tp1, tp, sizeof(Transform));
 595         return val;
 596 }
 597
 598 /*
 599  *  call-seq:
 600  *     new
 601  *     new(array)
 602  *     new(matrix)
 603  *
 604  *  Returns a new Transform object.
 605  *
 606  *  In the first form, an identity transform is returned.
 607  *
 608  *  In the second form, the array must be either of the following
 609  *  two forms:
 610  *  (1) [[a11 a21 a31] [a12 a22 a32] [a13 a23 a33] [a14 a24 a34]]
 611  *  (2) [a11 a21 a31 a12 a22 a32 a13 a23 a33 a14 a24 a34]
 612  *  where [a11..a33] denotes the rotation part and [a14 a24 a34] denotes
 613  *  the translation part. All vectors in (1) are column vectors.
 614  *
 615  *  In the third form, a new transform is built from a 3x4 matrix. The argument
 616  *  +matrix+ must respond to a method call <tt>matrix[col, row]</tt>
 617  *  where <tt>row</tt> is in <tt>0..2</tt> and <tt>col</tt> in <tt>0..3</tt>.
 618  */
 619 static VALUE
 620 s_Transform_Initialize(int argc, VALUE *argv, VALUE self)
 621 {
 622         VALUE val;
 623         Transform *tp;
 624         Data_Get_Struct(self, Transform, tp);
 625         rb_scan_args(argc, argv, "01", &val);
 626         if (!NIL_P(val))
 627                 TransformFromValue(val, tp);
 628         return Qnil;
 629 }
 630
 631 static VALUE
 632 s_Transform_NewFromTransform(Transform *tp)
 633 {
 634         Transform *tp1;
 635         VALUE retval = s_Transform_Alloc(rb_cTransform);
 636         Data_Get_Struct(retval, Transform, tp1);
 637         memmove(tp1, tp, sizeof(Transform));
 638         return retval;
 639 }
 640
 641 /*
 642  *  call-seq:
 643  *     from_columns(c1, c2, c3, c4)
 644  *
 645  *  Returns a new Transform object built from four column vectors. The arguments
 646  *  <tt>c1..c4</tt> are vectors of (at least) three-dimension. This is equivalent
 647  *  to <tt>Transform.new([c1, c2, c3, c4])</tt>.
 648  */
 649 static VALUE
 650 s_Transform_NewFromColumns(VALUE klass, VALUE val)
 651 {
 652         Transform tr;
 653         int i, j, n;
 654         VALUE *valp;
 655         static ID to_a_mid = 0;
 656         if (to_a_mid == 0)
 657                 to_a_mid = rb_intern("to_a");
 658         if (TYPE(val) != T_ARRAY)
 659                 val = rb_funcall(val, to_a_mid, 0);
 660         memset(tr, 0, sizeof(tr));
 661         n = RARRAY_LEN(val);
 662         valp = RARRAY_PTR(val);
 663         for (i = 0; i < 4; i++) {
 664                 int nn;
 665                 VALUE *valpp;
 666                 double w[3];
 667                 w[0] = w[1] = w[2] = 0.0;
 668                 if (i < n) {
 669                         val = valp[i];
 670                         if (TYPE(val) != T_ARRAY)
 671                                 val = rb_funcall(val, to_a_mid, 0);
 672                         nn = RARRAY_LEN(val);
 673                         valpp = RARRAY_PTR(val);
 674                         for (j = 0; j < 3 && j < nn; j++)
 675                                 w[j] = NUM2DBL(rb_Float(valpp[j]));
 676                 }
 677                 for (j = 0; j < 3; j++)
 678                         tr[i * 3 + j] = w[j];
 679         }
 680         return s_Transform_NewFromTransform(&tr);
 681 }
 682
 683 /*
 684  *  call-seq:
 685  *     from_rows(r1, r2, r3)
 686  *
 687  *  Returns a new Transform object built from three row vectors. The arguments
 688  *  <tt>r1, r2, r3</tt> are vectors of (at least) four-dimension.
 689  */
 690 static VALUE
 691 s_Transform_NewFromRows(VALUE klass, VALUE val)
 692 {
 693         Transform tr;
 694         int i, j, n;
 695         VALUE *valp;
 696         static ID to_a_mid = 0;
 697         if (to_a_mid == 0)
 698                 to_a_mid = rb_intern("to_a");
 699         if (TYPE(val) != T_ARRAY)
 700                 val = rb_funcall(val, to_a_mid, 0);
 701         memset(tr, 0, sizeof(tr));
 702         n = RARRAY_LEN(val);
 703         valp = RARRAY_PTR(val);
 704         for (i = 0; i < 3; i++) {
 705                 int nn;
 706                 VALUE *valpp;
 707                 double w[4];
 708                 w[0] = w[1] = w[2] = w[3] = 0.0;
 709                 if (i < n) {
 710                         val = valp[i];
 711                         if (TYPE(val) != T_ARRAY)
 712                                 val = rb_funcall(val, to_a_mid, 0);
 713                         nn = RARRAY_LEN(val);
 714                         valpp = RARRAY_PTR(val);
 715                         for (j = 0; j < 4 && j < nn; j++)
 716                                 w[j] = NUM2DBL(rb_Float(valpp[j]));
 717                 }
 718                 for (j = 0; j < 4; j++)
 719                         tr[j * 3 + i] = w[j];
 720         }
 721         return s_Transform_NewFromTransform(&tr);
 722 }
 723
 724 /*
 725  *  call-seq:
 726  *     self[i, j]  -> Float
 727  *
 728  *  Get the element (+i+,+j+) of the transform matrix, i.e. column +i+, row +j+.
 729  *  Be careful about the order of the arguments. It follows convention of multi-dimensional arrays
 730  *  rather than mathematical notation.
 731  */
 732 static VALUE
 733 s_Transform_ElementAtIndex(VALUE self, VALUE val1, VALUE val2)
 734 {
 735         Transform *tp;
 736         double w;
 737         int n1 = NUM2INT(val1);
 738         int n2 = NUM2INT(val2);
 739         Data_Get_Struct(self, Transform, tp);
 740         if (n1 < 0 || n1 >= 4 || n2 < 0 || n2 >= 3)
 741                 rb_raise(rb_eMolbyError, "index to Transform out of range");
 742         w = (*tp)[n1 * 3 + n2];
 743         return rb_float_new(w);
 744 }
 745
 746 /*
 747  *  call-seq:
 748  *     self[i, j] = val
 749  *
 750  *  Set the element (+i+,+j+) of the transform matrix, i.e. column +i+, row +j+.
 751  *  Be careful about the order of the arguments. It follows convention of multi-dimensional arrays
 752  *  rather than mathematical notation.
 753  */
 754 static VALUE
 755 s_Transform_SetElementAtIndex(VALUE self, VALUE idx1, VALUE idx2, VALUE val)
 756 {
 757         Transform *tp;
 758         double w;
 759         int n1 = NUM2INT(idx1);
 760         int n2 = NUM2INT(idx2);
 761         Data_Get_Struct(self, Transform, tp);
 762         if (n1 < 0 || n1 >= 4 || n2 < 0 || n2 >= 3)
 763                 rb_raise(rb_eMolbyError, "index to Transform out of range");
 764         w = NUM2DBL(rb_Float(val));
 765         (*tp)[n1 * 3 + n2] = w;
 766         return rb_float_new(w);
 767 }
 768
 769 /*
 770  *  call-seq:
 771  *     self == val  -> bool
 772  *
 773  *  Returns +true+ if and only if all the corresponding elements are equal.
 774  *  Usual caution about the comparison of floating-point numbers should be paid.
 775  */
 776 static VALUE
 777 s_Transform_IsEqual(VALUE self, VALUE val)
 778 {
 779         Transform *tp1, tr;
 780         int i;
 781         Data_Get_Struct(self, Transform, tp1);
 782         TransformFromValue(val, &tr);
 783         for (i = 0; i < 12; i++) {
 784                 if ((*tp1)[i] != tr[i])
 785                         return Qfalse;
 786         }
 787         return Qtrue;
 788 }
 789
 790 /*
 791  *  call-seq:
 792  *     self + val  -> (new) Transform
 793  *
 794  *  Returns a new transform corresponding to the sum of the two transform matrix.
 795  */
 796 static VALUE
 797 s_Transform_Add(VALUE self, VALUE val)
 798 {
 799         Transform *tp1, tr;
 800         int i;
 801         Data_Get_Struct(self, Transform, tp1);
 802         TransformFromValue(val, &tr);
 803         for (i = 0; i < 12; i++)
 804                 tr[i] += (*tp1)[i];
 805         return s_Transform_NewFromTransform(&tr);
 806 }
 807
 808 /*
 809  *  call-seq:
 810  *     self - val  -> (new) Transform
 811  *
 812  *  Returns a new transform corresponding to the difference of the two transform matrix.
 813  */
 814 static VALUE
 815 s_Transform_Subtract(VALUE self, VALUE val)
 816 {
 817         Transform *tp1, tr;
 818         int i;
 819         Data_Get_Struct(self, Transform, tp1);
 820         TransformFromValue(val, &tr);
 821         for (i = 0; i < 12; i++)
 822                 tr[i] = (*tp1)[i] - tr[i];
 823         return s_Transform_NewFromTransform(&tr);
 824 }
 825
 826 /*
 827  *  call-seq:
 828  *     self * numeric          -> (new) Transform
 829  *     self * Vector3D         -> (new) Vector3D
 830  *     self * other_transform  -> (new) Transform
 831  *
 832  *  Perform the matrix multiplication. In the first form, a new matrix with scaled elements
 833  *  is returned. In the second, the transformed vector is returned. In the third form,
 834  *  the multiple of the two matrices is returned.
 835  */
 836 static VALUE
 837 s_Transform_Multiply(VALUE self, VALUE val)
 838 {
 839         Transform *tp1, tr;
 840         int i;
 841         Data_Get_Struct(self, Transform, tp1);
 842         if (rb_obj_is_kind_of(val, rb_cNumeric)) {
 843                 double w = NUM2DBL(rb_Float(val));
 844                 for (i = 0; i < 12; i++)
 845                         tr[i] = (*tp1)[i] * w;
 846                 return s_Transform_NewFromTransform(&tr);
 847         } else {
 848                 static ID size_mid = 0;
 849                 if (size_mid == 0)
 850                         size_mid = rb_intern("size");
 851                 if (rb_respond_to(val, size_mid) && NUM2INT(rb_funcall(val, size_mid, 0)) == 3) {
 852                         /*  A 3D vector  */
 853                         Vector v;
 854                         VectorFromValue(val, &v);
 855                         TransformVec(&v, *tp1, &v);
 856                         return ValueFromVector(&v);
 857                 } else {
 858                         /*  Transform  */
 859                         TransformFromValue(val, &tr);
 860                         TransformMul(tr, *tp1, tr);
 861                         return s_Transform_NewFromTransform(&tr);
 862                 }
 863         }
 864 }
 865
 866 /*
 867  *  call-seq:
 868  *     identity  -> Transform
 869  *
 870  *  Returns an identity transform, <tt>[[1,0,0], [0,1,0], [0,0,1], [0,0,0]]</tt>.
 871  */
 872 static VALUE
 873 s_Transform_Identity(VALUE klass)
 874 {
 875         Transform tr;
 876         memset(tr, 0, sizeof(tr));
 877         tr[0] = tr[4] = tr[8] = 1.0;
 878         return s_Transform_NewFromTransform(&tr);
 879 }
 880
 881 /*
 882  *  call-seq:
 883  *     zero  -> Transform
 884  *
 885  *  Returns a zero transform, <tt>[[0,0,0], [0,0,0], [0,0,0], [0,0,0]]</tt>.
 886  */
 887 static VALUE
 888 s_Transform_Zero(VALUE klass)
 889 {
 890         Transform tr;
 891         memset(tr, 0, sizeof(tr));
 892         return s_Transform_NewFromTransform(&tr);
 893 }
 894
 895 /*
 896  *  call-seq:
 897  *     diagonal(Array)
 898  *     diagonal(f1, f2 = nil, f3 = nil)
 899  *
 900  *  Returns a diagonal transform (the translational componets are all zero).
 901  *  In the first form, <tt>array[0], array[1], array[2]</tt> are for the
 902  *  x, y, z components, respectively. In the second form, <tt>f1, f2, f3</tt>
 903  *  are the x, y, z components. If <tt>f3</tt> is not given, the <tt>f2</tt>
 904  *  is used for the z components. If <tt>f2</tt> is not given, the <tt>f1</tt>
 905  *  is used for the y and z components.
 906  */
 907 static VALUE
 908 s_Transform_Diagonal(int argc, VALUE *argv, VALUE klass)
 909 {
 910         Transform tr;
 911         VALUE arg1, arg2, arg3;
 912         memset(tr, 0, sizeof(tr));
 913         rb_scan_args(argc, argv, "12", &arg1, &arg2, &arg3);
 914         if (TYPE(arg1) == T_ARRAY) {
 915                 /*  [d1, d2, d3]  */
 916                 VALUE *valp;
 917                 int len;
 918                 len = RARRAY_LEN(arg1);
 919                 valp = RARRAY_PTR(arg1);
 920                 if (len >= 1)
 921                         tr[0] = tr[4] = tr[8] = NUM2DBL(rb_Float(valp[0]));
 922                 if (len >= 2)
 923                         tr[4] = tr[8] = NUM2DBL(rb_Float(valp[1]));
 924                 if (len >= 3)
 925                         tr[8] = NUM2DBL(rb_Float(valp[2]));
 926         } else {
 927                 tr[0] = tr[4] = tr[8] = NUM2DBL(rb_Float(arg1));
 928                 if (!NIL_P(arg2)) {
 929                         tr[4] = tr[8] = NUM2DBL(rb_Float(arg2));
 930                         if (!NIL_P(arg3))
 931                                 tr[8] = NUM2DBL(rb_Float(arg3));
 932                 }
 933         }
 934         return s_Transform_NewFromTransform(&tr);
 935 }
 936
 937 /*
 938  *  call-seq:
 939  *     inverse  -> (new) Transform
 940  *
 941  *  Returns the inverse transform. If the matrix is not regular, an exception is raised.
 942  */
 943 static VALUE
 944 s_Transform_Inverse(VALUE self)
 945 {
 946         Transform *tp1, tr;
 947         Data_Get_Struct(self, Transform, tp1);
 948         if (TransformInvert(tr, *tp1))
 949                 rb_raise(rb_eMolbyError, "the transform matrix is not regular");
 950         return s_Transform_NewFromTransform(&tr);
 951 }
 952
 953 /*
 954  *  call-seq:
 955  *     self / val -> (new) Transform
 956  *
 957  *  Returns self * val.invert. If val is not a regular transform,
 958  *  an exception is raised.
 959  */
 960 static VALUE
 961 s_Transform_Divide(VALUE self, VALUE val)
 962 {
 963         Transform *tp1, tr;
 964         Data_Get_Struct(self, Transform, tp1);
 965         TransformFromValue(val, &tr);
 966         if (TransformInvert(tr, tr))
 967                 rb_raise(rb_eMolbyError, "the transform matrix is not regular");
 968         TransformMul(tr, *tp1, tr);
 969         return s_Transform_NewFromTransform(&tr);
 970 }
 971
 972 /*
 973  *  call-seq:
 974  *     transpose -> (new) Transform
 975  *
 976  *  Returns a new transform in which the rotation component is transposed from the original.
 977  */
 978 static VALUE
 979 s_Transform_Transpose(VALUE self)
 980 {
 981         Transform *tp1, tr;
 982         Data_Get_Struct(self, Transform, tp1);
 983         TransformTranspose(tr, *tp1);
 984         return s_Transform_NewFromTransform(&tr);
 985 }
 986
 987 /*
 988  *  call-seq:
 989  *     determinant -> Float
 990  *
 991  *  Returns the determinant of the transform.
 992  */
 993 static VALUE
 994 s_Transform_Determinant(VALUE self)
 995 {
 996         Transform *tp1;
 997         Data_Get_Struct(self, Transform, tp1);
 998         return rb_float_new(TransformDeterminant(*tp1));
 999 }
1000
1001 /*
1002  *  call-seq:
1003  *     trace -> Float
1004  *
1005  *  Returns the trace (sum of the diagonal elements) of the transform.
1006  */
1007 static VALUE
1008 s_Transform_Trace(VALUE self)
1009 {
1010         Transform *tp1;
1011         Data_Get_Struct(self, Transform, tp1);
1012         return rb_float_new((*tp1)[0] + (*tp1)[4] + (*tp1)[8]);
1013 }
1014
1015 /*
1016  *  call-seq:
1017  *     column(index) -> Vector3D
1018  *
1019  *  Returns the index-th (0..3) column vector.
1020  */
1021 static VALUE
1022 s_Transform_Column(VALUE self, VALUE val)
1023 {
1024         Transform *tp1;
1025         Vector v;
1026         int n = NUM2INT(val);
1027         Data_Get_Struct(self, Transform, tp1);
1028         if (n < 0 || n >= 4)
1029                 rb_raise(rb_eMolbyError, "row index out of range");
1030         v.x = (*tp1)[n * 3];
1031         v.y = (*tp1)[n * 3 + 1];
1032         v.z = (*tp1)[n * 3 + 2];
1033         return ValueFromVector(&v);
1034 }
1035
1036 /*
1037  *  call-seq:
1038  *     eigenvalues -> [[k1, k2, k3], v1, v2, v3]
1039  *
1040  *  Calculate the eigenvalues and eigenvectors. The matrix must be symmetric.
1041  */
1042 static VALUE
1043 s_Transform_Eigenvalues(VALUE self)
1044 {
1045         Transform *tp1;
1046         Vector v[3];
1047         Double d[3];
1048         int info;
1049         Data_Get_Struct(self, Transform, tp1);
1050         if ((info = MatrixSymDiagonalize(*((Mat33 *)tp1), d, v)) != 0)
1051                 rb_raise(rb_eMolbyError, "cannot diagonalize the given matrix: info = %d", info);
1052         return rb_ary_new3(4, rb_ary_new3(3, rb_float_new(d[0]), rb_float_new(d[1]), rb_float_new(d[2])), ValueFromVector(v), ValueFromVector(v + 1), ValueFromVector(v + 2));
1053 }
1054
1055 /*
1056  *  call-seq:
1057  *     Transform[*args] -> (new) Transform
1058  *
1059  *  Create a new transform. Equivalent to Transform.new(args).
1060  */
1061 static VALUE
1062 s_Transform_Create(VALUE klass, VALUE args)
1063 {
1064         VALUE val = s_Transform_Alloc(klass);
1065         s_Transform_Initialize(1, &args, val);
1066         return val;
1067 }
1068
1069 /*
1070  *  call-seq:
1071  *     to_a  -> Array
1072  *
1073  *  Convert a transform to an array of 12 float numbers.
1074  */
1075 static VALUE
1076 s_Transform_ToArray(VALUE self)
1077 {
1078         Transform *tp1;
1079         VALUE val[12];
1080         int i;
1081         Data_Get_Struct(self, Transform, tp1);
1082         for (i = 0; i < 12; i++)
1083                 val[i] = rb_float_new((*tp1)[i]);
1084         return rb_ary_new4(12, val);
1085 }
1086
1087 /*
1088  *  call-seq:
1089  *     inspect  -> String
1090  *
1091  *  Convert a transform to a string like
1092  *  "Transform[[a11,a21,a31],[a12,a22,a32],[a13,a23,a33],[a14,a24,a34]]".
1093  */
1094 static VALUE
1095 s_Transform_Inspect(VALUE self)
1096 {
1097         Transform *tp;
1098         int i, j;
1099 /*      VALUE klass = CLASS_OF(self); */
1100 /*      VALUE val = rb_funcall(klass, rb_intern("name"), 0); */
1101         VALUE val = rb_str_new2("Transform");
1102         ID mid = rb_intern("<<");
1103         ID mid2 = rb_intern("inspect");
1104         rb_str_cat(val, "[", 1);
1105         Data_Get_Struct(self, Transform, tp);
1106         for (i = 0; i < 4; i++) {
1107                 rb_str_cat(val, "[", 1);
1108                 for (j = 0; j < 3; j++) {
1109                         double f;
1110                         f = (*tp)[i * 3 + j];
1111                         rb_funcall(val, mid, 1, rb_funcall(rb_float_new(f), mid2, 0));
1112                         if (j < 2)
1113                                 rb_str_cat(val, ",", 1);
1114                 }
1115                 rb_str_cat(val, "]", 1);
1116                 if (i < 3)
1117                         rb_str_cat(val, ",", 1);
1118         }
1119         rb_str_cat(val, "]", 1);
1120         return val;
1121 }
1122
1123 /*
1124  *  call-seq:
1125  *     translation(vec) -> (new) Transform
1126  *
1127  *  Returns a transform corresponding to translation along the given vector. Equivalent
1128  *  to <code>Transform[[1,0,0],[0,1,0],[0,0,1],vec]</code>.
1129  */
1130 static VALUE
1131 s_Transform_Translation(VALUE klass, VALUE vec)
1132 {
1133         Transform *tp;
1134         Vector v;
1135         VALUE val = s_Transform_Alloc(klass);
1136         Data_Get_Struct(val, Transform, tp);
1137         VectorFromValue(vec, &v);
1138         (*tp)[9] = v.x;
1139         (*tp)[10] = v.y;
1140         (*tp)[11] = v.z;
1141         return val;
1142 }
1143
1144 /*
1145  *  call-seq:
1146  *     rotation(axis, angle, center = [0,0,0]) -> (new) Transform
1147  *
1148  *  Returns a transform corresponding to the rotation along the given axis and angle.
1149  *  Angle is given in degree. If center is also given, that point will be the center of rotation.
1150  */
1151 static VALUE
1152 s_Transform_Rotation(int argc, VALUE *argv, VALUE klass)
1153 {
1154         Transform tr;
1155         VALUE axis, angle, center;
1156         Vector av, cv;
1157         double ang;
1158         rb_scan_args(argc, argv, "21", &axis, &angle, &center);
1159         VectorFromValue(axis, &av);
1160         if (NIL_P(center))
1161                 cv.x = cv.y = cv.z = 0.0;
1162         else
1163                 VectorFromValue(center, &cv);
1164         ang = NUM2DBL(rb_Float(angle)) * kDeg2Rad;
1165         if (TransformForRotation(tr, &av, ang, &cv))
1166                 rb_raise(rb_eMolbyError, "rotation axis cannot be a zero vector");
1167         return ValueFromTransform(&tr);
1168 }
1169
1170 /*
1171  *  call-seq:
1172  *     reflection(axis, center = [0,0,0]) -> (new)Transform
1173  *
1174  *  Returns a transform corresponding to the reflection along the given axis. If
1175  *  center is also given, that point will be fixed.
1176  */
1177 static VALUE
1178 s_Transform_Reflection(int argc, VALUE *argv, VALUE klass)
1179 {
1180         VALUE axis, center;
1181         Vector av, cv;
1182         Transform tr;
1183         rb_scan_args(argc, argv, "11", &axis, &center);
1184         VectorFromValue(axis, &av);
1185         if (NIL_P(center))
1186                 cv.x = cv.y = cv.z = 0.0;
1187         else
1188                 VectorFromValue(center, &cv);
1189         if (TransformForReflection(tr, &av, &cv))
1190                 rb_raise(rb_eMolbyError, "reflection axis cannot be a zero vector");
1191         return ValueFromTransform(&tr);
1192 }
1193
1194 /*
1195  *  call-seq:
1196  *     inversion(center = [0,0,0]) -> (new) Transform
1197  *
1198  *  Returns a transform corresponding to the inversion along the given point.
1199  */
1200 static VALUE
1201 s_Transform_Inversion(int argc, VALUE *argv, VALUE klass)
1202 {
1203         VALUE center;
1204         Vector cv;
1205         Transform tr;
1206         rb_scan_args(argc, argv, "01", &center);
1207         if (NIL_P(center))
1208                 cv.x = cv.y = cv.z = 0.0;
1209         else
1210                 VectorFromValue(center, &cv);
1211         TransformForInversion(tr, &cv);
1212         return ValueFromTransform(&tr);
1213 }
1214
1215 #pragma mark ====== IntGroup Class ======
1216
1217 IntGroup *
1218 IntGroupFromValue(VALUE val)
1219 {
1220         IntGroup *ig;
1221         if (!rb_obj_is_kind_of(val, rb_cIntGroup))
1222                 val = rb_funcall(rb_cIntGroup, rb_intern("new"), 1, val);
1223         Data_Get_Struct(val, IntGroup, ig);
1224         IntGroupRetain(ig);
1225         return ig;
1226 }
1227
1228 VALUE
1229 ValueFromIntGroup(IntGroup *ig)
1230 {
1231         if (ig == NULL)
1232                 return Qnil;
1233         IntGroupRetain(ig);
1234     return Data_Wrap_Struct(rb_cIntGroup, 0, (void (*)(void *))IntGroupRelease, ig);
1235 }
1236
1237 void
1238 IntGroup_RaiseIfError(int err)
1239 {
1240         if (err != 0) {
1241                 const char *s;
1242                 switch (err) {
1243                         case kIntGroupStatusOutOfMemory: s = "out of memory"; break;
1244                         case kIntGroupStatusOutOfRange: s = "out of range"; break;
1245                         default: s = ""; break;
1246                 }
1247                 rb_raise(rb_eMolbyError, "%s error occurred during IntGroup operation", s);
1248         }
1249 }
1250
1251 /*  Allocator  */
1252 VALUE
1253 IntGroup_Alloc(VALUE klass)
1254 {
1255         IntGroup *ig = IntGroupNew();
1256     return Data_Wrap_Struct(klass, 0, (void (*)(void *))IntGroupRelease, ig);
1257 }
1258
1259 /*  Iterator block for initializer  */
1260 static VALUE
1261 s_IntGroup_Initialize_i(VALUE val, VALUE ig1)
1262 {
1263         IntGroup_RaiseIfError(IntGroupAdd((IntGroup *)ig1, NUM2INT(val), 1));
1264         return Qnil;
1265 }
1266
1267 /*
1268  *  call-seq:
1269  *     new(arg1, arg2,...)
1270  *     new(arg1, arg2,...) {|i| ...}
1271  *
1272  *  Create a new integer group. If no arguments are given, an empty group is returned.
1273  *  The arguments are either IntGroup, Range, Enumerable, or Numeric. In either case,
1274  *  the non-negative integers included in the arguments are added to the result.
1275  *  If a block is given, the block is called with the each integer in the given arguments,
1276  *  and the integer group consisting with the returned integers is returned.
1277  */
1278 static VALUE
1279 s_IntGroup_Initialize(int argc, VALUE *argv, VALUE self)
1280 {
1281         IntGroup *ig1;
1282         Data_Get_Struct(self, IntGroup, ig1);
1283         while (argc-- > 0) {
1284                 VALUE arg = *argv++;
1285                 int type = TYPE(arg);
1286                 if (rb_obj_is_kind_of(arg, rb_cIntGroup))
1287                         rb_funcall(rb_cIntGroup, rb_intern("merge"), 1, arg);
1288                 else if (rb_obj_is_kind_of(arg, rb_cRange)) {
1289                         int sp, ep;
1290                         sp = NUM2INT(rb_funcall(arg, rb_intern("begin"), 0));
1291                         ep = NUM2INT(rb_funcall(arg, rb_intern("end"), 0));
1292                         if (RTEST(rb_funcall(arg, rb_intern("exclude_end?"), 0)))
1293                                 ep--;
1294                         if (ep >= sp)
1295                                 IntGroup_RaiseIfError(IntGroupAdd(ig1, sp, ep - sp + 1));
1296                 } else if (rb_respond_to(arg, rb_intern("each")) && type != T_STRING)
1297                         rb_iterate(rb_each, arg, s_IntGroup_Initialize_i, (VALUE)ig1);
1298                 else
1299                         IntGroup_RaiseIfError(IntGroupAdd(ig1, NUM2INT(arg), 1));
1300         }
1301         if (rb_block_given_p()) {
1302                 IntGroup *ig2 = IntGroupNew();
1303                 int i, n;
1304                 for (i = 0; (n = IntGroupGetNthPoint(ig1, i)) >= 0; i++) {
1305                         n = NUM2INT(rb_yield(INT2NUM(n)));
1306                         if (n >= 0)
1307                                 IntGroup_RaiseIfError(IntGroupAdd(ig2, n, 1));
1308                 }
1309                 IntGroup_RaiseIfError(IntGroupCopy(ig1, ig2));
1310         }
1311         return Qnil;
1312 }
1313
1314 /*
1315  *  call-seq:
1316  *     clear -> self
1317  *
1318  *  Discard all integers included in self.
1319  */
1320 static VALUE
1321 s_IntGroup_Clear(VALUE self)
1322 {
1323         IntGroup *ig;
1324         Data_Get_Struct(self, IntGroup, ig);
1325         IntGroupClear(ig);
1326         return self;
1327 }
1328
1329 /*
1330  *  call-seq:
1331  *     dup(IntGroup) -> (new) IntGroup
1332  *
1333  *  (Deep) copy the given IntGroup.
1334  */
1335 static VALUE
1336 s_IntGroup_InitializeCopy(VALUE self, VALUE val)
1337 {
1338         IntGroup *ig1, *ig2;
1339         Data_Get_Struct(self, IntGroup, ig1);
1340         if (!rb_obj_is_kind_of(val, rb_cIntGroup))
1341                 rb_raise(rb_eMolbyError, "IntGroup instance is expected");
1342     Data_Get_Struct(val, IntGroup, ig2);
1343         IntGroupCopy(ig1, ig2);
1344         return self;
1345 }
1346
1347 /*
1348  *  call-seq:
1349  *     length -> Integer
1350  *     size -> Integer
1351  *
1352  *  Returns the number of integers included in self.
1353  */
1354 static VALUE
1355 s_IntGroup_Length(VALUE self)
1356 {
1357         IntGroup *ig;
1358         Data_Get_Struct(self, IntGroup, ig);
1359         return INT2NUM(IntGroupGetCount(ig));
1360 }
1361
1362 /*
1363  *  call-seq:
1364  *     member?(val) -> bool
1365  *     include?(val) -> bool
1366  *
1367  *  Check whether the val is included in self.
1368  */
1369 static VALUE
1370 s_IntGroup_MemberP(VALUE self, VALUE val)
1371 {
1372         IntGroup *ig;
1373         int n = NUM2INT(val);
1374         Data_Get_Struct(self, IntGroup, ig);
1375         return (IntGroupLookup(ig, n, NULL) ? Qtrue : Qfalse);
1376 }
1377
1378 /*
1379  *  call-seq:
1380  *     self[index] -> Integer or nil
1381  *
1382  *  Get the index-th point in self. If the index is out of range, nil is returned.
1383  */
1384 static VALUE
1385 s_IntGroup_ElementAtIndex(VALUE self, VALUE val)
1386 {
1387         IntGroup *ig;
1388         int n;
1389         int index = NUM2INT(rb_Integer(val));
1390         Data_Get_Struct(self, IntGroup, ig);
1391         n = IntGroupGetNthPoint(ig, index);
1392         return (n >= 0 ? INT2NUM(n) : Qnil);
1393 }
1394
1395 /*
1396  *  call-seq:
1397  *     each {|i| ...}
1398  *
1399  *  Call the block with each integer in self.
1400  */
1401 static VALUE
1402 s_IntGroup_Each(VALUE self)
1403 {
1404         IntGroup *ig;
1405         int i, j, sp, ep;
1406         Data_Get_Struct(self, IntGroup, ig);
1407         for (i = 0; (sp = IntGroupGetStartPoint(ig, i)) >= 0; i++) {
1408                 ep = IntGroupGetEndPoint(ig, i);
1409                 for (j = sp; j < ep; j++) {
1410                         rb_yield(INT2NUM(j));
1411                 }
1412         }
1413         return self;
1414 }
1415
1416 /*
1417  *  call-seq:
1418  *     add(IntGroup) -> self
1419  *
1420  *  Add the points in the given group.
1421  */
1422 static VALUE
1423 s_IntGroup_Add(VALUE self, VALUE val)
1424 {
1425         IntGroup *ig, *ig2;
1426     if (OBJ_FROZEN(self))
1427                 rb_error_frozen("IntGroup");
1428         Data_Get_Struct(self, IntGroup, ig);
1429         if (rb_obj_is_kind_of(val, rb_cNumeric)) {
1430                 int n = NUM2INT(rb_Integer(val));
1431                 if (n < 0)
1432                         rb_raise(rb_eMolbyError, "the integer group can contain only non-negative values");
1433                 IntGroupAdd(ig, n, 1);
1434         } else {
1435                 ig2 = IntGroupFromValue(val);
1436                 IntGroupAddIntGroup(ig, ig2);
1437                 IntGroupRelease(ig2);
1438         }
1439         return self;
1440 }
1441
1442 /*
1443  *  call-seq:
1444  *     delete(IntGroup) -> self
1445  *
1446  *  Remove the points in the given group.
1447  */
1448 static VALUE
1449 s_IntGroup_Delete(VALUE self, VALUE val)
1450 {
1451         IntGroup *ig, *ig2;
1452     if (OBJ_FROZEN(self))
1453                 rb_error_frozen("IntGroup");
1454         Data_Get_Struct(self, IntGroup, ig);
1455         if (rb_obj_is_kind_of(val, rb_cNumeric)) {
1456                 int n = NUM2INT(rb_Integer(val));
1457                 if (n >= 0 && IntGroupLookup(ig, n, NULL))
1458                         IntGroupRemove(ig, n, 1);
1459         } else {
1460                 ig2 = IntGroupFromValue(val);
1461                 IntGroupRemoveIntGroup(ig, ig2);
1462                 IntGroupRelease(ig2);
1463         }
1464         return self;
1465 }
1466
1467 static VALUE
1468 s_IntGroup_Binary(VALUE self, VALUE val, int (*func)(const IntGroup *, const IntGroup *, IntGroup *))
1469 {
1470         IntGroup *ig1, *ig2, *ig3;
1471         VALUE retval;
1472         Data_Get_Struct(self, IntGroup, ig1);
1473         ig2 = IntGroupFromValue(val);
1474         retval = IntGroup_Alloc(rb_cIntGroup);
1475         Data_Get_Struct(retval, IntGroup, ig3);
1476         IntGroup_RaiseIfError(func(ig1, ig2, ig3));
1477         IntGroupRelease(ig2);
1478         return retval;
1479 }
1480
1481 /*
1482  *  call-seq:
1483  *     union(val) -> (new)IntGroup
1484  *     self + val -> (new)IntGroup
1485  *     self | val -> (new)IntGroup
1486  *
1487  *  Returns a union group.
1488  */
1489 static VALUE
1490 s_IntGroup_Union(VALUE self, VALUE val)
1491 {
1492         return s_IntGroup_Binary(self, val, IntGroupUnion);
1493 }
1494
1495 /*
1496  *  call-seq:
1497  *     intersection(val) -> (new)IntGroup
1498  *     self & val -> (new)IntGroup
1499  *
1500  *  Returns an intersection group.
1501  */
1502 static VALUE
1503 s_IntGroup_Intersection(VALUE self, VALUE val)
1504 {
1505         return s_IntGroup_Binary(self, val, IntGroupIntersect);
1506 }
1507
1508 /*
1509  *  call-seq:
1510  *     difference(val) -> (new)IntGroup
1511  *     self - val -> (new)IntGroup
1512  *
1513  *  Returns a difference group.
1514  */
1515 static VALUE
1516 s_IntGroup_Difference(VALUE self, VALUE val)
1517 {
1518         return s_IntGroup_Binary(self, val, IntGroupDifference);
1519 }
1520
1521 /*
1522  *  call-seq:
1523  *     sym_difference(val) -> (new)IntGroup
1524  *     self ^ val -> (new)IntGroup
1525  *
1526  *  Returns a symmetric-difference group (i.e. a group containing elements that are included
1527  *  in either self or val but not both).
1528  */
1529 static VALUE
1530 s_IntGroup_SymDifference(VALUE self, VALUE val)
1531 {
1532         return s_IntGroup_Binary(self, val, IntGroupXor);
1533 }
1534
1535 /*
1536  *  call-seq:
1537  *     convolute(val) -> (new)IntGroup
1538  *
1539  *  For each element n in self, get the n-th point in val, and return the result as a new group.
1540  *  If n is out of range, then that point is ignored.
1541  *
1542  *  <b>See Also:</b> IntGroup#deconvolute. If all points in self are within the range of val, then
1543  *  an equation <tt>self.convolute(val).deconvolute(val) == self</tt> holds.
1544  */
1545 static VALUE
1546 s_IntGroup_Convolute(VALUE self, VALUE val)
1547 {
1548         return s_IntGroup_Binary(self, val, IntGroupConvolute);
1549 }
1550
1551 /*
1552  *  call-seq:
1553  *     deconvolute(val) -> (new)IntGroup
1554  *
1555  *  For each element n in self, find the point n in val, and return the found indices as a new group.
1556  *  If n is not found in val, then that point is ignored.
1557  *
1558  *  <b>See Also:</b> IntGroup#convolute. If all points in self are found in val, then
1559  *  an equation <tt>self.deconvolute(val).convolute(val) == self</tt> holds.
1560  */
1561 static VALUE
1562 s_IntGroup_Deconvolute(VALUE self, VALUE val)
1563 {
1564         return s_IntGroup_Binary(self, val, IntGroupDeconvolute);
1565 }
1566
1567 /*
1568  *  call-seq:
1569  *     range_at(val) -> Range
1570  *
1571  *  Split self into consecutive chunks of integers, and return the val-th chunk as a Range.
1572  *  This method is relatively efficient, because it directly uses the internal representation
1573  *  of IntGroup.
1574  */
1575 static VALUE
1576 s_IntGroup_RangeAt(VALUE self, VALUE val)
1577 {
1578         IntGroup *ig;
1579         int n = NUM2INT(val);
1580         int sp, ep;
1581         Data_Get_Struct(self, IntGroup, ig);
1582         sp = IntGroupGetStartPoint(ig, n);
1583         if (sp < 0)
1584                 return Qnil;
1585         ep = IntGroupGetEndPoint(ig, n) - 1;
1586         return rb_funcall(rb_cRange, rb_intern("new"), 2, INT2NUM(sp), INT2NUM(ep));
1587 }
1588
1589 /*
1590 static VALUE
1591 s_IntGroup_Merge(VALUE self, VALUE val)
1592 {
1593         IntGroup *ig1, *ig2;
1594         int i, sp, interval;
1595     if (OBJ_FROZEN(self))
1596                 rb_error_frozen("IntGroup");
1597         Data_Get_Struct(self, IntGroup, ig1);
1598         ig2 = IntGroupFromValue(val);
1599         for (i = 0; (sp = IntGroupGetStartPoint(ig2, i)) >= 0; i++) {
1600                 interval = IntGroupGetInterval(ig2, i);
1601                 IntGroup_RaiseIfError(IntGroupAdd(ig1, sp, interval));
1602         }
1603         IntGroupRelease(ig2);
1604         return self;
1605 }
1606
1607 static VALUE
1608 s_IntGroup_Subtract(VALUE self, VALUE val)
1609 {
1610         IntGroup *ig1, *ig2;
1611         int i, sp, interval;
1612     if (OBJ_FROZEN(self))
1613                 rb_error_frozen("IntGroup");
1614         Data_Get_Struct(self, IntGroup, ig1);
1615         ig2 = IntGroupFromValue(val);
1616         for (i = 0; (sp = IntGroupGetStartPoint(ig2, i)) >= 0; i++) {
1617                 interval = IntGroupGetInterval(ig2, i);
1618                 IntGroup_RaiseIfError(IntGroupRemove(ig1, sp, interval));
1619         }
1620         IntGroupRelease(ig2);
1621         return self;
1622 }
1623 */
1624
1625 /*
1626  *  call-seq:
1627  *     offset(val) -> self
1628  *
1629  *  Move all points by an integer value. A negative val is allowed, but it
1630  *  must be no smaller than -(self[0]), otherwise an exception is thrown.
1631  */
1632 static VALUE
1633 s_IntGroup_Offset(VALUE self, VALUE ofs)
1634 {
1635         IntGroup *ig1, *ig2;
1636         int iofs;
1637         VALUE val;
1638         Data_Get_Struct(self, IntGroup, ig1);
1639         ig2 = IntGroupNewFromIntGroup(ig1);
1640         if (ig2 == NULL)
1641                 rb_raise(rb_eMolbyError, "Cannot duplicate IntGroup");
1642         iofs = NUM2INT(ofs);
1643         if (IntGroupOffset(ig2, iofs) != 0)
1644                 rb_raise(rb_eMolbyError, "Bad offset %d", iofs);
1645         val = ValueFromIntGroup(ig2);
1646         IntGroupRelease(ig2);
1647         return val;
1648 }
1649
1650 static VALUE
1651 s_IntGroup_Create(int argc, VALUE *argv, VALUE klass)
1652 {
1653         VALUE val = IntGroup_Alloc(klass);
1654         s_IntGroup_Initialize(argc, argv, val);
1655         return val;
1656 }
1657
1658 /*
1659  *  call-seq:
1660  *     inspect -> String
1661  *
1662  *  Create a String in the form "IntGroup[...]".
1663  */
1664 static VALUE
1665 s_IntGroup_Inspect(VALUE self)
1666 {
1667         int i, sp, ep;
1668         IntGroup *ig;
1669         char buf[64];
1670 /*      VALUE klass = CLASS_OF(self);
1671         VALUE val = rb_funcall(klass, rb_intern("name"), 0); */
1672 /*      rb_str_cat(val, "[", 1); */
1673         VALUE val = rb_str_new2("IntGroup[");
1674         Data_Get_Struct(self, IntGroup, ig);
1675         for (i = 0; (sp = IntGroupGetStartPoint(ig, i)) >= 0; i++) {
1676                 if (i > 0)
1677                         rb_str_cat(val, ", ", 2);
1678                 ep = IntGroupGetEndPoint(ig, i);
1679                 if (ep > sp + 1)
1680                         snprintf(buf, sizeof buf, "%d..%d", sp, ep - 1);
1681                 else
1682                         snprintf(buf, sizeof buf, "%d", sp);
1683                 rb_str_cat(val, buf, strlen(buf));
1684         }
1685         rb_str_cat(val, "]", 1);
1686         return val;
1687 }
1688
1689 void
1690 Init_MolbyTypes(void)
1691 {
1692         /*  class Vector3D  */
1693         rb_cVector3D = rb_define_class_under(rb_mMolby, "Vector3D", rb_cObject);
1694         rb_define_alloc_func(rb_cVector3D, s_Vector3D_Alloc);
1695         rb_define_method(rb_cVector3D, "initialize", s_Vector3D_Initialize, -1);
1696         rb_define_method(rb_cVector3D, "size", s_Vector3D_Size, 0);
1697         rb_define_method(rb_cVector3D, "[]", s_Vector3D_ElementAtIndex, 1);
1698         rb_define_method(rb_cVector3D, "[]=", s_Vector3D_SetElementAtIndex, 2);
1699         rb_define_method(rb_cVector3D, "==", s_Vector3D_IsEqual, 1);
1700         rb_define_method(rb_cVector3D, "+", s_Vector3D_Add, 1);
1701         rb_define_method(rb_cVector3D, "-", s_Vector3D_Subtract, 1);
1702         rb_define_method(rb_cVector3D, "*", s_Vector3D_Multiply, 1);
1703         rb_define_method(rb_cVector3D, "/", s_Vector3D_Divide, 1);
1704         rb_define_method(rb_cVector3D, "dot", s_Vector3D_Dot, 1);
1705         rb_define_method(rb_cVector3D, "cross", s_Vector3D_Cross, 1);
1706         rb_define_method(rb_cVector3D, "-@", s_Vector3D_UnaryMinus, 0);
1707         rb_define_method(rb_cVector3D, "length", s_Vector3D_Length, 0);
1708         rb_define_alias(rb_cVector3D, "r", "length"); /* size and length are not synonym!  */
1709         rb_define_method(rb_cVector3D, "length2", s_Vector3D_Length2, 0);
1710         rb_define_alias(rb_cVector3D, "r2", "length2");
1711         rb_define_method(rb_cVector3D, "normalize", s_Vector3D_Normalize, 0);
1712         rb_define_method(rb_cVector3D, "to_a", s_Vector3D_ToArray, 0);
1713         rb_define_method(rb_cVector3D, "each", s_Vector3D_Each, 0);
1714         rb_define_method(rb_cVector3D, "x", s_Vector3D_GetX, 0);
1715         rb_define_method(rb_cVector3D, "y", s_Vector3D_GetY, 0);
1716         rb_define_method(rb_cVector3D, "z", s_Vector3D_GetZ, 0);
1717         rb_define_method(rb_cVector3D, "x=", s_Vector3D_SetX, 1);
1718         rb_define_method(rb_cVector3D, "y=", s_Vector3D_SetY, 1);
1719         rb_define_method(rb_cVector3D, "z=", s_Vector3D_SetZ, 1);
1720         rb_define_method(rb_cVector3D, "inspect", s_Vector3D_Inspect, 0);
1721         rb_define_alias(rb_cVector3D, "to_s", "inspect");
1722         rb_define_singleton_method(rb_cVector3D, "[]", s_Vector3D_Create, -2);
1723
1724         /*  class Transform  */
1725         rb_cTransform = rb_define_class_under(rb_mMolby, "Transform", rb_cObject);
1726         rb_define_alloc_func(rb_cTransform, s_Transform_Alloc);
1727         rb_define_method(rb_cTransform, "initialize", s_Transform_Initialize, -1);
1728         rb_define_method(rb_cTransform, "[]", s_Transform_ElementAtIndex, 2);
1729         rb_define_method(rb_cTransform, "[]=", s_Transform_SetElementAtIndex, 3);
1730         rb_define_method(rb_cTransform, "==", s_Transform_IsEqual, 1);
1731         rb_define_method(rb_cTransform, "+", s_Transform_Add, 1);
1732         rb_define_method(rb_cTransform, "-", s_Transform_Subtract, 1);
1733         rb_define_method(rb_cTransform, "*", s_Transform_Multiply, 1);
1734         rb_define_method(rb_cTransform, "/", s_Transform_Divide, 1);
1735         rb_define_method(rb_cTransform, "inverse", s_Transform_Inverse, 0);
1736         rb_define_method(rb_cTransform, "transpose", s_Transform_Transpose, 0);
1737         rb_define_method(rb_cTransform, "determinant", s_Transform_Determinant, 0);
1738         rb_define_method(rb_cTransform, "trace", s_Transform_Trace, 0);
1739         rb_define_method(rb_cTransform, "column", s_Transform_Column, 1);
1740         rb_define_method(rb_cTransform, "eigenvalues", s_Transform_Eigenvalues, 0);
1741         rb_define_method(rb_cTransform, "to_a", s_Transform_ToArray, 0);
1742         rb_define_method(rb_cTransform, "inspect", s_Transform_Inspect, 0);
1743         rb_define_alias(rb_cTransform, "to_s", "inspect");
1744         rb_define_singleton_method(rb_cTransform, "diagonal", s_Transform_Diagonal, -1);
1745         rb_define_singleton_method(rb_cTransform, "[]", s_Transform_Create, -2);
1746         rb_define_singleton_method(rb_cTransform, "from_columns", s_Transform_NewFromColumns, -2);
1747         rb_define_singleton_method(rb_cTransform, "from_rows", s_Transform_NewFromRows, -2);
1748         rb_define_singleton_method(rb_cTransform, "identity", s_Transform_Identity, 0);
1749         rb_define_singleton_method(rb_cTransform, "zero", s_Transform_Zero, 0);
1750         rb_define_singleton_method(rb_cTransform, "translation", s_Transform_Translation, 1);
1751         rb_define_singleton_method(rb_cTransform, "rotation", s_Transform_Rotation, -1);
1752         rb_define_singleton_method(rb_cTransform, "reflection", s_Transform_Reflection, -1);
1753         rb_define_singleton_method(rb_cTransform, "inversion", s_Transform_Inversion, -1);
1754
1755         /*  class IntGroup  */
1756         rb_cIntGroup = rb_define_class_under(rb_mMolby, "IntGroup", rb_cObject);
1757         rb_include_module(rb_cIntGroup, rb_mEnumerable);
1758         rb_define_alloc_func(rb_cIntGroup, IntGroup_Alloc);
1759         rb_define_method(rb_cIntGroup, "clear", s_IntGroup_Clear, 0);
1760         rb_define_method(rb_cIntGroup, "initialize", s_IntGroup_Initialize, -1);
1761         rb_define_method(rb_cIntGroup, "initialize_copy", s_IntGroup_InitializeCopy, 1);
1762         rb_define_method(rb_cIntGroup, "length", s_IntGroup_Length, 0);
1763         rb_define_alias(rb_cIntGroup, "size", "length");
1764         rb_define_method(rb_cIntGroup, "member?", s_IntGroup_MemberP, 1);
1765         rb_define_alias(rb_cIntGroup, "include?", "member?");
1766         rb_define_method(rb_cIntGroup, "each", s_IntGroup_Each, 0);
1767         rb_define_method(rb_cIntGroup, "[]", s_IntGroup_ElementAtIndex, 1);
1768         rb_define_method(rb_cIntGroup, "add", s_IntGroup_Add, 1);
1769         rb_define_alias(rb_cIntGroup, "<<", "add");
1770         rb_define_method(rb_cIntGroup, "delete", s_IntGroup_Delete, 1);
1771         rb_define_method(rb_cIntGroup, "union", s_IntGroup_Union, 1);
1772         rb_define_method(rb_cIntGroup, "difference", s_IntGroup_Difference, 1);
1773         rb_define_method(rb_cIntGroup, "intersection", s_IntGroup_Intersection, 1);
1774         rb_define_method(rb_cIntGroup, "sym_difference", s_IntGroup_SymDifference, 1);
1775         rb_define_method(rb_cIntGroup, "convolute", s_IntGroup_Convolute, 1);
1776         rb_define_method(rb_cIntGroup, "deconvolute", s_IntGroup_Deconvolute, 1);
1777         rb_define_method(rb_cIntGroup, "offset", s_IntGroup_Offset, 1);
1778         rb_define_alias(rb_cIntGroup, "+", "union");
1779         rb_define_alias(rb_cIntGroup, "|", "union");
1780         rb_define_alias(rb_cIntGroup, "-", "difference");
1781         rb_define_alias(rb_cIntGroup, "&", "intersection");
1782         rb_define_alias(rb_cIntGroup, "^", "sym_difference");
1783         rb_define_method(rb_cIntGroup, "range_at", s_IntGroup_RangeAt, 1);
1784 /*      rb_define_method(rb_cIntGroup, "merge", s_IntGroup_Merge, -1);
1785         rb_define_method(rb_cIntGroup, "subtract", s_IntGroup_Subtract, -1); */
1786         rb_define_method(rb_cIntGroup, "inspect", s_IntGroup_Inspect, 0);
1787         rb_define_alias(rb_cIntGroup, "to_s", "inspect");
1788         rb_define_singleton_method(rb_cIntGroup, "[]", s_IntGroup_Create, -1);
1789         {
1790                 VALUE igval = IntGroup_Alloc(rb_cIntGroup);
1791                 IntGroup *ig;
1792                 Data_Get_Struct(igval, IntGroup, ig);
1793                 IntGroupAdd(ig, 0, ATOMS_MAX_NUMBER);
1794                 rb_define_global_const("All", igval);
1795         }
1796 }