[BACK]Return to math.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Annotation of OpenXM_contrib2/asir2000/builtin/math.c, Revision 1.13

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.13    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/math.c,v 1.12 2015/08/14 13:51:54 fujimoto Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include <math.h>
                     52: #include "parse.h"
1.12      fujimoto   53: #if defined(VISUAL) || defined(__MINGW32__)
1.10      ohara      54: #include <float.h>
                     55: #endif
1.1       noro       56:
                     57: void Pdsqrt(),Pdsin(),Pdcos(),Pdtan(),Pdasin(),Pdacos(),Pdatan(),Pdlog(),Pdexp();
1.9       noro       58: void Pabs(),Pdfloor(),Pdceil(),Pdrint(),Pdisnan();
1.1       noro       59:
                     60: struct ftab math_tab[] = {
1.13    ! noro       61:   {"dsqrt",Pdsqrt,1},
        !            62:   {"dabs",Pabs,1},
        !            63:   {"dsin",Pdsin,1},
        !            64:   {"dcos",Pdcos,1},
        !            65:   {"dtan",Pdtan,1},
        !            66:   {"dlog",Pdlog,1},
        !            67:   {"dexp",Pdexp,1},
        !            68:   {"dasin",Pdasin,1},
        !            69:   {"dacos",Pdacos,1},
        !            70:   {"datan",Pdatan,1},
        !            71:   {"floor",Pdfloor,1},
        !            72:   {"dfloor",Pdfloor,1},
        !            73:   {"ceil",Pdceil,1},
        !            74:   {"dceil",Pdceil,1},
        !            75:   {"rint",Pdrint,1},
        !            76:   {"drint",Pdrint,1},
        !            77:   {"disnan",Pdisnan,1},
        !            78:   {0,0,0},
1.1       noro       79: };
                     80:
1.8       noro       81: void get_ri(Num z,double *r,double *i)
                     82: {
1.13    ! noro       83:   if ( !z ) {
        !            84:     *r = 0; *i = 0; return;
        !            85:   }
        !            86:   if ( OID(z) != O_N )
        !            87:     error("get_ri : invalid argument");
        !            88:   switch ( NID(z) ) {
        !            89:     case N_Q: case N_R: case N_B:
        !            90:       *r = ToReal(z); *i = 0;
        !            91:       break;
        !            92:     case N_C:
        !            93:       *r = ToReal(((C)z)->r);
        !            94:       *i = ToReal(((C)z)->i);
        !            95:       break;
        !            96:     default:
        !            97:       error("get_ri : invalid argument");
        !            98:       break;
        !            99:   }
1.8       noro      100: }
                    101:
1.1       noro      102: void Pabs(arg,rp)
                    103: NODE arg;
                    104: Real *rp;
                    105: {
1.13    ! noro      106:   double s,r,i;
1.1       noro      107:
1.13    ! noro      108:   if ( !ARG0(arg) ) {
        !           109:     *rp = 0; return;
        !           110:   }
        !           111:   get_ri((Num)ARG0(arg),&r,&i);
        !           112:   if ( i == 0 )
        !           113:     s = fabs(r);
        !           114:   else if ( r == 0 )
        !           115:     s = fabs(i);
        !           116:   else
        !           117:     s = sqrt(r*r+i*i);
        !           118:   MKReal(s,*rp);
1.1       noro      119: }
                    120:
                    121: void Pdsqrt(arg,rp)
                    122: NODE arg;
1.8       noro      123: Num *rp;
1.1       noro      124: {
1.13    ! noro      125:   double s,r,i,a;
        !           126:   C z;
        !           127:   Real real;
        !           128:
        !           129:   if ( !ARG0(arg) ) {
        !           130:     *rp = 0; return;
        !           131:   }
        !           132:   get_ri((Num)ARG0(arg),&r,&i);
        !           133:   if ( i == 0 )
        !           134:     if ( r > 0 ) {
        !           135:       s = sqrt(r);
        !           136:       MKReal(s,real);
        !           137:       *rp = (Num)real;
        !           138:     } else {
        !           139:       NEWC(z);
        !           140:       z->r = 0;
        !           141:       s = sqrt(-r); MKReal(s,real); z->i = (Num)real;
        !           142:       *rp = (Num)z;
        !           143:     }
        !           144:   else {
        !           145:     a = sqrt(r*r+i*i);
        !           146:     NEWC(z);
        !           147:     s = sqrt((r+a)/2); MKReal(s,real); z->r = (Num)real;
        !           148:     s = i>0?sqrt((-r+a)/2):-sqrt((-r+a)/2);
        !           149:     MKReal(s,real); z->i = (Num)real;
        !           150:     *rp = (Num)z;
        !           151:   }
1.1       noro      152: }
                    153:
                    154: void Pdsin(arg,rp)
                    155: NODE arg;
                    156: Real *rp;
                    157: {
1.13    ! noro      158:   double s;
1.1       noro      159:
1.13    ! noro      160:   s = sin(ToReal(ARG0(arg)));
        !           161:   MKReal(s,*rp);
1.1       noro      162: }
                    163:
                    164: void Pdcos(arg,rp)
                    165: NODE arg;
                    166: Real *rp;
                    167: {
1.13    ! noro      168:   double s;
1.1       noro      169:
1.13    ! noro      170:   s = cos(ToReal(ARG0(arg)));
        !           171:   MKReal(s,*rp);
1.1       noro      172: }
                    173:
                    174: void Pdtan(arg,rp)
                    175: NODE arg;
                    176: Real *rp;
                    177: {
1.13    ! noro      178:   double s;
1.1       noro      179:
1.13    ! noro      180:   s = tan(ToReal(ARG0(arg)));
        !           181:   MKReal(s,*rp);
1.1       noro      182: }
                    183:
                    184: void Pdasin(arg,rp)
                    185: NODE arg;
                    186: Real *rp;
                    187: {
1.13    ! noro      188:   double s;
1.1       noro      189:
1.13    ! noro      190:   s = asin(ToReal(ARG0(arg)));
        !           191:   MKReal(s,*rp);
1.1       noro      192: }
                    193:
                    194: void Pdacos(arg,rp)
                    195: NODE arg;
                    196: Real *rp;
                    197: {
1.13    ! noro      198:   double s;
1.1       noro      199:
1.13    ! noro      200:   s = acos(ToReal(ARG0(arg)));
        !           201:   MKReal(s,*rp);
1.1       noro      202: }
                    203:
                    204: void Pdatan(arg,rp)
                    205: NODE arg;
                    206: Real *rp;
                    207: {
1.13    ! noro      208:   double s;
1.1       noro      209:
1.13    ! noro      210:   s = atan(ToReal(ARG0(arg)));
        !           211:   MKReal(s,*rp);
1.1       noro      212: }
                    213:
                    214: void Pdlog(arg,rp)
                    215: NODE arg;
                    216: Real *rp;
                    217: {
1.13    ! noro      218:   double s;
1.1       noro      219:
1.13    ! noro      220:   s = log(ToReal(ARG0(arg)));
        !           221:   MKReal(s,*rp);
1.1       noro      222: }
                    223:
                    224: void Pdexp(arg,rp)
                    225: NODE arg;
                    226: Real *rp;
                    227: {
1.13    ! noro      228:   double s;
1.1       noro      229:
1.13    ! noro      230:   s = exp(ToReal(ARG0(arg)));
        !           231:   MKReal(s,*rp);
1.1       noro      232: }
                    233:
                    234: void Pdfloor(arg,rp)
                    235: NODE arg;
                    236: Q *rp;
                    237: {
1.13    ! noro      238:   L a;
        !           239:   unsigned int au,al;
        !           240:   int sgn;
        !           241:   Q q;
        !           242:   double d;
        !           243:
        !           244:   if ( !ARG0(arg) ) {
        !           245:     *rp = 0;
        !           246:     return;
        !           247:   }
        !           248:   d = floor(ToReal(ARG0(arg)));
        !           249:   if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
        !           250:     error("dfloor : OverFlow");
        !           251:   if ( !d ) {
        !           252:     *rp = 0;
        !           253:     return;
        !           254:   }
        !           255:   a = (L)d;
        !           256:   if ( a < 0 ) {
        !           257:     sgn = -1;
        !           258:     a = -a;
        !           259:   } else
        !           260:     sgn = 1;
1.12      fujimoto  261: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.13    ! noro      262:   au = ((unsigned int *)&a)[1];
        !           263:   al = ((unsigned int *)&a)[0];
1.1       noro      264: #else
1.13    ! noro      265:   al = ((unsigned int *)&a)[1];
        !           266:   au = ((unsigned int *)&a)[0];
1.1       noro      267: #endif
1.13    ! noro      268:   if ( au ) {
        !           269:     NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
        !           270:     PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
        !           271:   } else {
        !           272:     UTOQ(al,q); SGN(q) = sgn;
        !           273:   }
        !           274:   *rp = q;
1.1       noro      275: }
                    276:
                    277: void Pdceil(arg,rp)
                    278: NODE arg;
                    279: Q *rp;
                    280: {
1.13    ! noro      281:   L a;
        !           282:   unsigned int au,al;
        !           283:   int sgn;
        !           284:   Q q;
        !           285:   double d;
        !           286:
        !           287:   if ( !ARG0(arg) ) {
        !           288:     *rp = 0;
        !           289:     return;
        !           290:   }
        !           291:   d = ceil(ToReal(ARG0(arg)));
        !           292:   if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
        !           293:     error("dceil : OverFlow");
        !           294:   if ( !d ) {
        !           295:     *rp = 0;
        !           296:     return;
        !           297:   }
        !           298:   a = (L)d;
        !           299:   if ( a < 0 ) {
        !           300:     sgn = -1;
        !           301:     a = -a;
        !           302:   } else
        !           303:     sgn = 1;
1.12      fujimoto  304: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.13    ! noro      305:   au = ((unsigned int *)&a)[1];
        !           306:   al = ((unsigned int *)&a)[0];
1.1       noro      307: #else
1.13    ! noro      308:   al = ((unsigned int *)&a)[1];
        !           309:   au = ((unsigned int *)&a)[0];
1.1       noro      310: #endif
1.13    ! noro      311:   if ( au ) {
        !           312:     NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
        !           313:     PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
        !           314:   } else {
        !           315:     UTOQ(al,q); SGN(q) = sgn;
        !           316:   }
        !           317:   *rp = q;
1.1       noro      318: }
                    319:
                    320: void Pdrint(arg,rp)
                    321: NODE arg;
                    322: Q *rp;
                    323: {
1.13    ! noro      324:   L a;
        !           325:   unsigned int au,al;
        !           326:   int sgn;
        !           327:   Q q;
        !           328:   double d;
        !           329:
        !           330:   if ( !ARG0(arg) ) {
        !           331:     *rp = 0;
        !           332:     return;
        !           333:   }
1.12      fujimoto  334: #if defined(VISUAL) || defined(__MINGW32__)
1.13    ! noro      335:   d = ToReal(ARG0(arg));
        !           336:   if ( d > 0 )
        !           337:     d = floor(d+0.5);
        !           338:   else
        !           339:     d = ceil(d-0.5);
1.1       noro      340: #else
1.13    ! noro      341:   d = rint(ToReal(ARG0(arg)));
1.1       noro      342: #endif
1.13    ! noro      343:   if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
        !           344:     error("drint : OverFlow");
        !           345:   a = (L)d;
        !           346:   if ( a < 0 ) {
        !           347:     sgn = -1;
        !           348:     a = -a;
        !           349:   } else
        !           350:     sgn = 1;
1.12      fujimoto  351: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.13    ! noro      352:   au = ((unsigned int *)&a)[1];
        !           353:   al = ((unsigned int *)&a)[0];
1.1       noro      354: #else
1.13    ! noro      355:   al = ((unsigned int *)&a)[1];
        !           356:   au = ((unsigned int *)&a)[0];
1.1       noro      357: #endif
1.13    ! noro      358:   if ( au ) {
        !           359:     NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
        !           360:     PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
        !           361:   } else if ( al ) {
        !           362:     UTOQ(al,q); SGN(q) = sgn;
        !           363:   } else
        !           364:     q = 0;
        !           365:   *rp = q;
1.1       noro      366: }
1.9       noro      367:
                    368: void Pdisnan(NODE arg,Q *rp)
                    369: {
1.13    ! noro      370:   Real r;
        !           371:   double d;
1.12      fujimoto  372: #if defined(VISUAL) || defined(__MINGW32__)
1.10      ohara     373:     int c;
                    374: #endif
1.9       noro      375:
1.13    ! noro      376:   r = (Real)ARG0(arg);
        !           377:   if ( !r || !NUM(r) || !REAL(r) ) {
        !           378:     *rp = 0;
        !           379:     return;
        !           380:   }
        !           381:   d = ToReal(r);
1.12      fujimoto  382: #if defined(VISUAL) || defined(__MINGW32__)
1.13    ! noro      383:   c = _fpclass(d);
        !           384:   if ( c == _FPCLASS_SNAN || c == _FPCLASS_QNAN ) *rp = ONE;
        !           385:   else if ( c == _FPCLASS_PINF || c == _FPCLASS_NINF ) STOQ(2,*rp);
1.10      ohara     386: #else
1.13    ! noro      387:   if ( isnan(d) ) *rp = ONE;
        !           388:   else if ( isinf(d) ) STOQ(2,*rp);
1.10      ohara     389: #endif
1.13    ! noro      390:   else *rp = 0;
1.9       noro      391: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>