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

Annotation of OpenXM_contrib2/asir2000/builtin/parif.c, Revision 1.15

1.3       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.4       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3       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.15    ! ohara      48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/parif.c,v 1.14 2003/02/14 22:29:07 ohara Exp $
1.3       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "parse.h"
                     52:
1.14      ohara      53: #if defined(PARI)
1.1       noro       54: #include "genpari.h"
                     55:
                     56: extern long prec;
                     57:
                     58: void patori(GEN,Obj *);
                     59: void patori_i(GEN,N *);
                     60: void ritopa(Obj,GEN *);
                     61: void ritopa_i(N,int,GEN *);
                     62:
1.11      noro       63: void Ptodouble();
1.1       noro       64: void Peval(),Psetprec(),p_pi(),p_e(),p_mul(),p_gcd();
1.5       noro       65: void asir_cgiv(GEN);
1.1       noro       66:
1.7       saito      67: #if defined(INTERVAL) || 1
                     68: void Psetprecword();
                     69: #endif
                     70:
1.1       noro       71: struct ftab pari_tab[] = {
1.7       saito      72:        {"eval",Peval,-2},
                     73:        {"setprec",Psetprec,-1},
1.11      noro       74:        {"todouble",Ptodouble,1},
1.7       saito      75: #if defined(INTERVAL) || 1
                     76:        {"setprecword",Psetprecword,-1},
                     77: #endif
                     78:        {0,0,0},
1.1       noro       79: };
                     80:
                     81: #define MKPREC(a,i,b) (argc(a)==(i)?mkprec(QTOS((Q)(b))):prec)
                     82:
                     83: #define CALLPARI1(f,a,p,r)\
1.5       noro       84: ritopa((Obj)a,&_pt1_); _pt2_ = f(_pt1_,p); patori(_pt2_,r); asir_cgiv(_pt2_); asir_cgiv(_pt1_)
1.1       noro       85: #define CALLPARI2(f,a,b,p,r)\
1.5       noro       86: ritopa((Obj)a,&_pt1_); ritopa((Obj)b,&_pt2_); _pt3_ = f(_pt1_,_pt2_,p); patori(_pt3_,r); asir_cgiv(_pt3_); asir_cgiv(_pt2_); asir_cgiv(_pt1_)
1.1       noro       87:
                     88: #define PARIF1P(f,pf)\
                     89: void f(NODE,Obj *);\
                     90: void f(ar,rp) NODE ar; Obj *rp;\
                     91: { GEN _pt1_,_pt2_; CALLPARI1(pf,ARG0(ar),MKPREC(ar,2,ARG1(ar)),rp); }
                     92: #define PARIF2P(f,pf)\
                     93: void f(NODE,Obj *);\
                     94: void f(ar,rp) NODE ar; Obj *rp;\
                     95: { GEN _pt1_,_pt2_,_pt3_; CALLPARI2(pf,ARG0(ar),ARG1(ar),MKPREC(ar,3,ARG2(ar)),rp); }
                     96:
1.7       saito      97: #if defined(INTERVAL)
                     98: #define        PREC_CONV       pariK1
1.15    ! ohara      99: #elif SIZEOF_LONG == 4
1.1       noro      100: #define PREC_CONV              0.103810253
1.15    ! ohara     101: #elif SIZEOF_LONG == 8
1.1       noro      102: #define PREC_CONV              0.051905126
                    103: #endif
                    104:
1.5       noro      105: /* XXX : we should be more careful when we free PARI pointers. */
                    106:
                    107: void asir_cgiv(ptr)
                    108: GEN ptr;
                    109: {
                    110:        if ( ptr != gzero && ptr != gun
                    111:        && ptr != gdeux && ptr != ghalf
                    112:        && ptr != polvar && ptr != gi )
                    113:                cgiv(ptr);
                    114: }
                    115:
1.1       noro      116: mkprec(p)
                    117: int p;
                    118: {
1.10      noro      119:        if ( p <= 0 )
                    120:                p = 1;
                    121:        return (int)(p*PREC_CONV+3);
1.11      noro      122: }
                    123:
                    124: void Ptodouble(arg,rp)
                    125: NODE arg;
                    126: Num *rp;
                    127: {
                    128:        double r,i;
                    129:        Real real,imag;
                    130:        Num num;
                    131:
                    132:        asir_assert(ARG0(arg),O_N,"todouble");
                    133:        num = (Num)ARG0(arg);
                    134:        if ( !num ) {
                    135:                *rp = 0;
                    136:                return;
                    137:        }
                    138:        switch ( NID(num) ) {
                    139:                case N_R: case N_Q: case N_B:
                    140:                        r = ToReal(num);
                    141:                        MKReal(r,real);
                    142:                        *rp = (Num)real;
                    143:                        break;
                    144:                case N_C:
                    145:                        r = ToReal(((C)num)->r);
                    146:                        i = ToReal(((C)num)->i);
                    147:                        MKReal(r,real);
                    148:                        MKReal(i,imag);
                    149:                        reimtocplx((Num)real,(Num)imag,rp);
                    150:                        break;
                    151:                default:
                    152:                        *rp = num;
                    153:                        break;
                    154:        }
1.1       noro      155: }
                    156:
                    157: void Peval(arg,rp)
                    158: NODE arg;
                    159: Obj *rp;
                    160: {
                    161:        asir_assert(ARG0(arg),O_R,"eval");
                    162:        evalr(CO,(Obj)ARG0(arg),argc(arg)==2?QTOS((Q)ARG1(arg)):0,rp);
                    163: }
                    164:
                    165: void Psetprec(arg,rp)
                    166: NODE arg;
                    167: Obj *rp;
                    168: {
                    169:        int p;
                    170:        Q q;
                    171:
1.7       saito     172: #if defined(INTERVAL) || 1
                    173:        p = (int)((prec-2)/PREC_CONV); STOQ(p,q); *rp = (Obj)q;
                    174:        if ( arg ) {
                    175:                asir_assert(ARG0(arg),O_N,"setprec");
                    176:                p = QTOS((Q)ARG0(arg));
                    177:                if ( p > 0 )
                    178:                        prec = (long)(p*PREC_CONV+3);
                    179:        }
                    180: #else
1.1       noro      181:        p = (int)((prec-3)/PREC_CONV); STOQ(p,q); *rp = (Obj)q;
                    182:        if ( arg ) {
                    183:                asir_assert(ARG0(arg),O_N,"setprec");
                    184:                prec = mkprec(QTOS((Q)ARG0(arg)));
                    185:        }
1.7       saito     186: #endif
1.1       noro      187: }
1.7       saito     188:
                    189: #if defined(INTERVAL) || 1
                    190: void Psetprecword(arg,rp)
                    191: NODE arg;
                    192: Obj *rp;
                    193: {
                    194:        int p;
                    195:        Q q;
                    196:
                    197:        p = (int)((prec-2)); STOQ(p,q); *rp = (Obj)q;
                    198:        if ( arg ) {
                    199:                asir_assert(ARG0(arg),O_N,"setprecword");
                    200:                p = QTOS((Q)ARG0(arg));
                    201:                if ( p > 0 ) {
                    202:                        prec = p + 2;
                    203:                }
                    204:        }
                    205: }
                    206: #endif
1.1       noro      207:
                    208: void p_pi(arg,rp)
                    209: NODE arg;
                    210: Obj *rp;
                    211: {
                    212:        GEN x;
                    213:
                    214:        x = mppi(MKPREC(arg,1,ARG0(arg)));
1.5       noro      215:        patori(x,rp); asir_cgiv(x);
1.1       noro      216: }
                    217:
                    218: void p_e(arg,rp)
                    219: NODE arg;
                    220: Obj *rp;
                    221: {
                    222:        GEN x;
                    223:
1.5       noro      224:        x = gexp(gun,MKPREC(arg,1,ARG0(arg))); patori(x,rp); asir_cgiv(x);
1.1       noro      225: }
                    226:
                    227: void p_mul(a,b,r)
                    228: Obj a,b,*r;
                    229: {
                    230:        GEN p1,p2,p3;
                    231:
                    232:        ritopa((Obj)a,&p1); ritopa((Obj)b,&p2);
                    233:        p3 = mulii(p1,p2);
1.5       noro      234:        patori(p3,r); asir_cgiv(p3); asir_cgiv(p2); asir_cgiv(p1);
1.1       noro      235: }
                    236:
                    237: void p_gcd(a,b,r)
                    238: N a,b,*r;
                    239: {
                    240:        GEN p1,p2,p3;
                    241:
                    242:        ritopa_i(a,1,&p1); ritopa_i(b,1,&p2);
                    243:        p3 = mppgcd(p1,p2);
1.5       noro      244:        patori_i(p3,r); asir_cgiv(p3); asir_cgiv(p2); asir_cgiv(p1);
1.1       noro      245: }
                    246:
                    247: PARIF1P(p_sin,gsin) PARIF1P(p_cos,gcos) PARIF1P(p_tan,gtan)
                    248: PARIF1P(p_asin,gasin) PARIF1P(p_acos,gacos) PARIF1P(p_atan,gatan)
                    249: PARIF1P(p_sinh,gsh) PARIF1P(p_cosh,gch) PARIF1P(p_tanh,gth)
                    250: PARIF1P(p_asinh,gash) PARIF1P(p_acosh,gach) PARIF1P(p_atanh,gath)
                    251: PARIF1P(p_exp,gexp) PARIF1P(p_log,glog)
                    252: PARIF1P(p_dilog,dilog) PARIF1P(p_erf,gerfc)
                    253: PARIF1P(p_eigen,eigen) PARIF1P(p_roots,roots)
                    254:
                    255: PARIF2P(p_pow,gpui)
                    256:
                    257: pointer evalparif(f,arg)
                    258: FUNC f;
                    259: NODE arg;
                    260: {
                    261:        GEN a,v;
                    262:        long ltop,lbot;
                    263:        pointer r;
1.9       noro      264:        int ac,opt,intarg,ret;
1.1       noro      265:        char buf[BUFSIZ];
1.9       noro      266:        Q q;
1.8       noro      267:        GEN (*dmy)();
1.1       noro      268:
                    269:        if ( !f->f.binf ) {
                    270:                sprintf(buf,"pari : %s undefined.",f->name);
                    271:                error(buf);
1.10      noro      272:                /* NOTREACHED */
                    273:                return 0;
1.1       noro      274:        }
                    275:        switch ( f->type ) {
1.9       noro      276:                case 0: /* in/out : integer */
                    277:                        ac = argc(arg);
                    278:                        if ( ac > 2 ) {
                    279:                                fprintf(stderr,"argument mismatch in %s()\n",NAME(f));
                    280:                                error("");
1.10      noro      281:                                /* NOTREACHED */
                    282:                                return 0;
1.9       noro      283:                        }
                    284:                        intarg = !ac ? 0 : QTOS((Q)ARG0(arg));
                    285:                        dmy = (GEN (*)())f->f.binf;
1.10      noro      286:                        ret = (int)(*dmy)(intarg);
1.9       noro      287:                        STOQ(ret,q);
                    288:                        return (pointer)q;
                    289:
1.1       noro      290:                case 1:
                    291:                        ac = argc(arg);
                    292:                        if ( !ac || ( ac > 2 ) ) {
                    293:                                fprintf(stderr,"argument mismatch in %s()\n",NAME(f));
                    294:                                error("");
1.10      noro      295:                                /* NOTREACHED */
                    296:                                return 0;
1.1       noro      297:                        }
                    298:                        ltop = avma;
                    299:                        ritopa((Obj)ARG0(arg),&a);
                    300:                        dmy = (GEN (*)())f->f.binf;
                    301:                        v = (*dmy)(a,MKPREC(arg,2,ARG1(arg)));
                    302:                        lbot = avma;
                    303:                        patori(v,(Obj *)&r); gerepile(ltop,lbot,0);
                    304:                        return r;
1.8       noro      305:
                    306:                case 2:
                    307:                        ac = argc(arg);
                    308:                        if ( !ac || ( ac > 2 ) ) {
                    309:                                fprintf(stderr,"argument mismatch in %s()\n",NAME(f));
                    310:                                error("");
1.10      noro      311:                                /* NOTREACHED */
                    312:                                return 0;
1.8       noro      313:                        }
                    314:                        if ( ac == 1 )
                    315:                                opt = 0;
                    316:                        else
                    317:                                opt = QTOS((Q)ARG1(arg));
                    318:                        ltop = avma;
                    319:                        ritopa((Obj)ARG0(arg),&a);
                    320:                        dmy = (GEN (*)())f->f.binf;
                    321:                        v = (*dmy)(a,opt);
                    322:                        lbot = avma;
                    323:                        patori(v,(Obj *)&r); gerepile(ltop,lbot,0);
                    324:                        return r;
                    325:
1.1       noro      326:                default:
                    327:                        error("evalparif : not implemented yet.");
1.10      noro      328:                        /* NOTREACHED */
                    329:                        return 0;
1.1       noro      330:        }
                    331: }
                    332:
                    333: struct pariftab {
                    334:        char *name;
                    335:        GEN (*f)();
                    336:        int type;
                    337: };
                    338:
1.8       noro      339: /*
                    340:  * type = 1 => argc = 1, second arg = precision
                    341:  * type = 2 => argc = 1, second arg = optional (long int)
                    342:  *
                    343:  */
                    344:
1.1       noro      345: struct pariftab pariftab[] = {
1.10      noro      346: {"allocatemem",(GEN(*)())allocatemoremem,0},
1.1       noro      347: {"abs",(GEN (*)())gabs,1},
                    348: {"adj",adj,1},
                    349: {"arg",garg,1},
                    350: {"bigomega",gbigomega,1},
                    351: {"binary",binaire,1},
                    352: {"ceil",gceil,1},
                    353: {"centerlift",centerlift,1},
                    354: {"cf",gcf,1},
                    355: {"classno",classno,1},
                    356: {"classno2",classno2,1},
                    357: {"conj",gconj,1},
                    358: {"content",content,1},
                    359: {"denom",denom,1},
                    360: {"det",det,1},
                    361: {"det2",det2,1},
                    362: {"dilog",dilog,1},
                    363: {"disc",discsr,1},
                    364: {"discf",discf,1},
                    365: {"divisors",divisors,1},
                    366: {"eigen",eigen,1},
                    367: {"eintg1",eint1,1},
                    368: {"erfc",gerfc,1},
                    369: {"eta",eta,1},
                    370: {"floor",gfloor,1},
                    371: {"frac",gfrac,1},
                    372: {"galois",galois,1},
                    373: {"galoisconj",galoisconj,1},
                    374: {"gamh",ggamd,1},
                    375: {"gamma",ggamma,1},
                    376: {"hclassno",classno3,1},
                    377: {"hermite",hnf,1},
                    378: {"hess",hess,1},
                    379: {"imag",gimag,1},
                    380: {"image",image,1},
                    381: {"image2",image2,1},
                    382: {"indexrank",indexrank,1},
                    383: {"indsort",indexsort,1},
                    384: {"initalg",initalg,1},
                    385: {"isfund",gisfundamental,1},
1.13      noro      386: {"isprime",gisprime,2},
1.1       noro      387: {"ispsp",gispsp,1},
                    388: {"isqrt",racine,1},
                    389: {"issqfree",gissquarefree,1},
                    390: {"issquare",gcarreparfait,1},
                    391: {"jacobi",jacobi,1},
                    392: {"jell",jell,1},
                    393: {"ker",ker,1},
                    394: {"keri",keri,1},
                    395: {"kerint",kerint,1},
                    396: {"kerintg1",kerint1,1},
                    397: {"length",(GEN(*)())glength,1},
                    398: {"lexsort",lexsort,1},
                    399: {"lift",lift,1},
                    400: {"lindep",lindep,1},
                    401: {"lll",lll,1},
                    402: {"lllgen",lllgen,1},
                    403: {"lllgram",lllgram,1},
                    404: {"lllgramgen",lllgramgen,1},
                    405: {"lllgramint",lllgramint,1},
                    406: {"lllgramkerim",lllgramkerim,1},
                    407: {"lllgramkerimgen",lllgramkerimgen,1},
                    408: {"lllint",lllint,1},
                    409: {"lllkerim",lllkerim,1},
                    410: {"lllkerimgen",lllkerimgen,1},
                    411: {"lngamma",glngamma,1},
                    412: {"logagm",glogagm,1},
                    413: {"mat",gtomat,1},
                    414: {"matrixqz2",matrixqz2,1},
                    415: {"matrixqz3",matrixqz3,1},
                    416: {"matsize",matsize,1},
                    417: {"modreverse",polymodrecip,1},
                    418: {"mu",gmu,1},
                    419: {"nextprime",nextprime,1},
                    420: {"norm",gnorm,1},
                    421: {"norml2",gnorml2,1},
                    422: {"numdiv",numbdiv,1},
                    423: {"numer",numer,1},
                    424: {"omega",gomega,1},
                    425: {"order",order,1},
                    426: {"ordred",ordred,1},
                    427: {"phi",phi,1},
                    428: {"pnqn",pnqn,1},
                    429: {"polred",polred,1},
                    430: {"polred2",polred2,1},
                    431: {"primroot",gener,1},
                    432: {"psi",gpsi,1},
                    433: {"quadgen",quadgen    ,1},
                    434: {"quadpoly",quadpoly    ,1},
                    435: {"real",greal,1},
                    436: {"recip",polrecip       ,1},
                    437: {"redreal",redreal       ,1},
                    438: {"regula",regula  ,1},
                    439: {"reorder",reorder  ,1},
                    440: {"reverse",recip  ,1},
                    441: {"rhoreal",rhoreal       ,1},
                    442: {"roots",roots,1},
                    443: {"round",ground,1},
                    444: {"sigma",sumdiv,1},
                    445: {"signat",signat,1},
                    446: {"simplify",simplify,1},
                    447: {"smalldiscf",smalldiscf,1},
                    448: {"smallfact",smallfact,1},
                    449: {"smallpolred",smallpolred,1},
                    450: {"smallpolred2",smallpolred2,1},
                    451: {"smith",smith,1},
                    452: {"smith2",smith2,1},
                    453: {"sort",sort,1},
                    454: {"sqr",gsqr,1},
                    455: {"sqred",sqred,1},
                    456: {"sqrt",gsqrt,1},
                    457: {"supplement",suppl,1},
                    458: {"trace",gtrace,1},
                    459: {"trans",gtrans,1},
                    460: {"trunc",gtrunc,1},
                    461: {"unit",fundunit,1},
                    462: {"vec",gtovec,1},
                    463: {"wf",wf,1},
                    464: {"wf2",wf2,1},
                    465: {"zeta",gzeta,1},
1.9       noro      466: {"factor",factor,1},
1.8       noro      467: {"factorint",factorint,2},
1.1       noro      468: {0,0,0},
                    469: };
                    470:
                    471: void parif_init() {
                    472:        int i;
                    473:
                    474:        for ( i = 0, parif = 0; pariftab[i].name; i++ )
                    475:                 appendparif(&parif,pariftab[i].name, (int (*)())pariftab[i].f,pariftab[i].type);
                    476: }
                    477: #else /* PARI */
                    478:
                    479: struct ftab pari_tab[] = {
                    480:        {0,0,0},
                    481: };
                    482:
                    483: void parif_init() {}
                    484:
1.2       noro      485: pointer evalparif(f,arg)
                    486: FUNC f;
                    487: NODE arg;
                    488: {
1.1       noro      489:        error("evalparif : PARI is not combined.");
                    490: }
                    491: #endif /*PARI */

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