===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/math.c,v
retrieving revision 1.2
retrieving revision 1.11
diff -u -p -r1.2 -r1.11
--- OpenXM_contrib2/asir2000/builtin/math.c	2000/08/21 08:31:20	1.2
+++ OpenXM_contrib2/asir2000/builtin/math.c	2015/08/06 10:01:52	1.11
@@ -23,7 +23,7 @@
  * shall be made on your publication or presentation in any form of the
  * results obtained by use of the SOFTWARE.
  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
- * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
+ * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
  * for such modification or the source code of the modified part of the
  * SOFTWARE.
  * 
@@ -45,14 +45,17 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/builtin/math.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/math.c,v 1.10 2011/12/15 16:53:26 ohara Exp $ 
 */
 #include "ca.h"
 #include <math.h>
 #include "parse.h"
+#if defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__)
+#include <float.h>
+#endif
 
 void Pdsqrt(),Pdsin(),Pdcos(),Pdtan(),Pdasin(),Pdacos(),Pdatan(),Pdlog(),Pdexp();
-void Pabs(),Pdfloor(),Pdceil(),Pdrint();
+void Pabs(),Pdfloor(),Pdceil(),Pdrint(),Pdisnan();
 
 struct ftab math_tab[] = {
 	{"dsqrt",Pdsqrt,1},
@@ -65,30 +68,87 @@ struct ftab math_tab[] = {
 	{"dasin",Pdasin,1},
 	{"dacos",Pdacos,1},
 	{"datan",Pdatan,1},
+	{"floor",Pdfloor,1},
 	{"dfloor",Pdfloor,1},
+	{"ceil",Pdceil,1},
 	{"dceil",Pdceil,1},
+	{"rint",Pdrint,1},
 	{"drint",Pdrint,1},
+	{"disnan",Pdisnan,1},
 	{0,0,0},
 };
 
+void get_ri(Num z,double *r,double *i)
+{
+	if ( !z ) {
+		*r = 0; *i = 0; return;
+	}
+	if ( OID(z) != O_N )
+		error("get_ri : invalid argument");
+	switch ( NID(z) ) {
+		case N_Q: case N_R: case N_B:
+			*r = ToReal(z); *i = 0;
+			break;
+		case N_C:
+			*r = ToReal(((C)z)->r);
+			*i = ToReal(((C)z)->i);
+			break;
+		default:
+			error("get_ri : invalid argument");
+			break;
+	}
+}
+
 void Pabs(arg,rp)
 NODE arg;
 Real *rp;
 {
-	double s;
+	double s,r,i;
 
-	s = fabs(ToReal(ARG0(arg)));
+	if ( !ARG0(arg) ) {
+		*rp = 0; return;
+	}
+	get_ri((Num)ARG0(arg),&r,&i);
+	if ( i == 0 )
+		s = fabs(r);
+	else if ( r == 0 )
+		s = fabs(i);
+	else
+		s = sqrt(r*r+i*i);
 	MKReal(s,*rp);
 }
 
 void Pdsqrt(arg,rp)
 NODE arg;
-Real *rp;
+Num *rp;
 {
-	double s;
+	double s,r,i,a;
+	C z;
+	Real real;
 
-	s = sqrt(ToReal(ARG0(arg)));
-	MKReal(s,*rp);
+	if ( !ARG0(arg) ) {
+		*rp = 0; return;
+	}
+	get_ri((Num)ARG0(arg),&r,&i);
+	if ( i == 0 )
+		if ( r > 0 ) {
+			s = sqrt(r);
+			MKReal(s,real);
+			*rp = (Num)real;
+		} else {
+			NEWC(z);
+			z->r = 0;
+			s = sqrt(-r); MKReal(s,real); z->i = (Num)real;
+			*rp = (Num)z;
+		}
+	else {
+		a = sqrt(r*r+i*i);
+		NEWC(z);
+		s = sqrt((r+a)/2); MKReal(s,real); z->r = (Num)real;
+		s = i>0?sqrt((-r+a)/2):-sqrt((-r+a)/2);
+		MKReal(s,real); z->i = (Num)real;
+		*rp = (Num)z;
+	}
 }
 
 void Pdsin(arg,rp)
@@ -188,13 +248,17 @@ Q *rp;
 	d = floor(ToReal(ARG0(arg)));
 	if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
 		error("dfloor : OverFlow");
+	if ( !d ) {
+		*rp = 0;
+		return;
+	}
 	a = (L)d;
 	if ( a < 0 ) {
 		sgn = -1;
 		a = -a;
 	} else
 		sgn = 1;
-#if defined(i386) || defined(__alpha) || defined(VISUAL)
+#if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__) || defined(__x86_64)
 	au = ((unsigned int *)&a)[1];
 	al = ((unsigned int *)&a)[0];
 #else
@@ -227,13 +291,17 @@ Q *rp;
 	d = ceil(ToReal(ARG0(arg)));
 	if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
 		error("dceil : OverFlow");
+	if ( !d ) {
+		*rp = 0;
+		return;
+	}
 	a = (L)d;
 	if ( a < 0 ) {
 		sgn = -1;
 		a = -a;
 	} else
 		sgn = 1;
-#if defined(i386) || defined(__alpha) || defined(VISUAL)
+#if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__) || defined(__x86_64)
 	au = ((unsigned int *)&a)[1];
 	al = ((unsigned int *)&a)[0];
 #else
@@ -263,7 +331,7 @@ Q *rp;
 		*rp = 0;
 		return;
 	}
-#if defined(VISUAL)
+#if defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__)
 	d = ToReal(ARG0(arg));
 	if ( d > 0 )
 		d = floor(d+0.5);
@@ -280,7 +348,7 @@ Q *rp;
 		a = -a;
 	} else
 		sgn = 1;
-#if defined(i386) || defined(__alpha) || defined(VISUAL)
+#if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__) || defined(__x86_64)
 	au = ((unsigned int *)&a)[1];
 	al = ((unsigned int *)&a)[0];
 #else
@@ -290,8 +358,34 @@ Q *rp;
 	if ( au ) {
 		NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
 		PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
-	} else {
+	} else if ( al ) {
 		UTOQ(al,q); SGN(q) = sgn;
-	}
+	} else
+		q = 0;
 	*rp = q;
+}
+
+void Pdisnan(NODE arg,Q *rp)
+{
+	Real r;
+	double d;
+#if defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__)
+    int c;
+#endif
+
+	r = (Real)ARG0(arg);
+	if ( !r || !NUM(r) || !REAL(r) ) {
+		*rp = 0;
+		return;
+	}
+	d = ToReal(r);
+#if defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__)
+	c = _fpclass(d);
+	if ( c == _FPCLASS_SNAN || c == _FPCLASS_QNAN ) *rp = ONE;
+	else if ( c == _FPCLASS_PINF || c == _FPCLASS_NINF ) STOQ(2,*rp);
+#else
+	if ( isnan(d) ) *rp = ONE;
+	else if ( isinf(d) ) STOQ(2,*rp);
+#endif
+	else *rp = 0;
 }