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

Annotation of OpenXM_contrib2/asir2018/builtin/parif.c, Revision 1.5

1.5     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2018/builtin/parif.c,v 1.4 2020/08/26 06:40:36 noro Exp $ */
1.1       noro        2: #include "ca.h"
                      3: #include "parse.h"
                      4: #include "ox.h"
                      5:
                      6: Q ox_pari_stream;
                      7: int ox_pari_stream_initialized = 0;
                      8: int ox_get_pari_result = 0;
                      9: P ox_pari_starting_function = 0;
                     10:
                     11: typedef void (*mpfr_func)(NODE,Obj *);
                     12:
                     13: void Pmpfr_ai();
                     14: void Pmpfr_eint(), Pmpfr_erf(),Pmpfr_li2();
                     15: void Pmpfr_zeta();
                     16: void Pmpfr_j0(), Pmpfr_j1();
                     17: void Pmpfr_y0(), Pmpfr_y1();
                     18: void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma();
                     19: void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil();
                     20:
                     21: void Pox_shutdown(NODE arg,Q *rp);
                     22: void Pox_launch_nox(NODE arg,Obj *rp);
1.4       noro       23: void Pox_launch(NODE arg,Obj *rp);
1.1       noro       24: void Pox_push_cmo(NODE arg,Obj *rp);
                     25: void Pox_pop_cmo(NODE arg,Obj *rp);
                     26: void Pox_execute_function(NODE arg,Obj *rp);
                     27:
1.4       noro       28: int debug_pari;
                     29:
1.1       noro       30: struct mpfr_tab_rec {
                     31:   char *name;
                     32:   mpfr_func func;
                     33: } mpfr_tab[] = {
                     34:   {"ai",Pmpfr_ai},
                     35:   {"zeta",Pmpfr_zeta},
                     36:   {"j0",Pmpfr_j0},
                     37:   {"j1",Pmpfr_j1},
                     38:   {"y0",Pmpfr_y0},
                     39:   {"y1",Pmpfr_y1},
                     40:   {"eint",Pmpfr_eint},
                     41:   {"erf",Pmpfr_erf},
                     42:   {"li2",Pmpfr_li2},
                     43:   {"gamma",Pmpfr_gamma},
                     44:   {"lngamma",Pmpfr_gamma},
                     45:   {"digamma",Pmpfr_gamma},
                     46:   {"floor",Pmpfr_floor},
                     47:   {"ceil",Pmpfr_ceil},
                     48:   {"round",Pmpfr_round},
                     49: };
                     50:
                     51: mpfr_func mpfr_search(char *name)
                     52: {
                     53:   int i,n;
                     54:
                     55:   n = sizeof(mpfr_tab)/sizeof(struct mpfr_tab_rec);
                     56:   for ( i = 0; i < n; i++ )
                     57:     if ( !strcmp(name,mpfr_tab[i].name) )
                     58:       return mpfr_tab[i].func;
                     59:   return 0;
                     60: }
                     61:
                     62: Obj list_to_vect(Obj a)
                     63: {
                     64:   int len,i;
                     65:   VECT v;
                     66:   NODE nd;
                     67:
                     68:   if ( !a || OID(a) != O_LIST ) return a;
                     69:   len = length(BDY((LIST)a));
                     70:   MKVECT(v,len);
                     71:   for ( i = 0, nd = BDY((LIST)a); nd; nd = NEXT(nd), i++ )
                     72:      v->body[i] = (pointer)list_to_vect((Obj)BDY(nd));
                     73:   return (Obj)v;
                     74: }
                     75:
                     76: Obj vect_to_mat(VECT v)
                     77: {
                     78:   MAT m;
                     79:   int len,col,i,j;
                     80:
                     81:   len = v->len;
                     82:   if ( v->body[0] && OID((Obj)v->body[0]) == O_VECT ) {
                     83:     col = ((VECT)v->body[0])->len;
                     84:   for ( i = 1; i < len; i++ )
                     85:     if ( !v->body[i] || OID((Obj)v->body[i]) != O_VECT
                     86:      || ((VECT)v->body[i])->len != col )
                     87:     break;
                     88:     if ( i == len ) {
                     89:     /* convert to a matrix */
                     90:     MKMAT(m,len,col);
                     91:     for ( i = 0; i < len; i++ )
                     92:       for ( j = 0; j < col; j++ )
                     93:       m->body[i][j] = ((VECT)v->body[i])->body[j];
                     94:     return (Obj)m;
                     95:   }
                     96:   }
                     97:   return (Obj)v;
                     98: }
                     99:
                    100: void reset_ox_pari()
                    101: {
                    102:   NODE nd;
                    103:   Q r;
                    104:
                    105:   if ( ox_get_pari_result ) {
                    106:   nd = mknode(1,ox_pari_stream);
                    107:   Pox_shutdown(nd,&r);
                    108:     ox_get_pari_result = 0;
                    109:   ox_pari_stream_initialized = 0;
                    110:   }
                    111: }
                    112:
1.5     ! noro      113: void pari_setprec(long n)
        !           114: {
        !           115:   struct oFUNC f;
        !           116:   Z z;
        !           117:   NODE arg;
        !           118:
        !           119:   f.name = f.fullname = "pari_setprec";
        !           120:   f.id = A_PARI;
        !           121:   f.argc = 1;
        !           122:   f.f.binf = 0;
        !           123:   STOZ(n,z);
        !           124:   arg = mknode(1,z);
        !           125:   evalparif(&f,arg);
        !           126: }
        !           127:
        !           128: /* decimal precision */
        !           129: long mpfrprec = 0;
        !           130:
1.1       noro      131: pointer evalparif(FUNC f,NODE arg)
                    132: {
1.5     ! noro      133:   long prec;
        !           134:   int ac,intarg,opt;
1.1       noro      135:   Q q,r,cmd;
                    136:   Z narg;
                    137:   Real sec;
                    138:   NODE nd,oxarg,t,t1,n;
                    139:   STRING name;
                    140:   USINT ui;
                    141:   LIST list;
                    142:   Obj ret,dmy;
                    143:   mpfr_func mpfr_function;
                    144:   V v;
                    145:
1.3       noro      146:   if ( arg && (!ARG0(arg) || (NUM(ARG0(arg)) && NID((Num)ARG0(arg)) != N_C))
1.1       noro      147:     && (mpfr_function = mpfr_search(f->name)) ) {
                    148:     (*mpfr_function)(arg,&ret);
                    149:     return (pointer) ret;
                    150:   }
                    151:
                    152:   if ( !ox_pari_stream_initialized ) {
                    153:   if ( ox_pari_starting_function && OID(ox_pari_starting_function) == O_P ) {
                    154:     v = VR(ox_pari_starting_function);
                    155:     if ( (long)v->attr != V_SR ) {
                    156:       error("pari : no handler.");
                    157:     }
                    158:     MKNODE(nd,0,0);
                    159:     r = (Q)bevalf((FUNC)v->priv,0);
                    160:   }else {
                    161: #if !defined(VISUAL)
                    162:   MKSTR(name,"ox_pari");
                    163:   nd = mknode(2,NULL,name);
1.4       noro      164:   if ( debug_pari )
                    165:     Pox_launch(nd,(Obj *)&r);
                    166:   else
                    167:     Pox_launch_nox(nd,(Obj *)&r);
1.1       noro      168: #else
                    169:   error("Please load names.rr from latest asir-contrib library before using pari functions.");
                    170: #endif
                    171:   }
                    172:   ox_pari_stream = r;
                    173:     ox_pari_stream_initialized = 1;
                    174:   }
                    175:
1.5     ! noro      176:   prec = mpfr_get_default_prec()*0.30103+1;
        !           177:   if ( prec != mpfrprec ) {
        !           178:     mpfrprec = prec;
        !           179:     pari_setprec(prec);
        !           180:   }
1.1       noro      181:   ac = argc(arg);
                    182:   /* reverse the arg list */
                    183:   for ( n = arg, t = 0; n; n = NEXT(n) ) {
                    184:     MKNODE(t1,BDY(n),t); t = t1;
                    185:   }
                    186:   /* push the reversed arg list */
                    187:   for ( ; t; t = NEXT(t) ) {
                    188:     oxarg = mknode(2,ox_pari_stream,BDY(t));
                    189:     Pox_push_cmo(oxarg,&dmy);
                    190:   }
                    191:   MKSTR(name,f->name);
1.2       noro      192:   STOZ(ac,narg);
1.1       noro      193:   oxarg = mknode(3,ox_pari_stream,name,narg);
                    194:   Pox_execute_function(oxarg,&dmy);
                    195:   ox_get_pari_result = 1;
                    196: #if defined(VISUAL) || defined(__MINGW32__)
                    197: #define SM_popCMO 262
1.2       noro      198:   STOZ(SM_popCMO,cmd);
1.1       noro      199:   oxarg = mknode(2,ox_pari_stream,cmd);
                    200:   Pox_push_cmd(oxarg,&dmy);
                    201:   nd = mknode(1,ox_pari_stream);
                    202:   MKLIST(list,nd);
                    203:   MKReal(1.0/8,sec);
                    204:   oxarg = mknode(2,list,sec);
                    205:   ret=0;
                    206:   do {
                    207:     check_intr();
                    208:     Pox_select(oxarg,&list);
                    209:     oxarg = mknode(1,list);
                    210:     Plength(oxarg,&ret);
                    211:   }while (!ret);
                    212:   oxarg = mknode(1,ox_pari_stream);
                    213:   Pox_get(oxarg,&ret);
                    214: #else
                    215:   oxarg = mknode(1,ox_pari_stream);
                    216:   Pox_pop_cmo(oxarg,&ret);
                    217: #endif
                    218:   ox_get_pari_result = 0;
                    219:   if ( ret && OID(ret) == O_ERR ) {
                    220:     char buf[BUFSIZ];
                    221:     soutput_init(buf);
                    222:     sprintexpr(CO,((ERR)ret)->body);
                    223:     error(buf);
                    224:   }
                    225:   if ( ret && OID(ret) == O_LIST ) {
                    226:     ret = list_to_vect(ret);
                    227:   ret = vect_to_mat((VECT)ret);
                    228:   }
                    229:   return ret;
                    230: }
                    231:
                    232: struct pariftab {
                    233:   char *name;
                    234:   int dmy;
                    235:   int type;
                    236: };
                    237:
                    238: /*
                    239:  * type = 1 => argc = 1, second arg = precision
                    240:  * type = 2 => argc = 1, second arg = (long int)0
                    241:  *
                    242:  */
                    243: /*
                    244: {"abs",0,1},
                    245: {"adj",0,1},
                    246: */
                    247:
                    248: struct pariftab pariftab[] = {
                    249: {"arg",0,1},
                    250: {"bigomega",0,1},
                    251: {"binary",0,1},
                    252: {"ceil",0,1},
                    253: {"centerlift",0,1},
                    254: {"cf",0,1},
                    255: {"classno",0,1},
                    256: {"classno2",0,1},
                    257: {"conj",0,1},
                    258: {"content",0,1},
                    259: {"denom",0,1},
                    260: {"det",0,1},
                    261: {"det2",0,1},
                    262: {"dilog",0,1},
                    263: {"disc",0,1},
                    264: {"discf",0,1},
                    265: {"divisors",0,1},
                    266: {"eigen",0,1},
                    267: {"eintg1",0,1},
                    268: {"erfc",0,1},
                    269: {"eta",0,1},
                    270: {"floor",0,1},
                    271: {"frac",0,1},
                    272: {"galois",0,1},
                    273: {"galoisconj",0,1},
                    274: {"gamh",0,1},
                    275: {"gamma",0,1},
                    276: {"hclassno",0,1},
                    277: {"hermite",0,1},
                    278: {"hess",0,1},
                    279: {"imag",0,1},
                    280: {"image",0,1},
                    281: {"image2",0,1},
                    282: {"indexrank",0,1},
                    283: {"indsort",0,1},
                    284: {"initalg",0,1},
                    285: {"isfund",0,1},
                    286: {"ispsp",0,1},
                    287: {"isqrt",0,1},
                    288: {"issqfree",0,1},
                    289: {"issquare",0,1},
                    290: {"jacobi",0,1},
                    291: {"jell",0,1},
                    292: {"ker",0,1},
                    293: {"keri",0,1},
                    294: {"kerint",0,1},
                    295: {"kerintg1",0,1},
                    296: {"length",0,1},
                    297: {"lexsort",0,1},
                    298: {"lift",0,1},
                    299: {"lindep",0,1},
                    300: {"lll",0,1},
                    301: {"lllgen",0,1},
                    302: {"lllgram",0,1},
                    303: {"lllgramgen",0,1},
                    304: {"lllgramint",0,1},
                    305: {"lllgramkerim",0,1},
                    306: {"lllgramkerimgen",0,1},
                    307: {"lllint",0,1},
                    308: {"lllkerim",0,1},
                    309: {"lllkerimgen",0,1},
                    310: {"lngamma",0,1},
                    311: {"logagm",0,1},
                    312: {"mat",0,1},
                    313: {"matrixqz2",0,1},
                    314: {"matrixqz3",0,1},
                    315: {"matsize",0,1},
                    316: {"modreverse",0,1},
                    317: {"mu",0,1},
                    318: {"nextprime",0,1},
                    319: {"norm",0,1},
                    320: {"norml2",0,1},
                    321: {"numdiv",0,1},
                    322: {"numer",0,1},
                    323: {"omega",0,1},
                    324: {"order",0,1},
                    325: {"ordred",0,1},
                    326: {"phi",0,1},
                    327: {"pnqn",0,1},
                    328: {"polred",0,1},
                    329: {"polred2",0,1},
                    330: {"primroot",0,1},
                    331: {"psi",0,1},
                    332: {"quadgen",0,1},
                    333: {"quadpoly",0,1},
                    334: {"real",0,1},
                    335: {"recip",0,1},
                    336: {"redreal",0,1},
                    337: {"regula",0,1},
                    338: {"reorder",0,1},
                    339: {"reverse",0,1},
                    340: {"rhoreal",0,1},
                    341: {"roots",0,1},
                    342: {"round",0,1},
                    343: {"sigma",0,1},
                    344: {"signat",0,1},
                    345: {"simplify",0,1},
                    346: {"smalldiscf",0,1},
                    347: {"smallfact",0,1},
                    348: {"smallpolred",0,1},
                    349: {"smallpolred2",0,1},
                    350: {"smith",0,1},
                    351: {"smith2",0,1},
                    352: {"sort",0,1},
                    353: {"sqr",0,1},
                    354: {"sqred",0,1},
                    355: {"sqrt",0,1},
                    356: {"supplement",0,1},
                    357: {"trace",0,1},
                    358: {"trans",0,1},
                    359: {"trunc",0,1},
                    360: {"unit",0,1},
                    361: {"vec",0,1},
                    362: {"wf",0,1},
                    363: {"wf2",0,1},
                    364: {"zeta",0,1},
                    365: {"factor",0,1},
                    366:
                    367: {"allocatemem",0,0},
                    368:
                    369: {"factpol",0,2},
                    370: {"isprime",0,2},
                    371: {"factorint",0,2},
                    372: {0,0,0},
                    373: };
                    374:
                    375: void parif_init() {
                    376:   int i;
                    377:
                    378:   for ( i = 0, parif = 0; pariftab[i].name; i++ )
                    379:      appendparif(&parif,pariftab[i].name, 0,pariftab[i].type);
                    380: }

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