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

Annotation of OpenXM_contrib2/asir2000/engine/pari.c, Revision 1.5

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.5     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/pari.c,v 1.4 2000/12/05 01:24:52 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #if PARI
                     52: #include "base.h"
                     53: #include <math.h>
                     54: #include "genpari.h"
                     55:
                     56: void patori(GEN,Obj *);
                     57: void patori_i(GEN,N *);
                     58: void ritopa(Obj,GEN *);
                     59: void ritopa_i(N,int,GEN *);
                     60:
                     61: extern long prec;
                     62: extern int paristack;
                     63:
                     64: void risa_pari_init() {
                     65:        pari_init(paristack,2);
                     66:        prec = 4;
                     67: }
                     68:
                     69: void create_pari_variable(index)
                     70: int index;
                     71: {
                     72:        static int max_varn;
                     73:        int i;
                     74:        char name[BUFSIZ];
                     75:
                     76:        if ( index > max_varn ) {
                     77:                for ( i = max_varn+1; i <= index; i++ ) {
                     78:                        sprintf(name,"x%d",i);
                     79:                        fetch_named_var(name,0);
                     80:                }
                     81:                max_varn = index;
                     82:        }
                     83: }
                     84:
                     85: int get_lg(a)
                     86: GEN a;
                     87: {
                     88:        return lg(a);
                     89: }
                     90:
                     91: void ritopa(a,rp)
                     92: Obj a;
                     93: GEN *rp;
                     94: {
                     95:        long ltop;
                     96:        GEN pnm,z,w;
                     97:        DCP dc;
                     98:        int i,j,l,row,col;
                     99:        VL vl;
                    100:        V v;
                    101:
                    102:        if ( !a ) {
                    103:                *rp = gzero; return;
                    104:        }
                    105:        switch ( OID(a) ) {
                    106:                case O_N:
                    107:                        switch ( NID(a) ) {
                    108:                                case N_Q:
                    109:                                        ltop = avma; ritopa_i(NM((Q)a),SGN((Q)a),&pnm);
                    110:                                        if ( INT((Q)a) )
                    111:                                                *rp = pnm;
                    112:                                        else {
                    113:                                                *rp = z = cgetg(3,4); z[1] = (long)pnm;
                    114:                                                ritopa_i(DN((Q)a),1,(GEN *)&z[2]);
                    115:                                        }
                    116:                                        break;
                    117:                                case N_R:
                    118:                                        *rp = dbltor(BDY((Real)a)); break;
                    119:                                case N_B:
                    120:                                        *rp = gcopy((GEN)BDY(((BF)a))); break;
                    121:                                case N_C:
                    122:                                        z = cgetg(3,6);
                    123:                                        ritopa((Obj)((C)a)->r,(GEN *)&z[1]); ritopa((Obj)((C)a)->i,(GEN *)&z[2]);
                    124:                                        *rp = z;
                    125:                                        break;
                    126:                                default:
                    127:                                        error("ritopa : not implemented"); break;
                    128:                        }
                    129:                        break;
                    130:                case O_P:
                    131:                        l = UDEG((P)a); *rp = z = cgetg(l+3,10);
                    132:                        setsigne(z,1);
                    133:                        for ( i = 0, vl = CO, v = VR((P)a); vl->v != v;
                    134:                                vl = NEXT(vl), i++ );
                    135:                        create_pari_variable(i);
                    136:                        setvarn(z,i);
                    137:                        setlgef(z,l+3);
                    138:                        for ( i = l+2; i >= 2; i-- )
                    139:                                z[i] = (long)gzero;
                    140:                        for ( dc = DC((P)a); dc; dc = NEXT(dc) )
                    141:                                ritopa((Obj)COEF(dc),(GEN *)&z[QTOS(DEG(dc))+2]);
                    142:                        break;
                    143:                case O_VECT:
                    144:                        l = ((VECT)a)->len; z = cgetg(l+1,17);
                    145:                        for ( i = 0; i < l; i++ )
                    146:                                ritopa((Obj)BDY((VECT)a)[i],(GEN *)&z[i+1]);
                    147:                        *rp = z;
                    148:                        break;
                    149:                case O_MAT:
                    150:                        row = ((MAT)a)->row; col = ((MAT)a)->col; z = cgetg(col+1,19);
                    151:                        for ( j = 0; j < col; j++ ) {
                    152:                                w = cgetg(row+1,18);
                    153:                                for ( i = 0; i < row; i++ )
                    154:                                        ritopa((Obj)BDY((MAT)a)[i][j],(GEN *)&w[i+1]);
                    155:                                z[j+1] = (long)w;
                    156:                        }
                    157:                        *rp = z;
                    158:                        break;
                    159:                default:
                    160:                        error("ritopa : not implemented"); break;
                    161:        }
                    162: }
                    163:
                    164: void patori(a,rp)
                    165: GEN a;
                    166: Obj *rp;
                    167: {
                    168:        Q q,qnm,qdn;
                    169:        BF r;
                    170:        C c;
                    171:        VECT v;
                    172:        MAT m;
                    173:        N n,nm,dn;
                    174:        DCP dc0,dc;
                    175:        P t;
                    176:        int s,i,j,l,row,col;
                    177:        GEN b;
                    178:        VL vl;
                    179:
                    180:        if ( gcmp0(a) )
                    181:                *rp = 0;
                    182:        else {
                    183:                switch ( typ(a) ) {
                    184:                        case 1: /* integer */
                    185:                                patori_i(a,&n); NTOQ(n,(char)signe(a),q); *rp = (Obj)q;
                    186:                                break;
                    187:                        case 2: /* real */
                    188:                                NEWBF(r,lg(a)); bcopy((char *)a,(char *)BDY(r),lg(a)*sizeof(long));
                    189:                                *rp = (Obj)r;
                    190:                                break;
                    191:                        case 4: /* rational; reduced */
                    192:                                patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
                    193:                                s = signe(a[1])*signe(a[2]);
                    194:                                if ( UNIN(dn) )
                    195:                                        NTOQ(nm,s,q);
                    196:                                else
                    197:                                        NDTOQ(nm,dn,s,q);
                    198:                                *rp = (Obj)q;
                    199:                                break;
                    200:                        case 5: /* rational; not necessarily reduced */
                    201:                                patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
                    202:                                s = signe(a[1])*signe(a[2]);
                    203:                                NTOQ(nm,s,qnm); NTOQ(dn,1,qdn); divq(qnm,qdn,(Q *)rp);
                    204:                                break;
                    205:                        case 6: /* complex */
                    206:                                if ( gcmp0((GEN)a[2]) )
                    207:                                        patori((GEN)a[1],rp);
                    208:                                else {
                    209:                                        NEWC(c); patori((GEN)a[1],(Obj *)&c->r); patori((GEN)a[2],(Obj *)&c->i);
                    210:                                        *rp = (Obj)c;
                    211:                                }
                    212:                                break;
                    213:                        case 10: /* polynomial */
                    214:                                for ( i = lgef(a)-1, dc0 = 0; i >= 2; i-- )
                    215:                                        if ( !gcmp0((GEN)a[i]) ) {
                    216:                                                NEXTDC(dc0,dc);
                    217:                                                patori((GEN)a[i],(Obj *)&COEF(dc)); STOQ(i-2,DEG(dc));
                    218:                                        }
                    219:                                if ( !dc0 )
                    220:                                        *rp = 0;
                    221:                                else {
                    222:                                        /* assuming that CO has not changed. */
                    223:                                        for ( s = varn(a), i = 0, vl = CO; i != s;
                    224:                                                i++, vl = NEXT(vl) );
                    225:                                        NEXT(dc) = 0; MKP(vl->v,dc0,t); *rp = (Obj)t;
                    226:                                }
                    227:                                break;
                    228:                        case 17: /* row vector */
                    229:                        case 18: /* column vector */
                    230:                                l = lg(a)-1; MKVECT(v,l);
                    231:                                for ( i = 0; i < l; i++ )
                    232:                                        patori((GEN)a[i+1],(Obj *)&BDY(v)[i]);
                    233:                                *rp = (Obj)v;
                    234:                                break;
                    235:                        case 19: /* matrix */
                    236:                                col = lg(a)-1; row = lg(a[1])-1; MKMAT(m,row,col);
                    237:                                for ( j = 0; j < col; j++ )
                    238:                                        for ( i = 0, b = (GEN)a[j+1]; i < row; i++ )
                    239:                                                patori((GEN)b[i+1],(Obj *)&BDY(m)[i][j]);
                    240:                                *rp = (Obj)m;
                    241:                                break;
                    242:                        default:
                    243:                                error("patori : not implemented");
                    244:                                break;
                    245:                }
                    246:        }
                    247: }
                    248:
                    249: #if defined(LONG_IS_32BIT)
                    250: void ritopa_i(a,s,rp)
                    251: N a;
                    252: int s;
                    253: GEN *rp;
                    254: {
                    255:        int j,l;
                    256:        unsigned int *b;
                    257:        GEN z;
                    258:
                    259:        l = PL(a); b = (unsigned int *)BD(a);
                    260:        z = cgeti(l+2);
                    261:        bzero((char *)&z[2],l*sizeof(int));
                    262:        for ( j = 0; j < l; j++ )
                    263:                z[l-j+1] = b[j];
                    264:        s = s>0?1:-1;
                    265:        setsigne(z,s);
                    266:        setlgefint(z,lg(z));
                    267:        *rp = z;
                    268: }
                    269:
                    270: void patori_i(g,rp)
                    271: GEN g;
                    272: N *rp;
                    273: {
                    274:        int j,l;
                    275:        unsigned int *a,*b;
                    276:        N z;
                    277:
                    278:        l = lgef(g)-2;
                    279:        a = (unsigned int *)g;
                    280:        *rp = z = NALLOC(l); PL(z) = l;
                    281:        for ( j = 0, b = (unsigned int *)BD(z); j < l; j++ )
                    282:                b[l-j-1] = ((unsigned int *)g)[j+2];
                    283: }
                    284: #endif
                    285:
                    286: #if defined(LONG_IS_64BIT)
                    287: void ritopa_i(a,s,rp)
                    288: N a;
                    289: int s;
                    290: GEN *rp;
                    291: {
                    292:        int j,l,words;
                    293:        unsigned int *b;
                    294:        GEN z;
                    295:
                    296:        l = PL(a); b = BD(a);
                    297:        words = (l+1)/2;
                    298:        z = cgeti(words+2);
                    299:        bzero((char *)&z[2],words*sizeof(long));
                    300:        for ( j = 0; j < words; j++ )
                    301:                z[words-j+1] = ((unsigned long)b[2*j])
                    302:                        |(((unsigned long)(2*j+1<l?b[2*j+1]:0))<<32);
                    303:        s = s>0?1:-1;
                    304:        setsigne(z,s);
                    305:        setlgefint(z,lg(z));
                    306:        *rp = z;
                    307: }
                    308:
                    309: void patori_i(g,rp)
                    310: GEN g;
                    311: N *rp;
                    312: {
                    313:        int j,l,words;
                    314:        unsigned long t;
                    315:        unsigned long *a;
                    316:        unsigned int *b;
                    317:        N z;
                    318:
                    319:        words = lgef(g)-2;
                    320:        l = 2*words;
                    321:        a = (unsigned long *)g;
                    322:        *rp = z = NALLOC(l); PL(z) = l;
                    323:        for ( j = 0, b = BD(z); j < words; j++ ) {
                    324:                t = a[words+1-j];
                    325:                b[2*j] = t&0xffffffff;
                    326:                b[2*j+1] = t>>32;
                    327:        }
                    328:        PL(z) = b[l-1] ? l : l-1;
                    329: }
                    330: #endif
                    331:
                    332: void strtobf(s,rp)
                    333: char *s;
                    334: BF *rp;
                    335: {
                    336:        GEN z;
                    337:
                    338:        z = lisexpr(s); patori(z,(Obj *)rp); cgiv(z);
                    339: }
                    340: #endif

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