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

Diff for /OpenXM_contrib2/asir2000/engine/pari.c between version 1.1 and 1.13

version 1.1, 1999/12/03 07:39:08 version 1.13, 2018/03/29 01:32:52
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/asir99/engine/pari.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */  /*
    * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
    * All rights reserved.
    *
    * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
    * non-exclusive and royalty-free license to use, copy, modify and
    * redistribute, solely for non-commercial and non-profit purposes, the
    * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
    * conditions of this Agreement. For the avoidance of doubt, you acquire
    * only a limited right to use the SOFTWARE hereunder, and FLL or any
    * third party developer retains all rights, including but not limited to
    * copyrights, in and to the SOFTWARE.
    *
    * (1) FLL does not grant you a license in any way for commercial
    * purposes. You may use the SOFTWARE only for non-commercial and
    * non-profit purposes only, such as academic, research and internal
    * business use.
    * (2) The SOFTWARE is protected by the Copyright Law of Japan and
    * international copyright treaties. If you make copies of the SOFTWARE,
    * with or without modification, as permitted hereunder, you shall affix
    * to all such copies of the SOFTWARE the above copyright notice.
    * (3) An explicit reference to this SOFTWARE and its copyright owner
    * shall be made on your publication or presentation in any form of the
    * results obtained by use of the SOFTWARE.
    * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
    * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
    * for such modification or the source code of the modified part of the
    * SOFTWARE.
    *
    * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
    * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
    * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
    * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
    * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
    * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
    * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
    * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
    * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
    * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
    * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
    * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
    * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
    * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
    * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
    * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
    * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
    *
    * $OpenXM: OpenXM_contrib2/asir2000/engine/pari.c,v 1.12 2014/07/31 08:01:29 noro Exp $
   */
 #include "ca.h"  #include "ca.h"
 #if PARI  
   #if defined(PARI)
 #include "base.h"  #include "base.h"
 #include <math.h>  #include <math.h>
 #include "genpari.h"  #include "genpari.h"
   
 #if defined(THINK_C)  
 void patori(GEN,Obj *);  void patori(GEN,Obj *);
 void patori_i(GEN,N *);  void patori_i(GEN,N *);
 void ritopa(Obj,GEN *);  void ritopa(Obj,GEN *);
 void ritopa_i(N,int,GEN *);  void ritopa_i(N,int,GEN *);
 #else  
 void patori();  
 void patori_i();  
 void ritopa();  
 void ritopa_i();  
 #endif  
   
   // PARI_VERSION(2,2,12) == 131596
   #if PARI_VERSION_CODE >= 131596
   #define prec precreal
   #endif
 extern long prec;  extern long prec;
 extern int paristack;  extern int paristack;
   
 void risa_pari_init() {  long get_pariprec() {
         char buf[BUFSIZ];    return prec;
         int i;  }
   void set_pariprec(long p) {
     prec = p;
   }
   
         pari_init(paristack,2);  void risa_pari_init() {
         prec = 4;    pari_init(paristack,2);
     set_pariprec(4);
 }  }
   
 void create_pari_variable(index)  void create_pari_variable(index)
 int index;  int index;
 {  {
         static int max_varn;    static int max_varn;
         int i;    int i;
         char name[BUFSIZ];    char name[BUFSIZ];
   
         if ( index > max_varn ) {    if ( index > max_varn ) {
                 for ( i = max_varn+1; i <= index; i++ ) {      for ( i = max_varn+1; i <= index; i++ ) {
                         sprintf(name,"x%d",i);        sprintf(name,"x%d",i);
                         fetch_named_var(name,0);  #if (PARI_VERSION_CODE < 131594)
                 }        fetch_named_var(name,0);
                 max_varn = index;  #else
         }        fetch_named_var(name);
   #endif
       }
       max_varn = index;
     }
 }  }
   
 int get_lg(a)  int get_lg(a)
 GEN a;  GEN a;
 {  {
         return lg(a);    return lg(a);
 }  }
   
   void gpui_ri(Obj a, Obj e, Obj *c)
   {
     GEN pa,pe,z;
     long ltop,lbot;
   
     ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
     z = gerepile(ltop,lbot,gpui(pa,pe,get_pariprec()));
     patori(z,c); cgiv(z);
   }
   
 void ritopa(a,rp)  void ritopa(a,rp)
 Obj a;  Obj a;
 GEN *rp;  GEN *rp;
 {  {
         long ltop;    long ltop;
         GEN pnm,z,w;    GEN pnm,z,w,u;
         DCP dc;    DCP dc;
         int i,j,l,row,col;    int i,j,l,row,col;
         VL vl;    VL vl;
         V v;    V v;
   
         if ( !a ) {    if ( !a ) {
                 *rp = gzero; return;      *rp = gzero; return;
         }    }
         switch ( OID(a) ) {    switch ( OID(a) ) {
                 case O_N:      case O_N:
                         switch ( NID(a) ) {        switch ( NID(a) ) {
                                 case N_Q:          case N_Q:
                                         ltop = avma; ritopa_i(NM((Q)a),SGN((Q)a),&pnm);            ltop = avma; ritopa_i(NM((Q)a),SGN((Q)a),&pnm);
                                         if ( INT((Q)a) )            if ( INT((Q)a) )
                                                 *rp = pnm;              *rp = pnm;
                                         else {            else {
                                                 *rp = z = cgetg(3,4); z[1] = (long)pnm;              *rp = z = cgetg(3,4); z[1] = (long)pnm;
                                                 ritopa_i(DN((Q)a),1,(GEN *)&z[2]);              ritopa_i(DN((Q)a),1,&u); z[2] = u;
                                         }            }
                                         break;            break;
                                 case N_R:          case N_R:
                                         *rp = dbltor(BDY((Real)a)); break;            *rp = dbltor(BDY((Real)a)); break;
                                 case N_B:          case N_B:
                                         *rp = gcopy((GEN)BDY(((BF)a))); break;            *rp = gcopy((GEN)BDY(((BF)a))); break;
                                 case N_C:          case N_C:
                                         z = cgetg(3,6);            z = cgetg(3,6);
                                         ritopa((Obj)((C)a)->r,(GEN *)&z[1]); ritopa((Obj)((C)a)->i,(GEN *)&z[2]);            ritopa((Obj)((C)a)->r,&u); z[1] = u;
                                         *rp = z;            ritopa((Obj)((C)a)->i,&u); z[2] = u;
                                         break;            *rp = z;
                                 default:            break;
                                         error("ritopa : not implemented"); break;          default:
                         }            error("ritopa : not implemented"); break;
                         break;        }
                 case O_P:        break;
                         l = UDEG((P)a); *rp = z = cgetg(l+3,10);      case O_P:
                         setsigne(z,1);        l = UDEG((P)a); *rp = z = cgetg(l+3,10);
                         for ( i = 0, vl = CO, v = VR((P)a); vl->v != v;        setsigne(z,1);
                                 vl = NEXT(vl), i++ );        for ( i = 0, vl = CO, v = VR((P)a); vl->v != v;
                         create_pari_variable(i);          vl = NEXT(vl), i++ );
                         setvarn(z,i);        create_pari_variable(i);
                         setlgef(z,l+3);        setvarn(z,i);
                         for ( i = l+2; i >= 2; i-- )        setlgef(z,l+3);
                                 z[i] = (long)gzero;        for ( i = l+2; i >= 2; i-- )
                         for ( dc = DC((P)a); dc; dc = NEXT(dc) )          z[i] = (long)gzero;
                                 ritopa((Obj)COEF(dc),(GEN *)&z[QTOS(DEG(dc))+2]);        for ( dc = DC((P)a); dc; dc = NEXT(dc) ) {
                         break;          ritopa((Obj)COEF(dc),&u); z[QTOS(DEG(dc))+2] = u;
                 case O_VECT:        }
                         l = ((VECT)a)->len; z = cgetg(l+1,17);        break;
                         for ( i = 0; i < l; i++ )      case O_VECT:
                                 ritopa((Obj)BDY((VECT)a)[i],(GEN *)&z[i+1]);        l = ((VECT)a)->len; z = cgetg(l+1,17);
                         *rp = z;        for ( i = 0; i < l; i++ ) {
                         break;          ritopa((Obj)BDY((VECT)a)[i],&u); z[i+1] = u;
                 case O_MAT:        }
                         row = ((MAT)a)->row; col = ((MAT)a)->col; z = cgetg(col+1,19);        *rp = z;
                         for ( j = 0; j < col; j++ ) {        break;
                                 w = cgetg(row+1,18);      case O_MAT:
                                 for ( i = 0; i < row; i++ )        row = ((MAT)a)->row; col = ((MAT)a)->col; z = cgetg(col+1,19);
                                         ritopa((Obj)BDY((MAT)a)[i][j],(GEN *)&w[i+1]);        for ( j = 0; j < col; j++ ) {
                                 z[j+1] = (long)w;          w = cgetg(row+1,18);
                         }          for ( i = 0; i < row; i++ ) {
                         *rp = z;            ritopa((Obj)BDY((MAT)a)[i][j],&u); w[i+1] = u;
                         break;          }
                 default:          z[j+1] = (long)w;
                         error("ritopa : not implemented"); break;        }
         }        *rp = z;
         break;
       default:
         error("ritopa : not implemented"); break;
     }
 }  }
   
 void patori(a,rp)  void patori(a,rp)
 GEN a;  GEN a;
 Obj *rp;  Obj *rp;
 {  {
         Q q,qnm,qdn;    Q q,qnm,qdn;
         BF r;    BF r;
         C c;    C c;
         VECT v;    VECT v;
         MAT m;    MAT m;
         N n,nm,dn;    N n,nm,dn;
         DCP dc0,dc;    DCP dc0,dc;
         P t;    P t;
         int s,i,j,l,row,col;    int s,i,j,l,row,col,ty;
         GEN b;    GEN b;
         VL vl;    VL vl;
   
         if ( gcmp0(a) )    if ( gcmp0(a) )
                 *rp = 0;      *rp = 0;
         else {    else {
                 switch ( typ(a) ) {      switch ( ty = typ(a) ) {
                         case 1: /* integer */        case 1: /* integer */
                                 patori_i(a,&n); NTOQ(n,(char)signe(a),q); *rp = (Obj)q;          patori_i(a,&n); NTOQ(n,(char)signe(a),q); *rp = (Obj)q;
                                 break;          break;
                         case 2: /* real */        case 2: /* real */
                                 NEWBF(r,lg(a)); bcopy((char *)a,(char *)BDY(r),lg(a)*sizeof(long));          NEWBF(r,lg(a)); bcopy((char *)a,(char *)BDY(r),lg(a)*sizeof(long));
                                 *rp = (Obj)r;          *rp = (Obj)r;
                                 break;          break;
                         case 4: /* rational; reduced */        case 3: /* integermod */
                                 patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);          MKVECT(v,2);
                                 s = signe(a[1])*signe(a[2]);          patori((GEN)a[1],(Obj *)&BDY(v)[0]);
                                 if ( UNIN(dn) )          patori((GEN)a[2],(Obj *)&BDY(v)[1]);
                                         NTOQ(nm,s,q);          *rp = (Obj)v;
                                 else          break;
                                         NDTOQ(nm,dn,s,q);        case 4: /* rational; reduced */
                                 *rp = (Obj)q;          patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
                                 break;          s = signe(a[1])*signe(a[2]);
                         case 5: /* rational; not necessarily reduced */          if ( UNIN(dn) )
                                 patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);            NTOQ(nm,s,q);
                                 s = signe(a[1])*signe(a[2]);          else
                                 NTOQ(nm,s,qnm); NTOQ(dn,1,qdn); divq(qnm,qdn,(Q *)rp);            NDTOQ(nm,dn,s,q);
                                 break;          *rp = (Obj)q;
                         case 6: /* complex */          break;
                                 if ( gcmp0((GEN)a[2]) )        case 5: /* rational; not necessarily reduced */
                                         patori((GEN)a[1],rp);          patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
                                 else {          s = signe(a[1])*signe(a[2]);
                                         NEWC(c); patori((GEN)a[1],(Obj *)&c->r); patori((GEN)a[2],(Obj *)&c->i);          NTOQ(nm,s,qnm); NTOQ(dn,1,qdn); divq(qnm,qdn,(Q *)rp);
                                         *rp = (Obj)c;          break;
                                 }        case 6: /* complex */
                                 break;          if ( gcmp0((GEN)a[2]) )
                         case 10: /* polynomial */            patori((GEN)a[1],rp);
                                 for ( i = lgef(a)-1, dc0 = 0; i >= 2; i-- )          else {
                                         if ( !gcmp0((GEN)a[i]) ) {            NEWC(c); patori((GEN)a[1],(Obj *)&c->r); patori((GEN)a[2],(Obj *)&c->i);
                                                 NEXTDC(dc0,dc);            *rp = (Obj)c;
                                                 patori((GEN)a[i],(Obj *)&COEF(dc)); STOQ(i-2,DEG(dc));          }
                                         }          break;
                                 if ( !dc0 )        case 10: /* polynomial */
                                         *rp = 0;          for ( i = lgef(a)-1, dc0 = 0; i >= 2; i-- )
                                 else {            if ( !gcmp0((GEN)a[i]) ) {
                                         /* assuming that CO has not changed. */              NEXTDC(dc0,dc);
                                         for ( s = varn(a), i = 0, vl = CO; i != s;              patori((GEN)a[i],(Obj *)&COEF(dc)); STOQ(i-2,DEG(dc));
                                                 i++, vl = NEXT(vl) );            }
                                         NEXT(dc) = 0; MKP(vl->v,dc0,t); *rp = (Obj)t;          if ( !dc0 )
                                 }            *rp = 0;
                                 break;          else {
                         case 17: /* row vector */            /* assuming that CO has not changed. */
                         case 18: /* column vector */            for ( s = varn(a), i = 0, vl = CO; i != s;
                                 l = lg(a)-1; MKVECT(v,l);              i++, vl = NEXT(vl) );
                                 for ( i = 0; i < l; i++ )            NEXT(dc) = 0; MKP(vl->v,dc0,t); *rp = (Obj)t;
                                         patori((GEN)a[i+1],(Obj *)&BDY(v)[i]);          }
                                 *rp = (Obj)v;          break;
                                 break;        case 17: /* row vector */
                         case 19: /* matrix */        case 18: /* column vector */
                                 col = lg(a)-1; row = lg(a[1])-1; MKMAT(m,row,col);          l = lg(a)-1; MKVECT(v,l);
                                 for ( j = 0; j < col; j++ )          for ( i = 0; i < l; i++ )
                                         for ( i = 0, b = (GEN)a[j+1]; i < row; i++ )            patori((GEN)a[i+1],(Obj *)&BDY(v)[i]);
                                                 patori((GEN)b[i+1],(Obj *)&BDY(m)[i][j]);          *rp = (Obj)v;
                                 *rp = (Obj)m;          break;
                                 break;        case 19: /* matrix */
                         default:          col = lg(a)-1; row = lg(a[1])-1; MKMAT(m,row,col);
                                 error("patori : not implemented");          for ( j = 0; j < col; j++ )
                                 break;            for ( i = 0, b = (GEN)a[j+1]; i < row; i++ )
                 }              patori((GEN)b[i+1],(Obj *)&BDY(m)[i][j]);
         }          *rp = (Obj)m;
           break;
         default:
           error("patori : not implemented");
           break;
       }
     }
 }  }
   
 #if defined(LONG_IS_32BIT)  #if SIZEOF_LONG == 4
 void ritopa_i(a,s,rp)  void ritopa_i(a,s,rp)
 N a;  N a;
 int s;  int s;
 GEN *rp;  GEN *rp;
 {  {
         int j,l;    int j,l;
         unsigned int *b;    unsigned int *b;
         GEN z;    GEN z;
   
         l = PL(a); b = (unsigned int *)BD(a);    l = PL(a); b = (unsigned int *)BD(a);
         z = cgeti(l+2);    z = cgeti(l+2);
         bzero((char *)&z[2],l*sizeof(int));    bzero((char *)&z[2],l*sizeof(int));
         for ( j = 0; j < l; j++ )    for ( j = 0; j < l; j++ )
                 z[l-j+1] = b[j];      z[l-j+1] = b[j];
         s = s>0?1:-1;    s = s>0?1:-1;
         setsigne(z,s);    setsigne(z,s);
         setlgefint(z,lg(z));    setlgefint(z,lg(z));
         *rp = z;    *rp = z;
 }  }
   
 void patori_i(g,rp)  void patori_i(g,rp)
 GEN g;  GEN g;
 N *rp;  N *rp;
 {  {
         int j,l;    int j,l;
         unsigned int *a,*b;    unsigned int *a,*b;
         N z;    N z;
   
         l = lgef(g)-2;    l = lgef(g)-2;
         a = (unsigned int *)g;    a = (unsigned int *)g;
         *rp = z = NALLOC(l); PL(z) = l;    *rp = z = NALLOC(l); PL(z) = l;
         for ( j = 0, b = (unsigned int *)BD(z); j < l; j++ )    for ( j = 0, b = (unsigned int *)BD(z); j < l; j++ )
                 b[l-j-1] = ((unsigned int *)g)[j+2];      b[l-j-1] = ((unsigned int *)g)[j+2];
 }  }
 #endif  
   
 #if defined(LONG_IS_64BIT)  #elif SIZEOF_LONG == 8
 void ritopa_i(a,s,rp)  void ritopa_i(a,s,rp)
 N a;  N a;
 int s;  int s;
 GEN *rp;  GEN *rp;
 {  {
         int j,l,words;    int j,l,words;
         unsigned int *b;    unsigned int *b;
         GEN z;    GEN z;
   
         l = PL(a); b = BD(a);    l = PL(a); b = BD(a);
         words = (l+1)/2;    words = (l+1)/2;
         z = cgeti(words+2);    z = cgeti(words+2);
         bzero((char *)&z[2],words*sizeof(long));    bzero((char *)&z[2],words*sizeof(long));
         for ( j = 0; j < words; j++ )    for ( j = 0; j < words; j++ )
                 z[words-j+1] = ((unsigned long)b[2*j])      z[words-j+1] = ((unsigned long)b[2*j])
                         |(((unsigned long)(2*j+1<l?b[2*j+1]:0))<<32);        |(((unsigned long)(2*j+1<l?b[2*j+1]:0))<<32);
         s = s>0?1:-1;    s = s>0?1:-1;
         setsigne(z,s);    setsigne(z,s);
         setlgefint(z,lg(z));    setlgefint(z,lg(z));
         *rp = z;    *rp = z;
 }  }
   
 void patori_i(g,rp)  void patori_i(g,rp)
 GEN g;  GEN g;
 N *rp;  N *rp;
 {  {
         int j,l,words;    int j,l,words;
         unsigned long t;    unsigned long t;
         unsigned long *a;    unsigned long *a;
         unsigned int *b;    unsigned int *b;
         N z;    N z;
   
         words = lgef(g)-2;    words = lgef(g)-2;
         l = 2*words;    l = 2*words;
         a = (unsigned long *)g;    a = (unsigned long *)g;
         *rp = z = NALLOC(l); PL(z) = l;    *rp = z = NALLOC(l); PL(z) = l;
         for ( j = 0, b = BD(z); j < words; j++ ) {    for ( j = 0, b = BD(z); j < words; j++ ) {
                 t = a[words+1-j];      t = a[words+1-j];
                 b[2*j] = t&0xffffffff;      b[2*j] = t&0xffffffff;
                 b[2*j+1] = t>>32;      b[2*j+1] = t>>32;
         }    }
         PL(z) = b[l-1] ? l : l-1;    PL(z) = b[l-1] ? l : l-1;
 }  }
 #endif  #endif
   
Line 295  void strtobf(s,rp)
Line 368  void strtobf(s,rp)
 char *s;  char *s;
 BF *rp;  BF *rp;
 {  {
         GEN z;    GEN z;
   
         z = lisexpr(s); patori(z,(Obj *)rp); cgiv(z);    z = lisexpr(s); patori(z,(Obj *)rp); cgiv(z);
 }  }
 #endif  #endif
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.13
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>