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

Diff for /OpenXM_contrib2/asir2000/engine/bf.c between version 1.5 and 1.12

version 1.5, 2009/03/27 14:42:29 version 1.12, 2015/08/20 08:42:07
Line 1 
Line 1 
 /*  /*
  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED   * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.11 2015/08/06 23:41:52 noro Exp $
  * 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/bf.c,v 1.4 2003/02/14 22:29:08 ohara Exp $  
 */  
 #include "ca.h"  #include "ca.h"
 #if defined(PARI)  
 #include "base.h"  #include "base.h"
 #include <math.h>  #include <math.h>
 #include "genpari.h"  
   
 void ritopa(Obj,GEN *);  extern int mpfr_roundmode;
 void patori(GEN,Obj *);  
   
 void addbf(a,b,c)  Num tobf(Num,int);
 Num a,b;  
 Num *c;  #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)
   
   void strtobf(char *s,BF *p)
 {  {
         GEN pa,pb,z;    BF r;
         long ltop,lbot;    NEWBF(r);
     mpfr_init(r->body);
     mpfr_set_str(r->body,s,10,mpfr_roundmode);
     *p = r;
   }
   
         if ( !a )  double mpfrtodbl(mpfr_t a)
                 *c = b;  {
         else if ( !b )    return mpfr_get_d(a,mpfr_roundmode);
                 *c = a;  
         else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )  
                 (*addnumt[MIN(NID(a),NID(b))])(a,b,c);  
         else {  
                 ltop = avma; ritopa((Obj)a,&pa);  
                 ritopa((Obj)b,&pb); lbot = avma;  
                 z = gerepile(ltop,lbot,gadd(pa,pb));  
                 patori(z,(Obj *)c); cgiv(z);  
         }  
 }  }
   
 void subbf(a,b,c)  Num tobf(Num a,int prec)
 Num a,b;  
 Num *c;  
 {  {
         GEN pa,pb,z;    mpfr_t r;
         long ltop,lbot;    mpz_t z;
     mpq_t q;
     BF d;
     N nm,dn;
     C c;
     Num re,im;
     int sgn;
   
         if ( !a )    if ( !a ) {
                 (*chsgnnumt[NID(b)])(b,c);      prec ? mpfr_init2(r,prec) : mpfr_init(r);
         else if ( !b )      mpfr_set_zero(r,1);
                 *c = a;      MPFRTOBF(r,d);
         else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )      return (Num)d;
                 (*subnumt[MIN(NID(a),NID(b))])(a,b,c);    } else {
         else {      switch ( NID(a) ) {
                 ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma;      case N_B:
                 z = gerepile(ltop,lbot,gsub(pa,pb));        return a;
                 patori(z,(Obj *)c); cgiv(z);        break;
         }      case N_R:
         prec ? mpfr_init2(r,prec) : mpfr_init(r);
         mpfr_init_set_d(r,((Real)a)->body,mpfr_roundmode);
         MPFRTOBF(r,d);
         return (Num)d;
         break;
       case N_Q:
         nm = NM((Q)a); dn = DN((Q)a); sgn = SGN((Q)a);
         if ( INT((Q)a) ) {
           mpz_init(z);
           mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
           if ( sgn < 0 ) mpz_neg(z,z);
           mpfr_init_set_z(r,z,mpfr_roundmode);
         } else {
           mpq_init(q);
           mpz_import(mpq_numref(q),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
           mpz_import(mpq_denref(q),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
           if ( sgn < 0 ) mpq_neg(q,q);
           mpfr_init_set_q(r,q,mpfr_roundmode);
         }
         MPFRTOBF(r,d);
         return (Num)d;
         break;
       case N_C:
         re = tobf(((C)a)->r,prec); im = tobf(((C)a)->i,prec);
         NEWC(c); c->r = re; c->i = im;
         return (Num)c;
         break;
       default:
         error("tobf : invalid argument");
         break;
       }
     }
 }  }
   
 void mulbf(a,b,c)  void addbf(Num a,Num b,Num *c)
 Num a,b;  
 Num *c;  
 {  {
         GEN pa,pb,z;    mpfr_t r;
         long ltop,lbot;    BF d;
     GZ z;
     GQ q;
   
         if ( !a || !b )    if ( !a )
                 *c = 0;      *c = b;
         else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( !b )
                 (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);      *c = a;
         else {    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                 ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma;      (*addnumt[MIN(NID(a),NID(b))])(a,b,c);
                 z = gerepile(ltop,lbot,gmul(pa,pb));    else if ( NID(a) == N_B ) {
                 patori(z,(Obj *)c); cgiv(z);      switch ( NID(b) ) {
         }      case N_Q:
         mpfr_init2(r,BFPREC(a));
         if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_add_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_add_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_add_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_add(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else { /* NID(b)==N_B */
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           mpfr_add_z(r,((BF)b)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)a);
           mpfr_add_q(r,((BF)b)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_add_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     }
     if ( !cmpbf(*c,0) ) *c = 0;
 }  }
   
 void divbf(a,b,c)  void subbf(Num a,Num b,Num *c)
 Num a,b;  
 Num *c;  
 {  {
         GEN pa,pb,z;    mpfr_t r,s;
         long ltop,lbot;    GZ z;
     GQ q;
     BF d;
   
         if ( !b )    if ( !a )
                 error("divbf : division by 0");      (*chsgnnumt[NID(b)])(b,c);
         else if ( !a )    else if ( !b )
                 *c = 0;      *c = a;
         else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                 (*divnumt[MIN(NID(a),NID(b))])(a,b,c);      (*subnumt[MIN(NID(a),NID(b))])(a,b,c);
         else {    else if ( NID(a) == N_B ) {
                 ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)b,&pb); lbot = avma;      switch ( NID(b) ) {
                 z = gerepile(ltop,lbot,gdiv(pa,pb));      case N_Q:
                 patori(z,(Obj *)c); cgiv(z);        mpfr_init2(r,BFPREC(a));
         }        if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_sub_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_sub_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_sub(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else { /* NID(b)==N_B */
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           mpfr_sub_z(r,((BF)b)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)a);
           mpfr_sub_q(r,((BF)b)->body,q->body,mpfr_roundmode);
         }
         mpfr_neg(r,r,mpfr_roundmode);
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     }
     if ( !cmpbf(*c,0) ) *c = 0;
 }  }
   
 void pwrbf(a,e,c)  void mulbf(Num a,Num b,Num *c)
 Num a,e;  
 Num *c;  
 {  {
         if ( !e )    mpfr_t r;
                 *c = (Num)ONE;    GZ z;
         else if ( !a )    GQ q;
                 *c = 0;    BF d;
         else {    int prec;
                 gpui_ri((Obj)a,(Obj)c,(Obj *)c);  
         }    if ( !a || !b )
       *c = 0;
     else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
       (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);
     else if ( NID(a) == N_B ) {
       switch ( NID(b) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(a));
         if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_mul_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_mul_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_mul(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else { /* NID(b)==N_B */
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           mpfr_mul_z(r,((BF)b)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)a);
           mpfr_mul_q(r,((BF)b)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     }
     if ( !cmpbf(*c,0) ) *c = 0;
 }  }
   
 void chsgnbf(a,c)  void divbf(Num a,Num b,Num *c)
 Num a,*c;  
 {  {
         BF t;    mpfr_t s,r;
         GEN z;    GZ z;
         int s;    GQ q;
     BF d;
   
         if ( !a )    if ( !b )
                 *c = 0;      error("divbf : division by 0");
         else if ( NID(a) <= N_A )    else if ( !a )
                 (*chsgnnumt[NID(a)])(a,c);      *c = 0;
         else {    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                 z = (GEN)((BF)a)->body; s = lg(z); NEWBF(t,s);      (*divnumt[MIN(NID(a),NID(b))])(a,b,c);
                 bcopy((char *)a,(char *)t,sizeof(struct oBF)+((s-1)*sizeof(long)));    else if ( NID(a) == N_B ) {
                 z = (GEN)((BF)t)->body; setsigne(z,-signe(z));      switch ( NID(b) ) {
                 *c = (Num)t;      case N_Q:
         }        mpfr_init2(r,BFPREC(a));
         if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_div_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_div_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_div_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else { /* NID(b)==N_B */
       switch ( NID(a) ) {
       case N_Q:
         /* XXX : mpfr_z_div and mpfr_q_div are not implemented */
         a = tobf(a,BFPREC(b));
         mpfr_init2(r,BFPREC(b));
         mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_d_div(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     }
     if ( !cmpbf(*c,0) ) *c = 0;
 }  }
   
 int cmpbf(a,b)  void pwrbf(Num a,Num b,Num *c)
 Num a,b;  
 {  {
         GEN pa,pb;    int prec;
         int s;    mpfr_t r;
     GZ z;
     BF d;
   
         if ( !a ) {    if ( !b )
                 if ( !b || (NID(b)<=N_A) )      *c = (Num)ONE;
                         return (*cmpnumt[NID(b)])(a,b);    else if ( !a )
                 else      *c = 0;
                         return -signe(((BF)b)->body);    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
         } else if ( !b ) {      (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);
                 if ( !a || (NID(a)<=N_A) )    else if ( NID(a) == N_B ) {
                         return (*cmpnumt[NID(a)])(a,b);      switch ( NID(b) ) {
                 else      case N_Q:
                         return signe(((BF)a)->body);        mpfr_init2(r,BFPREC(a));
         } else {        if ( INT((Q)b) ) {
                 ritopa((Obj)a,&pa); ritopa((Obj)b,&pb);          z = ztogz((Q)b);
                 s = gcmp(pa,pb); cgiv(pb); cgiv(pa);          mpfr_pow_z(r,((BF)a)->body,z->body,mpfr_roundmode);
                 return s;        } else {
         }          b = tobf(b,BFPREC(a));
           mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         prec = MAX(BFPREC(a),53);
         mpfr_init2(r,prec);
         b = tobf(b,prec);
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else { /* NID(b)==N_B */
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         a = tobf(a,BFPREC(b));
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       case N_R:
         /* double precision = 53 */
         prec = MAX(BFPREC(a),53);
         mpfr_init2(r,prec);
         a = tobf(a,prec);
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     }
     if ( !cmpbf(*c,0) ) *c = 0;
 }  }
 #endif /* PARI */  
   void chsgnbf(Num a,Num *c)
   {
     mpfr_t r;
     BF d;
   
     if ( !a )
       *c = 0;
     else if ( NID(a) <= N_A )
       (*chsgnnumt[NID(a)])(a,c);
     else {
       mpfr_init2(r,BFPREC(a));
       mpfr_neg(r,((BF)a)->body,mpfr_roundmode);
       MPFRTOBF(r,d);
       *c = (Num)d;
     }
   }
   
   int cmpbf(Num a,Num b)
   {
     int ret;
     GZ z;
     GQ q;
   
     if ( !a ) {
       if ( !b ) return 0;
       else if ((NID(b)<=N_A) )
         return (*cmpnumt[NID(b)])(a,b);
       else
         return -mpfr_sgn(((BF)a)->body);
     } else if ( !b ) {
       if ( !a || (NID(a)<=N_A) )
         return (*cmpnumt[NID(a)])(a,b);
       else
         return mpfr_sgn(((BF)a)->body);
     } else if ( NID(a) == N_B ) {
       switch ( NID(b) ) {
       case N_Q:
         if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           ret = mpfr_cmp_z(((BF)a)->body,z->body);
         } else {
           q = qtogq((Q)b);
           ret = mpfr_cmp_q(((BF)a)->body,q->body);
         }
         break;
       case N_R:
         /* double precision = 53 */
         ret = mpfr_cmp_d(((BF)a)->body,((Real)b)->body);
         break;
       case N_B:
         ret = mpfr_cmp(((BF)a)->body,((BF)b)->body);
         break;
       }
       return ret;
     } else { /* NID(b)==N_B */
       switch ( NID(a) ) {
       case N_Q:
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           ret = mpfr_cmp_z(((BF)b)->body,z->body);
         } else {
           q = qtogq((Q)a);
           ret = mpfr_cmp_q(((BF)b)->body,q->body);
         }
         break;
       case N_R:
         /* double precision = 53 */
         ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body);
         break;
       }
       return -ret;
     }
   }
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.12
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>