=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/math.c,v retrieving revision 1.1.1.1 retrieving revision 1.10 diff -u -p -r1.1.1.1 -r1.10 --- OpenXM_contrib2/asir2000/builtin/math.c 1999/12/03 07:39:07 1.1.1.1 +++ OpenXM_contrib2/asir2000/builtin/math.c 2011/12/15 16:53:26 1.10 @@ -1,10 +1,61 @@ -/* $OpenXM: OpenXM/src/asir99/builtin/math.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */ +/* + * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED + * All rights reserved. + * + * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, + * non-exclusive and royalty-free license to use, copy, modify and + * redistribute, solely for non-commercial and non-profit purposes, the + * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and + * conditions of this Agreement. For the avoidance of doubt, you acquire + * only a limited right to use the SOFTWARE hereunder, and FLL or any + * third party developer retains all rights, including but not limited to + * copyrights, in and to the SOFTWARE. + * + * (1) FLL does not grant you a license in any way for commercial + * purposes. You may use the SOFTWARE only for non-commercial and + * non-profit purposes only, such as academic, research and internal + * business use. + * (2) The SOFTWARE is protected by the Copyright Law of Japan and + * international copyright treaties. If you make copies of the SOFTWARE, + * with or without modification, as permitted hereunder, you shall affix + * to all such copies of the SOFTWARE the above copyright notice. + * (3) An explicit reference to this SOFTWARE and its copyright owner + * 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@sec.flab.fujitsu.co.jp of the detailed specification + * for such modification or the source code of the modified part of the + * SOFTWARE. + * + * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL + * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND + * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' + * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY + * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. + * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, + * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY + * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL + * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES + * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES + * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY + * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF + * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART + * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY + * 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.9 2011/08/24 07:20:09 noro Exp $ +*/ #include "ca.h" #include #include "parse.h" +#if defined(VISUAL) +#include +#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}, @@ -17,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) @@ -140,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(__x86_64) au = ((unsigned int *)&a)[1]; al = ((unsigned int *)&a)[0]; #else @@ -179,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(__x86_64) au = ((unsigned int *)&a)[1]; al = ((unsigned int *)&a)[0]; #else @@ -232,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(__x86_64) au = ((unsigned int *)&a)[1]; al = ((unsigned int *)&a)[0]; #else @@ -242,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) + int c; +#endif + + r = (Real)ARG0(arg); + if ( !r || !NUM(r) || !REAL(r) ) { + *rp = 0; + return; + } + d = ToReal(r); +#if defined(VISUAL) + 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; }