[BACK]Return to ox_pari.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / ox_pari

Annotation of OpenXM/src/ox_pari/ox_pari.c, Revision 1.3

1.3     ! noro        1: /*     $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.2 2015/08/06 09:15:32 noro Exp $      */
1.1       noro        2:
                      3: #include <stdio.h>
                      4: #include <stdlib.h>
                      5: #include <string.h>
1.2       noro        6: #include "pari/pari.h"
1.3     ! noro        7: #include "pari/paripriv.h"
1.1       noro        8: #include "gmp.h"
1.2       noro        9: #include "gmp-impl.h"
1.3     ! noro       10: #include "mpfr.h"
1.1       noro       11: #include "ox_toolkit.h"
                     12: OXFILE *fd_rw;
                     13:
1.3     ! noro       14: #define MPFR_PREC(x)      ((x)->_mpfr_prec)
        !            15: #define MPFR_EXP(x)       ((x)->_mpfr_exp)
        !            16: #define MPFR_MANT(x)      ((x)->_mpfr_d)
        !            17: #define MPFR_LAST_LIMB(x) ((MPFR_PREC (x) - 1) / GMP_NUMB_BITS)
        !            18: #define MPFR_LIMB_SIZE(x) (MPFR_LAST_LIMB (x) + 1)
        !            19:
1.1       noro       20: static int stack_size = 0;
                     21: static int stack_pointer = 0;
                     22: static cmo **stack = NULL;
                     23: extern int debug_print;
1.2       noro       24: long paristack=10000000;
1.1       noro       25:
                     26: void init_pari(void);
1.2       noro       27: cmo *GEN_to_cmo(GEN z);
                     28: cmo_zz *GEN_to_cmo_zz(GEN z);
1.3     ! noro       29: cmo_bf *GEN_to_cmo_bf(GEN z);
1.2       noro       30: cmo_list *GEN_to_cmo_list(GEN z);
                     31: GEN cmo_to_GEN(cmo *c);
                     32: GEN cmo_zz_to_GEN(cmo_zz *c);
1.3     ! noro       33: GEN cmo_bf_to_GEN(cmo_bf *c);
1.1       noro       34:
                     35: #define INIT_S_SIZE 2048
                     36: #define EXT_S_SIZE  2048
                     37:
1.3     ! noro       38: void *gc_realloc(void *p,size_t osize,size_t nsize)
        !            39: {
        !            40:   return (void *)GC_realloc(p,nsize);
        !            41: }
        !            42:
        !            43: void gc_free(void *p,size_t size)
        !            44: {
        !            45:   GC_free(p);
        !            46: }
        !            47:
        !            48: void init_gc()
        !            49: {
        !            50:        GC_INIT();
        !            51:   mp_set_memory_functions(GC_malloc,gc_realloc,gc_free);
        !            52: }
        !            53:
1.1       noro       54: void init_pari()
                     55: {
1.2       noro       56:   pari_init(paristack,2);
1.1       noro       57: }
                     58:
                     59: int initialize_stack()
                     60: {
                     61:        stack_pointer = 0;
                     62:        stack_size = INIT_S_SIZE;
                     63:        stack = MALLOC(stack_size*sizeof(cmo*));
                     64:        return 0;
                     65: }
                     66:
                     67: static int extend_stack()
                     68: {
                     69:        int size2 = stack_size + EXT_S_SIZE;
                     70:        cmo **stack2 = MALLOC(size2*sizeof(cmo*));
                     71:        memcpy(stack2, stack, stack_size*sizeof(cmo *));
                     72:        free(stack);
                     73:        stack = stack2;
                     74:        stack_size = size2;
                     75:        return 0;
                     76: }
                     77:
                     78: int push(cmo* m)
                     79: {
                     80:        stack[stack_pointer] = m;
                     81:        stack_pointer++;
                     82:        if(stack_pointer >= stack_size) {
                     83:                extend_stack();
                     84:        }
                     85:        return 0;
                     86: }
                     87:
                     88: cmo* pop()
                     89: {
                     90:        if(stack_pointer > 0) {
                     91:                stack_pointer--;
                     92:                return stack[stack_pointer];
                     93:        }
                     94:        return new_cmo_null();
                     95: }
                     96:
                     97: void pops(int n)
                     98: {
                     99:        stack_pointer -= n;
                    100:        if(stack_pointer < 0) {
                    101:                stack_pointer = 0;
                    102:        }
                    103: }
                    104:
                    105: #define OX_PARI_VERSION 20150731
                    106: #define ID_STRING  "2015/07/31 15:00:00"
                    107:
                    108: int sm_mathcap()
                    109: {
1.2       noro      110:        mathcap_init(OX_PARI_VERSION, ID_STRING, "ox_pari", NULL, NULL);
1.1       noro      111:        push((cmo*)oxf_cmo_mathcap(fd_rw));
                    112:        return 0;
                    113: }
                    114:
                    115: int sm_popCMO()
                    116: {
                    117:        cmo* m = pop();
                    118:
                    119:        if(m != NULL) {
                    120:                send_ox_cmo(fd_rw, m);
                    121:                return 0;
                    122:        }
                    123:        return SM_popCMO;
                    124: }
                    125:
                    126: cmo_error2 *make_error2(int code)
                    127: {
1.2       noro      128:        return (cmo_error2 *) new_cmo_int32(code);
1.1       noro      129: }
                    130:
                    131: int get_i()
                    132: {
                    133:        cmo *c = pop();
                    134:        if(c->tag == CMO_INT32) {
                    135:                return ((cmo_int32 *)c)->i;
                    136:        }else if(c->tag == CMO_ZZ) {
                    137:                return mpz_get_si(((cmo_zz *)c)->mpz);
                    138:        }
                    139:        make_error2(-1);
                    140:        return 0;
                    141: }
                    142:
                    143: char *get_str()
                    144: {
                    145:        cmo *c = pop();
                    146:        if(c->tag == CMO_STRING) {
                    147:                return ((cmo_string *)c)->s;
                    148:        }
                    149:        make_error2(-1);
                    150:        return "";
                    151: }
                    152:
                    153: int cmo2int(cmo *c)
                    154: {
                    155:        if(c->tag == CMO_INT32) {
                    156:                return ((cmo_int32 *)c)->i;
                    157:        }else if(c->tag == CMO_ZZ) {
                    158:                return mpz_get_si(((cmo_zz *)c)->mpz);
                    159:        } else if(c->tag == CMO_NULL){
                    160:                return 0;
                    161:        }
                    162:
                    163:        return 0;
                    164: }
                    165:
1.2       noro      166: GEN cmo_zz_to_GEN(cmo_zz *c)
                    167: {
                    168:   mpz_ptr mpz;
                    169:   GEN z;
                    170:   long *ptr;
                    171:   int j,sgn,len;
                    172:
                    173:   mpz = c->mpz;
                    174:   sgn = mpz_sgn(mpz);
                    175:   len = ABSIZ(mpz);
                    176:   ptr = (long *)PTR(mpz);
                    177:   z = cgeti(len+2);
                    178:   for ( j = 0; j < len; j++ )
                    179:     z[len-j+1] = ptr[j];
                    180:   setsigne(z,sgn);
                    181:   setlgefint(z,lg(z));
                    182:   return z;
                    183: }
                    184:
1.3     ! noro      185: GEN cmo_bf_to_GEN(cmo_bf *c)
        !           186: {
        !           187:   mpfr_ptr mpfr;
        !           188:   GEN z;
        !           189:   int sgn,len,j;
        !           190:   long exp;
        !           191:   long *ptr;
        !           192:
        !           193:   mpfr = c->mpfr;
        !           194:   sgn = MPFR_SIGN(mpfr);
        !           195:   exp = MPFR_EXP(mpfr)-1;
        !           196:   len = MPFR_LIMB_SIZE(mpfr);
        !           197:   ptr = (long *)MPFR_MANT(mpfr);
        !           198:   z = cgetr(len+2);
        !           199:   for ( j = 0; j < len; j++ )
        !           200:     z[len-j+1] = ptr[j];
        !           201:   z[1] = evalsigne(sgn)|evalexpo(exp);
        !           202:   setsigne(z,sgn);
        !           203:   return z;
        !           204: }
        !           205:
1.2       noro      206: cmo_zz *GEN_to_cmo_zz(GEN z)
                    207: {
                    208:   cmo_zz *c;
                    209:
                    210:   c = new_cmo_zz();
                    211:   mpz_import(c->mpz,lgef(z)-2,1,sizeof(long),0,0,&z[2]);
                    212:   if ( signe(z) < 0 )
                    213:     mpz_neg(c->mpz,c->mpz);
                    214:   return c;
                    215: }
                    216:
1.3     ! noro      217: cmo_bf *GEN_to_cmo_bf(GEN z)
        !           218: {
        !           219:   cmo_bf *c;
        !           220:   int len,prec,j;
        !           221:   long *ptr;
        !           222:
        !           223:   c = new_cmo_bf();
        !           224:   len = lg(z)-2;
        !           225:   prec = len*sizeof(long)*8;
        !           226:   mpfr_init2(c->mpfr,prec);
        !           227:   ptr = (long *)MPFR_MANT(c->mpfr);
        !           228:   for ( j = 0; j < len; j++ )
        !           229:     ptr[j] = z[len-j+1];
        !           230:   MPFR_EXP(c->mpfr) = (long long)(expo(z)+1);
        !           231:   MPFR_SIGN(c->mpfr) = gsigne(z);
        !           232:   return c;
        !           233: }
        !           234:
        !           235:
1.2       noro      236: cmo_list *GEN_to_cmo_list(GEN z)
                    237: {
                    238:   cmo_list *c;
                    239:   cmo *ob;
                    240:   int i,len;
                    241:
                    242:   c = new_cmo_list();
                    243:   len = lg(z)-1;
                    244:   for ( i = 1; i <= len; i++ ) {
                    245:     ob = GEN_to_cmo((GEN)z[i]);
                    246:     c = list_append(c,ob);
                    247:   }
                    248:   return c;
                    249: }
                    250:
                    251:
                    252: GEN cmo_to_GEN(cmo *c)
                    253: {
                    254:   switch ( c->tag ) {
                    255:   case CMO_ZERO:
1.3     ! noro      256:     return gen_0;
1.2       noro      257:   case CMO_ZZ: /* int */
                    258:     return cmo_zz_to_GEN((cmo_zz *)c);
1.3     ! noro      259:   case CMO_BIGFLOAT: /* bigfloat */
        !           260:     return cmo_bf_to_GEN((cmo_bf *)c);
1.2       noro      261:   default:
                    262:     return 0;
                    263:   }
                    264: }
                    265:
                    266: cmo *GEN_to_cmo(GEN z)
                    267: {
                    268:   if ( gcmp0(z) )
                    269:     return new_cmo_zero();
                    270:   switch ( typ(z) ) {
                    271:   case 1: /* int */
                    272:     return (cmo *)GEN_to_cmo_zz(z);
1.3     ! noro      273:   case 2: /* bigfloat */
        !           274:     return (cmo *)GEN_to_cmo_bf(z);
1.2       noro      275:   case 17: case 18: /* vector */
                    276:     return (cmo *)GEN_to_cmo_list(z);
                    277:   case 19: /* matrix */
                    278:     return (cmo *)GEN_to_cmo_list(shallowtrans(z));
                    279:   default:
                    280:     return (cmo *)make_error2(typ(z));
                    281:   }
                    282: }
                    283:
1.3     ! noro      284: struct parif {
        !           285:   char *name;
        !           286:   int type;
        !           287:   GEN (*f)();
        !           288: } parif_tab[] = {
        !           289:   {"erfc",1,gerfc},
        !           290:   {"factor",1,Z_factor},
        !           291:   {"isqrt",1,racine},
        !           292:   {"nextprime",1,nextprime},
        !           293:   {"det",1,det},
        !           294:   {"allocatemem",0,(GEN (*)())allocatemoremem},
        !           295: };
1.2       noro      296:
                    297: #define PARI_MAX_AC 64
                    298:
1.3     ! noro      299: struct parif *search_parif(char *name)
        !           300: {
        !           301:   int tablen,i;
        !           302:
        !           303:   tablen = sizeof(parif_tab)/sizeof(struct parif);
        !           304:   for ( i = 0; i < tablen; i++ ) {
        !           305:     if ( !strcmp(parif_tab[i].name,name) )
        !           306:       return &parif_tab[i];
        !           307:   }
        !           308:   return 0;
        !           309: }
        !           310:
1.1       noro      311: int sm_executeFunction()
                    312: {
1.3     ! noro      313:   long ltop,lbot;
1.2       noro      314:   int ac,i;
                    315:   cmo_int32 *c;
                    316:   cmo *av[PARI_MAX_AC];
                    317:   cmo *ret;
                    318:   GEN z,m;
1.3     ! noro      319:   struct parif *parif;
1.2       noro      320:
1.3     ! noro      321:   if ( setjmp(GP_DATA->env) ) {
        !           322:                printf("sm_executeFunction : an error occured.\n");fflush(stdout);
        !           323:                push((cmo*)make_error2(0));
        !           324:                return -1;
        !           325:   }
1.1       noro      326:        cmo_string *func = (cmo_string *)pop();
                    327:        if(func->tag != CMO_STRING) {
                    328:                printf("sm_executeFunction : func->tag is not CMO_STRING");fflush(stdout);
                    329:                push((cmo*)make_error2(0));
                    330:                return -1;
                    331:        }
                    332:
1.2       noro      333:        c = (cmo_int32 *)pop();
                    334:   ac = c->i;
                    335:   if ( ac > PARI_MAX_AC ) {
                    336:                push((cmo*)make_error2(0));
                    337:                return -1;
                    338:   }
                    339:   for ( i = 0; i < ac; i++ ) {
                    340:     av[i] = (cmo *)pop();
                    341:     fprintf(stderr,"arg%d:",i);
                    342:     print_cmo(av[i]);
                    343:     fprintf(stderr,"\n");
                    344:   }
1.3     ! noro      345:        if( strcmp( func->s, "exit" ) == 0 )
        !           346:                exit(0);
        !           347:
        !           348:   parif =search_parif(func->s);
        !           349:   if ( !parif ) {
        !           350:                push((cmo*)make_error2(0));
        !           351:                return -1;
        !           352:  } else if ( parif->type == 0 ) {
        !           353:     /* one long int variable */
        !           354:     int a = cmo_to_int(av[0]);
        !           355:     a = (int)(*parif->f)(a);
        !           356:     ret = (cmo *)new_cmo_int32(a);
1.2       noro      357:     push(ret);
                    358:                return 0;
1.3     ! noro      359:   } else if ( parif->type == 1 ) {
        !           360:     /* one variable possibly with prec */
        !           361:     unsigned long prec;
        !           362:
        !           363:     ltop = avma;
1.2       noro      364:     z = cmo_to_GEN(av[0]);
1.3     ! noro      365:     if ( ac == 2 ) {
        !           366:       prec = cmo_to_int(av[1])*3.32193/32+3;
        !           367:     } else
        !           368:       prec = precreal;
        !           369:     m = (*parif->f)(z,prec);
        !           370:     lbot = avma;
1.2       noro      371:     ret = GEN_to_cmo(m);
1.3     ! noro      372:    // gerepile(ltop,lbot,0);
1.2       noro      373:     push(ret);
1.1       noro      374:                return 0;
1.3     ! noro      375:   } else {
1.1       noro      376:                push((cmo*)make_error2(0));
                    377:                return -1;
1.3     ! noro      378:   }
1.1       noro      379: }
                    380:
                    381: int receive_and_execute_sm_command()
                    382: {
                    383:        int code = receive_int32(fd_rw);
                    384:        switch(code) {
                    385:        case SM_popCMO:
                    386:                sm_popCMO();
                    387:                break;
                    388:        case SM_executeFunction:
                    389:                sm_executeFunction();
                    390:                break;
                    391:        case SM_mathcap:
                    392:                sm_mathcap();
                    393:                break;
                    394:        case SM_setMathCap:
                    395:                pop();
                    396:                break;
                    397:        default:
                    398:                printf("receive_and_execute_sm_command : code=%d\n",code);fflush(stdout);
                    399:                break;
                    400:        }
                    401:        return 0;
                    402: }
                    403:
                    404: int receive()
                    405: {
                    406:        int tag;
                    407:
                    408:        tag = receive_ox_tag(fd_rw);
                    409:        switch(tag) {
                    410:        case OX_DATA:
1.2       noro      411:                printf("receive : ox_data %d\n",tag);fflush(stdout);
1.1       noro      412:                push(receive_cmo(fd_rw));
                    413:                break;
                    414:        case OX_COMMAND:
1.2       noro      415:                printf("receive : ox_command %d\n",tag);fflush(stdout);
1.1       noro      416:                receive_and_execute_sm_command();
                    417:                break;
                    418:        default:
                    419:                printf("receive : tag=%d\n",tag);fflush(stdout);
                    420:        }
                    421:        return 0;
                    422: }
                    423:
                    424: int main()
                    425: {
1.3     ! noro      426:   init_gc();
1.1       noro      427:        ox_stderr_init(stderr);
                    428:        initialize_stack();
                    429:        init_pari();
                    430:
                    431:        fprintf(stderr,"ox_pari\n");
                    432:
                    433:        fd_rw = oxf_open(3);
                    434:        oxf_determine_byteorder_server(fd_rw);
                    435:
                    436:        while(1){
                    437:                receive();
                    438:        }
                    439: }

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