Copied from rough/bigdecimal,documents & some sample programs added.

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@3625 b2dd03c8-39d4-4d8f-98ff-823fe69b080e

17 files changed, 6325 insertions, 0 deletions

diff --git a/ext/bigdecimal/MANIFEST b/ext/bigdecimal/MANIFEST
new file mode 100644
index 0000000000..375f16c2ab
--- /dev/null
+++ b/ext/bigdecimal/MANIFEST
@@ -0,0 +1,16 @@
+MANIFEST
+README
+bigdecimal.def
+bigdecimal.c
+bigdecimal.h
+depend
+extconf.rb
+bigdecimal_en.html
+bigdecimal_ja.html
+lib/nlsolve.rb
+lib/jacobian.rb
+lib/newton.rb
+lib/linear.rb
+lib/ludcmp.rb
+lib/pai.rb
+lib/bigdecimal-rational.rb
diff --git a/ext/bigdecimal/README b/ext/bigdecimal/README
new file mode 100644
index 0000000000..a233f47f64
--- /dev/null
+++ b/ext/bigdecimal/README
@@ -0,0 +1,60 @@
+
+  Ruby BIGDECIMAL(Variable Precision) extension library.
+  Copyright (C) 1999  by Shigeo Kobayashi([email protected])
+
+BigDecimal is copyrighted free software by Shigeo Kobayashi <[email protected]>.
+You can redistribute it and/or modify it under either the terms of the GPL
+(see COPYING file), or the conditions below:
+
+  1. You may make and give away verbatim copies of the source form of the
+     software without restriction, provided that you duplicate all of the
+     original copyright notices and associated disclaimers.
+
+  2. You may modify your copy of the software in any way, provided that
+     you do at least ONE of the following:
+
+       a) place your modifications in the Public Domain or otherwise
+          make them Freely Available, such as by posting said
+	  modifications to Usenet or an equivalent medium, or by allowing
+	  the author to include your modifications in the software.
+
+       b) use the modified software only within your corporation or
+          organization.
+
+       c) rename any non-standard executables so the names do not conflict
+	  with standard executables, which must also be provided.
+
+       d) make other distribution arrangements with the author.
+
+  3. You may distribute the software in object code or executable
+     form, provided that you do at least ONE of the following:
+
+       a) distribute the executables and library files of the software,
+	  together with instructions (in the manual page or equivalent)
+	  on where to get the original distribution.
+
+       b) accompany the distribution with the machine-readable source of
+	  the software.
+
+       c) give non-standard executables non-standard names, with
+          instructions on where to get the original software distribution.
+
+       d) make other distribution arrangements with the author.
+
+  4. You may modify and include the part of the software into any other
+     software (possibly commercial).
+
+  5. THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
+     IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
+     WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
+     PURPOSE.
+
+* The Author
+
+Feel free to send comments and bug reports to the author.  Here is the 
+author's latest mail address:
+
+  [email protected]
+
+-------------------------------------------------------
+created at: Thu Dec 22 1999
diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c
new file mode 100644
index 0000000000..3dd50921eb
--- /dev/null
+++ b/ext/bigdecimal/bigdecimal.c
@@ -0,0 +1,4198 @@
+/*
+ *
+ * Ruby BigDecimal(Variable decimal precision) extension library.
+ *
+ * Copyright(C) 2002 by Shigeo Kobayashi([email protected])
+ *
+ * You may distribute under the terms of either the GNU General Public
+ * License or the Artistic License, as specified in the README file
+ * of this BigDecimal distribution.
+ *
+ *  BigDecimal was originally copied from BigFloat(http://www.tinyforest.gr.jp/ruby/bigfloat_en.html),
+ *  and has been modified so that it is fit in Ruby style.
+ *  (Apr 2,2002).
+ *
+ *  2003-01-03
+ *    assign instead of asign(by knu),use string.h functions(by t.saito).
+ *  2002-12-06
+ *    The sqrt() bug reported by Bret Jolly fixed.
+ *  2002-5-6
+ *    The bug reported by Sako Hiroshi <[email protected]>(ruby-list:34988)
+ *    in to_i fixed.
+ *  2002-4-17
+ *    methods prec and double_fig(class method) added(S.K).
+ *
+ *  Version 1.0.0(2002/04/04)
+ *    Copied from BigFloat 1.1.9 and
+ *      hash method changed according to the suggestion from Akinori MUSHA <[email protected]>.
+ *      All ! class methods deactivated(but not actually removed).
+ *      to_s and to_s2 merged to one to_s[(n)].
+ *
+ */
+
+#include "ruby.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#ifdef NT
+#include <malloc.h>
+#ifdef _MSC_VER
+#include <float.h>
+#define isnan(x) _isnan(x)
+#define isinf(x) (!(_finite(x)))
+#endif /* _MSC_VER */
+#endif /* defined NT */
+#include "ruby.h"
+#include "math.h"
+#include "version.h"
+
+/* #define USE_MUTABLE_METHOD */
+
+VALUE rb_cBigDecimal;
+
+#include "bigdecimal.h"
+
+/* MACRO's to guard objects from GC by keeping it in stack */
+#define ENTER(n) volatile VALUE vStack[n];int iStack=0
+#define PUSH(x)  vStack[iStack++] = (unsigned long)(x);
+#define SAVE(p)  PUSH(p->obj);
+#define GUARD_OBJ(p,y) {p=y;SAVE(p);}
+
+/* ETC */
+#define MemCmp(x,y,z) memcmp(x,y,z)
+#define StrCmp(x,y)   strcmp(x,y)
+
+static int VpIsDefOP(Real *c,Real *a,Real *b,int sw);
+static int AddExponent(Real *a,S_INT n);
+static unsigned short VpGetException(void);
+static void  VpSetException(unsigned short f);
+static int VpAddAbs(Real *a,Real *b,Real *c);
+static int VpSubAbs(Real *a,Real *b,Real *c);
+static U_LONG VpSetPTR(Real *a,Real *b,Real *c,U_LONG *a_pos,U_LONG *b_pos,U_LONG *c_pos,U_LONG *av,U_LONG *bv);
+static int VpNmlz(Real *a);
+static void VpFormatSt(char *psz,S_INT fFmt);
+static int VpRdup(Real *m);
+static U_LONG SkipWhiteChar(char *szVal);
+
+/*
+ * ================== Ruby Interface part ==========================
+ */
+static ID coerce;
+
+/* Following functions borrowed from numeric.c */
+static VALUE
+coerce_body(VALUE *x)
+{
+    return rb_funcall(x[1], coerce, 1, x[0]);
+}
+
+static VALUE
+coerce_rescue(VALUE *x)
+{
+    rb_raise(rb_eTypeError, "%s can't be coerced into %s",
+        rb_special_const_p(x[1])?
+        STR2CSTR(rb_inspect(x[1])):
+        rb_class2name(CLASS_OF(x[1])),
+        rb_class2name(CLASS_OF(x[0])));
+    return (VALUE)0;
+}
+
+static void
+do_coerce(VALUE *x, VALUE *y)
+{
+    VALUE ary;
+    VALUE a[2];
+
+    a[0] = *x; a[1] = *y;
+    ary = rb_rescue(coerce_body, (VALUE)a, coerce_rescue, (VALUE)a);
+    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];
+}
+
+static VALUE
+DoSomeOne(VALUE x, VALUE y)
+{
+    do_coerce(&x, &y);
+    return rb_funcall(x, rb_frame_last_func(), 1, y);
+}
+
+static void
+BigDecimal_delete(Real *pv)
+{
+    VpFree(pv);
+}
+
+static VALUE
+ToValue(Real *p)
+{
+    if(VpIsNaN(p)) {
+        VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'(Not a Number)",0);
+    } else if(VpIsPosInf(p)) {
+        VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",0);
+    } else if(VpIsNegInf(p)) {
+        VpException(VP_EXCEPTION_INFINITY,"Computation results to '-Infinity'",0);
+    }
+    return p->obj;
+}
+
+static Real *
+GetVpValue(VALUE v, int must)
+{
+    double dv;
+    Real *pv;
+    VALUE bg;
+    char szD[128];
+
+    switch(TYPE(v))
+    {
+    case T_DATA:
+        if(RDATA(v)->dfree ==(void *) BigDecimal_delete) {
+            Data_Get_Struct(v, Real, pv);
+            return pv;
+        } else {
+            goto SomeOneMayDoIt;
+        }
+        break;
+    case T_FIXNUM:
+        sprintf(szD, "%d", NUM2INT(v));
+        return VpCreateRbObject(VpBaseFig() * 2 + 1, szD);
+    case T_FLOAT:
+        pv = VpCreateRbObject(VpDblFig()*2,"0");
+        dv = RFLOAT(v)->value;
+        /* From float */
+        if (isinf(dv)) {
+            VpException(VP_EXCEPTION_INFINITY,"Computation including infinity",0);
+            if(dv==VpGetDoublePosInf()) {
+                VpSetPosInf(pv);
+            } else {
+                VpSetNegInf(pv);
+            }
+        } else
+        if (isnan(dv)) {
+            VpException(VP_EXCEPTION_NaN,"Computation including NaN(Not a number)",0);
+            VpSetNaN(pv);
+        } else {
+            if (VpIsNegDoubleZero(dv)) {
+                VpSetNegZero(pv);
+            } else if(dv==0.0) {
+                VpSetPosZero(pv);
+            } else if(dv==1.0) {
+                VpSetOne(pv);
+            } else if(dv==-1.0) {
+                VpSetOne(pv);
+                pv->sign = -pv->sign;
+            } else {
+                VpDtoV(pv,dv);
+            }
+        }
+        return pv;
+    case T_STRING:
+        Check_SafeStr(v);
+        return VpCreateRbObject(strlen(RSTRING(v)->ptr) + VpBaseFig() + 1,
+                                RSTRING(v)->ptr);
+    case T_BIGNUM:
+        bg = rb_big2str(v, 10);
+        return VpCreateRbObject(strlen(RSTRING(bg)->ptr) + VpBaseFig() + 1,
+                                RSTRING(bg)->ptr);
+    default:
+        goto SomeOneMayDoIt;
+    }
+
+SomeOneMayDoIt:
+    if(must) {
+        rb_raise(rb_eTypeError, "%s can't be coerced into BigDecimal",
+                    rb_special_const_p(v)?
+                    STR2CSTR(rb_inspect(v)):
+                    rb_class2name(CLASS_OF(v))
+                );
+    }
+    return NULL; /* NULL means to coerce */
+}
+
+static VALUE
+BigDecimal_double_fig(VALUE self)
+{
+    return INT2FIX(VpDblFig());
+}
+
+static VALUE
+BigDecimal_prec(VALUE self)
+{
+    ENTER(1);
+    Real *p;
+    VALUE obj;
+
+    GUARD_OBJ(p,GetVpValue(self,1));
+    obj = rb_ary_new();
+    obj = rb_ary_push(obj,INT2NUM(p->Prec*VpBaseFig()));
+    obj = rb_ary_push(obj,INT2NUM(p->MaxPrec*VpBaseFig()));
+    return obj;
+}
+
+static VALUE
+BigDecimal_hash(VALUE self)
+{
+    ENTER(1);
+    Real *p;
+    U_LONG hash,i;
+
+    GUARD_OBJ(p,GetVpValue(self,1));
+    hash = (U_LONG)p->sign;
+    /* hash!=2: the case for 0(1),NaN(0) or +-Infinity(3) is sign itself */
+    if(hash==2) {
+        for(i = 0; i < p->Prec;i++) {
+            hash = 31 * hash + p->frac[i];
+            hash ^= p->frac[i];
+        }
+        hash += p->exponent;
+    }
+    return INT2FIX(hash);
+}
+
+static VALUE
+BigDecimal_dump(int argc, VALUE *argv, VALUE self)
+{
+    ENTER(5);
+    char sz[50];
+    Real *vp;
+    char *psz;
+    VALUE dummy;
+    rb_scan_args(argc, argv, "01", &dummy);
+    GUARD_OBJ(vp,GetVpValue(self,1));
+    sprintf(sz,"%d:",VpMaxPrec(vp)*VpBaseFig());
+    psz = ALLOCA_N(char,(unsigned int)VpNumOfChars(vp)+strlen(sz));
+    sprintf(psz,"%s",sz);
+    VpToString(vp, psz+strlen(psz), 0);
+    return rb_str_new2(psz);
+}
+
+static VALUE
+BigDecimal_load(VALUE self, VALUE str)
+{
+    ENTER(2);
+    Real *pv;
+    long len;
+    unsigned char *pch;
+    unsigned char ch;
+    unsigned long m=0;
+
+    Check_SafeStr(str);
+    pch = rb_str2cstr(str, &len);
+    /* First get max prec */
+    while((*pch)!=(unsigned char)'\0' && (ch=*pch++)!=(unsigned char)':') {
+        if(ch<'0' || ch>'9') {
+            rb_raise(rb_eTypeError, "Load failed: invalid character in the marshaled string");
+        }
+        m = m*10 + (unsigned long)(ch-'0');
+    }
+    if(m>VpBaseFig()) m -= VpBaseFig();
+    GUARD_OBJ(pv,VpNewRbClass(m,pch,self));
+    m /= VpBaseFig();
+    if(m && pv->MaxPrec>m) pv->MaxPrec = m+1;
+    return ToValue(pv);
+}
+
+static VALUE
+BigDecimal_mode(VALUE self, VALUE which, VALUE val)
+{
+    unsigned short fo = VpGetException();
+    unsigned short f;
+
+    if(TYPE(which)!=T_FIXNUM)  return INT2FIX(fo);
+    if(val!=Qfalse && val!=Qtrue) return INT2FIX(fo);
+
+    f = (unsigned short)NUM2INT(which);
+    if(f&VP_EXCEPTION_INFINITY) {
+        fo = VpGetException();
+        VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_INFINITY):
+                         (fo&(~VP_EXCEPTION_INFINITY))));
+    }
+    if(f&VP_EXCEPTION_NaN) {
+        fo = VpGetException();
+        VpSetException((unsigned short)((val==Qtrue)?(fo|VP_EXCEPTION_NaN):
+                         (fo&(~VP_EXCEPTION_NaN))));
+    }
+    fo = VpGetException();
+    return INT2FIX(fo);
+}
+
+static U_LONG
+GetAddSubPrec(Real *a, Real *b)
+{
+    U_LONG mxs;
+    U_LONG mx = a->Prec;
+    S_INT d;
+
+    if(!VpIsDef(a) || !VpIsDef(b)) return (-1L);
+    if(mx < b->Prec) mx = b->Prec;
+    if(a->exponent!=b->exponent) {
+        mxs = mx;
+        d = a->exponent - b->exponent;
+        if(d<0) d = -d;
+        mx = mx+(U_LONG)d;
+        if(mx<mxs) {
+            return VpException(VP_EXCEPTION_INFINITY,"Exponent overflow",0);
+        }
+    }
+    return mx;
+}
+
+static S_INT
+GetPositiveInt(VALUE v)
+{
+    S_INT n;
+    Check_Type(v, T_FIXNUM);
+    n = NUM2INT(v);
+    if(n <= 0) {
+        rb_fatal("Zero or negative argument not permitted.");
+    }
+    return n;
+}
+
+VP_EXPORT Real *
+VpNewRbClass(U_LONG mx, char *str, VALUE klass)
+{
+    Real *pv = VpAlloc(mx,str);
+    pv->obj = (VALUE)Data_Wrap_Struct(klass, 0, BigDecimal_delete, pv);
+    return pv;
+}
+
+VP_EXPORT Real *
+VpCreateRbObject(U_LONG mx, char *str)
+{
+    Real *pv = VpAlloc(mx,str);
+    pv->obj = (VALUE)Data_Wrap_Struct(rb_cBigDecimal, 0, BigDecimal_delete, pv);
+    return pv;
+}
+
+
+static VALUE
+BigDecimal_IsNaN(VALUE self)
+{
+    Real *p = GetVpValue(self,1);
+    if(VpIsNaN(p))  return Qtrue;
+    return Qfalse;
+}
+
+static VALUE
+BigDecimal_IsInfinite(VALUE self)
+{
+    Real *p = GetVpValue(self,1);
+    if(VpIsInf(p)) return Qtrue;
+    return Qfalse;
+}
+
+static VALUE
+BigDecimal_IsFinite(VALUE self)
+{
+    Real *p = GetVpValue(self,1);
+    if(VpIsNaN(p)) return Qfalse;
+    if(VpIsInf(p)) return Qfalse;
+    return Qtrue;
+}
+
+static VALUE
+BigDecimal_to_i(VALUE self)
+{
+    ENTER(5);
+    int e,n,i,nf;
+    U_LONG v,b,j;
+    char *psz,*pch;
+    Real *p;
+
+    GUARD_OBJ(p,GetVpValue(self,1));
+
+    if(!VpIsDef(p)) return Qnil; /* Infinity or NaN not converted. */
+
+    e = VpExponent10(p);
+    if(e<=0) return INT2FIX(0);
+    nf = VpBaseFig();
+    if(e<=nf) {
+        e = VpGetSign(p)*p->frac[0];
+        return INT2FIX(e);
+    }
+    psz = ALLOCA_N(char,(unsigned int)(e+nf+2));
+
+    n = (e+nf-1)/nf;
+    pch = psz;
+    if(VpGetSign(p)<0) *pch++ = '-';
+    for(i=0;i<n;++i) {
+        b = VpBaseVal()/10;
+        if(i>=(int)p->Prec) {
+            while(b) {
+                *pch++ = '0';
+                b /= 10;
+            }
+            continue;
+        }
+        v = p->frac[i];
+        while(b) {
+            j = v/b;
+            *pch++ = (char)(j + '0');
+            v -= j*b;
+            b /= 10;
+        }
+    }
+    *pch++ = 0;
+    return rb_cstr2inum(psz,10);
+}
+
+static VALUE
+BigDecimal_induced_from(VALUE self, VALUE x)
+{
+    Real *p = GetVpValue(x,1);
+    return p->obj;
+}
+
+static VALUE
+BigDecimal_coerce(VALUE self, VALUE other)
+{
+    ENTER(2);
+    VALUE obj;
+    Real *b;
+    GUARD_OBJ(b,GetVpValue(other,1));
+    obj = rb_ary_new();
+    obj = rb_ary_push(obj, b->obj);
+    obj = rb_ary_push(obj, self);
+    return obj;
+}
+
+static VALUE
+BigDecimal_uplus(VALUE self)
+{
+    return self;
+}
+
+static VALUE
+BigDecimal_add(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *c, *a, *b;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return DoSomeOne(self,r);
+    SAVE(b);
+
+    if(VpIsNaN(b)) return b->obj;
+    if(VpIsNaN(a)) return a->obj;
+    mx = GetAddSubPrec(a,b);
+    if(mx==(-1L)) {
+        GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0"));
+        VpAddSub(c, a, b, 1);
+    } else {
+        GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0"));
+        if(!mx) {
+            VpSetInf(c,VpGetSign(a));
+        } else {
+            VpAddSub(c, a, b, 1);
+        }
+    }
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_sub(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *c, *a, *b;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return DoSomeOne(self,r);
+    SAVE(b);
+
+    if(VpIsNaN(b)) return b->obj;
+    if(VpIsNaN(a)) return a->obj;
+
+    mx = GetAddSubPrec(a,b);
+    if(mx==(-1L)) {
+        GUARD_OBJ(c,VpCreateRbObject(VpBaseFig() + 1, "0"));
+        VpAddSub(c, a, b, -1);
+    } else {
+        GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0"));
+        if(!mx) {
+            VpSetInf(c,VpGetSign(a));
+        } else {
+            VpAddSub(c, a, b, -1);
+        }
+    }
+    return ToValue(c);
+}
+
+static S_INT
+BigDecimalCmp(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *a, *b;
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return DoSomeOne(self,r);
+    SAVE(b);
+    return VpComp(a, b);
+}
+
+static VALUE
+BigDecimal_zero(VALUE self)
+{
+    Real *a = GetVpValue(self,1);
+    return VpIsZero(a) ? Qtrue : Qfalse;
+}
+
+static VALUE
+BigDecimal_nonzero(VALUE self)
+{
+    Real *a = GetVpValue(self,1);
+    return VpIsZero(a) ? Qfalse : self;
+}
+
+static VALUE
+BigDecimal_comp(VALUE self, VALUE r)
+{
+    S_INT e;
+    e = BigDecimalCmp(self, r);
+    if(e==999) return rb_float_new(VpGetDoubleNaN());
+    return INT2FIX(e);
+}
+
+static VALUE
+BigDecimal_eq(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *a, *b;
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return Qfalse; /* Not comparable */
+    SAVE(b);
+    return VpComp(a, b)? Qfalse:Qtrue;
+}
+
+static VALUE
+BigDecimal_ne(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *a, *b;
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return Qtrue; /* Not comparable */
+    SAVE(b);
+    return VpComp(a, b) ? Qtrue : Qfalse;
+}
+
+static VALUE
+BigDecimal_lt(VALUE self, VALUE r)
+{
+    S_INT e;
+    e = BigDecimalCmp(self, r);
+    if(e==999) return Qfalse;
+    return(e < 0) ? Qtrue : Qfalse;
+}
+
+static VALUE
+BigDecimal_le(VALUE self, VALUE r)
+{
+    S_INT e;
+    e = BigDecimalCmp(self, r);
+    if(e==999) return Qfalse;
+    return(e <= 0) ? Qtrue : Qfalse;
+}
+
+static VALUE
+BigDecimal_gt(VALUE self, VALUE r)
+{
+    S_INT e;
+    e = BigDecimalCmp(self, r);
+    if(e==999) return Qfalse;
+    return(e > 0) ? Qtrue : Qfalse;
+}
+
+static VALUE
+BigDecimal_ge(VALUE self, VALUE r)
+{
+    S_INT e;
+    e = BigDecimalCmp(self, r);
+    if(e==999) return Qfalse;
+    return(e >= 0) ? Qtrue : Qfalse;
+}
+
+static VALUE
+BigDecimal_neg(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *c, *a;
+    GUARD_OBJ(a,GetVpValue(self,1));
+    GUARD_OBJ(c,VpCreateRbObject(a->Prec *(VpBaseFig() + 1), "0"));
+    VpAsgn(c, a, -1);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_mult(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *c, *a, *b;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return DoSomeOne(self,r);
+    SAVE(b);
+
+    mx = a->Prec + b->Prec;
+    GUARD_OBJ(c,VpCreateRbObject(mx *(VpBaseFig() + 1), "0"));
+    VpMult(c, a, b);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_divide(Real **c, Real **res, Real **div, VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *a, *b;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return DoSomeOne(self,r);
+    SAVE(b);
+    *div = b;
+    mx =(a->MaxPrec + b->MaxPrec) *VpBaseFig();
+    GUARD_OBJ((*c),VpCreateRbObject(mx, "0"));
+    GUARD_OBJ((*res),VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
+    VpDivd(*c, *res, a, b);
+    return (VALUE)0;
+}
+
+static VALUE
+BigDecimal_div(VALUE self, VALUE r)
+{
+    ENTER(5);
+    Real *c=NULL, *res=NULL, *div = NULL;
+    r = BigDecimal_divide(&c, &res, &div, self, r);
+    SAVE(c);SAVE(res);SAVE(div);
+    if(r!=(VALUE)0) return r; /* coerced by other */
+    if(res->frac[0]*2>=div->frac[0]) {
+        /* Round up */
+        VpRdup(c);
+    }
+    return ToValue(c);
+}
+
+/*
+ * %: mod = a%b = a - (a.to_f/b).floor * b
+ * div = (a.to_f/b).floor
+ */
+static VALUE
+BigDecimal_DoDivmod(VALUE self, VALUE r, Real **div, Real **mod)
+{
+    ENTER(8);
+    Real *c=NULL, *d=NULL, *res=NULL;
+    Real *a, *b;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return DoSomeOne(self,r);
+    SAVE(b);
+
+    mx = a->Prec;
+    if(mx<b->Prec) mx = b->Prec;
+    mx =(mx + 1) * VpBaseFig();
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    GUARD_OBJ(res,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
+    VpDivd(c, res, a, b);
+    mx = c->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(d,VpCreateRbObject(mx, "0"));
+    VpRound(d,c,1,3,0);
+    VpMult(res,d,b);
+    VpAddSub(c,a,res,-1);
+    *div = d;
+    *mod = c;
+    return (VALUE)0;
+}
+
+static VALUE
+BigDecimal_mod(VALUE self, VALUE r) /* %: a%b = a - (a.to_f/b).floor * b */
+{
+    ENTER(3);
+    VALUE obj;
+    Real *div=NULL, *mod=NULL;
+
+    obj = BigDecimal_DoDivmod(self,r,&div,&mod);
+    SAVE(div);SAVE(mod);
+    if(obj!=(VALUE)0) return obj;
+    return ToValue(mod);
+}
+
+static VALUE
+BigDecimal_divremain(VALUE self, VALUE r, Real **dv, Real **rv)
+{
+    ENTER(10);
+    U_LONG mx;
+    Real *a=NULL, *b=NULL, *c=NULL, *res=NULL, *d=NULL, *rr=NULL, *ff=NULL;
+    Real *f=NULL;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    b = GetVpValue(r,0);
+    if(!b) return DoSomeOne(self,r);
+    SAVE(b);
+
+    mx  =(a->MaxPrec + b->MaxPrec) *VpBaseFig();
+    GUARD_OBJ(c  ,VpCreateRbObject(mx, "0"));
+    GUARD_OBJ(res,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
+    GUARD_OBJ(rr ,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
+    GUARD_OBJ(ff ,VpCreateRbObject((mx+1) * 2 +(VpBaseFig() + 1), "#0"));
+
+    VpDivd(c, res, a, b);
+
+    mx = c->Prec *(VpBaseFig() + 1);
+
+    GUARD_OBJ(d,VpCreateRbObject(mx, "0"));
+    GUARD_OBJ(f,VpCreateRbObject(mx, "0"));
+
+    VpRound(d,c,1,1,0); /* 1: round off */
+
+    VpFrac(f, c);
+    VpMult(rr,f,b);
+    VpAddSub(ff,res,rr,1);
+
+    *dv = d;
+    *rv = ff;
+    return (VALUE)0;
+}
+
+static VALUE
+BigDecimal_remainder(VALUE self, VALUE r) /* remainder */
+{
+    VALUE  f;
+    Real  *d,*rv;
+    f = BigDecimal_divremain(self,r,&d,&rv);
+    if(f!=(VALUE)0) return f;
+    return ToValue(rv);
+}
+
+static VALUE
+BigDecimal_divmod(VALUE self, VALUE r)
+{
+    ENTER(5);
+    VALUE obj;
+    Real *div=NULL, *mod=NULL;
+
+    obj = BigDecimal_DoDivmod(self,r,&div,&mod);
+    if(obj!=(VALUE)0) return obj;
+    SAVE(div);SAVE(mod);
+    obj = rb_ary_new();
+    rb_ary_push(obj, ToValue(div));
+    rb_ary_push(obj, ToValue(mod));
+    return obj;
+}
+
+static VALUE
+BigDecimal_divmod2(VALUE self, VALUE b, VALUE n)
+{
+    ENTER(10);
+    VALUE obj;
+    Real *res=NULL;
+    Real *av=NULL, *bv=NULL, *cv=NULL;
+    U_LONG mx = (U_LONG)GetPositiveInt(n)+VpBaseFig();
+
+    obj = rb_ary_new();
+    GUARD_OBJ(cv,VpCreateRbObject(mx,"0"));
+    GUARD_OBJ(av,GetVpValue(self,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    mx = cv->MaxPrec+1;
+    GUARD_OBJ(res,VpCreateRbObject((mx * 2 + 1)*VpBaseFig(), "#0"));
+    VpDivd(cv,res,av,bv);
+    obj = rb_ary_push(obj, ToValue(cv));
+    obj = rb_ary_push(obj, ToValue(res));
+
+    return obj;
+}
+
+static VALUE
+BigDecimal_assign2(VALUE self, VALUE n, VALUE f)
+{
+    ENTER(5);
+    Real *cv;
+    Real *av;
+    U_LONG mx = (U_LONG)GetPositiveInt(n);
+    Check_Type(f, T_FIXNUM);
+    GUARD_OBJ(cv,VpCreateRbObject(mx,"0"));
+    GUARD_OBJ(av,GetVpValue(self,1));
+    VpAsgn(cv,av,NUM2INT(f));
+    return ToValue(cv);
+}
+
+static VALUE
+BigDecimal_add2(VALUE self, VALUE b, VALUE n)
+{
+    ENTER(5);
+    Real *av;
+    Real *bv;
+    Real *cv;
+    U_LONG mx = (U_LONG)GetPositiveInt(n);
+    GUARD_OBJ(cv,VpCreateRbObject(mx,"0"));
+    GUARD_OBJ(av,GetVpValue(self,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    VpAddSub(cv,av,bv,1);
+    return ToValue(cv);
+}
+
+static VALUE
+BigDecimal_sub2(VALUE self, VALUE b, VALUE n)
+{
+    ENTER(5);
+    Real *av;
+    Real *bv;
+    Real *cv;
+    U_LONG mx = (U_LONG)GetPositiveInt(n);
+    GUARD_OBJ(cv,VpCreateRbObject(mx,"0"));
+    GUARD_OBJ(av,GetVpValue(self,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    VpAddSub(cv,av,bv,-1);
+    return ToValue(cv);
+}
+
+static VALUE
+BigDecimal_mult2(VALUE self, VALUE b, VALUE n)
+{
+    ENTER(5);
+    Real *av;
+    Real *bv;
+    Real *cv;
+    U_LONG mx = (U_LONG)GetPositiveInt(n);
+    GUARD_OBJ(cv,VpCreateRbObject(mx,"0"));
+    GUARD_OBJ(av,GetVpValue(self,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    VpMult(cv,av,bv);
+    return ToValue(cv);
+}
+
+static VALUE
+BigDecimal_dup(VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    U_LONG mx;
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpAsgn(c, a, 1);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_abs(VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpAsgn(c, a, 1);
+    VpChangeSign(c,(S_INT)1);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_sqrt(VALUE self, VALUE nFig)
+{
+    ENTER(5);
+    Real *c, *a;
+    S_INT mx, n;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    mx *= 2;
+
+    n = GetPositiveInt(nFig) + VpBaseFig() + 1;
+    if(mx <= n) mx = n;
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpSqrt(c, a);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_fix(VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpRound(c,a,1,1,0); /* 1: round off */
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_round(int argc, VALUE *argv, VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    int iLoc;
+    int sw;
+    U_LONG mx;
+    VALUE vLoc;
+
+    if(rb_scan_args(argc,argv,"01",&vLoc)==0) {
+        iLoc = 0;
+    } else {
+        Check_Type(vLoc, T_FIXNUM);
+        iLoc = NUM2INT(vLoc);
+    }
+    sw = 2;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpRound(c,a,sw,1,iLoc);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_truncate(int argc, VALUE *argv, VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    int iLoc;
+    int sw;
+    U_LONG mx;
+    VALUE vLoc;
+
+    if(rb_scan_args(argc,argv,"01",&vLoc)==0) {
+        iLoc = 0;
+    } else {
+        Check_Type(vLoc, T_FIXNUM);
+        iLoc = NUM2INT(vLoc);
+    }
+    sw = 1; /* truncate */
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpRound(c,a,sw,1,iLoc);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_frac(VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    U_LONG mx;
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpFrac(c, a);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_floor(int argc, VALUE *argv, VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    U_LONG mx;
+    int iLoc;
+    VALUE vLoc;
+
+    if(rb_scan_args(argc,argv,"01",&vLoc)==0) {
+        iLoc = 0;
+    } else {
+        Check_Type(vLoc, T_FIXNUM);
+        iLoc = NUM2INT(vLoc);
+    }
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpRound(c,a,1,3,iLoc);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_ceil(int argc, VALUE *argv, VALUE self)
+{
+    ENTER(5);
+    Real *c, *a;
+    U_LONG mx;
+    int iLoc;
+    VALUE vLoc;
+
+    if(rb_scan_args(argc,argv,"01",&vLoc)==0) {
+        iLoc = 0;
+    } else {
+        Check_Type(vLoc, T_FIXNUM);
+        iLoc = NUM2INT(vLoc);
+    }
+
+    GUARD_OBJ(a,GetVpValue(self,1));
+    mx = a->Prec *(VpBaseFig() + 1);
+    GUARD_OBJ(c,VpCreateRbObject(mx, "0"));
+    VpRound(c,a,1,2,iLoc);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_to_s(int argc, VALUE *argv, VALUE self)
+{
+    ENTER(5);
+    Real *vp;
+    char *psz;
+    U_LONG nc;
+    S_INT mc = 0;
+    VALUE f;
+
+    GUARD_OBJ(vp,GetVpValue(self,1));
+    nc = VpNumOfChars(vp)+1;
+    if(rb_scan_args(argc,argv,"01",&f)==1) {
+        Check_Type(f, T_FIXNUM);
+        mc  = GetPositiveInt(f);
+        nc += (nc + mc - 1) / mc + 1;
+    }
+    psz = ALLOCA_N(char,(unsigned int)nc);
+    VpToString(vp, psz, mc);
+    return rb_str_new2(psz);
+}
+
+static VALUE
+BigDecimal_to_parts(VALUE self)
+{
+    ENTER(5);
+    Real *vp;
+    VALUE obj,obj1;
+    S_LONG e;
+    S_LONG s;
+    char *psz1;
+
+    GUARD_OBJ(vp,GetVpValue(self,1));
+    psz1 = ALLOCA_N(char,(unsigned int)VpNumOfChars(vp));
+    VpSzMantissa(vp,psz1);
+    s = 1;
+    if(psz1[0]=='-') {
+        int i=0;
+        s = -1;
+        while(psz1[i]=psz1[i+1]) i++ ;
+    }
+    if(psz1[0]=='N') s=0; /* NaN */
+    e = VpExponent10(vp);
+    obj1 = rb_str_new2(psz1);
+    obj  = rb_ary_new();
+    rb_ary_push(obj, INT2FIX(s));
+    rb_ary_push(obj, obj1);
+    rb_ary_push(obj, INT2FIX(10));
+    rb_ary_push(obj, INT2NUM(e));
+    return obj;
+}
+
+static VALUE
+BigDecimal_exponent(VALUE self)
+{
+    S_LONG e = VpExponent10(GetVpValue(self,1));
+    return INT2NUM(e);
+}
+
+static VALUE
+BigDecimal_inspect(VALUE self)
+{
+    ENTER(5);
+    Real *vp;
+    VALUE obj;
+    unsigned int nc;
+    char *psz1;
+    char *pszAll;
+
+    GUARD_OBJ(vp,GetVpValue(self,1));
+    nc = VpNumOfChars(vp);
+    nc +=(nc + 9) / 10;
+
+    psz1   = ALLOCA_N(char,nc);
+    pszAll = ALLOCA_N(char,nc+256);
+    VpToString(vp, psz1, 10);
+    sprintf(pszAll,"#<BigDecimal:%x,'%s',%u(%u)>",self,psz1,VpPrec(vp)*VpBaseFig(),VpMaxPrec(vp)*VpBaseFig());
+
+    obj = rb_str_new2(pszAll);
+    return obj;
+}
+
+static VALUE
+BigDecimal_power(VALUE self, VALUE p)
+{
+    ENTER(5);
+    Real *x, *y;
+    S_LONG mp, ma, n;
+
+    Check_Type(p, T_FIXNUM);
+    n = NUM2INT(p);
+    ma = n;
+    if(ma < 0)  ma = -ma;
+    if(ma == 0) ma = 1;
+
+    GUARD_OBJ(x,GetVpValue(self,1));
+    if(VpIsDef(x)) {
+        mp = x->Prec *(VpBaseFig() + 1);
+        GUARD_OBJ(y,VpCreateRbObject(mp *(ma + 1), "0"));
+    } else {
+        GUARD_OBJ(y,VpCreateRbObject(1, "0"));
+    }
+    VpPower(y, x, n);
+    return ToValue(y);
+}
+
+static VALUE
+BigDecimal_new(int argc, VALUE *argv, VALUE self)
+{
+    ENTER(5);
+    Real *pv;
+    S_LONG mf;
+    VALUE  nFig;
+    VALUE  iniValue;
+
+    if(rb_scan_args(argc,argv,"11",&iniValue,&nFig)==1) {
+        mf = 0;
+    } else {
+        mf = GetPositiveInt(nFig);
+    }
+    Check_SafeStr(iniValue);
+    GUARD_OBJ(pv,VpNewRbClass(mf, RSTRING(iniValue)->ptr,self));
+    return ToValue(pv);
+}
+
+static VALUE
+BigDecimal_limit(int argc, VALUE *argv, VALUE self)
+{
+    VALUE  nFig;
+    VALUE  nCur = INT2NUM(VpGetPrecLimit());
+
+    if(rb_scan_args(argc,argv,"01",&nFig)==1) {
+        Check_Type(nFig, T_FIXNUM);
+        VpSetPrecLimit(NUM2INT(nFig));
+    }
+    return nCur;
+}
+
+static VALUE
+BigDecimal_e(VALUE self, VALUE nFig)
+{
+    ENTER(5);
+    Real *pv;
+    S_LONG mf;
+
+    mf = GetPositiveInt(nFig);
+    GUARD_OBJ(pv,VpCreateRbObject(mf, "0"));
+    VpExp1(pv);
+    return ToValue(pv);
+}
+
+static VALUE
+BigDecimal_pai(VALUE self, VALUE nFig)
+{
+    ENTER(5);
+    Real *pv;
+    S_LONG mf;
+
+    mf = GetPositiveInt(nFig);
+    GUARD_OBJ(pv,VpCreateRbObject(mf, "0"));
+    VpPai(pv);
+    return ToValue(pv);
+}
+
+static VALUE
+BigDecimal_exp(VALUE self, VALUE nFig)
+{
+    ENTER(5);
+    Real *c, *y;
+    S_LONG mf;
+
+    GUARD_OBJ(y,GetVpValue(self,1));
+    mf = GetPositiveInt(nFig);
+    GUARD_OBJ(c,VpCreateRbObject(mf, "0"));
+    VpExp(c, y);
+    return ToValue(c);
+}
+
+static VALUE
+BigDecimal_sign(VALUE self)
+{ /* sign */
+    int s = GetVpValue(self,1)->sign;
+    return INT2FIX(s);
+}
+
+static VALUE
+BigDecimal_sincos(VALUE self, VALUE nFig)
+{
+    ENTER(5);
+    VALUE obj;
+    VALUE objSin;
+    VALUE objCos;
+    Real *pcos, *psin, *y;
+    S_LONG mf;
+
+    obj = rb_ary_new();
+    GUARD_OBJ(y,GetVpValue(self,1));
+    mf = GetPositiveInt(nFig);
+    GUARD_OBJ(pcos,VpCreateRbObject(mf, "0"));
+    GUARD_OBJ(psin,VpCreateRbObject(mf, "0"));
+    VpSinCos(psin, pcos, y);
+
+    objSin = ToValue(psin);
+    objCos = ToValue(pcos);
+    rb_ary_push(obj, objSin);
+    rb_ary_push(obj, objCos);
+    return obj;
+}
+
+
+#ifdef USE_MUTABLE_METHOD
+/**** Following methods are all MUTABLE and not currently activated. ****/
+static void
+CheckAssign(VALUE x, VALUE y)
+{
+    if(x==y)
+    rb_fatal("Bad assignment(the same object appears on both LHS and RHS).");
+}
+
+static VALUE
+BigDecimal_divmod4(VALUE self, VALUE c, VALUE r, VALUE a, VALUE b)
+{
+    ENTER(10);
+    U_LONG f;
+    Real *res=NULL;
+    Real *av=NULL, *bv=NULL, *cv=NULL;
+    CheckAssign(c,a);
+    CheckAssign(c,b);
+    CheckAssign(r,a);
+    CheckAssign(r,b);
+    CheckAssign(r,c);
+    GUARD_OBJ(cv,GetVpValue(c,1));
+    GUARD_OBJ(av,GetVpValue(a,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    GUARD_OBJ(res,GetVpValue(r,1));
+    f = VpDivd(cv,res,av,bv);
+    return INT2FIX(f);
+}
+
+static VALUE
+BigDecimal_assign(VALUE self, VALUE c, VALUE a, VALUE f)
+{
+    ENTER(5);
+    int v;
+    Real *av;
+    Real *cv;
+    CheckAssign(c,a);
+    Check_Type(f, T_FIXNUM);
+    GUARD_OBJ(cv,GetVpValue(c,1));
+    GUARD_OBJ(av,GetVpValue(a,1));
+    v = VpAsgn(cv,av,NUM2INT(f));
+    return INT2NUM(v);
+}
+
+static VALUE
+BigDecimal_add3(VALUE self, VALUE c, VALUE a, VALUE b)
+{
+    ENTER(5);
+    Real *av;
+    Real *bv;
+    Real *cv;
+    U_LONG f;
+    CheckAssign(c,a);
+    CheckAssign(c,b);
+    GUARD_OBJ(cv,GetVpValue(c,1));
+    GUARD_OBJ(av,GetVpValue(a,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    f = VpAddSub(cv,av,bv,1);
+    return INT2NUM(f);
+}
+
+static VALUE
+BigDecimal_sub3(VALUE self, VALUE c, VALUE a, VALUE b)
+{
+    ENTER(5);
+    Real *av;
+    Real *bv;
+    Real *cv;
+    U_LONG f;
+    CheckAssign(c,a);
+    CheckAssign(c,b);
+    GUARD_OBJ(cv,GetVpValue(c,1));
+    GUARD_OBJ(av,GetVpValue(a,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    f = VpAddSub(cv,av,bv,-1);
+    return INT2NUM(f);
+}
+
+static VALUE
+BigDecimal_mult3(VALUE self, VALUE c, VALUE a, VALUE b)
+{
+    ENTER(5);
+    Real *av;
+    Real *bv;
+    Real *cv;
+    U_LONG f;
+    CheckAssign(c,a);
+    CheckAssign(c,b);
+    GUARD_OBJ(cv,GetVpValue(c,1));
+    GUARD_OBJ(av,GetVpValue(a,1));
+    GUARD_OBJ(bv,GetVpValue(b,1));
+    f = VpMult(cv,av,bv);
+    return INT2NUM(f);
+}
+#endif /* USE_MUTABLE_METHOD */
+
+void
+Init_bigdecimal(void)
+{
+    /* Initialize VP routines */
+    VpInit((U_LONG)0);
+    coerce = rb_intern("coerce");
+    /* Class and method registration */
+    rb_cBigDecimal = rb_define_class("BigDecimal",rb_cNumeric);
+
+    /* Class methods */
+    rb_define_singleton_method(rb_cBigDecimal, "mode", BigDecimal_mode, 2);
+    rb_define_singleton_method(rb_cBigDecimal, "new", BigDecimal_new, -1);
+    rb_define_singleton_method(rb_cBigDecimal, "limit", BigDecimal_limit, -1);
+    rb_define_singleton_method(rb_cBigDecimal, "E", BigDecimal_e, 1);
+    rb_define_singleton_method(rb_cBigDecimal, "double_fig", BigDecimal_double_fig, 0);
+    rb_define_singleton_method(rb_cBigDecimal, "PI", BigDecimal_pai, 1);
+    rb_define_singleton_method(rb_cBigDecimal, "induced_from",BigDecimal_induced_from, 1);
+    rb_define_singleton_method(rb_cBigDecimal, "_load", BigDecimal_load, 1);
+
+    /* Constants */
+    rb_define_const(rb_cBigDecimal, "BASE", INT2FIX((S_INT)VpBaseVal()));
+    rb_define_const(rb_cBigDecimal, "EXCEPTION_ALL",INT2FIX(VP_EXCEPTION_ALL));
+    rb_define_const(rb_cBigDecimal, "EXCEPTION_NaN",INT2FIX(VP_EXCEPTION_NaN));
+    rb_define_const(rb_cBigDecimal, "EXCEPTION_INFINITY",INT2FIX(VP_EXCEPTION_INFINITY));
+    rb_define_const(rb_cBigDecimal, "EXCEPTION_UNDERFLOW",INT2FIX(VP_EXCEPTION_UNDERFLOW));
+    rb_define_const(rb_cBigDecimal, "EXCEPTION_OVERFLOW",INT2FIX(VP_EXCEPTION_OVERFLOW));
+    rb_define_const(rb_cBigDecimal, "EXCEPTION_ZERODIVIDE",INT2FIX(VP_EXCEPTION_ZERODIVIDE));
+
+    /* Constants for sign value */
+    rb_define_const(rb_cBigDecimal, "SIGN_NaN",INT2FIX(VP_SIGN_NaN));
+    rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_ZERO",INT2FIX(VP_SIGN_POSITIVE_ZERO));
+    rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_ZERO",INT2FIX(VP_SIGN_NEGATIVE_ZERO));
+    rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_FINITE",INT2FIX(VP_SIGN_POSITIVE_FINITE));
+    rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_FINITE",INT2FIX(VP_SIGN_NEGATIVE_FINITE));
+    rb_define_const(rb_cBigDecimal, "SIGN_POSITIVE_INFINITE",INT2FIX(VP_SIGN_POSITIVE_INFINITE));
+    rb_define_const(rb_cBigDecimal, "SIGN_NEGATIVE_INFINITE",INT2FIX(VP_SIGN_NEGATIVE_INFINITE));
+
+    /* instance methods */
+    rb_define_method(rb_cBigDecimal, "prec", BigDecimal_prec, 0);
+    rb_define_method(rb_cBigDecimal, "assign", BigDecimal_assign2, 2);
+    rb_define_method(rb_cBigDecimal, "add", BigDecimal_add2, 2);
+    rb_define_method(rb_cBigDecimal, "sub", BigDecimal_sub2, 2);
+    rb_define_method(rb_cBigDecimal, "mult", BigDecimal_mult2, 2);
+    rb_define_method(rb_cBigDecimal, "div",BigDecimal_divmod2, 2);
+    rb_define_method(rb_cBigDecimal, "hash", BigDecimal_hash, 0);
+    rb_define_method(rb_cBigDecimal, "to_s", BigDecimal_to_s, -1);
+    rb_define_method(rb_cBigDecimal, "to_i", BigDecimal_to_i, 0);
+    rb_define_method(rb_cBigDecimal, "to_parts", BigDecimal_to_parts, 0);
+    rb_define_method(rb_cBigDecimal, "+", BigDecimal_add, 1);
+    rb_define_method(rb_cBigDecimal, "-", BigDecimal_sub, 1);
+    rb_define_method(rb_cBigDecimal, "+@", BigDecimal_uplus, 0);
+    rb_define_method(rb_cBigDecimal, "-@", BigDecimal_neg, 0);
+    rb_define_method(rb_cBigDecimal, "*", BigDecimal_mult, 1);
+    rb_define_method(rb_cBigDecimal, "/", BigDecimal_div, 1);
+    rb_define_method(rb_cBigDecimal, "%", BigDecimal_mod, 1);
+    rb_define_method(rb_cBigDecimal, "modulo", BigDecimal_mod, 1);
+    rb_define_method(rb_cBigDecimal, "remainder", BigDecimal_remainder, 1);
+    rb_define_method(rb_cBigDecimal, "divmod", BigDecimal_divmod, 1);
+    rb_define_method(rb_cBigDecimal, "dup", BigDecimal_dup, 0);
+    rb_define_method(rb_cBigDecimal, "to_f", BigDecimal_dup, 0); /* to_f === dup */
+    rb_define_method(rb_cBigDecimal, "abs", BigDecimal_abs, 0);
+    rb_define_method(rb_cBigDecimal, "sqrt", BigDecimal_sqrt, 1);
+    rb_define_method(rb_cBigDecimal, "fix", BigDecimal_fix, 0);
+    rb_define_method(rb_cBigDecimal, "round", BigDecimal_round, -1);
+    rb_define_method(rb_cBigDecimal, "frac", BigDecimal_frac, 0);
+    rb_define_method(rb_cBigDecimal, "floor", BigDecimal_floor, -1);
+    rb_define_method(rb_cBigDecimal, "ceil", BigDecimal_ceil, -1);
+    rb_define_method(rb_cBigDecimal, "power", BigDecimal_power, 1);
+    rb_define_method(rb_cBigDecimal, "exp", BigDecimal_exp, 1);
+    rb_define_method(rb_cBigDecimal, "sincos", BigDecimal_sincos, 1);
+    rb_define_method(rb_cBigDecimal, "<=>", BigDecimal_comp, 1);
+    rb_define_method(rb_cBigDecimal, "==", BigDecimal_eq, 1);
+    rb_define_method(rb_cBigDecimal, "===", BigDecimal_eq, 1);
+    rb_define_method(rb_cBigDecimal, "eql?", BigDecimal_eq, 1);
+    rb_define_method(rb_cBigDecimal, "!=", BigDecimal_ne, 1);
+    rb_define_method(rb_cBigDecimal, "<", BigDecimal_lt, 1);
+    rb_define_method(rb_cBigDecimal, "<=", BigDecimal_le, 1);
+    rb_define_method(rb_cBigDecimal, ">", BigDecimal_gt, 1);
+    rb_define_method(rb_cBigDecimal, ">=", BigDecimal_ge, 1);
+    rb_define_method(rb_cBigDecimal, "zero?", BigDecimal_zero, 0);
+    rb_define_method(rb_cBigDecimal, "nonzero?", BigDecimal_nonzero, 0);
+    rb_define_method(rb_cBigDecimal, "coerce", BigDecimal_coerce, 1);
+    rb_define_method(rb_cBigDecimal, "inspect", BigDecimal_inspect, 0);
+    rb_define_method(rb_cBigDecimal, "exponent", BigDecimal_exponent, 0);
+    rb_define_method(rb_cBigDecimal, "sign", BigDecimal_sign, 0);
+    rb_define_method(rb_cBigDecimal, "nan?",      BigDecimal_IsNaN, 0);
+    rb_define_method(rb_cBigDecimal, "infinite?", BigDecimal_IsInfinite, 0);
+    rb_define_method(rb_cBigDecimal, "finite?",   BigDecimal_IsFinite, 0);
+    rb_define_method(rb_cBigDecimal, "truncate",  BigDecimal_truncate, -1);
+    rb_define_method(rb_cBigDecimal, "_dump", BigDecimal_dump, -1);
+
+#ifdef USE_MUTABLE_METHOD
+    rb_define_singleton_method(rb_cBigDecimal, "assign!", BigDecimal_assign, 3);
+    rb_define_singleton_method(rb_cBigDecimal, "add!", BigDecimal_add3, 3);
+    rb_define_singleton_method(rb_cBigDecimal, "sub!", BigDecimal_sub3, 3);
+    rb_define_singleton_method(rb_cBigDecimal, "mult!", BigDecimal_mult3, 3);
+    rb_define_singleton_method(rb_cBigDecimal, "div!",BigDecimal_divmod4, 4);
+#endif /* USE_MUTABLE_METHOD */
+}
+
+/*
+ *
+ *  ============================================================================
+ *
+ *  vp_ routines begins here
+ *
+ *  ============================================================================
+ *
+ */
+#ifdef _DEBUG
+static int gfDebug = 0;         /* Debug switch */
+static int gfCheckVal = 1;      /* Value checking flag in VpNmlz()  */
+#endif /* _DEBUG */
+
+static U_LONG gnPrecLimit = 0;  /* Global upper limit of the precision newly allocated */
+static U_LONG BASE_FIG = 4;     /* =log10(BASE)  */
+static U_LONG BASE = 10000L;    /* Base value(value must be 10**BASE_FIG) */
+                /* The value of BASE**2 + BASE must be represented */
+                /* within one U_LONG. */
+static U_LONG HALF_BASE = 5000L;/* =BASE/2  */
+static S_LONG DBLE_FIG = 8;    /* figure of double */
+static U_LONG BASE1 = 1000L;    /* =BASE/10  */
+
+static Real *VpConstOne;    /* constant 1.0 */
+static Real *VpPt5;        /* constant 0.5 */
+static U_LONG maxnr = 100;    /* Maximum iterations for calcurating sqrt. */
+                /* used in VpSqrt() */
+
+#ifdef _DEBUG
+static int gnAlloc=0; /* Memory allocation counter */
+#endif /* _DEBUG */
+
+VP_EXPORT void *
+VpMemAlloc(U_LONG mb)
+{
+    void *p = xmalloc((unsigned int)mb);
+    if(!p) {
+        VpException(VP_EXCEPTION_MEMORY,"failed to allocate memory",1);
+    }
+    memset(p,0,mb);
+#ifdef _DEBUG
+    gnAlloc++; /* Count allocation call */
+#endif /* _DEBUG */
+    return p;
+}
+
+VP_EXPORT void
+VpFree(Real *pv)
+{
+    if(pv != NULL) {
+        xfree(pv);
+#ifdef _DEBUG
+        gnAlloc--; /* Decrement allocation count */
+        if(gnAlloc==0) {
+            printf(" *************** All memories allocated freed ****************");
+            getchar();
+        }
+        if(gnAlloc<0) {
+            printf(" ??????????? Too many memory free calls(%d) ?????????????\n",gnAlloc);
+            getchar();
+        }
+#endif /* _DEBUG */
+    }
+}
+
+/*
+ * EXCEPTION Handling.
+ */
+static unsigned short gfDoException = 0; /* Exception flag */
+
+static unsigned short
+VpGetException (void)
+{
+    return gfDoException;
+}
+
+static void
+VpSetException(unsigned short f)
+{
+    gfDoException = f;
+}
+
+/* These 2 functions added at v1.1.7 */
+VP_EXPORT U_LONG
+VpGetPrecLimit(void)
+{
+    return gnPrecLimit;
+}
+
+VP_EXPORT U_LONG
+VpSetPrecLimit(U_LONG n)
+{
+    U_LONG s = gnPrecLimit;
+    gnPrecLimit = n;
+    return s;
+}
+
+/*
+ *  0.0 & 1.0 generator
+ *    These gZero_..... and gOne_..... can be any name
+ *    referenced from nowhere except Zero() and One().
+ *    gZero_..... and gOne_..... must have global scope
+ *    (to let the compiler know they may be changed in outside
+ *    (... but not actually..)).
+ */
+double gZero_ABCED9B1_CE73__00400511F31D = 0.0;
+double gOne_ABCED9B4_CE73__00400511F31D  = 1.0;
+static double
+Zero(void)
+{
+    return gZero_ABCED9B1_CE73__00400511F31D;
+}
+
+static double
+One(void)
+{
+    return gOne_ABCED9B4_CE73__00400511F31D;
+}
+
+VP_EXPORT U_LONG
+VpBaseFig(void)
+{
+    return BASE_FIG;
+}
+
+VP_EXPORT U_LONG
+VpDblFig(void)
+{
+    return DBLE_FIG;
+}
+
+VP_EXPORT U_LONG
+VpBaseVal(void)
+{
+    return BASE;
+}
+
+/*
+  ----------------------------------------------------------------
+  Value of sign in Real structure is reserved for future use.
+  short sign;
+                    ==0 : NaN
+                      1 : Positive zero
+                     -1 : Negative zero
+                      2 : Positive number
+                     -2 : Negative number
+                      3 : Positive infinite number
+                     -3 : Negative infinite number
+  ----------------------------------------------------------------
+*/
+
+VP_EXPORT double
+VpGetDoubleNaN(void) /* Returns the value of NaN */
+{
+    static double fNaN = 0.0;
+    if(fNaN==0.0) fNaN = Zero()/Zero();
+    return fNaN;
+}
+
+VP_EXPORT double
+VpGetDoublePosInf(void) /* Returns the value of +Infinity */
+{
+    static double fInf = 0.0;
+    if(fInf==0.0) fInf = One()/Zero();
+    return fInf;
+}
+
+VP_EXPORT double
+VpGetDoubleNegInf(void) /* Returns the value of -Infinity */
+{
+    static double fInf = 0.0;
+    if(fInf==0.0) fInf = -(One()/Zero());
+    return fInf;
+}
+
+VP_EXPORT double
+VpGetDoubleNegZero(void) /* Returns the value of -0 */
+{
+    static double nzero = 1000.0;
+    if(nzero!=0.0) nzero = (One()/VpGetDoubleNegInf());
+    return nzero;
+}
+
+VP_EXPORT int
+VpIsNegDoubleZero(double v)
+{
+    double z = VpGetDoubleNegZero();
+    return MemCmp(&v,&z,sizeof(v))==0;
+}
+
+VP_EXPORT int
+VpException(unsigned short f,char *str,int always)
+{
+    VALUE exc;
+    int   fatal=0;
+
+    if(f==VP_EXCEPTION_OP || f==VP_EXCEPTION_MEMORY) always = 1;
+
+    if(always||(gfDoException&f)) {
+        switch(f)
+        {
+        /*
+        case VP_EXCEPTION_ZERODIVIDE:
+        case VP_EXCEPTION_OVERFLOW:
+        */
+        case VP_EXCEPTION_INFINITY:
+             exc = rb_eFloatDomainError;
+             goto raise;
+        case VP_EXCEPTION_NaN:
+             exc = rb_eFloatDomainError;
+             goto raise;
+        case VP_EXCEPTION_UNDERFLOW:
+             exc = rb_eFloatDomainError;
+             goto raise;
+        case VP_EXCEPTION_OP:
+             exc = rb_eFloatDomainError;
+             goto raise;
+        case VP_EXCEPTION_MEMORY:
+             fatal = 1;
+             goto raise;
+        default:
+             fatal = 1;
+             goto raise;
+        }
+    }
+    return 0; /* 0 Means VpException() raised no exception */
+
+raise:
+    if(fatal) rb_fatal(str);
+    else   rb_raise(exc,str);
+    return 0;
+}
+
+/* Throw exception or returns 0,when resulting c is Inf or NaN */
+/*  sw=1:+ 2:- 3:* 4:/ */
+static int
+VpIsDefOP(Real *c,Real *a,Real *b,int sw)
+{
+ if(VpIsNaN(a) || VpIsNaN(b)) {
+        /* at least a or b is NaN */
+        VpSetNaN(c);
+        goto NaN;
+ }
+
+ if(VpIsInf(a)) {
+        if(VpIsInf(b)) {
+            switch(sw)
+            {
+            case 1: /* + */
+                if(VpGetSign(a)==VpGetSign(b)) {
+                    VpSetInf(c,VpGetSign(a));
+                    goto Inf;
+                } else {
+                    VpSetNaN(c);
+                    goto NaN;
+                }
+            case 2: /* - */
+                if(VpGetSign(a)!=VpGetSign(b)) {
+                    VpSetInf(c,VpGetSign(a));
+                    goto Inf;
+                } else {
+                    VpSetNaN(c);
+                    goto NaN;
+                }
+                break;
+            case 3: /* * */
+                VpSetInf(c,VpGetSign(a)*VpGetSign(b));
+                goto Inf;
+                break;
+            case 4: /* / */
+                VpSetNaN(c);
+                goto NaN;
+            }
+            VpSetNaN(c);
+            goto NaN;
+        }
+        /* Inf op Finite */
+        switch(sw)
+        {
+        case 1: /* + */
+        case 2: /* - */
+                VpSetInf(c,VpGetSign(a));
+                break;
+        case 3: /* * */
+                if(VpIsZero(b)) {
+                    VpSetNaN(c);
+                    goto NaN;
+                }
+                VpSetInf(c,VpGetSign(a)*VpGetSign(b));
+                break;
+        case 4: /* / */
+                VpSetInf(c,VpGetSign(a)*VpGetSign(b));
+        }
+        goto Inf;
+ }
+
+    if(VpIsInf(b)) {
+        switch(sw)
+        {
+        case 1: /* + */
+                VpSetInf(c,VpGetSign(b));
+                break;
+        case 2: /* - */
+                VpSetInf(c,-VpGetSign(b));
+                break;
+        case 3: /* * */
+                if(VpIsZero(a)) {
+                    VpSetNaN(c);
+                    goto NaN;
+                }
+                VpSetInf(c,VpGetSign(a)*VpGetSign(b));
+                break;
+        case 4: /* / */
+                VpSetZero(c,VpGetSign(a)*VpGetSign(b));
+        }
+        goto Inf;
+    }
+    return 1; /* Results OK */
+
+Inf:
+    return VpException(VP_EXCEPTION_INFINITY,"Computation results to 'Infinity'",0);
+NaN:
+    return VpException(VP_EXCEPTION_NaN,"Computation results to 'NaN'",0);
+}
+
+/*
+  ----------------------------------------------------------------
+*/
+
+/*
+ *    returns number of chars needed to represent vp.
+ */
+VP_EXPORT U_LONG
+VpNumOfChars(Real *vp)
+{
+    if(vp == NULL)   return BASE_FIG*2+6;
+    if(!VpIsDef(vp)) return 32; /* not sure,may be OK */
+    return     BASE_FIG *(vp->Prec + 2)+6; /* 3: sign + exponent chars */
+}
+
+/*
+ * Initializer for Vp routines and constants used.
+ * [Input]
+ *   BaseVal: Base value(assigned to BASE) for Vp calculation.
+ *   It must be the form BaseVal=10**n.(n=1,2,3,...)
+ *   If Base <= 0L,then the BASE will be calcurated so
+ *   that BASE is as large as possible satisfying the
+ *   relation MaxVal <= BASE*(BASE+1). Where the value
+ *   MaxVal is the largest value which can be represented
+ *   by one U_LONG word(LONG) in the computer used.
+ *
+ * [Returns]
+ * DBLE_FIG   ... OK
+ */
+VP_EXPORT U_LONG
+VpInit(U_LONG BaseVal)
+{
+    U_LONG w;
+    double v;
+
+    /* Setup +/- Inf  NaN -0 */
+    VpGetDoubleNaN();
+    VpGetDoublePosInf();
+    VpGetDoubleNegInf();
+    VpGetDoubleNegZero();
+
+    if(BaseVal <= 0) {
+        /* Base <= 0, then determine Base by calcuration. */
+        BASE = 1;
+        while(
+               (BASE > 0) &&
+               ((w = BASE *(BASE + 1)) > BASE) &&((w / BASE) ==(BASE + 1))
+            ) {
+            BaseVal = BASE;
+            BASE = BaseVal * 10L;
+        }
+    }
+    /* Set Base Values */
+    BASE = BaseVal;
+    HALF_BASE = BASE / 2;
+    BASE1 = BASE / 10;
+    BASE_FIG = 0;
+    while(BaseVal /= 10) ++BASE_FIG;
+    /* Allocates Vp constants. */
+    VpConstOne = VpAlloc((U_LONG)1, "1");
+    VpPt5 = VpAlloc((U_LONG)1, ".5");
+
+#ifdef _DEBUG
+    gnAlloc = 0;
+#endif /* _DEBUG */
+
+    /* Determine # of digits available in one 'double'. */
+
+    v = 1.0;
+    DBLE_FIG = 0;
+    while(v + 1.0 > 1.0) {
+        ++DBLE_FIG;
+        v /= 10;
+    }
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        printf("VpInit: BaseVal   = %u\n", BaseVal);
+        printf("  BASE   = %u\n", BASE);
+        printf("  HALF_BASE = %u\n", HALF_BASE);
+        printf("  BASE1  = %u\n", BASE1);
+        printf("  BASE_FIG  = %u\n", BASE_FIG);
+        printf("  DBLE_FIG  = %u\n", DBLE_FIG);
+    }
+#endif /* _DEBUG */
+
+    return DBLE_FIG;
+}
+
+/* If exponent overflows,then raise exception or returns 0 */
+static int
+AddExponent(Real *a,S_INT n)
+{
+    S_INT e = a->exponent;
+    S_INT m = e+n;
+    S_INT eb,mb;
+    if(e>0) {
+        if(n>0) {
+            mb = m*BASE_FIG;
+            eb = e*BASE_FIG;
+            if(mb<eb) goto overflow;
+        }
+    } else if(n<0) {
+        mb = m*BASE_FIG;
+        eb = e*BASE_FIG;
+        if(mb>eb) goto underflow;
+    }
+    a->exponent = m;
+    return 1;
+
+/* Overflow/Underflow ==> Raise exception or returns 0 */
+underflow:
+    VpSetZero(a,VpGetSign(a));
+    return VpException(VP_EXCEPTION_UNDERFLOW,"Exponent underflow",0);
+
+overflow:
+    VpSetInf(a,VpGetSign(a));
+    return VpException(VP_EXCEPTION_OVERFLOW,"Exponent overflow",0);
+}
+
+/*
+ * Allocates variable.
+ * [Input]
+ *   mx ... allocation unit, if zero then mx is determined by szVal.
+ *    The mx is the number of effective digits can to be stored.
+ *   szVal ... value assigned(char). If szVal==NULL,then zero is assumed.
+ *            If szVal[0]=='#' then Max. Prec. will not be considered(1.1.7),
+ *            full precision specified by szVal is allocated.
+ *
+ * [Returns]
+ *   Pointer to the newly allocated variable, or
+ *   NULL be returned if memory allocation is failed,or any error.
+ */
+VP_EXPORT Real *
+VpAlloc(U_LONG mx, char *szVal)
+{
+    U_LONG i, ni, ipf, nf, ipe, ne, nalloc;
+    char v;
+    int  sign=1;
+    Real *vp = NULL;
+    U_LONG mf = VpGetPrecLimit();
+    mx = (mx + BASE_FIG - 1) / BASE_FIG + 1;    /* Determine allocation unit. */
+    if(szVal) {
+        if(*szVal!='#') {
+             if(mf) {
+                mf = (mf + BASE_FIG - 1) / BASE_FIG + 1;
+                 if(mx>mf) {
+                    mx = mf;
+                }
+            }
+        } else {
+            ++szVal;
+        }
+    }
+
+    /* necessary to be able to store */
+    /* at least mx digits. */
+    if(szVal == NULL) {
+        /* szVal==NULL ==> allocate zero value. */
+        vp = (Real *) VpMemAlloc(sizeof(Real) + mx * sizeof(U_LONG));
+        /* xmalloc() alway returns(or throw interruption) */
+        vp->MaxPrec = mx;    /* set max precision */
+        VpSetZero(vp,1);    /* initialize vp to zero. */
+        return vp;
+    }
+    /* Check on Inf & NaN */
+    if(StrCmp(szVal,"+Infinity")==0 ||
+       StrCmp(szVal, "Infinity")==0 ) {
+        vp = (Real *) VpMemAlloc(sizeof(Real) + sizeof(U_LONG));
+        vp->MaxPrec = 1;    /* set max precision */
+        VpSetPosInf(vp);
+        return vp;
+    }
+    if(StrCmp(szVal,"-Infinity")==0) {
+        vp = (Real *) VpMemAlloc(sizeof(Real) + sizeof(U_LONG));
+        vp->MaxPrec = 1;    /* set max precision */
+        VpSetNegInf(vp);
+        return vp;
+    }
+    if(StrCmp(szVal,"NaN")==0) {
+        vp = (Real *) VpMemAlloc(sizeof(Real) + sizeof(U_LONG));
+        vp->MaxPrec = 1;    /* set max precision */
+        VpSetNaN(vp);
+        return vp;
+    }
+
+    /* check on number szVal[] */
+    i = SkipWhiteChar(szVal);
+    if  (szVal[i] == '-') {sign=-1;++i;}
+    else if(szVal[i] == '+')  ++i;
+    /* Skip digits */
+    ni = 0;            /* digits in mantissa */
+    while(v = szVal[i]) {
+        if((v > '9') ||(v < '0')) break;
+        ++i;
+        ++ni;
+    }
+    nf = 0;
+    ipf = 0;
+    ipe = 0;
+    ne = 0;
+    if(v) {
+        /* other than digit nor \0 */
+        if(szVal[i] == '.') {    /* xxx. */
+            ++i;
+            ipf = i;
+            while(v = szVal[i]) {    /* get fraction part. */
+                if((v > '9') ||(v < '0')) break;
+                ++i;
+                ++nf;
+            }
+        }
+        ipe = 0;        /* Exponent */
+
+        switch(szVal[i]) {
+        case '\0': break;
+        case 'e':
+        case 'E':
+        case 'd':
+        case 'D':
+            ++i;
+            ipe = i;
+            v = szVal[i];
+            if((v == '-') ||(v == '+')) ++i;
+            while(szVal[i]) {
+                ++i;
+                ++ne;
+            }
+            break;
+        default:
+            break;
+        }
+    }
+    nalloc =(ni + nf + BASE_FIG - 1) / BASE_FIG + 1;    /* set effective allocation  */
+    /* units for szVal[]  */
+    if(mx <= 0) mx = 1;
+    nalloc = Max(nalloc, mx);
+    mx = nalloc;
+    vp =(Real *) VpMemAlloc(sizeof(Real) + mx * sizeof(U_LONG));
+    /* xmalloc() alway returns(or throw interruption) */
+    vp->MaxPrec = mx;        /* set max precision */
+    VpSetZero(vp,sign);
+    VpCtoV(vp, szVal, ni, &(szVal[ipf]), nf, &(szVal[ipe]), ne);
+    return vp;
+}
+
+/*
+ * Assignment(c=a).
+ * [Input]
+ *   a   ... RHSV
+ *   isw ... switch for assignment.
+ *    c = a  when isw =  1 or 2
+ *    c = -a when isw = -1 or -1
+ *    when |isw|==1
+ *    if c->MaxPrec < a->Prec,then round up
+ *    will not be performed.
+ * [Output]
+ *  c  ... LHSV
+ */
+VP_EXPORT int
+VpAsgn(Real *c, Real *a, int isw)
+{
+    U_LONG j, n;
+    if(VpIsNaN(a)) {
+        VpSetNaN(c);
+        return 0;
+    }
+    if(VpIsInf(a)) {
+        VpSetInf(c,isw);
+        return 0;
+    }
+
+    /* check if the RHS is zero */
+    if(!VpIsZero(a)) {
+        c->exponent = a->exponent;    /* store  exponent */
+        VpSetSign(c,(isw*VpGetSign(a)));    /* set sign */
+        n =(a->Prec < c->MaxPrec) ?(a->Prec) :(c->MaxPrec);
+        c->Prec = n;
+        for(j=0;j < n; ++j) c->frac[j] = a->frac[j];
+        if(isw < 0) isw = -isw;
+        if(isw == 2) {
+            if(a->MaxPrec>n) {
+                if((c->Prec < a->Prec) &&
+                   (a->frac[n] >= HALF_BASE)) VpRdup(c);    /* round up/off */
+            }
+        }
+    } else {
+        /* The value of 'a' is zero.  */
+        VpSetZero(c,isw*VpGetSign(a));
+        return 1;
+    }
+    VpNmlz(c);
+    return c->Prec*BASE_FIG;
+}
+
+/*
+ *   c = a + b  when operation =  1 or 2
+ *  = a - b  when operation = -1 or -2.
+ *   Returns number of significant digits of c
+ */
+VP_EXPORT int
+VpAddSub(Real *c, Real *a, Real *b, int operation)
+{
+    S_INT sw, isw;
+    Real *a_ptr, *b_ptr;
+    U_LONG n, na, nb, i;
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpAddSub(enter) a=% \n", a);
+        VPrint(stdout, "     b=% \n", b);
+        printf(" operation=%d\n", operation);
+    }
+#endif /* _DEBUG */
+
+    if(!VpIsDefOP(c,a,b,(operation>0)?1:2)) return 0; /* No significant digits */
+
+    /* check if a or b is zero  */
+    if(VpIsZero(a)) {
+        /* a is zero,then assign b to c */
+        if(!VpIsZero(b)) {
+            VpAsgn(c, b, operation);
+        } else {
+            /* Both a and b are zero. */
+            if(VpGetSign(a)<0 && operation*VpGetSign(b)<0) {
+                /* -0 -0 */
+                VpSetZero(c,-1);
+            } else {
+                VpSetZero(c,1);
+            }
+            return 1; /* 0: 1 significant digits */
+        }
+        return c->Prec*BASE_FIG;
+    }
+    if(VpIsZero(b)) {
+        /* b is zero,then assign a to c. */
+        VpAsgn(c, a, 1);
+        return c->Prec*BASE_FIG;
+    }
+
+    if(operation < 0) sw = -1;
+    else     sw =  1;
+
+    /* compare absolute value. As a result,|a_ptr|>=|b_ptr| */
+    if(a->exponent > b->exponent) {
+        a_ptr = a;
+        b_ptr = b;
+    }         /* |a|>|b| */
+    else if(a->exponent < b->exponent) {
+        a_ptr = b;
+        b_ptr = a;
+    }                /* |a|<|b| */
+    else {
+        /* Exponent part of a and b is the same,then compare fraction */
+        /* part */
+        na = a->Prec;
+        nb = b->Prec;
+        n = Min(na, nb);
+        for(i=0;i < n; ++i) {
+            if(a->frac[i] > b->frac[i]) {
+                a_ptr = a;
+                b_ptr = b;
+                goto end_if;
+            } else if(a->frac[i] < b->frac[i]) {
+                a_ptr = b;
+                b_ptr = a;
+                goto end_if;
+            }
+        }
+        if(na > nb) {
+         a_ptr = a;
+            b_ptr = b;
+            goto end_if;
+        } else if(na < nb) {
+            a_ptr = b;
+            b_ptr = a;
+            goto end_if;
+        }
+        /* |a| == |b| */
+        if(VpGetSign(a) + sw *VpGetSign(b) == 0) {
+            VpSetZero(c,1);        /* abs(a)=abs(b) and operation = '-'  */
+            return c->Prec*BASE_FIG;
+        }
+        a_ptr = a;
+        b_ptr = b;
+    }
+
+end_if:
+    isw = VpGetSign(a) + sw *VpGetSign(b);
+    /*
+     *  isw = 0 ...( 1)+(-1),( 1)-( 1),(-1)+(1),(-1)-(-1)
+     *   = 2 ...( 1)+( 1),( 1)-(-1)
+     *   =-2 ...(-1)+(-1),(-1)-( 1)
+     *   If isw==0, then c =(Sign a_ptr)(|a_ptr|-|b_ptr|)
+     *     else c =(Sign of isw)(|a_ptr|+|b_ptr|)
+    */
+    if(isw) {            /* addition */
+        VpSetSign(c,(S_INT)1);
+        VpAddAbs(a_ptr, b_ptr, c);
+        VpSetSign(c,isw / 2);
+    } else {            /* subtraction */
+        VpSetSign(c,(S_INT)1);
+        VpSubAbs(a_ptr, b_ptr, c);
+        if(a_ptr == a) {
+            VpSetSign(c,VpGetSign(a));
+        } else    {
+            VpSetSign(c,VpGetSign(a_ptr) * sw);
+        }
+    }
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpAddSub(result) c=% \n", c);
+        VPrint(stdout, "     a=% \n", a);
+        VPrint(stdout, "     b=% \n", b);
+        printf(" operation=%d\n", operation);
+    }
+#endif /* _DEBUG */
+    return c->Prec*BASE_FIG;
+}
+
+/*
+ * Addition of two variable precisional variables
+ * a and b assuming abs(a)>abs(b).
+ *   c = abs(a) + abs(b) ; where |a|>=|b|
+ */
+static int
+VpAddAbs(Real *a, Real *b, Real *c)
+{
+    U_LONG word_shift;
+    U_LONG round;
+    U_LONG carry;
+    U_LONG ap;
+    U_LONG bp;
+    U_LONG cp;
+    U_LONG a_pos;
+    U_LONG b_pos;
+    U_LONG c_pos;
+    U_LONG av, bv;
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpAddAbs called: a = %\n", a);
+        VPrint(stdout, "     b = %\n", b);
+    }
+#endif /* _DEBUG */
+
+    word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
+    a_pos = ap;
+    b_pos = bp;
+    c_pos = cp;
+    if(word_shift==-1L) return 0; /* Overflow */
+    if(b_pos == -1L) goto Assign_a;
+
+    round =((av + bv) >= HALF_BASE) ? 1 : 0;
+
+    /* Just assign the last few digits of b to c because a has no  */
+    /* corresponding digits to be added. */
+    while(b_pos + word_shift > a_pos) {
+        --c_pos;
+        if(b_pos > 0) {
+            --b_pos;
+            c->frac[c_pos] = b->frac[b_pos];
+        } else {
+            --word_shift;
+            c->frac[c_pos] = 0;
+        }
+    }
+
+    /* Just assign the last few digits of a to c because b has no */
+    /* corresponding digits to be added. */
+    bv = b_pos + word_shift;
+    while(a_pos > bv) {
+        --c_pos;
+        --a_pos;
+        c->frac[c_pos] = a->frac[a_pos];
+    }
+    carry = 0;    /* set first carry be zero */
+
+    /* Now perform addition until every digits of b will be */
+    /* exhausted. */
+    while(b_pos > 0) {
+        --a_pos;
+        --b_pos;
+        --c_pos;
+        c->frac[c_pos] = a->frac[a_pos] + b->frac[b_pos] + carry;
+        if(c->frac[c_pos] >= BASE) {
+            c->frac[c_pos] -= BASE;
+            carry = 1;
+        } else {
+            carry = 0;
+        }
+    }
+
+    /* Just assign the first few digits of a with considering */
+    /* the carry obtained so far because b has been exhausted. */
+    while(a_pos > 0) {
+        --a_pos;
+        --c_pos;
+        c->frac[c_pos] = a->frac[a_pos] + carry;
+        if(c->frac[c_pos] >= BASE) {
+            c->frac[c_pos] -= BASE;
+            carry = 1;
+        } else {
+            carry = 0;
+        }
+    }
+    if(c_pos) c->frac[c_pos - 1] += carry;
+
+    if(round) VpRdup(c);        /* Roundup and normalize. */
+    else       VpNmlz(c);        /* normalize the result */
+    goto Exit;
+
+Assign_a:
+    VpAsgn(c, a, 1);
+
+Exit:
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpAddAbs exit: c=% \n", c);
+    }
+#endif /* _DEBUG */
+    return 1;
+}
+
+/*
+ * c = abs(a) - abs(b)
+ */
+static int
+VpSubAbs(Real *a, Real *b, Real *c)
+{
+    U_LONG word_shift;
+    U_LONG round;
+    U_LONG borrow;
+    U_LONG ap;
+    U_LONG bp;
+    U_LONG cp;
+    U_LONG a_pos;
+    U_LONG b_pos;
+    U_LONG c_pos;
+    U_LONG av, bv;
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpSubAbs called: a = %\n", a);
+        VPrint(stdout, "     b = %\n", b);
+    }
+#endif /* _DEBUG */
+
+    word_shift = VpSetPTR(a, b, c, &ap, &bp, &cp, &av, &bv);
+    a_pos = ap;
+    b_pos = bp;
+    c_pos = cp;
+    if(word_shift==-1L) return 0; /* Overflow */
+    if(b_pos == -1L) goto Assign_a;
+
+    if(av >= bv) {
+        round =((av -= bv) >= HALF_BASE) ? 1 : 0;
+        borrow = 0;
+    } else {
+        round = 0;
+        borrow = 1;
+    }
+
+    /* Just assign the values which are the BASE subtracted by   */
+    /* each of the last few digits of the b because the a has no */
+    /* corresponding digits to be subtracted. */
+    if(b_pos + word_shift > a_pos) {
+        borrow = 1;
+        --c_pos;
+        --b_pos;
+        c->frac[c_pos] = BASE - b->frac[b_pos];
+        while(b_pos + word_shift > a_pos) {
+            --c_pos;
+            if(b_pos > 0) {
+                --b_pos;
+                c->frac[c_pos] = BASE - b->frac[b_pos] - borrow;
+            } else {
+                --word_shift;
+                c->frac[c_pos] = BASE - borrow;
+            }
+        }
+    }
+    /* Just assign the last few digits of a to c because b has no */
+    /* corresponding digits to subtract. */
+
+    bv = b_pos + word_shift;
+    while(a_pos > bv) {
+        --c_pos;
+        --a_pos;
+        c->frac[c_pos] = a->frac[a_pos];
+    }
+
+    /* Now perform subtraction until every digits of b will be */
+    /* exhausted. */
+    while(b_pos > 0) {
+        --a_pos;
+        --b_pos;
+        --c_pos;
+        if(a->frac[a_pos] < b->frac[b_pos] + borrow) {
+            c->frac[c_pos] = BASE + a->frac[a_pos] - b->frac[b_pos] - borrow;
+            borrow = 1;
+        } else {
+            c->frac[c_pos] = a->frac[a_pos] - b->frac[b_pos] - borrow;
+            borrow = 0;
+        }
+    }
+
+    /* Just assign the first few digits of a with considering */
+    /* the borrow obtained so far because b has been exhausted. */
+    while(a_pos > 0) {
+        --c_pos;
+        --a_pos;
+        if(a->frac[a_pos] < borrow) {
+            c->frac[c_pos] = BASE + a->frac[a_pos] - borrow;
+            borrow = 1;
+        } else {
+            c->frac[c_pos] = a->frac[a_pos] - borrow;
+            borrow = 0;
+        }
+    }
+    if(c_pos) c->frac[c_pos - 1] -= borrow;
+    if(round) VpRdup(c);        /* Round up and normalize */
+    else   VpNmlz(c);        /* normalize the result */
+    goto Exit;
+
+Assign_a:
+    VpAsgn(c, a, 1);
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpSubAbs exit: c=% \n", c);
+    }
+#endif /* _DEBUG */
+    return 1;
+}
+
+/*
+ * Note: If(av+bv)>= HALF_BASE,then 1 will be added to the least significant
+ *    digit of c(In case of addition).
+ * ------------------------- figure of output -----------------------------------
+ *      a =  xxxxxxxxxxx
+ *      b =    xxxxxxxxxx
+ *      c =xxxxxxxxxxxxxxx
+ *   word_shift =  |   |
+ *   right_word =  |    | (Total digits in RHSV)
+ *   left_word  = |   |   (Total digits in LHSV)
+ *   a_pos   =    |
+ *   b_pos   =     |
+ *   c_pos   =      |
+ */
+static U_LONG
+VpSetPTR(Real *a, Real *b, Real *c, U_LONG *a_pos, U_LONG *b_pos, U_LONG *c_pos, U_LONG *av, U_LONG *bv)
+{
+    U_LONG left_word, right_word, word_shift;
+    c->frac[0] = 0;
+    *av = *bv = 0;
+    word_shift =((a->exponent) -(b->exponent));
+    left_word = b->Prec + word_shift;
+    right_word = Max((a->Prec),left_word);
+    left_word =(c->MaxPrec) - 1;    /* -1 ... prepare for round up */
+    /*
+     * check if 'round off' is needed.
+     */
+    if(right_word > left_word) {    /* round off ? */
+        /*---------------------------------
+         *  Actual size of a = xxxxxxAxx
+         *  Actual size of b = xxxBxxxxx
+         *  Max. size of   c = xxxxxx
+         *  Round off  =   |-----|
+         *  c_pos   =   |
+         *  right_word    =   |
+         *  a_pos   =    |
+         */
+        *c_pos = right_word = left_word + 1;    /* Set resulting precision */
+        /* be equal to that of c */
+        if((a->Prec) >=(c->MaxPrec)) {
+            /*
+             *   a =  xxxxxxAxxx
+             *   c =  xxxxxx
+             *  a_pos =  |
+             */
+            *a_pos = left_word;
+            *av = a->frac[*a_pos];    /* av is 'A' shown in above. */
+        } else {
+            /*
+             *   a = xxxxxxx
+             *   c = xxxxxxxxxx
+             *  a_pos =     |
+             */
+            *a_pos = a->Prec;
+        }
+        if((b->Prec + word_shift) >= c->MaxPrec) {
+            /*
+             *   a = xxxxxxxxx
+             *   b =  xxxxxxxBxxx
+             *   c = xxxxxxxxxxx
+             *  b_pos =   |
+             */
+            if(c->MaxPrec >=(word_shift + 1)) {
+                *b_pos = c->MaxPrec - word_shift - 1;
+                *bv = b->frac[*b_pos];
+            } else {
+                *b_pos = -1L;
+            }
+        } else {
+            /*
+             *   a = xxxxxxxxxxxxxxxx
+             *   b =  xxxxxx
+             *   c = xxxxxxxxxxxxx
+             *  b_pos =     |
+             */
+            *b_pos = b->Prec;
+        }
+    } else {            /* The MaxPrec of c - 1 > The Prec of a + b  */
+        /*
+         *    a =   xxxxxxx
+         *    b =   xxxxxx
+         *    c = xxxxxxxxxxx
+         *   c_pos =   |
+         */
+        *b_pos = b->Prec;
+        *a_pos = a->Prec;
+        *c_pos = right_word + 1;
+    }
+    c->Prec = *c_pos;
+    c->exponent = a->exponent;
+    if(!AddExponent(c,(S_LONG)1)) return (-1L);
+    return word_shift;
+}
+
+/*
+ * Return number og significant digits
+ *       c = a * b , Where a = a0a1a2 ... an
+ *             b = b0b1b2 ... bm
+ *             c = c0c1c2 ... cl
+ *          a0 a1 ... an   * bm
+ *       a0 a1 ... an   * bm-1
+ *         .   .    .
+ *       .   .   .
+ *        a0 a1 .... an    * b0
+ *      +_____________________________
+ *     c0 c1 c2  ......  cl
+ *     nc      <---|
+ *     MaxAB |--------------------|
+ */
+VP_EXPORT int
+VpMult(Real *c, Real *a, Real *b)
+{
+    U_LONG MxIndA, MxIndB, MxIndAB, MxIndC;
+    U_LONG ind_c, i, nc;
+    U_LONG ind_as, ind_ae, ind_bs, ind_be;
+    U_LONG Carry, s;
+    Real *w;
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpMult(Enter): a=% \n", a);
+        VPrint(stdout, "      b=% \n", b);
+    }
+#endif /* _DEBUG */
+
+    if(!VpIsDefOP(c,a,b,3)) return 0; /* No significant digit */
+
+    if(VpIsZero(a) || VpIsZero(b)) {
+        /* at least a or b is zero */
+        VpSetZero(c,VpGetSign(a)*VpGetSign(b));
+        return 1; /* 0: 1 significant digit */
+    }
+
+    if((a->Prec == 1) &&(a->frac[0] == 1) &&(a->exponent == 1)) {
+        VpAsgn(c, b, VpGetSign(a));
+        goto Exit;
+    }
+    if((b->Prec == 1) &&(b->frac[0] == 1) &&(b->exponent == 1)) {
+        VpAsgn(c, a, VpGetSign(b));
+        goto Exit;
+    }
+    if((b->Prec) >(a->Prec)) {
+        /* Adjust so that digits(a)>digits(b) */
+        w = a;
+        a = b;
+        b = w;
+    }
+    w = NULL;
+    MxIndA = a->Prec - 1;
+    MxIndB = b->Prec - 1;
+    MxIndC = c->MaxPrec - 1;
+    MxIndAB = a->Prec + b->Prec - 1;
+
+    if(MxIndC < MxIndAB) {    /* The Max. prec. of c < Prec(a)+Prec(b) */
+        w = c;
+        c = VpAlloc((U_LONG)((MxIndAB + 1) * BASE_FIG), "#0");
+        MxIndC = MxIndAB;
+    }
+
+    /* set LHSV c info */
+
+    c->exponent = a->exponent;    /* set exponent */
+    if(!AddExponent(c,b->exponent)) return 0;
+    VpSetSign(c,VpGetSign(a)*VpGetSign(b));    /* set sign  */
+    Carry = 0;
+    nc = ind_c = MxIndAB;
+    for(i = 0; i <= nc; i++) c->frac[i] = 0;        /* Initialize c  */
+    c->Prec = nc + 1;        /* set precision */
+    for(nc = 0; nc < MxIndAB; ++nc, --ind_c) {
+        if(nc < MxIndB) {    /* The left triangle of the Fig. */
+            ind_as = MxIndA - nc;
+            ind_ae = MxIndA;
+            ind_bs = MxIndB;
+            ind_be = MxIndB - nc;
+        } else if(nc <= MxIndA) {    /* The middle rectangular of the Fig. */
+            ind_as = MxIndA - nc;
+            ind_ae = MxIndA -(nc - MxIndB);
+            ind_bs = MxIndB;
+            ind_be = 0;
+        } else if(nc > MxIndA) {    /*  The right triangle of the Fig. */
+            ind_as = 0;
+            ind_ae = MxIndAB - nc - 1;
+            ind_bs = MxIndB -(nc - MxIndA);
+            ind_be = 0;
+        }
+
+        s = 0L;
+        for(i = ind_as; i <= ind_ae; ++i) s +=((a->frac[i]) *(b->frac[ind_bs--]));
+        Carry = s / BASE;
+        s = s -(Carry * BASE);
+
+        c->frac[ind_c] += s;
+        if(c->frac[ind_c] >= BASE) {
+            s = c->frac[ind_c] / BASE;
+            Carry += s;
+            c->frac[ind_c] -=(s * BASE);
+        }
+        i = ind_c;
+        if(Carry) {
+            while((--i) >= 0) {
+                c->frac[i] += Carry;
+                if(c->frac[i] >= BASE) {
+                    Carry = c->frac[i] / BASE;
+                    c->frac[i] -=(Carry * BASE);
+                } else {
+                    break;
+                }
+            }
+        }
+    }
+
+    VpNmlz(c);            /* normalize the result */
+    if(w != NULL) {        /* free work variable */
+        VpAsgn(w, c, 2);
+        VpFree(c);
+        c = w;
+    }
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpMult(c=a*b): c=% \n", c);
+        VPrint(stdout, "      a=% \n", a);
+        VPrint(stdout, "      b=% \n", b);
+    }
+#endif /*_DEBUG */
+    return c->Prec*BASE_FIG;
+}
+
+/*
+ *   c = a / b,  remainder = r
+ */
+VP_EXPORT int
+VpDivd(Real *c, Real *r, Real *a, Real *b)
+{
+    U_LONG word_a, word_b, word_c, word_r;
+    U_LONG i, n, ind_a, ind_b, ind_c, ind_r;
+    U_LONG nLoop;
+    U_LONG q, b1, b1p1, b1b2, b1b2p1, r1r2;
+    U_LONG borrow, borrow1, borrow2, qb;
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, " VpDivd(c=a/b)  a=% \n", a);
+        VPrint(stdout, "    b=% \n", b);
+    }
+#endif /*_DEBUG */
+
+    VpSetNaN(r);
+    if(!VpIsDefOP(c,a,b,4)) goto Exit;
+    if(VpIsZero(a)&&VpIsZero(b)) {
+        VpSetNaN(c);
+        return VpException(VP_EXCEPTION_NaN,"(VpDivd) 0/0 not defined(NaN)",0);
+    }
+    if(VpIsZero(b)) {
+        VpSetInf(c,VpGetSign(a)*VpGetSign(b));
+        return VpException(VP_EXCEPTION_ZERODIVIDE,"(VpDivd) Divide by zero",0);
+    }
+    if(VpIsZero(a)) {
+        /* numerator a is zero  */
+        VpSetZero(c,VpGetSign(a)*VpGetSign(b));
+        VpSetZero(r,VpGetSign(a)*VpGetSign(b));
+        goto Exit;
+    }
+
+    if((b->Prec == 1) &&(b->frac[0] == 1) &&(b->exponent == 1)) {
+        /* divide by one  */
+        VpAsgn(c, a, VpGetSign(b));
+        VpSetZero(r,VpGetSign(a));
+        goto Exit;
+    }
+
+    word_a = a->Prec;
+    word_b = b->Prec;
+    word_c = c->MaxPrec;
+    word_r = r->MaxPrec;
+
+    ind_c = 0;
+    ind_r = 1;
+
+    if(word_a >= word_r) goto space_error;
+
+    r->frac[0] = 0;
+    while(ind_r <= word_a) {
+        r->frac[ind_r] = a->frac[ind_r - 1];
+        ++ind_r;
+    }
+
+    while(ind_r < word_r) r->frac[ind_r++] = 0;
+    while(ind_c < word_c) c->frac[ind_c++] = 0;
+
+    /* initial procedure */
+    b1 = b1p1 = b->frac[0];
+    if(b->Prec <= 1) {
+        b1b2p1 = b1b2 = b1p1 * BASE;
+    } else {
+        b1p1 = b1 + 1;
+        b1b2p1 = b1b2 = b1 * BASE + b->frac[1];
+        if(b->Prec > 2) ++b1b2p1;
+    }
+
+    /* */
+    /* loop start */
+    ind_c = word_r - 1;
+    nLoop = Min(word_c,ind_c);
+    ind_c = 1;
+    while(ind_c < nLoop) {
+        if(r->frac[ind_c] == 0) {
+            ++ind_c;
+            continue;
+        }
+        r1r2 = r->frac[ind_c] * BASE + r->frac[ind_c + 1];
+        if(r1r2 == b1b2) {
+            /* The first two word digits is the same */
+            ind_b = 2;
+            ind_a = ind_c + 2;
+            while(ind_b < word_b) {
+                if(r->frac[ind_a] < b->frac[ind_b]) goto div_b1p1;
+                if(r->frac[ind_a] > b->frac[ind_b]) break;
+                ++ind_a;
+                ++ind_b;
+            }
+            /* The first few word digits of r and b is the same and */
+            /* the first different word digit of w is greater than that */
+            /* of b, so quotinet is 1 and just subtract b from r. */
+            borrow = 0;        /* quotient=1, then just r-b */
+            ind_b = b->Prec - 1;
+            ind_r = ind_c + ind_b;
+            if(ind_r >= word_r) goto space_error;
+            n = ind_b;
+            for(i = 0; i <= n; ++i) {
+                if(r->frac[ind_r] < b->frac[ind_b] + borrow) {
+                    r->frac[ind_r] +=(BASE -(b->frac[ind_b] + borrow));
+                    borrow = 1;
+                } else {
+                    r->frac[ind_r] = r->frac[ind_r] - b->frac[ind_b] - borrow;
+                    borrow = 0;
+                }
+                --ind_r;
+                --ind_b;
+            }
+            ++(c->frac[ind_c]);
+            goto carry;
+        }
+        /* The first two word digits is not the same, */
+        /* then compare magnitude, and divide actually. */
+        if(r1r2 >= b1b2p1) {
+            q = r1r2 / b1b2p1;
+            c->frac[ind_c] += q;
+            ind_r = b->Prec + ind_c - 1;
+            goto sub_mult;
+        }
+
+div_b1p1:
+        if(ind_c + 1 >= word_c) goto out_side;
+        q = r1r2 / b1p1;
+        c->frac[ind_c + 1] += q;
+        ind_r = b->Prec + ind_c;
+
+sub_mult:
+        borrow1 = borrow2 = 0;
+        ind_b = word_b - 1;
+        if(ind_r >= word_r) goto space_error;
+        n = ind_b;
+        for(i = 0; i <= n; ++i) {
+            /* now, perform r = r - q * b */
+            qb = q *(b->frac[ind_b]);
+            if(qb < BASE) borrow1 = 0;
+            else {
+                borrow1 = qb / BASE;
+                qb = qb - borrow1 * BASE;
+            }
+            if(r->frac[ind_r] < qb) {
+                r->frac[ind_r] +=(BASE - qb);
+                borrow2 = borrow2 + borrow1 + 1;
+            } else {
+                r->frac[ind_r] -= qb;
+                borrow2 += borrow1;
+            }
+            if(borrow2) {
+                if(r->frac[ind_r - 1] < borrow2) {
+                    r->frac[ind_r - 1] +=(BASE - borrow2);
+                    borrow2 = 1;
+                } else {
+                    r->frac[ind_r - 1] -= borrow2;
+                    borrow2 = 0;
+                }
+            }
+            --ind_r;
+            --ind_b;
+        }
+
+        r->frac[ind_r] -= borrow2;
+carry:
+        ind_r = ind_c;
+        while(c->frac[ind_r] >= BASE) {
+            c->frac[ind_r] -= BASE;
+            --ind_r;
+            ++(c->frac[ind_r]);
+        }
+    }
+    /* End of operation, now final arrangement */
+out_side:
+    c->Prec = word_c;
+    c->exponent = a->exponent;
+    if(!AddExponent(c,(S_LONG)2))   return 0;
+    if(!AddExponent(c,-(b->exponent))) return 0;
+
+    VpSetSign(c,VpGetSign(a)*VpGetSign(b));
+    VpNmlz(c);            /* normalize c */
+    r->Prec = word_r;
+    r->exponent = a->exponent;
+    if(!AddExponent(r,(S_LONG)1)) return 0;
+    VpSetSign(r,VpGetSign(a));
+    VpNmlz(r);            /* normalize r(remainder) */
+    goto Exit;
+
+space_error:
+    rb_fatal("ERROR(VpDivd): space for remainder too small.\n");
+#ifdef _DEBUG
+    if(gfDebug) {
+        printf("   word_a=%d\n", word_a);
+        printf("   word_b=%d\n", word_b);
+        printf("   word_c=%d\n", word_c);
+        printf("   word_r=%d\n", word_r);
+        printf("   ind_r =%d\n", ind_r);
+    }
+#endif /* _DEBUG */
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, " VpDivd(c=a/b), c=% \n", c);
+        VPrint(stdout, "    r=% \n", r);
+    }
+#endif /* _DEBUG */
+    return c->Prec*BASE_FIG;
+}
+
+/*
+ *  Input  a = 00000xxxxxxxx En(5 preceeding zeros)
+ *  Output a = xxxxxxxx En-5
+ */
+static int
+VpNmlz(Real *a)
+{
+    U_LONG ind_a, i, j;
+
+    if(VpIsZero(a)) {
+        VpSetZero(a,VpGetSign(a));
+        return 1;
+    }
+    ind_a = a->Prec;
+    while(ind_a--) {
+        if(a->frac[ind_a]) {
+            a->Prec = ind_a + 1;
+            i = j = 0;
+            while(a->frac[i] == 0) ++i;        /* skip the first few zeros */
+            if(i) {
+                a->Prec -= i;
+                if(!AddExponent(a,-((S_INT)i))) return 0;
+                while(i <= ind_a) {
+                    a->frac[j] = a->frac[i];
+                    ++i;
+                    ++j;
+                }
+            }
+#ifdef _DEBUG
+            if(gfCheckVal)    VpVarCheck(a);
+#endif /* _DEBUG */
+            return 1;
+        }
+    }
+    /* a is zero(no non-zero digit) */
+    VpSetZero(a,VpGetSign(a));
+    return 1;
+}
+
+/*
+ *  VpComp = 0  ... if a=b,
+ *   Pos  ... a>b,
+ *   Neg  ... a<b.
+ *   999  ... result undefined(NaN)
+ */
+VP_EXPORT int
+VpComp(Real *a, Real *b)
+{
+    int val;
+    U_LONG mx, ind;
+    int e;
+    val = 0;
+    if(VpIsNaN(a)||VpIsNaN(b)) return 999;
+    if(!VpIsDef(a)) {
+        if(!VpIsDef(b)) e = a->sign - b->sign;
+        else             e = a->sign;
+        if(e>0)   return  1;
+        else if(e<0) return -1;
+        else   return  0;
+    }
+    if(!VpIsDef(b)) {
+        e = -b->sign;
+        if(e>0) return  1;
+        else return -1;
+    }
+    /* Zero check */
+    if(VpIsZero(a)) {
+        if(VpIsZero(b))      return 0; /* both zero */
+        val = -VpGetSign(b);
+        goto Exit;
+    }
+    if(VpIsZero(b)) {
+        val = VpGetSign(a);
+        goto Exit;
+    }
+
+    /* compare sign */
+    if(VpGetSign(a) > VpGetSign(b)) {
+        val = 1;        /* a>b */
+        goto Exit;
+    }
+    if(VpGetSign(a) < VpGetSign(b)) {
+        val = -1;        /* a<b */
+        goto Exit;
+    }
+
+    /* a and b have same sign, && signe!=0,then compare exponent */
+    if((a->exponent) >(b->exponent)) {
+        val = VpGetSign(a);
+        goto Exit;
+    }
+    if((a->exponent) <(b->exponent)) {
+        val = -VpGetSign(b);
+        goto Exit;
+    }
+
+    /* a and b have same exponent, then compare significand. */
+    mx =((a->Prec) <(b->Prec)) ?(a->Prec) :(b->Prec);
+    ind = 0;
+    while(ind < mx) {
+        if((a->frac[ind]) >(b->frac[ind])) {
+            val = VpGetSign(a);
+         goto Exit;
+        }
+        if((a->frac[ind]) <(b->frac[ind])) {
+            val = -VpGetSign(b);
+            goto Exit;
+        }
+        ++ind;
+    }
+    if((a->Prec) >(b->Prec)) {
+        val = VpGetSign(a);
+    } else if((a->Prec) <(b->Prec)) {
+        val = -VpGetSign(b);
+    }
+
+Exit:
+    if  (val> 1) val =  1;
+    else if(val<-1) val = -1;
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, " VpComp a=%\n", a);
+        VPrint(stdout, "  b=%\n", b);
+        printf("  ans=%d\n", val);
+    }
+#endif /* _DEBUG */
+    return (int)val;
+}
+
+/*
+ *    cntl_chr ... ASCIIZ Character, print control characters
+ *     Available control codes:
+ *      %  ... VP variable. To print '%', use '%%'.
+ *      \n ... new line
+ *      \b ... backspace
+ *           ... tab
+ *     Note: % must must not appear more than once
+ *    a  ... VP variable to be printed
+ */
+VP_EXPORT int
+VPrint(FILE *fp, char *cntl_chr, Real *a)
+{
+    U_LONG i, j, nc, nd, ZeroSup;
+    U_LONG n, m, e, nn;
+
+    /* Check if NaN & Inf. */
+    if(VpIsNaN(a)) {
+        fprintf(fp,"NaN");
+        return 8;
+    }
+    if(VpIsPosInf(a)) {
+        fprintf(fp,"Infinity");
+        return 8;
+    }
+    if(VpIsNegInf(a)) {
+        fprintf(fp,"-Infinity");
+        return 9;
+    }
+    if(VpIsZero(a)) {
+        fprintf(fp,"0.0");
+        return 3;
+    }
+
+    j = 0;
+    nd = nc = 0;        /*  nd : number of digits in fraction part(every 10 digits, */
+    /*    nd<=10). */
+    /*  nc : number of caracters printed  */
+    ZeroSup = 1;        /* Flag not to print the leading zeros as 0.00xxxxEnn */
+    while(*(cntl_chr + j)) {
+        if((*(cntl_chr + j) == '%') &&(*(cntl_chr + j + 1) != '%')) {
+         nc = 0;
+         if(!VpIsZero(a)) {
+                if(VpGetSign(a) < 0) {
+                    fprintf(fp, "-");
+                    ++nc;
+                }
+                nc += fprintf(fp, "0.");
+                n = a->Prec;
+                for(i=0;i < n;++i) {
+                 m = BASE1;
+                    e = a->frac[i];
+                    while(m) {
+                        nn = e / m;
+                        if((!ZeroSup) || nn) {
+                            nc += fprintf(fp, "%u", nn);    /* The reading zero(s) */
+                            /* as 0.00xx will not */
+                            /* be printed. */
+                            ++nd;
+                            ZeroSup = 0;    /* Set to print succeeding zeros */
+                        }
+                        if(nd >= 10) {    /* print ' ' after every 10 digits */
+                            nd = 0;
+                            nc += fprintf(fp, " ");
+                        }
+                        e = e - nn * m;
+                        m /= 10;
+                    }
+                }
+                nc += fprintf(fp, "E%d", VpExponent10(a));
+            } else {
+                nc += fprintf(fp, "0.0");
+            }
+        } else {
+            ++nc;
+            if(*(cntl_chr + j) == '\\') {
+                switch(*(cntl_chr + j + 1)) {
+                case 'n':
+                    fprintf(fp, "\n");
+                    ++j;
+                    break;
+                case 't':
+                    fprintf(fp, "\t");
+                    ++j;
+                 break;
+                case 'b':
+                    fprintf(fp, "\n");
+                    ++j;
+                    break;
+                default:
+                    fprintf(fp, "%c", *(cntl_chr + j));
+                    break;
+                }
+            } else {
+                fprintf(fp, "%c", *(cntl_chr + j));
+                if(*(cntl_chr + j) == '%') ++j;
+            }
+        }
+        j++;
+    }
+    return (int)nc;
+}
+
+static void
+VpFormatSt(char *psz,S_INT fFmt)
+{
+    U_LONG ie;
+    U_LONG i, j;
+    S_INT nf;
+    char ch;
+    int fDot = 0;
+
+    ie = strlen(psz);
+    for(i = 0; i < ie; ++i) {
+        ch = psz[i];
+        if(!ch) break;
+        if(ch == '.') {
+            nf = 0;
+            fDot = 1;
+            continue;
+        }
+        if(!fDot)  continue;
+        if(ch == 'E') break;
+        nf++;
+        if(nf > fFmt) {
+            for(j = ie; j >= i; --j)
+            psz[j + 1] = psz[j];
+            ++ie;
+            nf = 0;
+            psz[i] = ' ';
+        }
+    }
+}
+
+VP_EXPORT S_LONG
+VpExponent10(Real *a)
+{
+    S_LONG ex;
+    U_LONG n;
+
+    if(!VpIsDef(a)) return 0;
+    if(VpIsZero(a)) return 0;
+
+    ex =(a->exponent) * BASE_FIG;
+    n = BASE1;
+    while((a->frac[0] / n) == 0) {
+         --ex;
+         n /= 10;
+    }
+    return ex;
+}
+
+VP_EXPORT void
+VpSzMantissa(Real *a,char *psz)
+{
+    U_LONG i, ZeroSup;
+    U_LONG n, m, e, nn;
+
+    if(VpIsNaN(a)) {
+        sprintf(psz,"NaN");
+        return;
+    }
+    if(VpIsPosInf(a)) {
+        sprintf(psz,"Infinity");
+        return;
+    }
+    if(VpIsNegInf(a)) {
+        sprintf(psz,"-Infinity");
+        return;
+    }
+
+    ZeroSup = 1;        /* Flag not to print the leading zeros as 0.00xxxxEnn */
+    if(!VpIsZero(a)) {
+        if(VpGetSign(a) < 0) *psz++ = '-';
+        n = a->Prec;
+        for(i=0;i < n;++i) {
+            m = BASE1;
+            e = a->frac[i];
+            while(m) {
+                nn = e / m;
+                if((!ZeroSup) || nn) {
+                    sprintf(psz, "%u", nn);    /* The reading zero(s) */
+                    psz += strlen(psz);
+                    /* as 0.00xx will be ignored. */
+                    ZeroSup = 0;    /* Set to print succeeding zeros */
+                }
+                e = e - nn * m;
+                m /= 10;
+            }
+        }
+        *psz = 0;
+    } else {
+        if(VpIsPosZero(a)) sprintf(psz, "0");
+        else      sprintf(psz, "-0");
+    }
+}
+
+VP_EXPORT void
+VpToString(Real *a,char *psz,int fFmt)
+{
+    U_LONG i, ZeroSup;
+    U_LONG n, m, e, nn;
+    char *pszSav = psz;
+    S_LONG ex;
+
+    if(VpIsNaN(a)) {
+        sprintf(psz,"NaN");
+        return;
+    }
+    if(VpIsPosInf(a)) {
+        sprintf(psz,"Infinity");
+        return;
+    }
+    if(VpIsNegInf(a)) {
+        sprintf(psz,"-Infinity");
+        return;
+    }
+
+    ZeroSup = 1;    /* Flag not to print the leading zeros as 0.00xxxxEnn */
+    if(!VpIsZero(a)) {
+        if(VpGetSign(a) < 0) *psz++ = '-';
+        *psz++ = '0';
+        *psz++ = '.';
+        n = a->Prec;
+        for(i=0;i < n;++i) {
+            m = BASE1;
+            e = a->frac[i];
+            while(m) {
+                nn = e / m;
+                if((!ZeroSup) || nn) {
+                    sprintf(psz, "%u", nn);    /* The reading zero(s) */
+                    psz += strlen(psz);
+                    /* as 0.00xx will be ignored. */
+                    ZeroSup = 0;    /* Set to print succeeding zeros */
+                }
+                e = e - nn * m;
+                m /= 10;
+            }
+        }
+        ex =(a->exponent) * BASE_FIG;
+        n = BASE1;
+        while((a->frac[0] / n) == 0) {
+            --ex;
+            n /= 10;
+        }
+        sprintf(psz, "E%d", ex);
+    } else {
+        if(VpIsPosZero(a)) sprintf(psz, "0.0");
+        else      sprintf(psz, "-0.0");
+    }
+    if(fFmt) VpFormatSt(pszSav, fFmt);
+}
+
+/*
+ *  [Output]
+ *   a[]  ... variable to be assigned the value.
+ *  [Input]
+ *   int_chr[]  ... integer part(may include '+/-').
+ *   ni   ... number of characters in int_chr[],not including '+/-'.
+ *   frac[]  ... fraction part.
+ *   nf   ... number of characters in frac[].
+ *   exp_chr[]  ... exponent part(including '+/-').
+ *   ne   ... number of characters in exp_chr[],not including '+/-'.
+ */
+VP_EXPORT int
+VpCtoV(Real *a, char *int_chr, U_LONG ni, char *frac, U_LONG nf, char *exp_chr, U_LONG ne)
+{
+    U_LONG i, j, ind_a, ma, mi, me;
+    U_LONG loc;
+    S_INT  e,es, eb, ef;
+    S_INT  sign, signe;
+    /* get exponent part */
+    e = 0;
+    ma = a->MaxPrec;
+    mi = ni;
+    me = ne;
+    signe = 1;
+    for(i=0;i < ma;++i) a->frac[i] = 0;
+    if(ne > 0) {
+        i = 0;
+        if(exp_chr[0] == '-') {
+            signe = -1;
+            ++i;
+            ++me;
+        } else if(exp_chr[0] == '+') {
+            ++i;
+            ++me;
+        }
+        while(i < me) {
+            es = e*BASE_FIG;
+            e = e * 10 + exp_chr[i] - '0';
+            if(es>e*((S_INT)BASE_FIG)) {
+                return VpException(VP_EXCEPTION_INFINITY,"Exponent overflow",0);
+            }
+            ++i;
+        }
+    }
+
+    /* get integer part */
+    i = 0;
+    sign = 1;
+    if(ni > 0) {
+        if(int_chr[0] == '-') {
+            sign = -1;
+            ++i;
+            ++mi;
+        } else if(int_chr[0] == '+') {
+            ++i;
+            ++mi;
+        }
+    }
+
+    e = signe * e;        /* e: The value of exponent part. */
+    e = e + ni;        /* set actual exponent size. */
+
+    if(e > 0)    signe = 1;
+    else        signe = -1;
+
+    /* Adjust the exponent so that it is the multiple of BASE_FIG. */
+    j = 0;
+    ef = 1;
+    while(ef) {
+        if(e>=0) eb =  e;
+        else  eb = -e;
+        ef = eb / BASE_FIG;
+        ef = eb - ef * BASE_FIG;
+        if(ef) {
+            ++j;        /* Means to add one more preceeding zero */
+            ++e;
+        }
+    }
+
+    eb = e / ((S_INT)BASE_FIG);
+
+    ind_a = 0;
+    while(i < mi) {
+        a->frac[ind_a] = 0;
+        while((j < (U_LONG)BASE_FIG) &&(i < mi)) {
+            a->frac[ind_a] = a->frac[ind_a] * 10 + int_chr[i] - '0';
+            ++j;
+            ++i;
+        }
+        if(i < mi) {
+            ++ind_a;
+            if(ind_a >= ma) goto over_flow;
+            j = 0;
+        }
+    }
+    loc = 1;
+
+    /* get fraction part */
+
+    i = 0;
+    while(i < nf) {
+        while((j < (U_LONG)BASE_FIG) &&(i < nf)) {
+            a->frac[ind_a] = a->frac[ind_a] * 10 + frac[i] - '0';
+            ++j;
+            ++i;
+        }
+        if(i < nf) {
+            ++ind_a;
+            if(ind_a >= ma) goto over_flow;
+            j = 0;
+        }
+    }
+    goto Final;
+
+over_flow:
+    rb_warn("Conversion from String to BigDecimal overflow (last few digits discarded).");
+
+Final:
+    if(ind_a >= ma) ind_a = ma - 1;
+    while(j < (U_LONG)BASE_FIG) {
+        a->frac[ind_a] = a->frac[ind_a] * 10;
+        ++j;
+    }
+    a->Prec = ind_a + 1;
+    a->exponent = eb;
+    VpSetSign(a,sign);
+    VpNmlz(a);
+    return 1;
+}
+
+/*
+ * [Input]
+ *   *m  ... Real
+ * [Output]
+ *   *d  ... fraction part of m(d = 0.xxxxxxx). where # of 'x's is fig.
+ *   *e  ... U_LONG,exponent of m.
+ * DBLE_FIG ... Number of digits in a double variable.
+ *
+ *  m -> d*10**e, 0<d<BASE
+ */
+VP_EXPORT void
+VpVtoD(double *d, S_LONG *e, Real *m)
+{
+    U_LONG ind_m, mm, fig;
+    double div;
+
+    fig =(DBLE_FIG + BASE_FIG - 1) / BASE_FIG;
+    if(VpIsPosZero(m)) {
+        *d = 0.0;
+        *e = 0;
+        goto Exit;
+    } else
+    if(VpIsNegZero(m)) {
+        *d = VpGetDoubleNegZero();
+        *e = 0;
+        goto Exit;
+    }
+    ind_m = 0;
+    mm = Min(fig,(m->Prec));
+    *d = 0.0;
+    div = 1.;
+    while(ind_m < mm) {
+        div /=(double)((S_INT)BASE);
+        *d = *d +((double) ((S_INT)m->frac[ind_m++])) * div;
+    }
+    *e = m->exponent * ((S_INT)BASE_FIG);
+    *d *= VpGetSign(m);
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, " VpVtoD: m=%\n", m);
+        printf("   d=%e * 10 **%d\n", *d, *e);
+        printf("   DBLE_FIG = %d\n", DBLE_FIG);
+    }
+#endif /*_DEBUG */
+    return;
+}
+
+/*
+ * m <- d
+ */
+VP_EXPORT void
+VpDtoV(Real *m, double d)
+{
+    U_LONG i, ind_m, mm;
+    U_LONG ne;
+    double  val, val2;
+
+    if(isnan(d)) {
+        VpSetNaN(m);
+        goto Exit;
+    }
+    if(isinf(d)) {
+        if(d>0.0) VpSetPosInf(m);
+        else   VpSetNegInf(m);
+        goto Exit;
+    }
+
+    if(d == 0.0) {
+        VpSetZero(m,1);
+        goto Exit;
+    }
+    val =(d > 0.) ? d :(-d);
+    ne = 0;
+    if(val >= 1.0) {
+        while(val >= 1.0) {
+            val /=(double)((S_INT)BASE);
+            ++ne;
+        }
+    } else {
+        val2 = 1.0 /(double)((S_INT)BASE);
+        while(val < val2) {
+            val *=(double)((S_INT)BASE);
+            --ne;
+        }
+    }
+    /* Now val = 0.xxxxx*BASE**ne */
+
+    mm = m->MaxPrec;
+    for(ind_m = 0;ind_m < mm;ind_m++) m->frac[ind_m] = 0;
+    for(ind_m = 0;val > 0.0 && ind_m < mm;ind_m++) {
+        val *=(double)((S_INT)BASE);
+        i =(U_LONG) val;
+        val -=(double)((S_INT)i);
+        m->frac[ind_m] = i;
+    }
+    if(ind_m >= mm) ind_m = mm - 1;
+    if(d > 0.0) {
+        VpSetSign(m, (S_INT)1);
+    } else {
+        VpSetSign(m,-(S_INT)1);
+    }
+    m->Prec = ind_m + 1;
+    m->exponent = ne;
+    if(val*((double)((S_INT)BASE)) >=(double)((S_INT)HALF_BASE)) VpRdup(m);
+    VpNmlz(m);
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        printf("VpDtoV d=%30.30e\n", d);
+        VPrint(stdout, "  m=%\n", m);
+    }
+#endif /* _DEBUG */
+    return;
+}
+
+/*
+ *  m <- ival
+ */
+VP_EXPORT void
+VpItoV(Real *m, S_INT ival)
+{
+    U_LONG mm, ind_m;
+    U_LONG val, v1, v2, v;
+    int isign;
+    S_INT ne;
+
+    if(ival == 0) {
+        VpSetZero(m,1);
+        goto Exit;
+    }
+    isign = 1;
+    val = ival;
+    if(ival < 0) {
+        isign = -1;
+        val =(U_LONG)(-ival);
+    }
+    ne = 0;
+    ind_m = 0;
+    mm = m->MaxPrec;
+    while(ind_m < mm) {
+        m->frac[ind_m] = 0;
+        ++ind_m;
+    }
+    ind_m = 0;
+    while(val > 0) {
+        if(val) {
+         v1 = val;
+         v2 = 1;
+            while(v1 >= BASE) {
+                v1 /= BASE;
+                v2 *= BASE;
+            }
+            val = val - v2 * v1;
+            v = v1;
+        } else {
+            v = 0;
+        }
+        m->frac[ind_m] = v;
+        ++ind_m;
+        ++ne;
+    }
+    m->Prec = ind_m - 1;
+    m->exponent = ne;
+    VpSetSign(m,isign);
+    VpNmlz(m);
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        printf(" VpItoV i=%ld\n", ival);
+        VPrint(stdout, "  m=%\n", m);
+    }
+#endif /* _DEBUG */
+    return;
+}
+
+/*
+ * y = SQRT(x),  y*y - x =>0
+ */
+VP_EXPORT int
+VpSqrt(Real *y, Real *x)
+{
+    Real *f = NULL;
+    Real *r = NULL;
+    S_LONG y_prec, f_prec;
+	S_LONG n;
+	S_LONG e;
+    S_LONG prec;
+    S_LONG nr;
+    double val;
+
+    if(!VpIsDef(x)) {
+        VpAsgn(y,x,1);
+        goto Exit;
+    }
+
+    if(VpIsZero(x)) {
+        VpSetZero(y,VpGetSign(x));
+        goto Exit;
+    }
+
+    if(VpGetSign(x) < 0) {
+        VpSetZero(y,VpGetSign(x));
+        return VpException(VP_EXCEPTION_OP,"(VpSqrt) SQRT(negative valuw)",0);
+    }
+
+    n = (S_LONG)y->MaxPrec;
+    if((S_LONG)x->MaxPrec > n) n = (S_LONG)x->MaxPrec;
+    /* allocate temporally variables  */
+    f = VpAlloc(y->MaxPrec *(BASE_FIG + 2), "#1");
+    r = VpAlloc((n + n) *(BASE_FIG + 2), "#1");
+
+    nr = 0;
+    y_prec = (S_LONG)y->MaxPrec;
+    f_prec = (S_LONG)f->MaxPrec;
+
+    VpAsgn(y, x, 1);        /* assign initial guess. y <= x */
+    prec = x->exponent;
+    if(prec > 0)    ++prec;
+    else            --prec;
+    prec = prec / 2 - (S_LONG)y->MaxPrec;
+    /*
+     *  y  = 0.yyyy yyyy yyyy YYYY
+     *  BASE_FIG =   |  |
+     *  prec  =(0.YYYY*BASE-4)
+     */
+    VpVtoD(&val, &e, y);    /* val <- y  */
+    e /= ((S_LONG)BASE_FIG);
+    n = e / 2;
+    if(e - n * 2 != 0) {
+        val /=(double)((S_INT)BASE);
+        n =(e + 1) / 2;
+    }
+    VpDtoV(y, sqrt(val));    /* y <- sqrt(val) */
+    y->exponent += n;
+    n = (DBLE_FIG + BASE_FIG - 1) / BASE_FIG;
+    y->MaxPrec = (U_LONG)Min(n , y_prec);
+    f->MaxPrec = y->MaxPrec + 1;
+    n = y_prec*((S_LONG)BASE_FIG);
+    if((U_LONG)n<maxnr) n = (U_LONG)maxnr;
+    do {
+        y->MaxPrec *= 2;
+        if(y->MaxPrec > (U_LONG)y_prec) y->MaxPrec = (U_LONG)y_prec;
+        f->MaxPrec = y->MaxPrec;
+        VpDivd(f, r, x, y);    /* f = x/y  */
+        VpAddSub(r, y, f, 1);    /* r = y + x/y  */
+        VpMult(f, VpPt5, r);    /* f = 0.5*r  */
+        VpAddSub(r, f, y, -1);
+        if(VpIsZero(r))         goto converge;
+        if(r->exponent <= prec) goto converge;
+        VpAsgn(y, f, 1);
+    } while(++nr < n);
+    /* */
+#ifdef _DEBUG
+    if(gfDebug) {
+        printf("ERROR(VpSqrt): did not converge within %d iterations.\n",
+            nr);
+    }
+#endif /* _DEBUG */
+    y->MaxPrec = y_prec;
+    goto Exit;
+
+converge:
+    VpChangeSign(y,(S_INT)1);
+#ifdef _DEBUG
+    if(gfDebug) {
+        VpMult(r, y, y);
+        VpAddSub(f, x, r, -1);
+        printf("VpSqrt: iterations = %d\n", nr);
+        VPrint(stdout, "  y =% \n", y);
+        VPrint(stdout, "  x =% \n", x);
+        VPrint(stdout, "  x-y*y = % \n", f);
+    }
+#endif /* _DEBUG */
+    y->MaxPrec = y_prec;
+
+Exit:
+    VpFree(f);
+    VpFree(r);
+    return 1;
+}
+
+/*
+ *
+ * f = 1: round, 2:ceil, 3: floor
+ *
+ */
+VP_EXPORT void
+VpRound(Real *y, Real *x, int sw, int f, int nf)
+{
+    int n,i,j,ix,ioffset;
+    U_LONG v;
+    U_LONG div;
+
+    if(!VpIsDef(x)) {
+        VpAsgn(y,x,1);
+        goto Exit;
+    }
+
+    /* First,assign whole value */
+    VpAsgn(y, x, sw);
+    nf += y->exponent*((int)BASE_FIG);
+    /* ix: x->fraq[ix] contains round position */
+    ix = (nf + ((int)BASE_FIG))/((int)BASE_FIG)-1;
+    if(ix<0 || ((U_LONG)ix)>=y->Prec) goto Exit; /* Unable to round */
+    ioffset = nf - ix*((int)BASE_FIG);
+    for(j=ix+1;j<(int)y->Prec;++j) y->frac[j] = 0;
+    VpNmlz(y);
+    v = y->frac[ix];
+    /* drop digits after pointed digit */
+    n = BASE_FIG - ioffset - 1;
+    for(i=0;i<n;++i) v /= 10;
+    div = v/10;
+    v = v - div*10;
+    switch(f){
+    case 1: /* Round */
+        if(sw==2 && v>=5) {
+            ++div;
+        }
+        break;
+    case 2: /* ceil */
+        if(v) {
+            if(VpGetSign(x)>0) ++div;
+        }
+        break;
+    case 3: /* floor */
+        if(v) {
+            if(VpGetSign(x)<0) ++div;
+        }
+        break;
+    }
+    for(i=0;i<=n;++i) div *= 10;
+    if(div>=BASE) {
+        y->frac[ix] = 0;
+        if(ix) {
+            VpNmlz(y);
+            VpRdup(y);
+        } else {
+            VpSetOne(y);
+            VpSetSign(y,VpGetSign(x));
+        }
+    } else {
+        y->frac[ix] = div;
+        VpNmlz(y);
+    }
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpRound y=%\n", y);
+        VPrint(stdout, "  x=%\n", x);
+    }
+#endif /*_DEBUG */
+    return;
+}
+
+/*
+ *  Rounds up m(plus one to final digit of m).
+ */
+static int
+VpRdup(Real *m)
+{
+    U_LONG ind_m, carry;
+    ind_m = m->Prec;
+    carry = 1;
+    while(carry > 0 && ind_m) {
+        --ind_m;
+        m->frac[ind_m] += carry;
+        if(m->frac[ind_m] >= BASE) m->frac[ind_m] -= BASE;
+        else                       carry = 0;
+    }
+    if(carry > 0) {        /* Overflow,count exponent and set fraction part be 1  */
+        if(!AddExponent(m,(S_LONG)1)) return 0;
+        m->Prec = m->frac[0] = 1;
+    } else {
+        VpNmlz(m);
+    }
+    return 1;
+}
+
+/*
+ *  y = x - fix(x)
+ */
+VP_EXPORT void
+VpFrac(Real *y, Real *x)
+{
+    U_LONG my, ind_y, ind_x;
+
+    if(!VpIsDef(x) || VpIsZero(x)) {
+        VpAsgn(y,x,1);
+        goto Exit;
+    }
+
+    if(x->exponent > 0 && (U_LONG)x->exponent >= x->Prec) {
+        VpSetZero(y,VpGetSign(x));
+        goto Exit;
+    } else if(x->exponent <= 0) {
+        VpAsgn(y, x, 1);
+        goto Exit;
+    }
+    y->Prec = x->Prec -(U_LONG) x->exponent;
+    y->Prec = Min(y->Prec, y->MaxPrec);
+    y->exponent = 0;
+    VpSetSign(y,VpGetSign(x));
+    ind_y = 0;
+    my = y->Prec;
+    ind_x = x->exponent;
+    while(ind_y <= my) {
+        y->frac[ind_y] = x->frac[ind_x];
+        ++ind_y;
+        ++ind_x;
+    }
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpFrac y=%\n", y);
+        VPrint(stdout, "    x=%\n", x);
+    }
+#endif /* _DEBUG */
+    return;
+}
+
+/*
+ *   y = x ** n
+ */
+VP_EXPORT int
+VpPower(Real *y, Real *x, S_INT n)
+{
+    U_LONG s, ss;
+    S_LONG sign;
+    Real *w1 = NULL;
+    Real *w2 = NULL;
+
+    if(VpIsZero(x)) {
+        if(n<0) n = -n;
+        VpSetZero(y,(n%2)?VpGetSign(x):(-VpGetSign(x)));
+        goto Exit;
+    }
+    if(!VpIsDef(x)) {
+        VpSetNaN(y); /* Not sure !!! */
+        goto Exit;
+    }
+
+    if((x->exponent == 1) &&(x->Prec == 1) &&(x->frac[0] == 1)) {
+        /* abs(x) = 1 */
+        VpSetOne(y);
+        if(VpGetSign(x) > 0) goto Exit;
+        if((n % 2) == 0) goto Exit;
+        VpSetSign(y,-(S_INT)1);
+        goto Exit;
+    }
+
+    if(n > 0) sign = 1;
+    else if(n < 0) {
+        sign = -1;
+        n = -n;
+    } else {
+        VpSetOne(y);
+        goto Exit;
+    }
+
+    /* Allocate working variables  */
+
+    w1 = VpAlloc((x->Prec + 2) * BASE_FIG, "#0");
+    w2 = VpAlloc((w1->MaxPrec * 2 + 1) * BASE_FIG, "#0");
+    /* calculation start */
+
+    VpAsgn(y, x, 1);
+    --n;
+    while(n > 0) {
+        VpAsgn(w1, x, 1);
+        s = 1;
+loop1:        ss = s;
+        s += s;
+        if(s >(U_LONG) n) goto out_loop1;
+        VpMult(w2, w1, w1);
+        VpAsgn(w1, w2, 1);
+        goto loop1;
+out_loop1:
+        n -= ss;
+        VpMult(w2, y, w1);
+        VpAsgn(y, w2, 1);
+    }
+    if(sign < 0) {
+        VpDivd(w1, w2, VpConstOne, y);
+        VpAsgn(y, w1, 1);
+    }
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "VpPower y=%\n", y);
+        VPrint(stdout, "VpPower x=%\n", x);
+        printf("  n=%d\n", n);
+    }
+#endif /* _DEBUG */
+    VpFree(w2);
+    VpFree(w1);
+    return 1;
+}
+
+/*
+ * Calculates pai(=3.141592653589793238462........).
+ */
+VP_EXPORT void
+VpPai(Real *y)
+{
+    Real *n, *n25, *n956, *n57121;
+    Real *r, *f, *t;
+    U_LONG p;
+    U_LONG nc;
+    U_LONG i1,i2;
+
+    p = y->MaxPrec *(BASE_FIG + 2) + 2;
+    if(p<maxnr) nc = maxnr;
+    else  nc = p;
+
+    /* allocate temporally variables  */
+    r = VpAlloc(p * 2, "#0");
+    f = VpAlloc(p, "#0");
+    t = VpAlloc(p, "#-80");
+
+    n = VpAlloc((U_LONG)10, "1");
+    n25 = VpAlloc((U_LONG)2, "-0.04"); /*-25");*/
+    n956 = VpAlloc((U_LONG)3, "956");
+    n57121 = VpAlloc((U_LONG)5, "-57121");
+
+    VpSetZero(y,1);        /* y   = 0 */
+    i1 = 0;
+    do {
+        ++i1;
+        /* VpDivd(f, r, t, n25); */ /* f = t/(-25) */
+        VpMult(f,t,n25);
+        VpAsgn(t, f, 1);       /* t = f    */
+
+        VpDivd(f, r, t, n); /* f = t/n  */
+
+        VpAddSub(r, y, f, 1);  /* r = y + f   */
+        VpAsgn(y, r, 1);       /* y = r    */
+
+        VpRdup(n);            /* n   = n + 1 */
+        VpRdup(n);            /* n   = n + 1 */
+        if(VpIsZero(f)) break;
+    } while((f->exponent > 0 ||    ((U_LONG)(-(f->exponent)) < y->MaxPrec)) &&
+            i1<nc
+    );
+
+    VpSetOne(n);
+    VpAsgn(t, n956,1);
+    i2 = 0;
+    do {
+        ++i2;
+        VpDivd(f, r, t, n57121); /* f = t/(-57121) */
+        VpAsgn(t, f, 1);      /* t = f    */
+
+        VpDivd(f, r, t, n);   /* f = t/n  */
+        VpAddSub(r, y, f, 1); /* r = y + f   */
+
+        VpAsgn(y, r, 1);      /* y = r    */
+        VpRdup(n);           /* n   = n + 1  */
+        VpRdup(n);           /* n   = n + 1 */
+        if(VpIsZero(f)) break;
+    } while((f->exponent > 0 || ((U_LONG)(-(f->exponent)) < y->MaxPrec)) &&
+            i2<nc
+        );
+
+    VpFree(n);
+    VpFree(n25);
+    VpFree(n956);
+    VpFree(n57121);
+
+    VpFree(t);
+    VpFree(f);
+    VpFree(r);
+#ifdef _DEBUG
+    printf("VpPai: # of iterations=%d+%d\n",i1,i2);
+#endif /* _DEBUG */
+}
+
+/*
+ * Calculates the value of e(=2.18281828459........).
+ * [Output] *y ... Real , the value of e.
+ *
+ *     y = e
+ */
+VP_EXPORT void
+VpExp1(Real *y)
+{
+    Real *n, *r, *f, *add;
+    U_LONG p;
+    U_LONG nc;
+    U_LONG i;
+
+    p = y->MaxPrec*(BASE_FIG + 2) + 2;
+    if(p<maxnr) nc = maxnr;
+    else  nc = p;
+
+    /* allocate temporally variables  */
+
+    r = VpAlloc(p *(BASE_FIG + 2), "#0");
+    f = VpAlloc(p, "#1");
+    n = VpAlloc(p, "#1");    /* n   = 1 */
+    add = VpAlloc(p, "#1");    /* add = 1 */
+
+    VpSetOne(y);            /* y   = 1 */
+    VpRdup(y);            /* y   = y + 1  */
+    i = 0;
+    do {
+        ++i;
+        VpRdup(n);        /* n   = n + 1  */
+        VpDivd(f, r, add, n);    /* f   = add/n(=1/n!)  */
+        VpAsgn(add, f, 1);    /* add = 1/n!  */
+        VpAddSub(r, y, f, 1);
+        VpAsgn(y, r, 1);    /* y = y + 1/n! */
+    } while((f->exponent > 0 || ((U_LONG)(-(f->exponent)) <= y->MaxPrec)) &&
+            i<nc
+        );
+
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "vpexp e=%\n", y);
+        printf("   r=%d\n", f[3]);
+    }
+#endif /* _DEBUG */
+    VpFree(add);
+    VpFree(n);
+    VpFree(f);
+    VpFree(r);
+}
+
+/*
+ * Calculates y=e**x where e(=2.18281828459........).
+ */
+VP_EXPORT void
+VpExp(Real *y, Real *x)
+{
+    Real *z=NULL, *div=NULL, *n=NULL, *r=NULL, *c=NULL;
+    U_LONG p;
+    U_LONG nc;
+    U_LONG i;
+
+    if(!VpIsDef(x)) {
+        VpSetNaN(y); /* Not sure */
+        goto Exit;
+    }
+    if(VpIsZero(x)) {
+        VpSetOne(y);
+        goto Exit;
+    }
+    p = y->MaxPrec;
+    if(p < x->Prec) p = x->Prec;
+    p = p *(BASE_FIG + 2) + 2;
+    if(p<maxnr) nc = maxnr;
+    else  nc = p;
+
+    /* allocate temporally variables  */
+    z = VpAlloc(p, "#1");
+    div = VpAlloc(p, "#1");
+
+    r = VpAlloc(p * 2, "#0");
+    c = VpAlloc(p, "#0");
+    n = VpAlloc(p, "#1");    /* n   = 1 */
+
+    VpSetOne(r);        /* y = 1 */
+    VpAddSub(y, r, x, 1);    /* y = 1 + x/1 */
+    VpAsgn(z, x, 1);    /* z = x/1  */
+
+    i = 0;
+    do {
+        ++i;
+        VpRdup(n);        /* n   = n + 1  */
+        VpDivd(div, r, x, n);    /* div = x/n */
+        VpMult(c, z, div);    /* c   = x/(n-1)! * x/n */
+        VpAsgn(z, c, 1);    /* z   = x*n/n! */
+        VpAsgn(r, y, 1);    /* Save previous val. */
+        VpAddSub(div, y, z, 1);    /*  */
+        VpAddSub(c, div, r, -1);    /* y = y(new) - y(prev) */
+        VpAsgn(y, div, 1);    /* y = y(new) */
+    } while(((!VpIsZero(c)) &&(c->exponent >= 0 ||((U_LONG)(-c->exponent)) <= y->MaxPrec)) &&
+            i<nc
+        );
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "vpexp e=%\n", y);
+    }
+#endif /* _DEBUG */
+    VpFree(div);
+    VpFree(n);
+    VpFree(c);
+    VpFree(r);
+    VpFree(z);
+}
+
+VP_EXPORT void
+VpSinCos(Real *psin,Real *pcos,Real *x)
+/*
+ * Calculates sin(x) & cos(x)
+ *(Assumes psin->MaxPrec==pcos->MaxPrec)
+ */
+{
+    Real *z=NULL, *div=NULL, *n=NULL, *r=NULL, *c=NULL;
+    U_LONG p;
+    int fcos;
+    int fsin;
+    int which;
+    U_LONG nc;
+    U_LONG i;
+
+    if(!VpIsDef(x)) {
+        VpSetNaN(psin);
+        VpSetNaN(pcos);
+        goto Exit;
+    }
+
+    p = pcos->MaxPrec;
+    if(p < x->Prec) p = x->Prec;
+    p = p *(BASE_FIG + 2) + 2;
+    if(p<maxnr) nc = maxnr;
+    else  nc = p;
+
+    /* allocate temporally variables  */
+    z = VpAlloc(p, "#1");
+    div = VpAlloc(p, "#1");
+
+    r = VpAlloc(p * 2, "#0");
+    c = VpAlloc(p , "#0");
+    n = VpAlloc(p, "#1");    /* n   = 1 */
+
+    VpSetOne(pcos);        /* cos = 1 */
+    VpAsgn(psin, x, 1);    /* sin = x/1 */
+    VpAsgn(z, x, 1);        /* z = x/1  */
+    fcos = 1;
+    fsin = 1;
+    which = 1;
+    i = 0;
+    do {
+        ++i;
+        VpRdup(n);        /* n   = n + 1  */
+        VpDivd(div, r, x, n);    /* div = x/n */
+        VpMult(c, z, div);    /* c   = x/(n-1)! * x/n */
+        VpAsgn(z, c, 1);    /* z   = x*n/n! */
+        if(which) {
+            /* COS */
+            which = 0;
+            fcos *= -1;
+            VpAsgn(r, pcos, 1);    /* Save previous val. */
+            VpAddSub(div, pcos, z, fcos);    /*  */
+            VpAddSub(c, div, r, -1);    /* cos = cos(new) - cos(prev) */
+            VpAsgn(pcos, div, 1);    /* cos = cos(new) */
+        } else {
+            /* SIN */
+            which = 1;
+            fsin *= -1;
+            VpAsgn(r, psin, 1);    /* Save previous val. */
+            VpAddSub(div, psin, z, fsin);    /*  */
+            VpAddSub(c, div, r, -1);    /* sin = sin(new) - sin(prev) */
+            VpAsgn(psin, div, 1);    /* sin = sin(new) */
+        }
+    } while(((!VpIsZero(c)) &&(c->exponent >= 0 || ((U_LONG)(-c->exponent)) <= pcos->MaxPrec)) &&
+            i<nc
+    );
+
+Exit:
+#ifdef _DEBUG
+    if(gfDebug) {
+        VPrint(stdout, "cos=%\n", pcos);
+        VPrint(stdout, "sin=%\n", psin);
+    }
+#endif /* _DEBUG */
+    VpFree(div);
+    VpFree(n);
+    VpFree(c);
+    VpFree(r);
+    VpFree(z);
+}
+
+#ifdef _DEBUG
+int
+VpVarCheck(Real * v)
+/*
+ * Checks the validity of the Real variable v.
+ * [Input]
+ *   v ... Real *, variable to be checked.
+ * [Returns]
+ *   0  ... correct v.
+ *   other ... error
+ */
+{
+    U_LONG i;
+
+    if(v->MaxPrec <= 0) {
+        printf("ERROR(VpVarCheck): Illegal Max. Precision(=%u)\n",
+            v->MaxPrec);
+        return 1;
+    }
+    if((v->Prec <= 0) ||((v->Prec) >(v->MaxPrec))) {
+        printf("ERROR(VpVarCheck): Illegal Precision(=%u)\n", v->Prec);
+        printf("       Max. Prec.=%u\n", v->MaxPrec);
+        return 2;
+    }
+    for(i = 0; i < v->Prec; ++i) {
+        if((v->frac[i] >= BASE)) {
+            printf("ERROR(VpVarCheck): Illegal fraction\n");
+            printf("       Frac[%d]=%u\n", i, v->frac[i]);
+            printf("       Prec.   =%u\n", v->Prec);
+            printf("       Exp. =%d\n", v->exponent);
+            printf("       BASE =%u\n", BASE);
+            return 3;
+        }
+    }
+    return 0;
+}
+#endif /* _DEBUG */
+
+static U_LONG
+SkipWhiteChar(char *szVal)
+{
+    char ch;
+    U_LONG i = 0;
+    while(ch = szVal[i++]) {
+        if(IsWhiteChar(ch)) continue;
+        break;
+    }
+    return i - 1;
+}
diff --git a/ext/bigdecimal/bigdecimal.def b/ext/bigdecimal/bigdecimal.def
new file mode 100644
index 0000000000..8450e164e6
--- /dev/null
+++ b/ext/bigdecimal/bigdecimal.def
@@ -0,0 +1,2 @@
+EXPORTS
+Init_bigdecimal
diff --git a/ext/bigdecimal/bigdecimal.h b/ext/bigdecimal/bigdecimal.h
new file mode 100644
index 0000000000..b7d2474e0b
--- /dev/null
+++ b/ext/bigdecimal/bigdecimal.h
@@ -0,0 +1,190 @@
+/*
+ *
+ * Ruby BigDecimal(Variable decimal precision) extension library. 
+ *
+ * Copyright(C) 2002 by Shigeo Kobayashi([email protected]) 
+ *
+ * You may distribute under the terms of either the GNU General Public 
+ * License or the Artistic License, as specified in the README file 
+ * of this BigDecimal distribution. 
+ *
+ *  2002-12-06
+ *    The sqrt() bug reported by Bret Jolly fixed.
+ *
+ */
+
+#ifndef  ____BIG_DECIMAL__H____
+#define  ____BIG_DECIMAL__H____
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/*
+ *   #define VP_EXPORT other than static to let VP_ routines 
+ *   be called from outside of this module.
+ */
+#define VP_EXPORT static 
+
+#define U_LONG unsigned long
+#define S_LONG long
+#define U_INT  unsigned int
+#define S_INT  int
+
+/* Exception codes */
+#define VP_EXCEPTION_ALL        ((unsigned short)0xFFFF)
+#define VP_EXCEPTION_INFINITY   ((unsigned short)0x0001)
+#define VP_EXCEPTION_NaN        ((unsigned short)0x0002)
+#define VP_EXCEPTION_UNDERFLOW  ((unsigned short)0x0004)
+#define VP_EXCEPTION_OVERFLOW   ((unsigned short)0x0001) /* 0x0008) */
+#define VP_EXCEPTION_ZERODIVIDE ((unsigned short)0x0001) /* 0x0010) */
+
+/* Following 2 exceptions cann't controlled by user */
+#define VP_EXCEPTION_OP         ((unsigned short)0x0020)
+#define VP_EXCEPTION_MEMORY     ((unsigned short)0x0040)
+
+#define VP_SIGN_NaN                0 /* NaN                      */
+#define VP_SIGN_POSITIVE_ZERO      1 /* Positive zero            */
+#define VP_SIGN_NEGATIVE_ZERO     -1 /* Negative zero            */
+#define VP_SIGN_POSITIVE_FINITE    2 /* Positive finite number   */
+#define VP_SIGN_NEGATIVE_FINITE   -2 /* Negative finite number   */
+#define VP_SIGN_POSITIVE_INFINITE  3 /* Positive infinite number */
+#define VP_SIGN_NEGATIVE_INFINITE -3 /* Negative infinite number */
+
+/*
+ * VP representation
+ *  r = 0.xxxxxxxxx *BASE**exponent
+ */
+typedef struct {
+    VALUE  obj;     /* Back pointer(VALUE) for Ruby object.     */
+    U_LONG MaxPrec; /* Maximum precision size                   */
+                    /* This is the actual size of pfrac[]       */
+                    /*(frac[0] to frac[MaxPrec] are available). */
+    U_LONG Prec;    /* Current precision size.                  */
+                    /* This indicates how much the.             */
+                    /* the array frac[] is actually used.       */
+    S_INT  exponent;/* Exponent part.                           */
+    short  sign;    /* Attributes of the value.                 */
+                    /*
+                     *        ==0 : NaN
+                     *          1 : Positive zero
+                     *         -1 : Negative zero
+                     *          2 : Positive number
+                     *         -2 : Negative number
+                     *          3 : Positive infinite number
+                     *         -3 : Negative infinite number
+                     */
+    short  flag;    /* Not used in vp_routines,space for user.  */
+    U_LONG frac[1]; /* Pointer to array of fraction part.       */
+} Real;
+
+/*  
+ *  ------------------
+ *   EXPORTables.
+ *  ------------------
+ */
+
+VP_EXPORT  Real *
+VpNewRbClass(U_LONG mx,char *str,VALUE klass);
+
+VP_EXPORT  Real *VpCreateRbObject(U_LONG mx,char *str);
+
+VP_EXPORT U_LONG VpBaseFig(void);
+VP_EXPORT U_LONG VpDblFig(void);
+VP_EXPORT U_LONG VpBaseVal(void);
+
+/* Zero,Inf,NaN (isinf(),isnan() used to check) */
+VP_EXPORT double VpGetDoubleNaN(void);
+VP_EXPORT double VpGetDoublePosInf(void);
+VP_EXPORT double VpGetDoubleNegInf(void);
+VP_EXPORT double VpGetDoubleNegZero(void);
+
+/* These 2 functions added at v1.1.7 */
+VP_EXPORT U_LONG VpGetPrecLimit(void);
+VP_EXPORT U_LONG VpSetPrecLimit(U_LONG n);
+
+VP_EXPORT int VpException(unsigned short f,char *str,int always);
+VP_EXPORT int VpIsNegDoubleZero(double v);
+VP_EXPORT U_LONG VpNumOfChars(Real *vp);
+VP_EXPORT U_LONG VpInit(U_LONG BaseVal);
+VP_EXPORT void *VpMemAlloc(U_LONG mb);
+VP_EXPORT void VpFree(Real *pv);
+VP_EXPORT Real *VpAlloc(U_LONG mx, char *szVal);
+VP_EXPORT int VpAsgn(Real *c,Real *a,int isw);
+VP_EXPORT int VpAddSub(Real *c,Real *a,Real *b,int operation);
+VP_EXPORT int VpMult(Real *c,Real *a,Real *b);
+VP_EXPORT int VpDivd(Real *c,Real *r,Real *a,Real *b);
+VP_EXPORT int VpComp(Real *a,Real *b);
+VP_EXPORT S_LONG VpExponent10(Real *a);
+VP_EXPORT void VpSzMantissa(Real *a,char *psz);
+VP_EXPORT void VpToString(Real *a,char *psz,int fFmt);
+VP_EXPORT int VpCtoV(Real *a,char *int_chr,U_LONG ni,char *frac,U_LONG nf,char *exp_chr,U_LONG ne);
+VP_EXPORT void VpVtoD(double *d,S_LONG *e,Real *m);
+VP_EXPORT void VpDtoV(Real *m,double d);
+VP_EXPORT void VpItoV(Real *m,S_INT ival);
+VP_EXPORT int VpSqrt(Real *y,Real *x);
+VP_EXPORT void VpRound(Real *y,Real *x,int sw,int f,int il);
+VP_EXPORT void VpFrac(Real *y,Real *x);
+VP_EXPORT int VpPower(Real *y,Real *x,S_INT n);
+VP_EXPORT void VpPai(Real *y);
+VP_EXPORT void VpExp1(Real *y);
+VP_EXPORT void VpExp(Real *y,Real *x);
+VP_EXPORT void VpSinCos(Real *psin,Real *pcos,Real *x);
+VP_EXPORT int VPrint(FILE *fp,char *cntl_chr,Real *a);
+
+/*  
+ *  ------------------
+ *  MACRO definitions.
+ *  ------------------
+ */
+#define Abs(a)     (((a)>= 0)?(a):(-(a)))
+#define Max(a, b)  (((a)>(b))?(a):(b))
+#define Min(a, b)  (((a)>(b))?(b):(a))
+
+#define IsWhiteChar(ch) (((ch==' ')||(ch=='\n')||(ch=='\t')||(ch=='\b'))?1:0)
+
+#define VpMaxPrec(a)   ((a)->MaxPrec)
+#define VpPrec(a)      ((a)->Prec)
+#define VpGetFlag(a)   ((a)->flag)
+
+/* Sign */
+
+/* VpGetSign(a) returns 1,-1 if a>0,a<0 respectively */
+#define VpGetSign(a) (((a)->sign>0)?1:(-1))
+/* Change sign of a to a>0,a<0 if s = 1,-1 respectively */
+#define VpChangeSign(a,s) {if((s)>0) (a)->sign=(short)Abs((S_LONG)(a)->sign);else (a)->sign=-(short)Abs((S_LONG)(a)->sign);}
+/* Sets sign of a to a>0,a<0 if s = 1,-1 respectively */
+#define VpSetSign(a,s)    {if((s)>0) (a)->sign=(short)VP_SIGN_POSITIVE_FINITE;else (a)->sign=(short)VP_SIGN_NEGATIVE_FINITE;}
+
+/* 1 */
+#define VpSetOne(a)       {(a)->frac[0]=(a)->exponent=(a)->Prec=1;(a)->sign=VP_SIGN_POSITIVE_FINITE;}
+
+/* ZEROs */
+#define VpIsPosZero(a)  ((a)->sign==VP_SIGN_POSITIVE_ZERO)
+#define VpIsNegZero(a)  ((a)->sign==VP_SIGN_NEGATIVE_ZERO)
+#define VpIsZero(a)     (VpIsPosZero(a) || VpIsNegZero(a))
+#define VpSetPosZero(a) ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_POSITIVE_ZERO)
+#define VpSetNegZero(a) ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_NEGATIVE_ZERO)
+#define VpSetZero(a,s)  ( ((s)>0)?VpSetPosZero(a):VpSetNegZero(a) )
+
+/* NaN */
+#define VpIsNaN(a)      ((a)->sign==VP_SIGN_NaN)
+#define VpSetNaN(a)     ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_NaN)
+
+/* Infinity */
+#define VpIsPosInf(a)   ((a)->sign==VP_SIGN_POSITIVE_INFINITE)
+#define VpIsNegInf(a)   ((a)->sign==VP_SIGN_NEGATIVE_INFINITE)
+#define VpIsInf(a)      (VpIsPosInf(a) || VpIsNegInf(a))
+#define VpIsDef(a)      ( !(VpIsNaN(a)||VpIsInf(a)) )
+#define VpSetPosInf(a)  ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_POSITIVE_INFINITE)
+#define VpSetNegInf(a)  ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_NEGATIVE_INFINITE)
+#define VpSetInf(a,s)   ( ((s)>0)?VpSetPosInf(a):VpSetNegInf(a) )
+
+#ifdef _DEBUG
+int VpVarCheck(Real * v);
+#endif /* _DEBUG */
+
+#if defined(__cplusplus)
+}  /* extern "C" { */
+#endif
+#endif //____BIG_DECIMAL__H____
diff --git a/ext/bigdecimal/bigdecimal_en.html b/ext/bigdecimal/bigdecimal_en.html
new file mode 100644
index 0000000000..c87cd5ccd8
--- /dev/null
+++ b/ext/bigdecimal/bigdecimal_en.html
@@ -0,0 +1,767 @@
+<HTML>
+<HEAD>
+<META HTTP-EQUIV="Content-Type" CONTENT="text/html">
+<style type="text/css"><!--
+body {
+  color: #3f0f0f;
+  background: #fefeff;
+  margin-left: 2em; margin-right: 2em;
+}
+h1 {
+  color: #ffffff;
+  background-color: #3939AD;
+  border-color: #FF00FF;
+  width: 100%;
+  border-style: solid;
+  border-top-width: 0.1em;
+  border-bottom-width: 0.1em;
+  border-right: none;
+  border-left: none;
+  padding: 0.1em;
+  font-weight: bold;
+  font-size: 160%;
+  text-align: center;
+}
+h2 {
+  color: #00007f;
+  background-color: #e7e7ff;
+  border-color: #000094;
+  width: 100%;
+  border-style: solid;
+  border-left: none;
+  border-right: none;
+  border-top-width: 0.1em;
+  border-bottom-width: 0.1em;
+  padding: 0.1em;
+  font-weight: bold;
+  font-size: 110%;
+}
+h3 {
+  color: #00007f;
+  padding: 0.2em;
+  font-size: 110%;
+}
+h4, h5 {
+  color: #000000;
+  padding: 0.2em;
+  font-size: 100%;
+}
+table {
+  margin-top: 0.2em; margin-bottom: 0.2em;
+  margin-left: 2em; margin-right: 2em;
+}
+caption {
+  color: #7f0000;
+  font-weight: bold;
+}
+th {
+  background: #e7e7ff;
+  padding-left: 0.2em; padding-right: 0.2em;
+}
+td {
+  background: #f3f7ff;
+  padding-left: 0.2em; padding-right: 0.2em;
+}
+code {
+  color: #0000df;
+}
+dt {
+  margin-top: 0.2em;
+}
+li {
+  margin-top: 0.2em;
+}
+pre
+{
+    BACKGROUND-COLOR: #d0d0d0;
+    BORDER-BOTTOM: medium none;
+    BORDER-LEFT: medium none;
+    BORDER-RIGHT: medium none;
+    BORDER-TOP: medium none;
+    LINE-HEIGHT: 100%;
+    MARGIN: 12px 12px 12px 12px;
+    PADDING-BOTTOM: 12px;
+    PADDING-LEFT: 12px;
+    PADDING-RIGHT: 12px;
+    PADDING-TOP: 12px;
+    WHITE-SPACE: pre;
+    WIDTH: 100%
+}
+--></style>
+
+<TITLE>BigDecimal:An extension library for Ruby</TITLE>
+</HEAD>
+<BODY BGCOLOR=#FFFFE0>
+<H1>BigDecimal(Variable Precision Floating Library for Ruby)</H1>
+<DIV align="right"><A HREF="./bigdecimal_ja.html">Japanese</A></DIV><BR>
+BigDecimal is an extension library for the Ruby interpreter. 
+Using BigDecimal class, you can obtain any number of significant digits in computation. 
+For the details about Ruby see:<BR>
+<UL>
+<LI><A HREF="http://www.ruby-lang.org/en/">http://www.ruby-lang.org/en/</A>:Official Ruby page(English).</LI>
+<LI><A HREF="http://ruby.freak.ne.jp/">http://ruby.freak.ne.jp/</A>:Ruby informations(Japanese).</LI>
+<LI><A HREF="http://kahori.com/ruby/ring/">http://kahori.com/ruby/ring/</A>:Mutually linked pages relating to Ruby(Japanese).
+</LI>
+</UL> 
+NOTE:<BR>
+ This software is provided "AS IS" and without any express or
+ implied warranties,including,without limitation,the implied
+ warranties of merchantibility and fitness for a particular
+ purpose. For the details,see COPYING and README included in this
+ distribution.
+<BR>
+<hr>
+
+<H2>Contents</H2>
+<UL>
+<LI><A HREF="#INTRO">Introduction</LI>
+<LI><A HREF="#SPEC">Usage and methods</A></LI>
+<LI><A HREF="#UNDEF">Infinity,NaN,Zero</A></LI>
+<LI><A HREF="#STRUCT">Internal structure</A></LI>
+<LI><A HREF="#BASE">Binary or decimal number representation</A></LI>
+<LI><A HREF="#PREC">Resulting number of significant digits</A></LI>
+</UL>
+<HR>
+
+<A NAME="#INTRO">
+<H2>Introduction</H2>
+Ruby already has builtin (variable length integer number) class Bignum. Using Bignum class,you can obtain
+ any integer value in magnitude. But, variable length floating number class is not yet built in. 
+This is why I made variable length floating class BigDecimal.
+Feel free to send any comments or bug reports to me.
+<A HREF="mailto:[email protected]">[email protected]</A>
+I will try(but can't promise) to fix bugs reported. 
+<hr>
+<H2>Installation</H2>
+The Ruby latest version can be downloaded from <A HREF="http://www.ruby-lang.org/en/">Official Ruby page</A>.
+Once decompress the downloaded Ruby archive,follow the normal installation procedures according to the 
+documents included.
+
+<A NAME="#SPEC">
+<H2>Usage and methods</H2>
+Suppose you already know Ruby programming,
+to create BigDecimal objects,the program would like:<BR>
+
+<CODE><PRE>
+   require 'bigdecimal'
+   a=BigDecimal::new("0.123456789123456789")
+   b=BigDecimal::new("123456.78912345678",40)
+   c=a+b
+</PRE></CODE>
+
+<H3>List of methods</H3>
+In the following explanations,n specifies the minimum number of resulting significant digits,
+not exactly but slightly excess memories will be allocated to newly created object.
+In 32 bits integer system,every 4 digits(in decimal) are computed simultaneously.
+This means the number of significant digits in BigDecimal is always a multiple of 4.
+
+<UL>
+<LI><B>new</B></LI><BR>
+"new" method creates a new BigDecimal object.<BR>
+a=BigDecimal::new(s[,n])<BR>
+where:<BR>
+s: Initial value string.<BR>
+n: Maximum number of significant digits of a. n must be a Fixnum object.
+If n is omitted or is equal to 0,then the maximum number of significant digits of a is determined from the length of s.
+
+<LI><B>double_fig</B></LI><BR>
+double_fig is a class method which returns the number of digits 
+the Float class can have.
+<CODE><PRE>
+  p BigDecimal::double_fig  # ==> 20 (depends on the CPU etc.)
+</PRE></CODE>
+The equivalent C programs which calculates the value of
+double_fig is:
+<CODE><PRE>
+ double v          = 1.0;
+ int    double_fig = 0;
+ while(v + 1.0 > 1.0) {
+    ++double_fig;
+    v /= 10;
+ }
+</PRE></CODE>
+
+<LI><B>prec</B></LI><BR>
+r,m = a.prec<BR>
+where r is the number of significant digits of a,
+m is the maximum number of significant digits a can hold.<br>
+<CODE><PRE>
+  require "bigdecimal"
+  a = BigDecimal.new("0.12345")
+  p a.prec                              # ==> [8, 12]
+  b = BigDecimal.new("0.1234500000000")
+  p b.prec                              # ==> [8, 20]
+  c = BigDecimal.new("0.12345",20)
+  p c.prec                              # ==> [8, 24]
+</PRE></CODE>
+r and m are always the multiple of log10(BigDecimal::BASE).
+
+<LI><B>+</B></LI><BR>
+addition(c = a + b)<BR>
+For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
+
+<LI><B>-</B></LI><BR>
+subtraction (c = a - b) or negation (c = -a)<BR>
+For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
+
+<LI><B>*</B></LI><BR>
+multiplication(c = a * b)<BR>
+For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
+
+<LI><B>/</B></LI><BR>
+division(c = a / b)<BR>
+For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
+
+<LI><B>assign</B></LI><BR>
+c = a.assign(n,f)<BR>
+assigns the value of a to c.<BR>
+n is the number of significant digits of resulting c.<BR>
+If f > 0,then a is assigned to c.<BR>
+If f < 0,then -a is assigned to c.<BR>
+The absolute value of f (|f|) must be 1 or 2.
+If |f|=2,then proper round operation over c is performed,when the maximum 
+number of significant digits of c is less than current 
+number of significant digits of a.
+If |f|=1 then extra digits are discarded when the maximum 
+number of significant digits of c is less than current 
+number of significant digits of a.
+
+<LI><B>add</B></LI><BR>
+c = a.add(b,n)<BR>
+c = a.add(b,n) performs c = a + b.
+If n is less than the actual significant digits of a + b,
+then c is rounded properly.
+
+<LI><B>sub</B></LI><BR>
+c = a.sub(b,n)<BR>
+c = a.sub(b,n) performs c = a - b.
+If n is less than the actual significant digits of a - b,
+then c is rounded properly.
+
+<LI><B>mult</B></LI><BR>
+c = a.mult(b,n)<BR>
+c = a.mult(b,n) performs c = a * b.
+If n is less than the actual significant digits of a * b,
+then c is rounded properly.
+
+<LI><B>div</B></LI><BR>
+c,r = a.div(b,n)<BR>
+c,r = a.div(b,n) performs c = a / b, r is the residue of a / b.
+If necessary,the divide operation continues to n digits which c
+ can hold.
+Unlike the divmod method,c is not always an integer.
+c is never rounded,and the equation a = c*b + r is always 
+valid unless c is NaN or Infinity.
+
+<LI><B>%</B></LI><BR>
+r = a%b <BR>
+is the same as:<BR>
+r = a-((a/b).floor)*b<BR>
+
+<LI><B>fix</B></LI><BR>
+c = a.fix<BR>
+returns integer part of a.<BR>
+
+<LI><B>frac</B></LI><BR>
+c = a.frac<BR>
+returns fraction part of a.<BR>
+
+<LI><B>floor[(n)]</B></LI><BR>
+c = a.floor<BR>
+returns the maximum integer value (in BigDecimal) which is less than or equal to a.<BR>
+As shown in the following example,an optional integer argument (n) specifying the position 
+of 'floor'ed digit can be given.
+If n> 0,then the (n+1)th digit counted from the decimal point in fraction part is 'floor'ed.
+If n<0,then the n-th digit counted from the decimal point in integer part is 'floor'ed.<BR>
+
+c = BigDecimal::new("1.23456")<BR>
+d = c.floor(4) # d = 1.2345<BR>
+c = BigDecimal::new("15.23456")<BR>
+d = c.floor(-1) # d = 10.0<BR>
+
+<LI><B>ceil[(n)]</B></LI><BR>
+c = a.ceil<BR>
+returns the minimum integer value (in BigDecimal) which is greater than or equal to a.<BR>
+As shown in the following example,an optional integer argument (n) specifying the position 
+of 'ceil'ed digit can be given.
+If n>0,then the (n+1)th digit counted from the decimal point in fraction part is 'ceil'ed.
+If n<0,then the n-th digit counted from the decimal point in integer part is 'ceil'ed.<BR>
+
+c = BigDecimal::new("1.23456")<BR>
+d = c.ceil(4) # d = 1.2346<BR>
+c = BigDecimal::new("15.23456")<BR>
+d = c.ceil(-1) # d = 20.0<BR>
+
+<LI><B>round[(n)]</B></LI><BR>
+c = a.round<BR>
+round off a to the nearest 1�D<BR>
+As shown in the following example,an optional integer argument (n) specifying the position 
+of rounded digit can be given.
+If n>0,then the (n+1)th digit counted from the decimal point in fraction part is rounded.
+If n<0,then the n-th digit counted from the decimal point in integer part is rounded.<BR>
+
+c = BigDecimal::new("1.23456")<BR>
+d = c.round(4)  # d = 1.235 <BR>
+c = BigDecimal::new("15.23456")<BR>
+d = c.round(-1) # d = 20.0<BR>
+
+<LI><B>truncate[(n)]</B></LI><BR>
+c = a.truncate<BR>
+truncate a to the nearest 1�D<BR>
+As shown in the following example,an optional integer argument (n) specifying the position 
+of truncated digit can be given.
+If n>0,then the (n+1)th digit counted from the decimal point in fraction part is truncated.
+If n<0,then the n-th digit counted from the decimal point in integer part is truncated.<BR>
+
+c = BigDecimal::new("1.23456")<BR>
+d = c.truncate(4)  # d = 1.2345<BR>
+c = BigDecimal::new("15.23456")<BR>
+d = c.truncate(-1) # d = 10.0<BR>
+
+<LI><B>divmod</B></LI><BR>
+c,r = a.divmod(b) # a = c*b + r<BR>
+returns the quotient and remainder of a/b.<BR>
+a = c * b + r is always satisfied.<BR>
+where c is the integer sutisfying
+c = (a/b).floor <BR>
+and,therefore 
+r = a - c*b<BR>
+
+<LI><B>remainder</B></LI><BR>
+r=a.remainder(b)<BR>
+returns the remainder of a/b.<BR>
+where c is the integer sutisfying 
+c = (a/b).fix <BR>
+and,therefore: 
+r = a - c*b<BR>
+
+<LI><B>abs</B></LI><BR>
+c = a.abs<BR>
+returns an absolute value of a.<BR>
+
+<LI><B>to_i</B></LI><BR>
+changes a to an integer.<BR>
+i = a.to_i<BR>
+i becomes to Fixnum or Bignum.
+IF a is Infinity or NaN,then i becomes to nil.
+
+<LI><B>to_s[(n)]</B></LI><BR>
+converts to string(results look like "0.xxxxxEn").<BR>
+s = a.to_s<BR>
+If n is given,then a space is inserted after every n digits for readability.<BR>
+s = a.to_s(n)
+
+<LI><B>exponent</B></LI><BR>
+returns an integer holding exponent value of a.<BR>
+n = a.exponent <BR>
+means a = 0.xxxxxxx*10**n.
+
+<LI><B>to_f</B></LI><BR>
+same as dup method.
+creates a new BigDecimal object having same value.
+
+<LI><B>E</B></LI><BR>
+e = BigDecimal::E(n)<BR>
+where e(=2.718281828....) is the base value of natural logarithm.<BR>
+n specifies the length of significant digits of e.
+
+<LI><B>PI</B></LI><BR>
+e = BigDecimal::PI(n)<BR>
+returns at least n digits of the ratio of the circumference of a circle to its dirmeter
+(pi=3.14159265358979....) using J.Machin's formula.<BR>
+
+<LI><B>BASE</B></LI><BR>
+Base value used in the BigDecimal calculation.
+On 32 bit integer system,the value of BASE is 10000.<BR>
+b = BigDecimal::BASE<BR>
+
+<LI><B>mode</B></LI><BR>
+mode method controls BigDecimal computation.
+Following usage are defined.<BR>
+
+f = BigDecimal::mode(BigDecimal::EXCEPTION_NaN,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_INFINITY,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_UNDERFLOW,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_OVERFLOW,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_ZERODIVIDE,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_ALL,flag)<BR>
+
+EXCEPTION_NaN controls the execution once computation results to NaN.
+EXCEPTION_INFINITY controls the execution once computation results to Infinity(�}Infinity).
+EXCEPTION_UNDERFLOW controls the execution once computation underflows.
+EXCEPTION_OVERFLOW controls the execution once computation overflows.
+EXCEPTION_ZERODIVIDE controls the execution once zero-division occures.
+EXCEPTION_ALL controls the execution for any exception defined occures.
+If the flag is true,then the relating exception is thrown.
+No exception is thrown when the flag is false(default) and computation 
+continues with the result:<BR>
+EXCEPTION_NaN results to NaN<BR>
+EXCEPTION_INFINITY results to +Infinity or -Infinity<BR>
+EXCEPTION_UNDERFLOW results to 0.<BR>
+EXCEPTION_OVERFLOW results to +Infinity or -Infinity<BR>
+EXCEPTION_ZERODIVIDE results to +Infinity or -Infinity<BR>
+EXCEPTION_INFINITY,EXCEPTION_OVERFLOW, and EXCEPTION_ZERODIVIDE are
+ currently the same.<BR>
+The return value of mode method is the value set.
+Suppose the return value of the mode method is f,then 
+ f & BigDecimal::EXCEPTION_NaN !=0 means EXCEPTION_NaN is set to on.
+If the value of the argument flag is other than nil,true nor false then 
+current mode status is returned.
+
+<LI><B>limit[(n)]</B></LI><BR>
+Limits the maximum digits that the newly created BigDecimal objects can hold
+never exceed n. Returns maximum value before set.
+Zero,the default value,means no upper limit.<BR>
+mf = BigDecimal::limit(n)<BR>
+
+<LI><B>sign</B></LI><BR>
+returns the 'attribute'.
+n = a.sign <BR>
+where the value of n means that a is:<BR>
+n = BigDecimal::SIGN_NaN(0) : a is NaN<BR>
+n = BigDecimal::SIGN_POSITIVE_ZERO(1) : a is +0<BR>
+n = BigDecimal::SIGN_NEGATIVE_ZERO(-1) : a is -0<BR>
+n = BigDecimal::SIGN_POSITIVE_FINITE(2) : a is positive<BR>
+n = BigDecimal::SIGN_NEGATIVE_FINITE(-2) : a is negative<BR>
+n = BigDecimal::SIGN_POSITIVE_INFINITE(3) : a is +Infinity<BR>
+n = BigDecimal::SIGN_NEGATIVE_INFINITE(-3) : a is -Infinity<BR>
+The value in () is the actual value,see (<A HREF="#STRUCT">Internal structure</A>.<BR>
+
+<LI><B>nan?</B></LI><BR>
+a.nan? returns True when a is NaN.
+
+<LI><B>infinite?</B></LI><BR>
+a.infinite? returns True when a is +�� or -��.
+
+<LI><B>finite?</B></LI><BR>
+a.finite? returns True when a is neither �� nor NaN.
+
+<LI><B>to_parts</B></LI><BR>
+decomposes a BigDecimal value to 4 parts.
+All 4 parts are returned as an array.<BR>
+Parts consist of a sign(0 when the value is NaN,+1 for positive and
+ -1 for negative value), a string representing fraction part,base value(always 10 currently),and an integer(Fixnum) for exponent respectively.
+a=BigDecimal::new("3.14159265",10)<BR>
+f,x,y,z = a.to_parts<BR>
+where f=+1,x="314159265",y=10 and z=1<BR>
+therefore,you can translate BigDecimal value to Float as:<BR>
+s = "0."+x<BR>
+b = f*(s.to_f)*(y**z)<BR>
+
+<LI><B>inspect</B></LI><BR>
+is used for debugging output.<BR>
+p a=BigDecimal::new("3.14",10)<BR>
+should produce output like "#&lt;0x112344:'0.314E1',4(12)%gt;".
+where "0x112344" is the address,
+'0.314E1' is the value,4 is the number of the significant digits,
+and 12 is the maximum number of the significant digits 
+the object can hold.
+
+<LI><B>dup</B></LI><BR>
+creates a new BigDecimal object having same value.
+
+<LI><B>sqrt</B></LI><BR>
+c = a.sqrt(n)<BR>
+computes square root value of a with significant digit number n at least.<BR>
+
+<LI><B>sincos</B></LI><BR>
+computes and returns sine and cosine value of a with significant digit number n at least.<BR>
+sin,cos = a.sincos(n)<BR>
+
+<LI><B>exp</B></LI><BR>
+c = a.exp(n)<BR>
+computes the base of natural logarithm value(e=2.718281828....) powered by a
+with significant digit number n at least.<BR>
+
+<LI><B>power</B></LI><BR>
+c = a.power(n)<BR>
+returns the value of a powered by n(c=a**n).
+n must be an integer.<BR>
+
+<LI><B>zero?</B></LI><BR>
+c = a.zero?<BR>
+returns true if a is equal to 0,otherwise returns false<BR>
+
+<LI><B>nonzero?</B></LI><BR>
+c = a.nonzero?<BR>
+returns false if a is 0,otherwise returns a itself.<BR>
+
+<LI><B>&lt;=&gt;</B></LI><BR>
+c = a &lt;=&gt; b <BR>
+returns 0 if a==b,1 if a &gt b,and returns -1 if a &lt b.<BR>
+</UL>
+
+Following methods need no explanation.<BR>
+<UL>
+<LI>==</LI>
+<LI>===</LI>
+same as ==,used in case statement.
+<LI>!=</LI>
+<LI>&lt;</LI>
+<LI>&lt;=</LI>
+<LI>&gt;</LI>
+<LI>&gt;=</LI>
+</UL>
+<HR>
+
+<H3>About 'coerce'</H3>
+<B>For the binary operation like A op B:</B>
+<DL>
+<DT> 1.Both A and B are BigDecimal objects</DT>
+<DD> A op B is normally performed.</DD>
+<DT> 2.A is the BigDecimal object but B is other than BigDecimal object</DT>
+<DD> Operation is performed,after B is translated to correcponding BigDecimal object(because BigDecimal supports coerce method).</DD>
+<DT> 3.A is not the BigDecimal object but B is BigDecimal object</DT>
+<DD>If A has coerce mthod,then B will translate A to corresponding 
+BigDecimal object and the operation is performed,otherwise an error occures.</DD>
+</DL>
+
+Attention must be paid when a String is to be translated to BigDecimal.
+Translation stops without error at the character representing non digit.
+For instance,"10XX" is translated to 10,"XXXX" is translated to 0.<BR>
+String representing zero or infinity such as "Infinity","+Infinity","-Infinity",and "NaN" can also be translated to BigDecimal unless false is specified by mode method.<BR>
+
+BigDecimal class supports coerce method(for the details about coerce method,see Ruby documentations). This means the most binary operation can be performed if the BigDecimal object is at the left hand side of the operation.<BR><BR>
+
+ For example:
+<CODE><PRE>
+  a = BigDecimal.E(20)
+  c = a * "0.123456789123456789123456789" # A String is changed to BigDecimal object.
+</PRE></CODE>
+is performed normally.<BR>
+ But,because String does not have coerce method,the following example can not be performed.<BR>
+
+<CODE><PRE>
+  a = BigDecimal.E(20)
+  c = "0.123456789123456789123456789" * a # ERROR
+</PRE></CODE>
+
+If you actually have any inconvenience about the error above.
+You can define a new class derived from String class,
+and define coerce method within the new class.<BR>
+
+<hr>
+<A NAME="#UNDEF">
+<H2>Infinity,Not a Number(NaN),Zero</H2>
+Infinite numbers and NaN can be represented by string writing "+Infinity"(or "Infinity"),"-Infinity",and "NaN" respectively in your program.
+Infinite numbers can be obtained by 1.0/0.0(=Infinity) or -1.0/0.0(=-Infinity).
+<BR><BR>
+NaN(Not a number) can be obtained by undefined computation like 0.0/0.0 
+or Infinity-Infinity.
+Any computation including NaN results to NaN.
+Comparisons with NaN never become true,including comparison with NaN itself.
+<BR><BR>
+Zero has two different variations as +0.0 and -0.0.
+But,still, +0.0==-0.0 is true.
+<BR><BR>
+Computation results including Infinity,NaN,+0.0 or -0.0 become complicated.
+Run following program and comfirm the results.
+Send me any incorrect result if you find.
+
+<PRE><CODE>
+ require "bigdecimal"
+ aa  = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
+ ba  = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
+ opa = %w(+ - * / <=> > >=  < == != <=)
+ for a in aa
+  for b in ba
+    for op in opa
+      x = BigDecimal::new(a)
+      y = BigDecimal::new(b)
+      eval("ans= x #{op} y;print a,' ',op,' ',b,' ==> ',ans.to_s,\"\n\"")
+    end
+  end
+ end
+</CODE></PRE>
+<hr>
+
+<A NAME="#STRUCT">
+<H2>Internal structure</H2>
+BigDecimal number is defined by the structure Real in BigDecimal.h.
+Digits representing a float number are kept in the array frac[] defined in the structure.
+In the program,any floating number(BigDecimal number) is represented as:<BR>
+ <BigDecimal number> = 0.xxxxxxxxx*BASE**n<BR><BR>
+where 'x' is any digit representing mantissa(kept in the array frac[]),
+BASE is base value(=10000 in 32 bit integer system),
+and n is the exponent value.<BR>
+Larger BASE value enables smaller size of the array frac[],and increases computation speed.
+The value of BASE is defined ind VpInit(). In 32 bit integer system,this value is 
+10000. In 64 bit integer system,the value becomes larger.
+BigDecimal has not yet been compiled and tested on 64 bit integer system.
+It will be very nice if anyone try to run BigDecimal on 64 bit system and
+ inform me the results.
+When BASE is 10000,an element of the array frac[] can have vale of from 0 to 9999.
+(up to 4 digits).<BR>
+The structure Real is defined in bigdecimal.h as:<BR>
+<CODE><PRE>
+  typedef struct {
+     VALUE  obj;     /* Back pointer(VALUE) for Ruby object.         */
+     unsigned long MaxPrec; /* The size of the array frac[]          */
+     unsigned long Prec;    /* Current size of frac[] actually used. */
+     short    sign;         /* Attribute of the value.  */
+                            /*  ==0 : NaN               */
+                            /*    1 : +0                */
+                            /*   -1 : -0                */
+                            /*    2 : Positive number   */
+                            /*   -2 : Negative number   */
+                            /*    3 : +Infinity         */
+                            /*   -3 : -Infinity         */
+     unsigned short flag;   /* Control flag             */
+     int      exponent;     /* Exponent value(0.xxxx*BASE**exponent) */
+     unsigned long frac[1]; /* An araay holding mantissa(Variable)   */
+  } Real;
+</CODE></PRE>
+The decimal value 1234.56784321 is represented as(BASE=10000):<BR>
+<PRE>
+    0.1234 5678 4321*(10000)**1
+</PRE>
+where frac[0]=1234,frac[1]=5678,frac[2]=4321,
+Prec=3,sign=2,exponent=1. MaxPrec can be any value greater than or equal to 
+Prec.
+<hr>
+
+<A NAME="#BASE">
+<H2>Binary or decimal number representation</H2>
+I adopted decimal number representation for BigDecimal implementation.
+Of cource,binary number representation is common on the most computers.
+
+<H3>Advantages using decimal representation</H3>
+The reason why I adopted decimal number representation for BigDecimal is:<BR>
+<DL>
+<DT>Easy for debugging
+<DD>The floating number 1234.56784321 can be easily represented as:<BR>
+  frac[0]=1234,frac[1]=5678,frac[2]=4321,exponent=1,and sign=2.
+<DT>Exact representation
+<DD>Following program can add all numbers(in decimal) in a file
+ without any error(no round operation).<BR>
+
+<PRE><CODE>
+   file = File::open(....,"r")
+   s = BigDecimal::new("0")
+   while line = file.gets
+      s = s + line
+   end
+</CODE></PRE>
+
+If the internal representation is binary,translation from decimal to 
+binary is required and the translation error is inevitable.
+For example, 0.1 can not exactly be represented in binary.<BR>
+0.1 => b1*2**(-1)+b1*2**(-2)+b3*2**(-3)+b4*2**(-4)....<BR>
+where b1=0,b2=0,b3=0,b4=1...<BR>
+bn(n=1,2,3,...) is infinite series of digit with value of 0 or 1,
+and rounding operation is necessary but where we should round the series ?
+Of cource,exact "0.1" is printed if the rouding operation is properly done,
+<DT>Significant digit we can have is automatically determined
+<DD>In binary representation,0.1 can not be represented in finite series of digit.
+
+But we only need one element(frac[0]=1) in decimal representation.
+This means that we can always determine the size of the array frac[] in Real 
+structure.
+</DL>
+
+<H3>Disadvantage of decimal representation</H3>
+Advantages stated so far can also be disadvantages if the input from outside is
+ represented in binary.
+Translation error from decimal to binary or vice versa is inevitable.
+So,translation from Float(binary) to BigDecimal(decimal) is not alway done exactly.
+
+<H4>Which is the first input?</H4>
+Because most people uses decimal notatin for numeric data representation,
+BigDecimal can handle numeric data without loss of translation error.
+<hr>
+
+<A NAME="#PREC">
+<H2>Resulting number of significant digits</H2>
+For the fundamental arithmetics such as addition,subtraction,
+multiplication,and division,I prepared 2 group of methods<BR>
+
+<H3>1. +,-,*,/</H3>
+For the operation + - * /,you can not specify the resulting 
+number of significant digits.<BR>
+Resulting number of significant digits are defined as:<BR>
+1.1 For * and /,resulting number of significant digits is the sum of the 
+significant digits of both side of the operator.<BR>
+1.2 For + and -,resulting number of significant digits is determined so that
+ no round operation is needed. <br>
+For example, c has more than 100 siginificant digits if c is computed as:<BR>
+c = 0.1+0.1*10**(-100)<br>
+<BR>
+As +,-,and * are always exact(no round operation is performed),
+which means more momories are required to keep computation results.
+
+As for the division as c = a/b,the significant digits of c is the same 
+as a*b. Division such as c=1.0/3.0 will be rounded.<BR>
+
+<H3>2. assign,add,sub,mult,div</H3>
+The length of the significant digits obtained from +,-,*,/ 
+is always defined by that of right and left side of the operator.
+To specify the length of the significant digits by your self,
+use methos assign,add,sub,mult,div, or limit(class method).
+Following example compute the ratio of the circumference of a circle to 
+its dirmeter(pi=3.14159265358979....) using J.Machin's formula.
+<BR><BR>
+<CODE><PRE>
+#!/usr/local/bin/ruby
+
+#
+# pai.rb
+#  USAGE: ruby pai.rb n
+#   where n is the number of digits required.
+#  EX.: ruby pai.rb 1000
+#
+
+require "bigdecimal"
+#
+# Calculates 3.1415.... using J. Machin's formula.
+#
+def pai(sig) # sig: Number of significant figures
+  exp    = -sig
+  pi     = BigDecimal::new("0")
+  two    = BigDecimal::new("2")
+  m25    = BigDecimal::new("-0.04")
+  m57121 = BigDecimal::new("-57121")
+
+  u = BigDecimal::new("1")
+  k = BigDecimal::new("1")
+  w = BigDecimal::new("1")
+  t = BigDecimal::new("-80")
+  while (u.exponent >= exp) 
+    t   = t*m25
+    u,r = t.div(k,sig)
+    pi  = pi + u
+    k   = k+two
+  end
+
+  u = BigDecimal::new("1")
+  k = BigDecimal::new("1")
+  w = BigDecimal::new("1")
+  t = BigDecimal::new("956")
+  while (u.exponent >= exp )
+    t,r = t.div(m57121,sig)
+    u,r = t.div(k,sig)
+    pi  = pi + u
+    k   = k+two
+  end
+  pi
+end
+
+if $0 == __FILE__
+  print "PAI("+ARGV[0]+"):\n"
+  p pai(ARGV[0].to_i)
+end
+
+</PRE></CODE>
+<HR>
+<FONT size=2>
+<I>
+<A HREF="http://www.tinyforest.gr.jp">
+Shigeo Kobayashi
+</A>
+(E-Mail:<A HREF="mailto:[email protected]">&lt;[email protected]&gt;</U></A>)
+</I>
+</FONT>
+</TD>
+</TR>
+</TABLE>
+</BODY>
+</HTML>
diff --git a/ext/bigdecimal/bigdecimal_ja.html b/ext/bigdecimal/bigdecimal_ja.html
new file mode 100644
index 0000000000..bfddabd3ba
--- /dev/null
+++ b/ext/bigdecimal/bigdecimal_ja.html
@@ -0,0 +1,706 @@
+<HTML>
+<HEAD>
+<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=Shift_JIS">
+<style type="text/css"><!--
+body {
+  color: #3f0f0f;
+  background: #fefeff;
+  margin-left: 2em; margin-right: 2em;
+}
+h1 {
+  color: #ffffff;
+  background-color: #3939AD;
+  border-color: #FF00FF;
+  width: 100%;
+  border-style: solid;
+  border-top-width: 0.1em;
+  border-bottom-width: 0.1em;
+  border-right: none;
+  border-left: none;
+  padding: 0.1em;
+  font-weight: bold;
+  font-size: 160%;
+  text-align: center;
+}
+h2 {
+  color: #00007f;
+  background-color: #e7e7ff;
+  border-color: #000094;
+  width: 100%;
+  border-style: solid;
+  border-left: none;
+  border-right: none;
+  border-top-width: 0.1em;
+  border-bottom-width: 0.1em;
+  padding: 0.1em;
+  font-weight: bold;
+  font-size: 110%;
+}
+h3 {
+  color: #00007f;
+  padding: 0.2em;
+  font-size: 110%;
+}
+h4, h5 {
+  color: #000000;
+  padding: 0.2em;
+  font-size: 100%;
+}
+table {
+  margin-top: 0.2em; margin-bottom: 0.2em;
+  margin-left: 2em; margin-right: 2em;
+}
+caption {
+  color: #7f0000;
+  font-weight: bold;
+}
+th {
+  background: #e7e7ff;
+  padding-left: 0.2em; padding-right: 0.2em;
+}
+td {
+  background: #f3f7ff;
+  padding-left: 0.2em; padding-right: 0.2em;
+}
+code {
+  color: #0000df;
+}
+dt {
+  margin-top: 0.2em;
+}
+li {
+  margin-top: 0.2em;
+}
+pre
+{
+    BACKGROUND-COLOR: #d0d0d0;
+    BORDER-BOTTOM: medium none;
+    BORDER-LEFT: medium none;
+    BORDER-RIGHT: medium none;
+    BORDER-TOP: medium none;
+    LINE-HEIGHT: 100%;
+    MARGIN: 12px 12px 12px 12px;
+    PADDING-BOTTOM: 12px;
+    PADDING-LEFT: 12px;
+    PADDING-RIGHT: 12px;
+    PADDING-TOP: 12px;
+    WHITE-SPACE: pre;
+    WIDTH: 100%
+}
+--></style>
+
+<TITLE>BigDecimal:An extension library for Ruby</TITLE>
+</HEAD>
+<BODY BGCOLOR=#FFFFE0>
+<H1>BigDecimal(�•ϒ����������_���Z�p�g�����C�u����)</H1>
+<DIV align="right"><A HREF="./bigdecimal_en.html">English</A></DIV><BR>
+BigDecimal �̓I�u�W�F�N�g�w���̋��͂ȃX�N���v�g����ł��� Ruby �ɉ•ϒ����������_
+�v�Z�@�\��lj����邽�߂̊g�����C�u�����ł��B
+Ruby �ɂ‚��Ă̏ڂ������e�͈ȉ���URL���Q�Ƃ��Ă��������B
+<UL>
+<LI><A HREF="http://www.ruby-lang.org/ja/">http://www.ruby-lang.org/ja/</A>�FRuby�����y�[�W</LI>
+<LI><A HREF="http://ruby.freak.ne.jp/">http://ruby.freak.ne.jp/</A>�FRuby�Ɋւ�����y�[�W</LI>
+<LI><A HREF="http://kahori.com/ruby/ring/">http://kahori.com/ruby/ring/</A>�FRuby�Ɋւ���y�[�W��H��܂�</LI>
+</UL> 
+<hr>
+<H2>�ڎ�</H2>
+<UL>
+<LI><A HREF="#INTRO">�͂��߂�</LI>
+<LI><A HREF="#SPEC">�g�p���@�ƃ��\�b�h�̈ꗗ</A></LI>
+<LI><A HREF="#UNDEF">�����A�񐔁A�[���̈���</A></LI>
+<LI><A HREF="#STRUCT">�����\��</A></LI>
+<LI><A HREF="#BASE">2�i��10�i</A></LI>
+<LI><A HREF="#PREC">�v�Z���x�ɂ‚���</A></LI>
+</UL>
+
+<HR>
+<A NAME="#INTRO">
+<H2>�͂��߂�</H2>
+Ruby �ɂ� Bignum �Ƃ����N���X������A���S���̐����ł��v�Z���邱�Ƃ��ł��܂��B
+�����A�C�ӌ��̕��������_���Z�p�N���X�������悤�ł��B�����ŁA
+�C�ӌ��̕��������_���Z�p�g�����C�u���� BigDecimal ���쐬���܂����B
+�s��⏕���E��Ă�����ꍇ�ǂ��ǂ��A
+<A HREF="mailto:[email protected]">[email protected]</A>
+�܂ł��m�点���������B�s��𒼂��C�͑傢�ɂ���܂��B�����A���ԂȂǂ̊֌W�Ŗ�
+�͂ł��܂���B�܂��A���ʂɂ‚��Ă��ۏ؂ł�����̂ł͂���܂���B
+�\�߁A���������������B
+<BR><BR>
+���̃v���O�����́A���R�ɔz�z�E���ς��č\���܂���B�������A���쌠�͕������Ă��܂���B
+�z�z�E���ϓ��̌����� Ruby �̂���ɏ����܂��B�ڂ����� README ��ǂ�ł��������B
+
+<hr>
+<H2>�C���X�g�[���ɂ‚���</H2>
+BigDecimal ���܂� Ruby �̍ŐV�ł�<A HREF="http://www.ruby-lang.org/ja/">Ruby�����y�[�W</A>����_�E�����[�h�ł��܂��B
+�_�E�����[�h�����ŐV�ł��𓀂�����A�ʏ�̃C���X�g�[���菇�����s���ĉ������B
+Ruby ���������C���X�g�[�������΁A������ BigDecimal �����p�ł���悤�ɂȂ�͂��ł��B
+�\�[�X�t�@�C���� 
+bigdecimal.c,bigdecimal.h 
+�̂Q�‚݂̂ł��B<BR>
+
+<hr>
+<A NAME="#SPEC">
+<H2>�g�p���@�ƃ��\�b�h�̈ꗗ</H2>
+�uRuby�͊��ɏ�����v�Ƃ����O��ŁA<br><br>
+<CODE>
+<PRE>
+require 'bigdecimal'
+a=BigDecimal::new("0.123456789123456789")
+b=BigDecimal::new("123456.78912345678",40)
+c=a+b
+</PRE>
+</CODE>
+<br>
+�Ƃ����悤�Ȋ����Ŏg�p���܂��B
+
+<H3>���\�b�h�ꗗ</H3>
+�ȉ��̂悤�ȃ��\�b�h�����p�”\�ł��B<BR>
+�L�q��ABigDecimal �I�u�W�F�N�g�� a,b,c,r�ŁAString�i������j�I�u�W�F�N�g��
+ s�A������ n �ŕ\�L���܂��B�܂��A�u�L�������v�Ƃ� BigDecimal �����x��ۏ؂���
+�����ł��B�҂�����ł͂���܂���A�኱�̗]�T�������Čv�Z����܂��B�܂��A
+�Ⴆ�΂R�Q�r�b�g�̃V�X�e���ł͂P�O�i�łS�����Ɍv�Z���܂��B�]���āA����ł́A
+�����́u�L�������v�͂S�̔{���ƂȂ��Ă��܂��B
+
+<UL>
+<LI>new</LI><BR>
+�V���� BigDecimal �I�u�W�F�N�g�𐶐����܂��B<BR>
+a=BigDecimal::new(s[,n])<BR>
+s �͏����l�𕶎���Ŏw�肵�܂��D
+n �͕K�v�ȗL�������ia �̍ő�L�������j�𐮐��Ŏw�肵�܂��B
+n �� 0 �܂��͏ȗ����ꂽ�Ƃ��́An �̒l�� s �̗L�������Ƃ݂Ȃ���܂��B
+s �̗L��������� n ���������Ƃ��� n=0 �̂Ƃ��Ɠ����ł��B
+a �̍ő�L�������� n ���኱�傢�l���̗p����܂��B
+<LI>+</LI><BR>
+���Z�ic = a + b�j<BR>
+c �̐��x�ɂ‚��Ắu<A HREF="#PREC">�v�Z���x�ɂ‚���</A>�v���Q�Ƃ��Ă��������B
+<LI>-</LI><BR>
+���Z�ic = a - b�j�A�܂��͕������]�ic = -a�j<BR>
+c �̐��x�ɂ‚��Ắu<A HREF="#PREC">�v�Z���x�ɂ‚���</A>�v���Q�Ƃ��Ă��������B
+<LI>*</LI><BR>
+��Z(c = a * b)<BR>
+c�̐��x��(a�̐��x)+(b�̐��x)���x�ł��B<br>
+�ڂ����́u<A HREF="#PREC">�v�Z���x�ɂ‚���</A>�v���Q�Ƃ��Ă��������B
+
+<LI>/</LI><BR>
+���Z(c = a / b)<BR>
+c �̐��x�ɂ‚��Ắu<A HREF="#PREC">�v�Z���x�ɂ‚���</A>�v���Q�Ƃ��Ă��������B
+
+<LI>assign</LI><BR>
+�ȉ��̂悤�Ɏg�p���܂��B<BR>
+c = a.assign(n,f)<BR>
+f > 0 �Ȃ�Aa �� c �ɁA���̂܂ܑ�����܂��B
+f < 0 �Ȃ�A-a �� c �ɑ�����܂��B
+f �̐�Βl�i|f|�j�͂P���Q���w�肵�Ă��������B
+|f|=2 �̂Ƃ��́Ac �̍ő吸�x�� a �̎����x��菬�����Ƃ��ɂ�
+�ۂ߂��܂��B|f|=1 �̂Ƃ��͐؂�̂Ă��܂��B
+n �� c �̗L�������ł��in ���ȏ�̐��x������ c ����������܂��j�B
+
+<LI>add</LI><BR>
+�ȉ��̂悤�Ɏg�p���܂��B<BR>
+c = a.add(b,n)<BR>
+c = a + b ���ő�� n ���܂Ōv�Z���܂��B
+a + b �̐��x�� n ���傫���Ƃ��͊ۂ߂��܂��B
+
+<LI>sub</LI><BR>
+�ȉ��̂悤�Ɏg�p���܂��B<BR>
+c = a.sub(b,n)<BR>
+c = a - b ���ő�� n ���܂Ōv�Z���܂��B
+a - b �̐��x�� n ���傫���Ƃ��͊ۂ߂��܂��B
+
+<LI>mult</LI><BR>
+�ȉ��̂悤�Ɏg�p���܂��B<BR>
+c = a.mult(b,n)<BR>
+c = a * b ���ő�� n ���܂Ōv�Z���܂��B
+a * b �̐��x�� n ���傫���Ƃ��͊ۂ߂��܂��B
+
+<LI>div</LI><BR>
+�ȉ��̂悤�Ɏg�p���܂��B<BR>
+c,r = a.div(b,n)<BR>
+c=a/b �̌v�Z�����܂��B r �ɂ͏�]���������܂��Ba/b��
+�K�v�Ȃ�n ���܂Ōv�Z����܂��Bdivmod ���\�b�h
+�ƈقȂ�Ac �͐����Ƃ͌���܂���B
+�܂��A c �͊ۂ߂��邱�Ƃ͂���܂���B
+ a = c*b + r �̊֌W�͐������܂��B
+<LI>%</LI><BR>
+r = a%b <BR>
+a/b �̗]����v�Z���܂��B�ȉ��̌v�Z�Ɠ������̂ł��B<BR>
+r = a-((a/b).floor)*b<BR>
+<LI>fix</LI><BR>
+a �̏����_�ȉ��̐؂�̂āB<BR>
+c = a.fix
+<LI>frac</LI><BR>
+a �̐��������̐؂�̂āB<BR>
+c = a.frac
+<LI>floor</LI><BR>
+a �ȉ��̍ő吮����\���l�iBigDecimal �l�j��Ԃ��܂��B<BR>
+c = a.floor<BR>
+�ȉ��̂悤�Ɉ�����^���āA�����_�ȉ� n+1 �ʂ̐����𑀍삷�邱�Ƃ��ł��܂�
+(�����_�ȉ����A�ő� n ���ɂ��܂�)�B<BR>
+c = BigDecimal("1.23456")<BR>
+d = c.floor(4) # d = 1.2345 �ɂȂ�܂��B<BR>
+ n �����̂Ƃ��͏����_�ȏ� n ���ڂ𑀍삵�܂��B<BR>
+c = BigDecimal("15.23456")<BR>
+d = c.floor(-1) # d = 10.0 �ɂȂ�܂��B<BR>
+
+<LI>ceil</LI><BR>
+a �ȏ�̐����̂����A�ł��������������v�Z���A���̒l�iBigDecimal �l�j��Ԃ��܂��B<BR>
+c = a.ceil<BR>
+�ȉ��̂悤�Ɉ�����^���āA�����_�ȉ� n+1 �ʂ̐����𑀍삷�邱�Ƃ��ł��܂�
+(�����_�ȉ����A�ő� n ���ɂ��܂�)�B<BR>
+c = BigDecimal::new("1.23456")<BR>
+d = c.ceil(4) # d = 1.2346 �ɂȂ�܂��B<BR>
+ n �����̂Ƃ��͏����_�ȏ� n ���ڂ��𑀍삵�܂��B<BR>
+c = BigDecimal::new("15.23456")<BR>
+d = c.ceil(-1) # d = 20.0 �ɂȂ�܂��B<BR>
+
+<LI>round</LI><BR>
+�����_�ȉ����ʂ̐����l�̌ܓ����Đ����iBigDecimal �l�j�ɂ��܂��B<BR>
+c = a.round<BR>
+�ȉ��̂悤�Ɉ�����^���āA�����_�ȉ� n+1 �ʂ̐����𑀍삷�邱�Ƃ��ł��܂�
+(�����_�ȉ����A�ő� n ���ɂ��܂�)�B<BR>
+n �����̎��́A�����_�ȉ� n+1 �ʂ̐������l�̌ܓ����܂��B
+c = BigDecimal::new("1.23456")<BR>
+d = c.round(4)  # d = 1.235 �ɂȂ�܂��B<BR>
+ n �����̂Ƃ��͏����_�ȏ� n ���ڂ��𑀍삵�܂��B<BR>
+c = BigDecimal::new("15.23456")<BR>
+d = c.round(-1) # d = 20.0 �ɂȂ�܂��B<BR>
+
+<LI>truncate</LI><BR>
+�����_�ȉ��̐���؂�̂ĂĐ����iBigDecimal �l�j�ɂ��܂��B<BR>
+c = a.truncate<BR>
+�ȉ��̂悤�Ɉ�����^���āA�����_�ȉ� n+1 �ʂ̐����𑀍삷�邱�Ƃ��ł��܂�
+(�����_�ȉ����A�ő� n ���ɂ��܂�)�B<BR>
+n �����̎��́A�����_�ȉ� n+1 �ʂ̐�����؂�̂Ă܂��B
+c = BigDecimal::new("1.23456")<BR>
+d = c.truncate(4)  # d = 1.2345 �ɂȂ�܂��B<BR>
+ n �����̂Ƃ��͏����_�ȏ� n ���ڂ��𑀍삵�܂��B<BR>
+c = BigDecimal::new("15.23456")<BR>
+d = c.truncate(-1) # d = 10.0 �ɂȂ�܂��B<BR>
+
+<LI>divmod</LI><BR>
+���Ə�]�̔z���Ԃ��܂��B<BR>
+c,r = a.divmod(b) # a = c*b + r<BR>
+divmod���\�b�h�� a = c * b + r �ƂȂ� a / b �̕��������_�^�̏� c �Ə�] r ��
+�v�Z���܂��B������ c �͐����i���������̂Ȃ������j�ɂȂ�܂��B<BR>
+c = (a/b).floor <BR>
+r = a - c*b<BR>
+�Ōv�Z����܂��B
+<LI>remainder</LI><BR>
+r=a.remainder(b)<BR>
+a/b �̏�] r ���v�Z���܂��B<BR>
+c = (a/b).fix <BR>
+r = a - c*b<BR>
+�Ōv�Z����܂��B
+
+<LI>abs</LI><BR>
+���̐�Βl<BR>
+c = a.abs<BR>
+
+<LI>to_i</LI><BR>
+�����_�ȉ���؂�̂ĂĐ����ɕϊ����܂��B<BR>
+i = a.to_i<BR>
+i �͒l�ɉ����� Fixnum �� Bignum �ɂȂ�܂��B
+a �� Infinity �� NaN �̂Ƃ��Ai �� nil �ɂȂ�܂��B
+<LI>to_f</LI><BR>
+dup �ƑS�������ł��B
+�����l�� BigDecimal �I�u�W�F�N�g�𐶐����܂��B
+<LI>to_s</LI><BR>
+������ɕϊ����܂�("0.xxxxxEn"�̌`�ɂȂ�܂��j�B<BR>
+s = a.to_s
+<LI>to_s2</LI><BR>
+������ɕϊ����܂��B���������� n �����ɋ󔒂ŋ�؂�܂��B<BR>
+s = a.to_s2(n)
+<LI>exponent</LI><BR>
+�w�����𐮐��l�ŕԂ��܂��B
+n = a.exponent <BR>
+�� a �̒l�� 0.xxxxxxx*10**n ���Ӗ����܂��B
+
+<LI>E</LI><BR>
+���R�ΐ��̒�e(=2.718281828....)���v�Z���܂��i�����Ƀe�C���[�W�J�Łj�B<BR>
+e = BigDecimal::E(n)<BR>
+n�͕K�v�ȗL�������𐮐��Ŏw�肵�܂��B
+<LI>PI</LI><BR>
+�~����(=3.14159265358979....)���v�Z���܂��iJ�DMachin�̌�����p���܂��j�B<BR>
+e = BigDecimal::PI(n)<BR>
+n �͕K�v�ȗL�������𐮐��Ŏw�肵�܂��B
+<LI>BASE</LI><BR>
+�����Ŏg�p������̒l�ł��B������ 32 �r�b�g�̏����n�ł�10000�ł��B<BR>
+b = BigDecimal::BASE<BR>
+<LI>mode</LI><BR>
+BigDecimal�̎��s���ʂ𐧌䂵�܂��B�ȉ��̎g�p���@����`����Ă��܂��B<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_NaN,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_INFINITY,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_UNDERFLOW,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_OVERFLOW,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_ZERODIVIDE,flag)<BR>
+f = BigDecimal::mode(BigDecimal::EXCEPTION_ALL,flag)<BR>
+
+EXCEPTION_NaN �͌��ʂ� NaN �ɂȂ����Ƃ��̎w��ł��B
+EXCEPTION_INFINITY �͌��ʂ�������(�}Infinity)
+�ɂȂ����Ƃ��̎w��ł��B
+EXCEPTION_UNDERFLOW �͎w�������A���_�[�t���[����Ƃ��̎w��ł��B
+EXCEPTION_OVERFLOW �͎w�������I�[�o�[�t���[����Ƃ��̎w��ł��B
+EXCEPTION_ZERODIVIDE �̓[���ɂ�銄��Z�����s�����Ƃ��̎w��ł��B
+EXCEPTION_ALL �́A�”\�ȑS�Ăɑ΂��Ĉꊇ���Đݒ肷��Ƃ���
+�g�p���܂��B
+flag �� true �̂Ƃ��́A�w�肵����ԂɂȂ����Ƃ��ɗ�O�𔭍s
+����悤�ɂȂ�܂��B
+flag �� false�i�f�t�H���g�j�Ȃ�A��O�͔��s����܂���B�v�Z���ʂ�
+�ȉ��̂悤�ɂȂ�܂��B<BR>
+EXCEPTION_NaN �̂Ƃ��A��(NaN)<BR>
+EXCEPTION_INFINITY �̂Ƃ��A����(+ or -Infinity)<BR>
+EXCEPTION_UNDERFLOW �̂Ƃ��A�[��<BR>
+EXCEPTION_OVERFLOW �̂Ƃ��A+Infinity �� -Infinity<BR>
+EXCEPTION_ZERODIVIDE �̂Ƃ��A+Infinity �� -Infinity<BR>
+EXCEPTION_INFINITY�AEXCEPTION_OVERFLOW�AEXCEPTION_ZERODIVIDE
+�͍��̂Ƃ��듯���ł��B<BR>
+�߂�l�́A�ݒ��̒l�ł��B�u�l�v�̈Ӗ��́A�Ⴆ��
+BigDecimal::EXCEPTION_NaN�Ɓu�l�v�� & �� �[���ȊO�Ȃ��
+EXCEPTION_NaN���ݒ肳��Ă���Ƃ����Ӗ��ł��B
+flag �� nil�A�܂��́Atrue �� false �ȊO�Ȃ猻�݂̐ݒ�l���Ԃ�܂��B
+
+<LI>limit([n])</LI><BR>
+���������BigDecimal�I�u�W�F�N�g�̍ő包����n���ɐ������܂��B�߂�l��
+�ݒ肷��O�̒l�ł��B�ݒ�l�̃f�t�H���g�l�͂O�ŁA�����������Ƃ����Ӗ��ł��B
+n���w�肵�Ȃ��ꍇ�́A����̍ő包�����Ԃ�܂��B<BR>
+mf = BigDecimal::limit(n)<BR>
+<LI>sign</LI><BR>
+�l�̑�����Ԃ��܂��B
+n = a.sign <BR>
+�Ƃ����Ƃ� n �̒l�� a ���ȉ��̂Ƃ����Ӗ����܂��B<BR>
+() �̒��̐����́A���ۂ̒l�ł�(<A HREF="#STRUCT">�u�����\���v</A>���Q��)�B<BR>
+n = BigDecimal::SIGN_NaN(0) : a �� NaN<BR>
+n = BigDecimal::SIGN_POSITIVE_ZERO(1) : a �� +0<BR>
+n = BigDecimal::SIGN_NEGATIVE_ZERO(-1) : a �� -0<BR>
+n = BigDecimal::SIGN_POSITIVE_FINITE(2) : a �͐��̒l<BR>
+n = BigDecimal::SIGN_NEGATIVE_FINITE(-2) : a �͕��̒l<BR>
+n = BigDecimal::SIGN_POSITIVE_INFINITE(3) : a ��+Infinity<BR>
+n = BigDecimal::SIGN_NEGATIVE_INFINITE(-3) : a ��-Infinity<BR>
+
+<LI>nan?</LI><BR>
+a.nan? �� a ��NaN�̂Ƃ��^��Ԃ��܂��B
+<LI>infinite?</LI><BR>
+a.infinite? �� a ��+���܂���-���̂Ƃ��^��Ԃ��܂��B
+<LI>finite?</LI><BR>
+a.finite? �� a �����܂��� NaN �łȂ��Ƃ��^��Ԃ��܂��B
+
+<LI>to_parts</LI><BR>
+BigDecimal �l�� 0.xxxxxxx*10**n �ƕ\�������Ƃ��ɁA�����iNaN�̂Ƃ���
+0�A����ȊO��+1��-1�ɂȂ�܂��j�A
+���������̕�����i"xxxxxxx"�j�ƁA��i10�j�A�X�Ɏw�� n ��z���
+�Ԃ��܂��B<BR>
+a=BigDecimal::new("3.14159265",10)<BR>
+f,x,y,z = a.to_parts<BR>
+�Ƃ���ƁAf=+1�Ax="314159265"�Ay=10�Az=1�ɂȂ�܂��B<BR>
+�]���āA<BR>
+s = "0."+x<BR>
+b = f*(s.to_f)*(y**z)<BR>
+�� Float �ɕϊ����邱�Ƃ��ł��܂��B
+<LI>inspect</LI><BR>
+�f�o�b�O�o�͂Ɏg�p����܂��B<BR>
+p a=BigDecimal::new("3.14",10)<BR>
+�Ƃ���ƁA[0x112344:'0.314E1',4(12)]�̂悤�ɏo�͂���܂��B
+�ŏ���16�i���̓I�u�W�F�N�g�̃A�h���X�A���� '0.314E1' �͒l�A
+����4�͌��݂̗L������(�\�����኱�傫�����Ƃ�����܂�)�A
+�Ō�̓I�u�W�F�N�g����蓾��ő包���ɂȂ�܂��B
+
+<LI>dup</LI><BR>
+�����l�� BigDecimal �I�u�W�F�N�g�𐶐����܂��B
+<LI>sqrt</LI><BR>
+a�̗L���� n ���̕������in �̕������ł͂���܂���j�B
+����܂��A�����Ƀj���[�g���@�Ōv�Z���܂��B<BR>
+c = a.sqrt(n)<BR>
+<LI>sincos</LI><BR>
+a �̗L���� n ���� sin �� cos �𓯎��Ɂi�e�C���[�W�J�Łj�v�Z���āA
+ sin �� cos �̔z���Ԃ��܂��B
+n �͕K�v�ȗL�������ł��i n �� sin �� cos ���v�Z����킯�ł͂���܂���j�B
+<BR>
+sin,cos = a.sincos(n)<BR>
+<LI>exp</LI><BR>
+���R�ΐ��̒�e(=2.718281828....)�� a ����v�Z���܂��B<BR>
+c = a.exp(n)<BR>
+n �͕K�v�ȗL�������ł��B
+<LI>power</LI><BR>
+a �� n ����v�Z���܂��B���͐����B<BR>
+c = a.power(n)<BR>
+���ʂƂ��� c �̗L������ a �� n �{�ȏ�ɂȂ�̂Œ��ӁB
+<LI>zero?</LI><BR>
+a �� 0 �Ȃ� true �ɂȂ�܂��B<BR>
+c = a.zero?<BR>
+<LI>nonzero?</LI><BR>
+a �� 0 �Ȃ� false�A0 �ȊO�Ȃ� a ���̂��̂��Ԃ�܂��B<BR>
+c = a.nonzero?<BR>
+<LI>&lt=&gt</LI><BR>
+a==b �Ȃ� 0�Aa &gt b �Ȃ� 1�Aa &lt b �Ȃ� -1 �ɂȂ�܂��B<BR>
+c = a &lt=&gt b <BR>
+</UL>
+��́A�ǂ�Ŏ��̔@���ł��B<BR>
+<UL>
+<LI>==</LI>
+<LI>===</LI>
+�u==�v�Ɠ����ł��� case ���Ŏg�p����܂��B
+<LI>!=</LI>
+<LI>&lt</LI>
+<LI>&lt=</LI>
+<LI>&gt</LI>
+<LI>&gt=</LI>
+</UL>
+<H3>coerce�ɂ‚���</H3>
+BigDecimal �I�u�W�F�N�g���Z�p���Z�q�̍��ɂ���Ƃ��́ABigDecimal �I�u�W�F�N�g��
+�E�ɂ���I�u�W�F�N�g��(�K�v�Ȃ�) BigDecimal �ɕϊ����Ă���v�Z���܂��B
+�]���āABigDecimal �I�u�W�F�N�g�ȊO�ł����l���Ӗ�������̂Ȃ�E�ɒu����
+���Z�͉”\�ł��B<BR><BR>
+������Ő��l��^����ꍇ�͒��ӂ��K�v�ł��B���l�ɕϊ��ł��Ȃ�����������ƁA
+�P�ɕϊ����~�߂邾���ŃG���[�ɂ͂Ȃ�܂���B"10XX"�Ȃ�P�O�A"XXXX"�͂O
+�ƈ����܂��B<BR>
+<CODE><PRE>
+   a = BigDecimal.E(20)
+   c = a * "0.123456789123456789123456789" # ������ BigDecimal �ɕϊ����Ă���v�Z
+</PRE></CODE>
+�������񐔂�\�������Ƃ��āA"Infinity"�A"+Infinity"�A"-Infinity"�A"NaN"
+���g�p�ł��܂�(�啶���E����������ʂ��܂�)�B�������Amode ���\�b�h�� false ��
+�w�肵���ꍇ�͗�O���������܂��B
+<BR>
+�܂��ABigDecimal�N���X�� coerce�iRuby�{�Q�Ɓj���T�|�[�g���Ă��܂��B
+�]���āABigDecimal �I�u�W�F�N�g���E�ɂ���ꍇ�����͑��v�ł��B
+�����A���݂� Ruby �C���^�v���^�̎d�l��A�����񂪍��ɂ���ƌv�Z�ł��܂���B<BR>
+<CODE><PRE>
+  a = BigDecimal.E(20)
+  c = "0.123456789123456789123456789" * a # �G���[
+</PRE></CODE>
+�K�v��������Ƃ͎v���܂��񂪁A�ǂ����Ă��ƌ����l��
+ String �I�u�W�F�N�g���p�������V���ȃN���X���쐬���Ă���A
+���̃N���X�� coerce ���T�|�[�g���Ă��������B
+
+<hr>
+<A NAME="#UNDEF">
+<H2>�����A�񐔁A�[���̈���</H2>
+�u�����v�Ƃ͕\���ł��Ȃ����炢�傫�Ȑ��ł��B���ʂɈ������߂�
+ +Infinity�i���̖�����j�� -Infinity�i���̖�����j�Ƃ���
+�悤�ɕ\�L����܂��B
+������ 1.0/0.0 �̂悤�Ƀ[���Ŋ���悤�Ȍv�Z�������Ƃ��ɐ�������܂��B
+<BR><BR>
+�u�񐔁v�� 0.0/0.0 �� Infinity-Infinity ���̌��ʂ���`�ł��Ȃ�
+�v�Z�������Ƃ��ɐ�������܂��B�񐔂� NaN�iNot a Number�j�ƕ\�L����܂��B
+NaN ���܂ތv�Z�͑S�� NaN �ɂȂ�܂��B�܂� NaN �͎������܂߂āA�ǂ�Ȑ�
+�Ƃ���v���܂���B
+<BR><BR>
+�[���� +0.0 �� -0.0 �����݂��܂��B�������A+0.0==-0.0 �� true �ł��B
+<BR><BR>
+Infinity�ANaN�A +0.0 �� -0.0 �����܂񂾌v�Z���ʂ͑g�ݍ��킹��
+��蕡�G�ł��B�����̂���l�́A�ȉ��̃v���O���������s���Č��ʂ�
+�m�F���Ă��������i���ʂɂ‚��āA�^���ԈႢ�𔭌����ꂽ����
+���m�点�肢�܂��j�B
+
+<PRE>
+<CODE>
+require "bigdecimal"
+
+aa  = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
+ba  = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
+opa = %w(+ - * / <=> > >=  < == != <=)
+
+for a in aa
+  for b in ba
+    for op in opa
+      x = BigDecimal::new(a)
+      y = BigDecimal::new(b)
+      eval("ans= x #{op} y;print a,' ',op,' ',b,' ==> ',ans.to_s,\"\n\"")
+    end
+  end
+end
+</CODE>
+</PRE>
+
+<hr>
+<A NAME="#STRUCT">
+<H2>�����\��</H2>
+BigDecimal�����ŕ��������_�͍\����(Real)�ŕ\������܂��B
+���̂����������� unsigned long �̔z��(�ȉ��̍\���̗v�ffrac)�ŊǗ�����܂��B
+�T�O�I�ɂ́A�ȉ��̂悤�ɂȂ�܂��B<BR><BR>
+ <���������_��> = 0.xxxxxxxxx*BASE**n<BR><BR>
+�����ŁAx�͉�������\�������ABASE�͊�i�P�O�i�Ȃ�P�O�j�An�͎w������\��
+�����l�ł��BBASE���傫���قǁA�傫�Ȑ��l���\���ł��܂��B�‚܂�A�z��̃T�C�Y��
+���Ȃ��ł��܂��BBASE�͑傫���قǓs�����悢�킯�ł����A�f�o�b�O�̂��₷���Ȃǂ�
+�l�����āA10000�ɂȂ��Ă��܂��iBASE��VpInit()�֐��Ŏ����I�Ɍv�Z���܂��j�B
+����́A32�r�b�g�����̏ꍇ�ł��B64�r�b�g�����̏ꍇ�͂����Ƒ傫�Ȓl�ɂȂ�܂��B
+�c�O�Ȃ���A64�r�b�g�����ł̃e�X�g�͂܂�����Ă��܂���i�����A���ꂽ���������
+���ʂ������Ă���������΂��肪�����ł��j�B
+BASE��10000�̂Ƃ��́A�ȉ��̉������̔z��(frac)�̊e�v�f�ɂ͍ő�łS����
+�������i�[����܂��B<BR><BR>
+���������_�\����(Real)�͈ȉ��̂悤�ɂȂ��Ă��܂��B
+<BR>
+<CODE><PRE>
+  typedef struct {
+     unsigned long MaxPrec; // �ő吸�x(frac[]�̔z��T�C�Y)
+     unsigned long Prec;    // ���x(frac[]�̎g�p�T�C�Y)
+     short    sign;         // �ȉ��̂悤�ɕ������̏�Ԃ��`���܂��B
+                            //  ==0 : NaN
+                            //    1 : +0
+                            //   -1 : -0
+                            //    2 : ���̒l
+                            //   -2 : ���̒l
+                            //    3 : +Infinity
+                            //   -3 : -Infinity
+     unsigned short flag;   // �e��̐���t���b�O
+     int      exponent;     // �w�����̒l(������*BASE**exponent)
+     unsigned long frac[1]; // �������̔z��(�•�)
+  } Real;
+</CODE></PRE>
+�Ⴆ�� 1234.56784321 �Ƃ���������(BASE=10000�Ȃ�)<BR>
+<PRE>
+    0.1234 5678 4321*(10000)**1
+</PRE>
+�ł����� frac[0]=1234�Afrac[1]=5678�Afrac[2]=4321�A
+Prec=3�Asign=2�Aexponent=1 �ƂȂ�܂��BMaxPrec��
+Prec ���傫����΂����‚ł����܂��܂���Bflag ��
+�g�p���@�͎����Ɉˑ����ē����Ŏg�p����܂��B
+
+<hr>
+<A NAME="#BASE">
+<H2>2�i��10�i</H2>
+BigDecimal �� <���������_��> = 0.xxxxxxxxx*10**n �Ƃ���10�i�`���Ő��l��ێ����܂��B
+�������A�v�Z�@�̕��������_���̓����\���́A�����܂ł��Ȃ� <���������_��> = 0.bbbbbbbb*2**n �Ƃ���
+2�i�`�������ʂł�(x �� 0 ���� 9 �܂ŁAb �� 0 �� 1 �̐���)�B
+BigDecimal ���Ȃ�10�i�̓����\���`�����̗p�����̂����ȉ��ɐ������܂��B
+<H4>10�i�̃����b�g</H4>
+<DL>
+<DT>�f�o�b�O�̂��₷��
+<DD>�܂��A�v���O�����쐬���y�ł��Bfrac[0]=1234�Afrac[1]=5678�Afrac[2]=4321�A
+exponent=1�Asign=2 �Ȃ琔�l�� 1234.56784321 �ł���̂͌���Β����ɕ�����܂��B
+
+<DT>10�i�\�L���ꂽ���l�Ȃ�m���ɓ����\���ɕϊ��ł���
+<DD>�Ⴆ�΁A�ȉ��̂悤�ȃv���O�����͑S���덷������
+�v�Z���邱�Ƃ��ł��܂��B�ȉ��̗�́A��s�Ɉ�‚̐��l
+�������Ă���t�@�C�� file �̍��v���l�����߂���̂ł��B
+<PRE><CODE>
+   file = File::open(....,"r")
+   s = BigDecimal::new("0")
+   while line = file.gets
+      s = s + line
+   end
+</CODE></PRE>
+���̗��2�i���ł��ƌ덷�����荞�މ”\��������܂��B
+�Ⴆ�� 0.1 ��2�i�ŕ\������� 0.1 = b1*2**(-1)+b1*2**(-2)+b3*2**(-3)+b4*2**(-4)....
+�Ɩ����ɑ����Ă��܂��܂�(b1=0,b2=0,b3=0,b4=1...)�B������ bn(n=1,2,3,...) ��
+2�i��\������ 0 �� 1 �̐�����ł��B�]���āA�ǂ����őł��؂�K�v������܂��B
+�����ŕϊ��덷������܂��B�������A������ēx10�i�\�L�ɂ��Ĉ������悤��
+�ꍇ�͓K�؂Ȋۂߑ���i�l�̌ܓ��j�ɂ���čĂ� "0.1" �ƕ\������܂��B�������A
+�����ł͐��m�� 0.1 �ł͂���܂���B
+
+<DT>�L�������͗L���ł���i�‚܂莩������ł���j
+<DD>0.1 ��\�����邽�߂̗̈�͂�������‚̔z��v�f�i frac[0]=1 �j�ōς݂܂��B
+�z��v�f�̐���10�i���l���玩���I�Ɍ���ł��܂��B����́A�•ϒ����������_���Z�ł�
+�厖�Ȃ��Ƃł��B�t�� 0.1 ��2�i�\�������Ƃ��ɂ�2�i�̗L�����������‚ɂ���̂� 0.1 ��
+���������ł͌���ł��܂���B
+</DL>
+
+<H3>10�i�̃f�����b�g</H3>
+���͍��܂ł̃����b�g�́A���̂܂܃f�����b�g�ɂ��Ȃ�܂��B
+���������A10�i��2�i�A2�i��10�i�ɕϊ�����悤�ȑ���͕ϊ��덷
+�𔺂��ꍇ��������邱�Ƃ͂ł��܂���B
+���Ɍv�Z�@�����Ɏ�荞�܂ꂽ2�i���l�� BigDecimal �̓����\����
+�ϊ�����Ƃ��ɂ͌덷���������Ȃ��ꍇ������܂��B
+
+<H3>�ŏ��͉����H</H3>
+�����Ōv�Z����Ƃ��ɂ킴�킴2�i�����g���l�͋ɂ߂Ă܂�ł��B
+�v�Z�@�Ƀf�[�^����͂���Ƃ����قƂ�ǂ̏ꍇ�A
+10�i���œ��͂��܂��B���̌��ʁAdouble ���̌v�Z�@����
+�\���͍ŏ�����덷�������Ă���ꍇ������܂��B
+BigDecimal �̓��[�U���͂��덷�����Ŏ�荞�ނ��Ƃ��ł��܂��B
+�f�o�b�O�����₷���̂ƁA�f�[�^�ǂ݂��ݎ��Ɍ덷������Ȃ�
+�Ƃ����̂����ۂ̃����b�g�ł��B
+
+<hr>
+<A NAME="#PREC">
+<H2>�v�Z���x�ɂ‚���</H2>
+c = a op b �Ƃ����v�Z(op �� + - * /)�������Ƃ��̓����
+�ȉ��̂悤�ɂȂ�܂��B<BR><BR>
+�P�D��Z�Ə��Z��(a �̗L������)+(a �̗L������)���̍ő包���i���ۂ́A�]�T�������āA
+���������傫���Ȃ�܂��j�����•ϐ� c ��V���ɐ������܂��B
+�����Z�̏ꍇ�́A�덷���o�Ȃ������̐��x������ c �𐶐����܂��B�Ⴆ��
+ c = 0.1+0.1*10**(-100) �̂悤�ȏꍇ�Ac �̐��x�͂P�O�O���ȏ�̐��x��
+���‚悤�ɂȂ�܂��B
+<BR>
+�Q�D���� c = a op b �̌v�Z�����s���܂��B<BR><BR>
+���̂悤�ɁA�����Z�Ə�Z�ł� c �͕K���u�덷���o�Ȃ��v�����̐��x��
+�����Đ�������܂��B���Z��(a �̗L������)+(a �̗L������)���̍ő包��
+������ c ����������܂����Ac = 1.0/3.0 �̂悤�Ȍv�Z�Ŗ��炩�Ȃ悤�ɁA
+ c �̍ő吸�x�𒴂���Ƃ���Ōv�Z���ł��؂���ꍇ������܂��B<BR><BR>
+������ɂ���Ac �̍ő吸�x�� a �� b ���傫���Ȃ�܂��̂� c ���K�v�Ƃ���
+�������[�̈�͑傫���Ȃ邱�Ƃɒ��ӂ��ĉ������B
+<BR><BR>
+���ӁF�u+,-,*,/�v�ł͌��ʂ̐��x�i�L�������j�������Ŏw��ł��܂���B
+���x���R���g���[���������ꍇ�́A�ȉ��� add,sub ���̃��\�b�h
+���g�p���܂��B<BR>
+
+<H3>�����Ő��x���R���g���[���������ꍇ</H3>
+�����Ő��x(�L������)���R���g���[���������ꍇ�� assign�Aadd�Asub�Amult�Adiv ���̃��\�b�h
+���g�p�ł��܂��B
+�ȉ��̉~�������v�Z����v���O������̂悤�ɁA
+���߂錅���͎����Ŏw�肷�邱�Ƃ��ł��܂��B
+<BR><BR>
+<CODE><PRE>
+#!/usr/local/bin/ruby
+
+#
+# pai.rb
+#  USAGE: ruby pai.rb n
+#   where n is the number of digits required.
+#  EX.: ruby pai.rb 1000
+#
+
+require "bigdecimal"
+#
+# Calculates 3.1415.... using J. Machin's formula.
+#
+def pai(sig) # sig: Number of significant figures
+  exp    = -sig
+  pi     = BigDecimal::new("0")
+  two    = BigDecimal::new("2")
+  m25    = BigDecimal::new("-0.04")
+  m57121 = BigDecimal::new("-57121")
+
+  u = BigDecimal::new("1")
+  k = BigDecimal::new("1")
+  w = BigDecimal::new("1")
+  t = BigDecimal::new("-80")
+  while (u.exponent >= exp) 
+    t   = t*m25
+    u,r = t.div(k,sig)
+    pi  = pi + u
+    k   = k+two
+  end
+
+  u = BigDecimal::new("1")
+  k = BigDecimal::new("1")
+  w = BigDecimal::new("1")
+  t = BigDecimal::new("956")
+  while (u.exponent >= exp )
+    t,r = t.div(m57121,sig)
+    u,r = t.div(k,sig)
+    pi  = pi + u
+    k   = k+two
+  end
+  pi
+end
+
+if $0 == __FILE__
+  print "PAI("+ARGV[0]+"):\n"
+  p pai(ARGV[0].to_i)
+end
+
+</PRE></CODE>
+<HR>
+<FONT size=2>
+<I>
+<A HREF="http://www.tinyforest.gr.jp">
+���� �ΗY
+</A>
+(E-Mail:<A HREF="mailto:[email protected]">&[email protected]&gt</U></A>)
+</I>
+</FONT>
+</TD>
+</TR>
+</TABLE>
+</BODY>
+</HTML>
diff --git a/ext/bigdecimal/depend b/ext/bigdecimal/depend
new file mode 100644
index 0000000000..402cae95dd
--- /dev/null
+++ b/ext/bigdecimal/depend
@@ -0,0 +1 @@
+bigdecimal.o: bigdecimal.c bigdecimal.h $(hdrdir)/ruby.h
diff --git a/ext/bigdecimal/extconf.rb b/ext/bigdecimal/extconf.rb
new file mode 100644
index 0000000000..a68a656044
--- /dev/null
+++ b/ext/bigdecimal/extconf.rb
@@ -0,0 +1,2 @@
+require 'mkmf'
+create_makefile('bigdecimal')
diff --git a/ext/bigdecimal/lib/bigdecimal-rational.rb b/ext/bigdecimal/lib/bigdecimal-rational.rb
new file mode 100644
index 0000000000..de999d19cf
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal-rational.rb
@@ -0,0 +1,30 @@
+#
+# BigDecimal <-> Rational 
+#
+class BigDecimal
+    # Convert BigDecimal to Rational
+    def to_r 
+       sign,digits,base,power = self.to_parts
+       numerator = sign*digits.to_i
+       denomi_power = power - digits.size # base is always 10
+       if denomi_power < 0
+          denominator = base ** (-denomi_power)
+       else
+          denominator = base ** denomi_power
+       end
+       Rational.new(numerator,denominator)
+    end
+end
+
+class Rational
+  # Convert Rational to BigDecimal
+  # to_d returns an array [quotient,residue]
+  def to_d(nFig=0)
+     num = self.numerator.to_s
+     if nFig<=0
+        nFig = BigDecimal.double_fig*2+1
+     end
+     BigDecimal.new(num).div(self.denominator,nFig)
+  end
+end
+
diff --git a/ext/bigdecimal/lib/delcr b/ext/bigdecimal/lib/delcr
new file mode 100644
index 0000000000..df4483ce34
--- /dev/null
+++ b/ext/bigdecimal/lib/delcr
@@ -0,0 +1,7 @@
+tr -d '\r' < bigdecimal-rational.rb > unix/bigdecimal-rational.rb
+tr -d '\r' <jacobian.rb > unix/jacobian.rb
+tr -d '\r' <linear.rb > unix/linear.rb
+tr -d '\r' <ludcmp.rb > unix/ludcmp.rb
+tr -d '\r' <newton.rb > unix/newton.rb
+tr -d '\r' <nlsolve.rb > unix/nlsolve.rb
+tr -d '\r' <pai.rb > unix/pai.rb
diff --git a/ext/bigdecimal/lib/jacobian.rb b/ext/bigdecimal/lib/jacobian.rb
new file mode 100644
index 0000000000..34a60ae67a
--- /dev/null
+++ b/ext/bigdecimal/lib/jacobian.rb
@@ -0,0 +1,63 @@
+#
+# jacobian.rb
+#
+# Computes Jacobian matrix of f at x
+#
+module Jacobian
+  def isEqual(a,b,zero=0.0,e=1.0e-8)
+    aa = a.abs
+    bb = b.abs
+    if aa == zero &&  bb == zero then
+          true
+    else
+          if ((a-b)/(aa+bb)).abs < e then
+             true
+          else
+             false
+          end
+    end
+  end
+
+  def dfdxi(f,fx,x,i)
+    nRetry = 0
+    n = x.size
+    xSave = x[i]
+    ok = 0
+    ratio = f.ten*f.ten*f.ten
+    dx = x[i].abs/ratio
+    dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps)
+    dx = f.one/f.ten     if isEqual(dx,f.zero,f.zero,f.eps)
+    until ok>0 do
+      s = f.zero
+      deriv = []
+      if(nRetry>100) then
+         raize "Singular Jacobian matrix. No change at x[" + i.to_s + "]"
+      end
+      dx = dx*f.two
+      x[i] += dx
+      fxNew = f.values(x)
+      for j in 0...n do
+        if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then
+           ok += 1
+           deriv <<= (fxNew[j]-fx[j])/dx
+        else
+           deriv <<= f.zero
+        end
+      end
+      x[i] = xSave
+    end
+    deriv
+  end
+
+  def jacobian(f,fx,x)
+    n = x.size
+    dfdx = Array::new(n*n)
+    for i in 0...n do
+      df = dfdxi(f,fx,x,i)
+      for j in 0...n do
+         dfdx[j*n+i] = df[j]
+      end
+    end
+    dfdx
+  end
+end
diff --git a/ext/bigdecimal/lib/linear.rb b/ext/bigdecimal/lib/linear.rb
new file mode 100644
index 0000000000..f93404fb6f
--- /dev/null
+++ b/ext/bigdecimal/lib/linear.rb
@@ -0,0 +1,46 @@
+#!/usr/local/bin/ruby
+
+#
+# linear.rb
+#
+# Solves linear equation system(A*x = b) by LU decomposition method.
+#  where  A is a coefficient matrix,x is an answer vector,b is a constant vector.
+#
+require "bigdecimal"
+require "ludcmp"
+
+include LUSolve
+
+def rd_order
+   printf("Number of equations ?")
+   n = gets().chomp.to_i
+end
+
+zero = BigDecimal::new("0.0")
+one  = BigDecimal::new("1.0")
+
+while (n=rd_order())>0
+  a = []
+  as= []
+  b = []
+  printf("\nEnter coefficient matrix element A[i,j]\n");
+  for i in 0...n do
+    for j in 0...n do
+      printf("A[%d,%d]? ",i,j); s = gets
+      a <<=BigDecimal::new(s);
+      as<<=BigDecimal::new(s);
+    end
+    printf("Contatant vector element b[%d] ? ",i);b<<=BigDecimal::new(gets);
+  end
+  printf "ANS="
+  x = lusolve(a,b,ludecomp(a,n,zero,one),zero)
+  p x
+  printf "A*x-b\n"
+  for i in 0...n do
+    s = zero
+    for j in 0...n do
+       s = s + as[i*n+j]*x[j]
+    end
+    p s-b[i]
+  end
+end
diff --git a/ext/bigdecimal/lib/ludcmp.rb b/ext/bigdecimal/lib/ludcmp.rb
new file mode 100644
index 0000000000..c36f0dea5b
--- /dev/null
+++ b/ext/bigdecimal/lib/ludcmp.rb
@@ -0,0 +1,75 @@
+#
+# ludcmp.rb
+#
+module LUSolve
+  def ludecomp(a,n,zero=0.0,one=1.0)
+    ps     = []
+    scales = []
+    for i in 0...n do  # pick up largest(abs. val.) element in each row.
+      ps <<= i
+      nrmrow  = zero
+      ixn = i*n
+      for j in 0...n do
+         biggst = a[ixn+j].abs
+         nrmrow = biggst if biggst>nrmrow
+      end
+      if nrmrow>zero then
+         scales <<= one/nrmrow
+      else 
+         raise "Singular matrix"
+      end
+    end
+    n1          = n - 1
+    for k in 0...n1 do # Gaussian elimination with partial pivoting.
+      biggst  = zero;
+      for i in k...n do
+         size = a[ps[i]*n+k].abs*scales[ps[i]]
+         if size>biggst then
+            biggst = size
+            pividx  = i
+         end
+      end
+      raise "Singular matrix" if biggst<=zero
+      if pividx!=k then
+        j = ps[k]
+        ps[k] = ps[pividx]
+        ps[pividx] = j
+      end
+      pivot   = a[ps[k]*n+k]
+      for i in (k+1)...n do
+        psin = ps[i]*n
+        a[psin+k] = mult = a[psin+k]/pivot
+        if mult!=zero then
+           pskn = ps[k]*n
+           for j in (k+1)...n do
+             a[psin+j] -= mult*a[pskn+j]
+           end
+        end
+      end
+    end
+    raise "Singular matrix" if a[ps[n1]*n+n1] == zero
+    ps
+  end
+
+  def lusolve(a,b,ps,zero=0.0)
+    n = ps.size
+    x = []
+    for i in 0...n do
+      dot = zero
+      psin = ps[i]*n
+      for j in 0...i do
+        dot = a[psin+j]*x[j] + dot
+      end
+      x <<= b[ps[i]] - dot
+    end
+    (n-1).downto(0) do |i|
+       dot = zero
+       psin = ps[i]*n
+       for j in (i+1)...n do
+         dot = a[psin+j]*x[j] + dot
+       end
+       x[i]  = (x[i]-dot)/a[psin+i]
+    end
+    x
+  end
+end
diff --git a/ext/bigdecimal/lib/newton.rb b/ext/bigdecimal/lib/newton.rb
new file mode 100644
index 0000000000..67a92474ac
--- /dev/null
+++ b/ext/bigdecimal/lib/newton.rb
@@ -0,0 +1,75 @@
+#
+# newton.rb 
+#
+# Solves nonlinear algebraic equation system f = 0 by Newton's method.
+#  (This program is not dependent on BigDecimal)
+#
+# To call:
+#    n = nlsolve(f,x)
+#  where n is the number of iterations required.
+#        x is the solution vector.
+#        f is the object to be solved which must have following methods.
+#
+#   f ... Object to compute Jacobian matrix of the equation systems.
+#       [Methods required for f]
+#         f.values(x) returns values of all functions at x.
+#         f.zero      returns 0.0
+#         f.one       returns 1.0
+#         f.two       returns 1.0
+#         f.ten       returns 10.0
+#         f.eps       convergence criterion
+#   x ... initial values
+#
+require "ludcmp"
+require "jacobian"
+
+module Newton
+  include LUSolve
+  include Jacobian
+  
+  def norm(fv,zero=0.0)
+    s = zero
+    n = fv.size
+    for i in 0...n do
+      s += fv[i]*fv[i]
+    end
+    s
+  end
+
+  def nlsolve(f,x)
+    nRetry = 0
+    n = x.size
+
+    f0 = f.values(x)
+    zero = f.zero
+    one  = f.one
+    two  = f.two
+    p5 = one/two
+    d  = norm(f0,zero)
+    minfact = f.ten*f.ten*f.ten
+    minfact = one/minfact
+    e = f.eps
+    while d >= e do
+      nRetry += 1
+      # Not yet converged. => Compute Jacobian matrix
+      dfdx = jacobian(f,f0,x)
+      # Solve dfdx*dx = -f0 to estimate dx
+      dx = lusolve(dfdx,f0,ludecomp(dfdx,n,zero,one),zero)
+      fact = two
+      xs = x.dup
+      begin
+        fact *= p5
+        if fact < minfact then
+          raize "Failed to reduce function values."
+        end
+        for i in 0...n do
+          x[i] = xs[i] - dx[i]*fact
+        end
+        f0 = f.values(x)
+        dn = norm(f0,zero)
+      end while(dn>=d)
+      d = dn
+    end
+    nRetry
+  end
+end
diff --git a/ext/bigdecimal/lib/nlsolve.rb b/ext/bigdecimal/lib/nlsolve.rb
new file mode 100644
index 0000000000..08f17f9ecd
--- /dev/null
+++ b/ext/bigdecimal/lib/nlsolve.rb
@@ -0,0 +1,38 @@
+#!/usr/local/bin/ruby
+
+#
+# nlsolve.rb
+# An example for solving nonlinear algebraic equation system.
+#
+
+require "bigdecimal"
+require "newton"
+include Newton
+
+class Function
+  def initialize()
+    @zero = BigDecimal::new("0.0")
+    @one  = BigDecimal::new("1.0")
+    @two  = BigDecimal::new("2.0")
+    @ten  = BigDecimal::new("10.0")
+    @eps  = BigDecimal::new("1.0e-16")
+  end
+  def zero;@zero;end
+  def one ;@one ;end
+  def two ;@two ;end
+  def ten ;@ten ;end
+  def eps ;@eps ;end
+  def values(x) # <= defines functions solved
+    f = []
+    f1 = x[0]*x[0] + x[1]*x[1] - @two # f1 = x**2 + y**2 - 2 => 0
+    f2 = x[0] - x[1]                  # f2 = x    - y        => 0
+    f <<= f1
+    f <<= f2
+    f
+  end
+end
+ f = BigDecimal::limit(100)
+ f = Function.new
+ x = [f.zero,f.zero]      # Initial values
+ n = nlsolve(f,x)
+ p x
diff --git a/ext/bigdecimal/lib/pai.rb b/ext/bigdecimal/lib/pai.rb
new file mode 100644
index 0000000000..4e2dd7104a
--- /dev/null
+++ b/ext/bigdecimal/lib/pai.rb
@@ -0,0 +1,49 @@
+#!/usr/local/bin/ruby
+
+#
+# pai.rb
+#
+
+require "bigdecimal"
+#
+# Calculates 3.1415.... using J. Machin's formula.
+#
+def pai(sig) # sig: Number of significant figures
+  exp    = -sig
+  pi     = BigDecimal::new("0")
+  two    = BigDecimal::new("2")
+  m25    = BigDecimal::new("-0.04")
+  m57121 = BigDecimal::new("-57121")
+
+  u = BigDecimal::new("1")
+  k = BigDecimal::new("1")
+  w = BigDecimal::new("1")
+  t = BigDecimal::new("-80")
+  while (u.exponent >= exp) 
+    t   = t*m25
+    u,r = t.div(k,sig)
+    pi  = pi + u
+    k   = k+two
+  end
+
+  u = BigDecimal::new("1")
+  k = BigDecimal::new("1")
+  w = BigDecimal::new("1")
+  t = BigDecimal::new("956")
+  while (u.exponent >= exp )
+    t,r = t.div(m57121,sig)
+    u,r = t.div(k,sig)
+    pi  = pi + u
+    k   = k+two
+  end
+  pi
+end
+
+if $0 == __FILE__
+  if ARGV.size == 1
+    print "PAI("+ARGV[0]+"):\n"
+    p pai(ARGV[0].to_i)
+  else
+    print "TRY: ruby pai.rb 1000 \n"
+  end
+end