[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.9

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.9     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/math.c,v 1.8 2003/12/26 05:47:37 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include <math.h>
                     52: #include "parse.h"
                     53:
                     54: void Pdsqrt(),Pdsin(),Pdcos(),Pdtan(),Pdasin(),Pdacos(),Pdatan(),Pdlog(),Pdexp();
1.9     ! noro       55: void Pabs(),Pdfloor(),Pdceil(),Pdrint(),Pdisnan();
1.1       noro       56:
                     57: struct ftab math_tab[] = {
                     58:        {"dsqrt",Pdsqrt,1},
                     59:        {"dabs",Pabs,1},
                     60:        {"dsin",Pdsin,1},
                     61:        {"dcos",Pdcos,1},
                     62:        {"dtan",Pdtan,1},
                     63:        {"dlog",Pdlog,1},
                     64:        {"dexp",Pdexp,1},
                     65:        {"dasin",Pdasin,1},
                     66:        {"dacos",Pdacos,1},
                     67:        {"datan",Pdatan,1},
1.5       noro       68:        {"floor",Pdfloor,1},
1.1       noro       69:        {"dfloor",Pdfloor,1},
1.5       noro       70:        {"ceil",Pdceil,1},
1.1       noro       71:        {"dceil",Pdceil,1},
1.5       noro       72:        {"rint",Pdrint,1},
1.1       noro       73:        {"drint",Pdrint,1},
1.9     ! noro       74:        {"disnan",Pdisnan,1},
1.1       noro       75:        {0,0,0},
                     76: };
                     77:
1.8       noro       78: void get_ri(Num z,double *r,double *i)
                     79: {
                     80:        if ( !z ) {
                     81:                *r = 0; *i = 0; return;
                     82:        }
                     83:        if ( OID(z) != O_N )
                     84:                error("get_ri : invalid argument");
                     85:        switch ( NID(z) ) {
                     86:                case N_Q: case N_R: case N_B:
                     87:                        *r = ToReal(z); *i = 0;
                     88:                        break;
                     89:                case N_C:
                     90:                        *r = ToReal(((C)z)->r);
                     91:                        *i = ToReal(((C)z)->i);
                     92:                        break;
                     93:                default:
                     94:                        error("get_ri : invalid argument");
                     95:                        break;
                     96:        }
                     97: }
                     98:
1.1       noro       99: void Pabs(arg,rp)
                    100: NODE arg;
                    101: Real *rp;
                    102: {
1.8       noro      103:        double s,r,i;
1.1       noro      104:
1.8       noro      105:        if ( !ARG0(arg) ) {
                    106:                *rp = 0; return;
                    107:        }
                    108:        get_ri((Num)ARG0(arg),&r,&i);
                    109:        if ( i == 0 )
                    110:                s = fabs(r);
                    111:        else if ( r == 0 )
                    112:                s = fabs(i);
                    113:        else
                    114:                s = sqrt(r*r+i*i);
1.1       noro      115:        MKReal(s,*rp);
                    116: }
                    117:
                    118: void Pdsqrt(arg,rp)
                    119: NODE arg;
1.8       noro      120: Num *rp;
1.1       noro      121: {
1.8       noro      122:        double s,r,i,a;
                    123:        C z;
                    124:        Real real;
                    125:
                    126:        if ( !ARG0(arg) ) {
                    127:                *rp = 0; return;
                    128:        }
                    129:        get_ri((Num)ARG0(arg),&r,&i);
                    130:        if ( i == 0 )
                    131:                if ( r > 0 ) {
                    132:                        s = sqrt(r);
                    133:                        MKReal(s,real);
                    134:                        *rp = (Num)real;
                    135:                } else {
                    136:                        NEWC(z);
                    137:                        z->r = 0;
                    138:                        s = sqrt(-r); MKReal(s,real); z->i = (Num)real;
                    139:                        *rp = (Num)z;
                    140:                }
                    141:        else {
                    142:                a = sqrt(r*r+i*i);
                    143:                NEWC(z);
                    144:                s = sqrt((r+a)/2); MKReal(s,real); z->r = (Num)real;
                    145:                s = i>0?sqrt((-r+a)/2):-sqrt((-r+a)/2);
                    146:                MKReal(s,real); z->i = (Num)real;
                    147:                *rp = (Num)z;
                    148:        }
1.1       noro      149: }
                    150:
                    151: void Pdsin(arg,rp)
                    152: NODE arg;
                    153: Real *rp;
                    154: {
                    155:        double s;
                    156:
                    157:        s = sin(ToReal(ARG0(arg)));
                    158:        MKReal(s,*rp);
                    159: }
                    160:
                    161: void Pdcos(arg,rp)
                    162: NODE arg;
                    163: Real *rp;
                    164: {
                    165:        double s;
                    166:
                    167:        s = cos(ToReal(ARG0(arg)));
                    168:        MKReal(s,*rp);
                    169: }
                    170:
                    171: void Pdtan(arg,rp)
                    172: NODE arg;
                    173: Real *rp;
                    174: {
                    175:        double s;
                    176:
                    177:        s = tan(ToReal(ARG0(arg)));
                    178:        MKReal(s,*rp);
                    179: }
                    180:
                    181: void Pdasin(arg,rp)
                    182: NODE arg;
                    183: Real *rp;
                    184: {
                    185:        double s;
                    186:
                    187:        s = asin(ToReal(ARG0(arg)));
                    188:        MKReal(s,*rp);
                    189: }
                    190:
                    191: void Pdacos(arg,rp)
                    192: NODE arg;
                    193: Real *rp;
                    194: {
                    195:        double s;
                    196:
                    197:        s = acos(ToReal(ARG0(arg)));
                    198:        MKReal(s,*rp);
                    199: }
                    200:
                    201: void Pdatan(arg,rp)
                    202: NODE arg;
                    203: Real *rp;
                    204: {
                    205:        double s;
                    206:
                    207:        s = atan(ToReal(ARG0(arg)));
                    208:        MKReal(s,*rp);
                    209: }
                    210:
                    211: void Pdlog(arg,rp)
                    212: NODE arg;
                    213: Real *rp;
                    214: {
                    215:        double s;
                    216:
                    217:        s = log(ToReal(ARG0(arg)));
                    218:        MKReal(s,*rp);
                    219: }
                    220:
                    221: void Pdexp(arg,rp)
                    222: NODE arg;
                    223: Real *rp;
                    224: {
                    225:        double s;
                    226:
                    227:        s = exp(ToReal(ARG0(arg)));
                    228:        MKReal(s,*rp);
                    229: }
                    230:
                    231: void Pdfloor(arg,rp)
                    232: NODE arg;
                    233: Q *rp;
                    234: {
                    235:        L a;
                    236:        unsigned int au,al;
                    237:        int sgn;
                    238:        Q q;
                    239:        double d;
                    240:
                    241:        if ( !ARG0(arg) ) {
                    242:                *rp = 0;
                    243:                return;
                    244:        }
                    245:        d = floor(ToReal(ARG0(arg)));
                    246:        if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
                    247:                error("dfloor : OverFlow");
1.4       noro      248:        if ( !d ) {
                    249:                *rp = 0;
                    250:                return;
                    251:        }
1.1       noro      252:        a = (L)d;
                    253:        if ( a < 0 ) {
                    254:                sgn = -1;
                    255:                a = -a;
                    256:        } else
                    257:                sgn = 1;
1.7       noro      258: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__x86_64)
1.1       noro      259:        au = ((unsigned int *)&a)[1];
                    260:        al = ((unsigned int *)&a)[0];
                    261: #else
                    262:        al = ((unsigned int *)&a)[1];
                    263:        au = ((unsigned int *)&a)[0];
                    264: #endif
                    265:        if ( au ) {
                    266:                NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
                    267:                PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
                    268:        } else {
                    269:                UTOQ(al,q); SGN(q) = sgn;
                    270:        }
                    271:        *rp = q;
                    272: }
                    273:
                    274: void Pdceil(arg,rp)
                    275: NODE arg;
                    276: Q *rp;
                    277: {
                    278:        L a;
                    279:        unsigned int au,al;
                    280:        int sgn;
                    281:        Q q;
                    282:        double d;
                    283:
                    284:        if ( !ARG0(arg) ) {
                    285:                *rp = 0;
                    286:                return;
                    287:        }
                    288:        d = ceil(ToReal(ARG0(arg)));
                    289:        if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
                    290:                error("dceil : OverFlow");
1.4       noro      291:        if ( !d ) {
                    292:                *rp = 0;
                    293:                return;
                    294:        }
1.1       noro      295:        a = (L)d;
                    296:        if ( a < 0 ) {
                    297:                sgn = -1;
                    298:                a = -a;
                    299:        } else
                    300:                sgn = 1;
1.7       noro      301: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__x86_64)
1.1       noro      302:        au = ((unsigned int *)&a)[1];
                    303:        al = ((unsigned int *)&a)[0];
                    304: #else
                    305:        al = ((unsigned int *)&a)[1];
                    306:        au = ((unsigned int *)&a)[0];
                    307: #endif
                    308:        if ( au ) {
                    309:                NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
                    310:                PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
                    311:        } else {
                    312:                UTOQ(al,q); SGN(q) = sgn;
                    313:        }
                    314:        *rp = q;
                    315: }
                    316:
                    317: void Pdrint(arg,rp)
                    318: NODE arg;
                    319: Q *rp;
                    320: {
                    321:        L a;
                    322:        unsigned int au,al;
                    323:        int sgn;
                    324:        Q q;
                    325:        double d;
                    326:
                    327:        if ( !ARG0(arg) ) {
                    328:                *rp = 0;
                    329:                return;
                    330:        }
                    331: #if defined(VISUAL)
                    332:        d = ToReal(ARG0(arg));
                    333:        if ( d > 0 )
                    334:                d = floor(d+0.5);
                    335:        else
                    336:                d = ceil(d-0.5);
                    337: #else
                    338:        d = rint(ToReal(ARG0(arg)));
                    339: #endif
                    340:        if ( d < -9.223372036854775808e18 || d >= 9.223372036854775808e18 )
                    341:                error("drint : OverFlow");
                    342:        a = (L)d;
                    343:        if ( a < 0 ) {
                    344:                sgn = -1;
                    345:                a = -a;
                    346:        } else
                    347:                sgn = 1;
1.7       noro      348: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__x86_64)
1.1       noro      349:        au = ((unsigned int *)&a)[1];
                    350:        al = ((unsigned int *)&a)[0];
                    351: #else
                    352:        al = ((unsigned int *)&a)[1];
                    353:        au = ((unsigned int *)&a)[0];
                    354: #endif
                    355:        if ( au ) {
                    356:                NEWQ(q); SGN(q) = sgn; NM(q)=NALLOC(2); DN(q)=0;
                    357:                PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
1.6       noro      358:        } else if ( al ) {
1.1       noro      359:                UTOQ(al,q); SGN(q) = sgn;
1.6       noro      360:        } else
                    361:                q = 0;
1.1       noro      362:        *rp = q;
                    363: }
1.9     ! noro      364:
        !           365: void Pdisnan(NODE arg,Q *rp)
        !           366: {
        !           367:        Real r;
        !           368:        double d;
        !           369:
        !           370:        r = (Real)ARG0(arg);
        !           371:        if ( !r || !NUM(r) || !REAL(r) ) {
        !           372:                *rp = 0;
        !           373:                return;
        !           374:        }
        !           375:        d = ToReal(r);
        !           376:        if ( isnan(d) ) *rp = ONE;
        !           377:        else if ( isinf(d) ) STOQ(2,*rp);
        !           378:        else *rp = 0;
        !           379: }

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