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

1.4     ! noro        1: /*  $OpenXM: OpenXM/src/ox_pari/ox_pari.c,v 1.3 2015/08/17 05:18:35 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: {
1.4     ! noro       50:   GC_INIT();
1.3       noro       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: {
1.4     ! noro       61:   stack_pointer = 0;
        !            62:    stack_size = INIT_S_SIZE;
        !            63:   stack = MALLOC(stack_size*sizeof(cmo*));
        !            64:   return 0;
1.1       noro       65: }
                     66:
                     67: static int extend_stack()
                     68: {
1.4     ! noro       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;
1.1       noro       76: }
                     77:
                     78: int push(cmo* m)
                     79: {
1.4     ! noro       80:   stack[stack_pointer] = m;
        !            81:   stack_pointer++;
        !            82:   if(stack_pointer >= stack_size) {
        !            83:     extend_stack();
        !            84:   }
        !            85:   return 0;
1.1       noro       86: }
                     87:
                     88: cmo* pop()
                     89: {
1.4     ! noro       90:   if(stack_pointer > 0) {
        !            91:     stack_pointer--;
        !            92:     return stack[stack_pointer];
        !            93:   }
        !            94:   return new_cmo_null();
1.1       noro       95: }
                     96:
                     97: void pops(int n)
                     98: {
1.4     ! noro       99:   stack_pointer -= n;
        !           100:   if(stack_pointer < 0) {
        !           101:     stack_pointer = 0;
        !           102:   }
1.1       noro      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.4     ! noro      110:   mathcap_init(OX_PARI_VERSION, ID_STRING, "ox_pari", NULL, NULL);
        !           111:   push((cmo*)oxf_cmo_mathcap(fd_rw));
        !           112:   return 0;
1.1       noro      113: }
                    114:
                    115: int sm_popCMO()
                    116: {
1.4     ! noro      117:   cmo* m = pop();
1.1       noro      118:
1.4     ! noro      119:   if(m != NULL) {
        !           120:     send_ox_cmo(fd_rw, m);
        !           121:     return 0;
        !           122:   }
        !           123:   return SM_popCMO;
1.1       noro      124: }
                    125:
                    126: cmo_error2 *make_error2(int code)
                    127: {
1.4     ! noro      128:   return (cmo_error2 *) new_cmo_int32(code);
1.1       noro      129: }
                    130:
                    131: int get_i()
                    132: {
1.4     ! noro      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;
1.1       noro      141: }
                    142:
                    143: char *get_str()
                    144: {
1.4     ! noro      145:   cmo *c = pop();
        !           146:   if(c->tag == CMO_STRING) {
        !           147:     return ((cmo_string *)c)->s;
        !           148:   }
        !           149:   make_error2(-1);
        !           150:   return "";
1.1       noro      151: }
                    152:
                    153: int cmo2int(cmo *c)
                    154: {
1.4     ! noro      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:   }
1.1       noro      162:
1.4     ! noro      163:   return 0;
1.1       noro      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;
1.4     ! noro      286:   GEN (*f)();
1.3       noro      287:   int type;
                    288: } parif_tab[] = {
1.4     ! noro      289: /* (ulong)allocatemoremem(ulong) */
        !           290:   {"allocatemem",(GEN (*)())allocatemoremem,0},
        !           291: /* num/num */
        !           292:   {"abs",gabs,1},
        !           293:   {"erfc",gerfc,1},
        !           294:   {"arg",garg,1},
        !           295:   {"isqrt",racine,1},
        !           296:   {"gamma",ggamma,1},
        !           297:   {"zeta",gzeta,1},
        !           298:   {"floor",gfloor,1},
        !           299:   {"frac",gfrac,1},
        !           300:   {"imag",gimag,1},
        !           301:   {"conj",gconj,1},
        !           302:   {"ceil",gceil,1},
        !           303:   {"isprime",gisprime,2},
        !           304:   {"bigomega",gbigomega,1},
        !           305:   {"denom",denom,1},
        !           306:   {"numer",numer,1},
        !           307:   {"lngamma",glngamma,1},
        !           308:   {"logagm",glogagm,1},
        !           309:   {"classno",classno,1},
        !           310:   {"classno2",classno2,1},
        !           311:   {"dilog",dilog,1},
        !           312:   {"disc",discsr,1},
        !           313:   {"discf",discf,1},
        !           314:   {"nextprime",nextprime,1},
        !           315:   {"eintg1",eint1,1},
        !           316:   {"eta",eta,1},
        !           317:   {"issqfree",gissquarefree,1},
        !           318:   {"issquare",gcarreparfait,1},
        !           319:   {"gamh",ggamd,1},
        !           320:   {"hclassno",classno3,1},
        !           321:
        !           322:   /* num/array */
        !           323:   {"binary",binaire,1},
        !           324:   {"factorint",factorint,2},
        !           325:   {"factor",Z_factor,1},
        !           326:   {"cf",gcf,1},
        !           327:   {"divisors",divisors,1},
        !           328:   {"smallfact",smallfact,1},
        !           329:
        !           330:   /* mat/mat */
        !           331:   {"adj",adj,1},
        !           332:   {"lll",lll,1},
        !           333:   {"lllgen",lllgen,1},
        !           334:   {"lllgram",lllgram,1},
        !           335:   {"lllgramgen",lllgramgen,1},
        !           336:   {"lllgramint",lllgramint,1},
        !           337:   {"lllgramkerim",lllgramkerim,1},
        !           338:   {"lllgramkerimgen",lllgramkerimgen,1},
        !           339:   {"lllint",lllint,1},
        !           340:   {"lllkerim",lllkerim,1},
        !           341:   {"lllkerimgen",lllkerimgen,1},
        !           342:   {"trans",gtrans,1},
        !           343:   {"eigen",eigen,1},
        !           344:   {"hermite",hnf,1},
        !           345:   {"mat",gtomat,1},
        !           346:   {"matrixqz2",matrixqz2,1},
        !           347:   {"matrixqz3",matrixqz3,1},
        !           348:   {"hess",hess,1},
        !           349:
        !           350:   /* mat/poly */
        !           351:   {"det",det,1},
        !           352:   {"det2",det2,1},
        !           353:
        !           354:   /* poly/poly */
        !           355:   {"centerlift",centerlift,1},
        !           356:   {"content",content,1},
        !           357:
        !           358:   /* poly/array */
        !           359:   {"galois",galois,1},
        !           360:   {"roots",roots,1},
        !           361:
1.3       noro      362: };
1.2       noro      363:
                    364: #define PARI_MAX_AC 64
                    365:
1.3       noro      366: struct parif *search_parif(char *name)
                    367: {
                    368:   int tablen,i;
                    369:
                    370:   tablen = sizeof(parif_tab)/sizeof(struct parif);
                    371:   for ( i = 0; i < tablen; i++ ) {
                    372:     if ( !strcmp(parif_tab[i].name,name) )
                    373:       return &parif_tab[i];
                    374:   }
                    375:   return 0;
                    376: }
                    377:
1.1       noro      378: int sm_executeFunction()
                    379: {
1.3       noro      380:   long ltop,lbot;
1.2       noro      381:   int ac,i;
                    382:   cmo_int32 *c;
                    383:   cmo *av[PARI_MAX_AC];
                    384:   cmo *ret;
                    385:   GEN z,m;
1.3       noro      386:   struct parif *parif;
1.2       noro      387:
1.3       noro      388:   if ( setjmp(GP_DATA->env) ) {
1.4     ! noro      389:     printf("sm_executeFunction : an error occured.\n");fflush(stdout);
        !           390:     push((cmo*)make_error2(0));
        !           391:     return -1;
        !           392:   }
        !           393:   cmo_string *func = (cmo_string *)pop();
        !           394:   if(func->tag != CMO_STRING) {
        !           395:     printf("sm_executeFunction : func->tag is not CMO_STRING");fflush(stdout);
        !           396:     push((cmo*)make_error2(0));
        !           397:     return -1;
        !           398:   }
1.1       noro      399:
1.4     ! noro      400:   c = (cmo_int32 *)pop();
1.2       noro      401:   ac = c->i;
                    402:   if ( ac > PARI_MAX_AC ) {
1.4     ! noro      403:     push((cmo*)make_error2(0));
        !           404:     return -1;
1.2       noro      405:   }
                    406:   for ( i = 0; i < ac; i++ ) {
                    407:     av[i] = (cmo *)pop();
                    408:     fprintf(stderr,"arg%d:",i);
                    409:     print_cmo(av[i]);
                    410:     fprintf(stderr,"\n");
                    411:   }
1.4     ! noro      412:   if( strcmp( func->s, "exit" ) == 0 )
        !           413:     exit(0);
1.3       noro      414:
                    415:   parif =search_parif(func->s);
                    416:   if ( !parif ) {
1.4     ! noro      417:     push((cmo*)make_error2(0));
        !           418:     return -1;
1.3       noro      419:  } else if ( parif->type == 0 ) {
                    420:     /* one long int variable */
                    421:     int a = cmo_to_int(av[0]);
1.4     ! noro      422:     a = (int)(parif->f)(a);
1.3       noro      423:     ret = (cmo *)new_cmo_int32(a);
1.2       noro      424:     push(ret);
1.4     ! noro      425:     return 0;
1.3       noro      426:   } else if ( parif->type == 1 ) {
                    427:     /* one variable possibly with prec */
                    428:     unsigned long prec;
                    429:
                    430:     ltop = avma;
1.2       noro      431:     z = cmo_to_GEN(av[0]);
1.3       noro      432:     if ( ac == 2 ) {
                    433:       prec = cmo_to_int(av[1])*3.32193/32+3;
                    434:     } else
                    435:       prec = precreal;
                    436:     m = (*parif->f)(z,prec);
                    437:     lbot = avma;
1.2       noro      438:     ret = GEN_to_cmo(m);
1.4     ! noro      439: //    gerepile(ltop,lbot,0);
1.2       noro      440:     push(ret);
1.4     ! noro      441:     return 0;
1.3       noro      442:   } else {
1.4     ! noro      443:     push((cmo*)make_error2(0));
        !           444:     return -1;
1.3       noro      445:   }
1.1       noro      446: }
                    447:
                    448: int receive_and_execute_sm_command()
                    449: {
1.4     ! noro      450:   int code = receive_int32(fd_rw);
        !           451:   switch(code) {
        !           452:   case SM_popCMO:
        !           453:     sm_popCMO();
        !           454:     break;
        !           455:   case SM_executeFunction:
        !           456:     sm_executeFunction();
        !           457:     break;
        !           458:   case SM_mathcap:
        !           459:     sm_mathcap();
        !           460:     break;
        !           461:   case SM_setMathCap:
        !           462:     pop();
        !           463:     break;
        !           464:   default:
        !           465:     printf("receive_and_execute_sm_command : code=%d\n",code);fflush(stdout);
        !           466:     break;
        !           467:   }
        !           468:   return 0;
1.1       noro      469: }
                    470:
                    471: int receive()
                    472: {
1.4     ! noro      473:   int tag;
1.1       noro      474:
1.4     ! noro      475:   tag = receive_ox_tag(fd_rw);
        !           476:   switch(tag) {
        !           477:   case OX_DATA:
        !           478:     printf("receive : ox_data %d\n",tag);fflush(stdout);
        !           479:     push(receive_cmo(fd_rw));
        !           480:     break;
        !           481:   case OX_COMMAND:
        !           482:     printf("receive : ox_command %d\n",tag);fflush(stdout);
        !           483:     receive_and_execute_sm_command();
        !           484:     break;
        !           485:   default:
        !           486:     printf("receive : tag=%d\n",tag);fflush(stdout);
        !           487:   }
        !           488:   return 0;
1.1       noro      489: }
                    490:
                    491: int main()
                    492: {
1.3       noro      493:   init_gc();
1.4     ! noro      494:   ox_stderr_init(stderr);
        !           495:   initialize_stack();
        !           496:   init_pari();
        !           497:
        !           498:   fprintf(stderr,"ox_pari\n");
        !           499:
        !           500:   fd_rw = oxf_open(3);
        !           501:   oxf_determine_byteorder_server(fd_rw);
        !           502:
        !           503:   while(1){
        !           504:     receive();
        !           505:   }
1.1       noro      506: }

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