1: //========================================================================================
     2: //  QFP.java
     3: //    en:Quad-precision floating-point number
     4: //    ja:四倍精度浮動小数点数
     5: //  Copyright (C) 2003-2023 Makoto Kamada
     6: //
     7: //  This file is part of the XEiJ (X68000 Emulator in Java).
     8: //  You can use, modify and redistribute the XEiJ if the conditions are met.
     9: //  Read the XEiJ License for more details.
    10: //  https://stdkmd.net/xeij/
    11: //========================================================================================
    12: 
    13: //----------------------------------------------------------------------------------------
    14: //  Double-Doubleにフラグと指数部の下駄を加えたもの
    15: //  メモ
    16: //    Double-Doubleの弱点は表現できる数の範囲(指数部の範囲)がdoubleと同じであることだが、doubleと同じ範囲の計算ができるわけではない
    17: //    doubleではオーバーフローしない上限に近い数の計算がDouble-Doubleでは途中でオーバーフローしてしまい完遂できないことがある
    18: //----------------------------------------------------------------------------------------
    19: 
    20: package xeij;
    21: 
    22: import java.lang.*;  //Boolean,Character,Class,Comparable,Double,Exception,Float,IllegalArgumentException,Integer,Long,Math,Number,Object,Runnable,SecurityException,String,StringBuilder,System
    23: 
    24: public class QFP {
    25: 
    26:   //定数
    27:   //  flg
    28:   public static final int QFP_P  = 0x00000000;  //+
    29:   public static final int QFP_M  = 0x08000000;  //-
    30:   public static final int QFP_Z  = 0x04000000;  //±0
    31:   public static final int QFP_I  = 0x02000000;  //±Inf
    32:   public static final int QFP_N  = 0x01000000;  //NaN
    33:   public static final int QFP_ZIN = QFP_Z | QFP_I | QFP_N;  //±0,±Inf,NaN
    34:   public static final int QFP_MZIN = QFP_M | QFP_Z | QFP_I | QFP_N;  //-,±0,±Inf,NaN
    35:   //  fpsr
    36:   public static final int QFP_BS = 0x00008000;  //比較不能な分岐/セット
    37:   public static final int QFP_SN = 0x00004000;  //シグナリングNaN
    38:   public static final int QFP_OE = 0x00002000;  //オペランドエラー
    39:   public static final int QFP_OF = 0x00001000;  //オーバーフロー
    40:   public static final int QFP_UF = 0x00000800;  //アンダーフロー
    41:   public static final int QFP_DZ = 0x00000400;  //ゼロ除算
    42:   public static final int QFP_X2 = 0x00000200;  //不正確な結果
    43:   public static final int QFP_X1 = 0x00000100;  //不正確な10進数
    44:   public static final int QFP_EXC = 0x0000ff00;
    45:   //  指数部の下駄
    46:   //  加減算の軽量化のため、dvlの指数部の範囲は1.5倍+αしたときdoubleに収まらなければならない
    47:   //  dvlの指数部の範囲はfloatをカバーできるものとする
    48:   private static final int QFP_GETA_SIZE = 512;  //dvlの指数部の上限-下限+1
    49:   private static final int QFP_GETA_BASE = -QFP_GETA_SIZE >> 1;  //dvlの指数部の下限
    50: 
    51:   //  echo read("efp.gp");qfppub("QFP_ONE",eval("1"),"1") | gp-2.7 -q
    52:   public static final QFP QFP_ONE = new QFP (QFP_P, 0, 0x1.0p0, 0.0);  //=1=1
    53:   //  echo read("efp.gp");qfppub("QFP_TEN",eval("10"),"10") | gp-2.7 -q
    54:   public static final QFP QFP_TEN = new QFP (QFP_P, 0, 0x1.4p3, 0.0);  //=10=10
    55:   //  echo read("efp.gp");qfppub("QFP_TENTO4",eval("10^4"),"10^4") | gp-2.7 -q
    56:   public static final QFP QFP_TENTO4 = new QFP (QFP_P, 0, 0x1.388p13, 0.0);  //=10^4=10000
    57:   //  echo read("efp.gp");qfppub("QFP_PI",Pi,"pi") | gp-2.7 -q
    58:   public static final QFP QFP_PI = new QFP (QFP_P, 0, 0x1.921fb54442d18p1, 0x1.1a62633145c07p-53);  //>pi=3.14159265358979323846264338328
    59:   //  echo read("efp.gp");qfppub("QFP_LN_2_2",log(2)/2,"log(2)/2") | gp-2.7 -q
    60:   public static final QFP QFP_LN_2_2 = new QFP (QFP_P, 0, 0x1.62e42fefa39efp-2, 0x1.abc9e3b39803fp-57);  //<log(2)/2=0.34657359027997265470861606072909
    61: 
    62:   //クラスフィールド
    63:   public static int fpsr;
    64: 
    65:   //インスタンスフィールド
    66:   public int flg;  //フラグ。±0,±Inf,NaN以外のときMはdvlの符号のコピー
    67:   public int epp;  //指数部の下駄。指数部-QFP_GETA_BASE&-QFP_GETA_SIZE
    68:   public double dvl;  //上位。符号を含む。指数部の範囲はQFP_GETA_BASE..QFP_GETA_BASE+QFP_GETA_SIZE-1
    69:   public double cvl;  //下位。符号を含む
    70: 
    71:   //公開コンストラクタ
    72:   public QFP () {
    73:     this.set0 ();
    74:   }  //new QFP()
    75:   public QFP (double d) {
    76:     this.setd (d);
    77:   }  //new QFP(double)
    78:   public QFP (float f) {
    79:     this.setf (f);
    80:   }  //new QFP(float)
    81:   public QFP (int i) {
    82:     this.seti (i);
    83:   }  //new QFP(int)
    84:   public QFP (long l) {
    85:     this.setl (l);
    86:   }  //new QFP(long)
    87:   public QFP (QFP x) {
    88:     this.flg = x.flg;
    89:     this.epp = x.epp;
    90:     this.dvl = x.dvl;
    91:     this.cvl = x.cvl;
    92:   }  //new QFP(QFP)
    93: 
    94:   //内部コンストラクタ
    95:   private QFP (int f, int e, double d, double c) {
    96:     this.flg = f;
    97:     this.epp = e;
    98:     this.dvl = d;
    99:     this.cvl = c;
   100:   }  //new QFP(int,int,double,double)
   101: 
   102:   //qfp = qfp.abs ()
   103:   //qfp = qfp.abs (x)
   104:   //  絶対値
   105:   public QFP abs () {
   106:     if ((this.flg & QFP_M) != 0) {  //-x,-0,-Inf
   107:       this.flg ^= QFP_M;
   108:       this.dvl = -this.dvl;
   109:       this.cvl = -this.cvl;
   110:     }
   111:     return this;
   112:   }  //qfp.abs()
   113:   public QFP abs (QFP x) {
   114:     if ((x.flg & QFP_M) != 0) {  //-x,-0,-Inf
   115:       this.flg = x.flg ^ QFP_M;
   116:       this.epp = x.epp;
   117:       this.dvl = -x.dvl;
   118:       this.cvl = -x.cvl;
   119:     } else {
   120:       this.flg = x.flg;
   121:       this.epp = x.epp;
   122:       this.dvl = x.dvl;
   123:       this.cvl = x.cvl;
   124:     }
   125:     return this;
   126:   }  //qfp.abs(QFP)
   127: 
   128:   //qfp = qfp.add (y)
   129:   //qfp = qfp.add (x, y)
   130:   //  加算
   131:   //  フラグ
   132:   //    (NaN)+(NaN,+Inf,0<y,+0,-0,y<0,-Inf)=NaN
   133:   //    (NaN,+Inf,0<x,+0,-0,x<0,-Inf)+(NaN)=NaN
   134:   //    (+Inf)+(-Inf)=NaN,OE
   135:   //    (-Inf)+(+Inf)=NaN,OE
   136:   //    (+Inf)+(+Inf,0<y,+0,-0,y<0)=+Inf
   137:   //    (-Inf)+(0<y,+0,-0,y<0,-Inf)=-Inf
   138:   //    (+Inf,0<x,+0,-0,x<0)+(+Inf)=+Inf
   139:   //    (0<x,+0,-0,x<0,-Inf)+(-Inf)=-Inf
   140:   //    (+0)+(+0)=+0
   141:   //    (+0)+(-0)=+0
   142:   //    (-0)+(+0)=+0
   143:   //    (-0)+(-0)=-0
   144:   //    (+0,-0)+(0<y,y<0)=y
   145:   //    (0<x,x<0)+(+0,-0)=x
   146:   //    (0<x,x<0)+(0<y,y<0)=z
   147:   //    +------+------+------+------+------+------+------+
   148:   //    |      | +Inf |  0<y |   +0 |   -0 |  y<0 | -Inf |
   149:   //    +------+------+------+------+------+------+------+
   150:   //    | +Inf | +Inf | +Inf | +Inf | +Inf | +Inf |  NaN |
   151:   //    |  0<x | +Inf |    z |    x |    x |    z | -Inf |
   152:   //    |   +0 | +Inf |    y |   +0 |   +0 |    y | -Inf |
   153:   //    |   -0 | +Inf |    y |   +0 |   -0 |    y | -Inf |
   154:   //    |  x<0 | +Inf |    z |    x |    x |    z | -Inf |
   155:   //    | -Inf |  NaN | -Inf | -Inf | -Inf | -Inf | -Inf |
   156:   //    +------+------+------+------+------+------+------+
   157:   private static final int[] QFP_ADD_FLG = {
   158:     //  perl -e "use strict;my($M,$Z,$I,$N,$OE,$DZ)=(0x08000000,0x04000000,0x02000000,0x01000000,0x00002000,0x00000400);for my$a(0..15){my$x=$a<<24;for my$b(0..15){my$y=$b<<24;($b&7)==0 and print'      ';printf('0x%08x,',0b11101010>>($x>>24&7)&1||0b11101010>>($y>>24&7)&1?$N:($x&$y)&$I&&$x!=$y?$N|$OE:($x|$y)&$I?$x&$I?$x:$y:($x&$y)&$Z?($x&$y)&$M|$Z:$x&$Z?0xc0000000:$y&$Z?0x80000000:0);($b&7)==7 and print chr 10;}}"
   159:     0x00000000,0x01000000,0x02000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
   160:     0x00000000,0x01000000,0x0a000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
   161:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   162:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   163:     0x02000000,0x01000000,0x02000000,0x01000000,0x02000000,0x01000000,0x01000000,0x01000000,
   164:     0x02000000,0x01000000,0x01002000,0x01000000,0x02000000,0x01000000,0x01000000,0x01000000,
   165:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   166:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   167:     0xc0000000,0x01000000,0x02000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
   168:     0xc0000000,0x01000000,0x0a000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
   169:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   170:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   171:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   172:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   173:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   174:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   175:     0x00000000,0x01000000,0x02000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
   176:     0x00000000,0x01000000,0x0a000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
   177:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   178:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   179:     0x0a000000,0x01000000,0x01002000,0x01000000,0x0a000000,0x01000000,0x01000000,0x01000000,
   180:     0x0a000000,0x01000000,0x0a000000,0x01000000,0x0a000000,0x01000000,0x01000000,0x01000000,
   181:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   182:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   183:     0xc0000000,0x01000000,0x02000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
   184:     0xc0000000,0x01000000,0x0a000000,0x01000000,0x0c000000,0x01000000,0x01000000,0x01000000,
   185:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   186:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   187:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   188:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   189:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   190:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   191:   };
   192:   public QFP add (QFP y) {
   193:     int xf = this.flg;
   194:     int xe = this.epp;
   195:     double xd = this.dvl;
   196:     double xc = this.cvl;
   197:     int yf = y.flg;
   198:     int ye = y.epp;
   199:     double yd = y.dvl;
   200:     double yc = y.cvl;
   201:     //フラグを設定する
   202:     int zf = QFP_ADD_FLG[xf >>> 24 - 4 | yf >>> 24];
   203:     if ((zf & (0xc0000000 | QFP_ZIN)) != 0) {  //x,y,±0,±Inf,NaN
   204:       if (0 <= zf) {  //±0,±Inf,NaN
   205:         QFP.fpsr |= zf & QFP_EXC;
   206:         this.flg = zf & QFP_MZIN;
   207:       } else if (zf << 1 < 0) {  //y
   208:         this.flg = yf;
   209:         this.epp = ye;
   210:         this.dvl = yd;
   211:         this.cvl = yc;
   212:         //} else {  //x
   213:         //  this.flg = xf;
   214:         //  this.epp = xe;
   215:         //  this.dvl = xd;
   216:         //  this.cvl = xc;
   217:       }
   218:       return this;
   219:     }
   220:     //両方±0,±Inf,NaN以外
   221:     //指数部の下駄を合わせる
   222:     if (xe != ye) {
   223:       if (xe + 128 == ye) {
   224:         yd = Math.scalb (yd, 128);
   225:         yc = Math.scalb (yc, 128);
   226:       } else if (xe - 128 == ye) {
   227:         yd = Math.scalb (yd, -128);
   228:         yc = Math.scalb (yc, -128);
   229:       } else {
   230:         if (xe < ye) {  //yが大きすぎるときy
   231:           this.flg = yf;
   232:           this.epp = ye;
   233:           this.dvl = yd;
   234:           this.cvl = yc;
   235:           //} else {  //xが大きすぎるときx
   236:           //  this.flg = xf;
   237:           //  this.epp = xe;
   238:           //  this.dvl = xd;
   239:           //  this.cvl = xc;
   240:         }
   241:         return this;
   242:       }
   243:     }
   244:     //加減算を行う
   245:     {
   246:       double t1 = xd + yd;
   247:       double t2 = xd - t1;
   248:       t2 = (((xd - (t1 + t2)) + (t2 + yd)) + xc) + yc;
   249:       xd = t1 + t2;
   250:       xc = (t1 - xd) + t2;
   251:     }
   252:     //フラグを設定する
   253:     //  両方±0,±Inf,NaN以外の加減算で結果がNaNになることはない
   254:     //  両方絶対値が2^192未満の加減算で結果が±Infになることはない
   255:     if (xd == 0.0) {  //±0
   256:       this.flg = QFP_P | QFP_Z;  //+0
   257:       return this;
   258:     }
   259:     this.flg = (int) (Double.doubleToLongBits (xd) >>> 36) & QFP_M;
   260:     //指数部の下駄を設定する
   261:     int e = Math.getExponent (xd) - QFP_GETA_BASE & -QFP_GETA_SIZE;
   262:     if (e != 0) {
   263:       xe += e;
   264:       xd = Math.scalb (xd, -e);
   265:       xc = Math.scalb (xc, -e);
   266:     }
   267:     //結果を格納する
   268:     this.epp = xe;
   269:     this.dvl = xd;
   270:     this.cvl = xc;
   271:     return this;
   272:   }  //qfp.add(QFP)
   273:   public QFP add (QFP x, QFP y) {
   274:     int xf = x.flg;
   275:     int xe = x.epp;
   276:     double xd = x.dvl;
   277:     double xc = x.cvl;
   278:     int yf = y.flg;
   279:     int ye = y.epp;
   280:     double yd = y.dvl;
   281:     double yc = y.cvl;
   282:     //フラグを設定する
   283:     int zf = QFP_ADD_FLG[xf >>> 24 - 4 | yf >>> 24];
   284:     if ((zf & (0xc0000000 | QFP_ZIN)) != 0) {  //x,y,±0,±Inf,NaN
   285:       if (0 <= zf) {  //±0,±Inf,NaN
   286:         QFP.fpsr |= zf & QFP_EXC;
   287:         this.flg = zf & QFP_MZIN;
   288:       } else if (zf << 1 < 0) {  //y
   289:         this.flg = yf;
   290:         this.epp = ye;
   291:         this.dvl = yd;
   292:         this.cvl = yc;
   293:       } else {  //x
   294:         this.flg = xf;
   295:         this.epp = xe;
   296:         this.dvl = xd;
   297:         this.cvl = xc;
   298:       }
   299:       return this;
   300:     }
   301:     //両方±0,±Inf,NaN以外
   302:     //指数部の下駄を合わせる
   303:     if (xe != ye) {
   304:       if (xe + 128 == ye) {
   305:         yd = Math.scalb (yd, 128);
   306:         yc = Math.scalb (yc, 128);
   307:       } else if (xe - 128 == ye) {
   308:         yd = Math.scalb (yd, -128);
   309:         yc = Math.scalb (yc, -128);
   310:       } else {
   311:         if (xe < ye) {  //yが大きすぎるときy
   312:           this.flg = yf;
   313:           this.epp = ye;
   314:           this.dvl = yd;
   315:           this.cvl = yc;
   316:         } else {  //xが大きすぎるときx
   317:           this.flg = xf;
   318:           this.epp = xe;
   319:           this.dvl = xd;
   320:           this.cvl = xc;
   321:         }
   322:         return this;
   323:       }
   324:     }
   325:     //加減算を行う
   326:     {
   327:       double t1 = xd + yd;
   328:       double t2 = xd - t1;
   329:       t2 = (((xd - (t1 + t2)) + (t2 + yd)) + xc) + yc;
   330:       xd = t1 + t2;
   331:       xc = (t1 - xd) + t2;
   332:     }
   333:     //フラグを設定する
   334:     //  両方±0,±Inf,NaN以外の加減算で結果がNaNになることはない
   335:     //  両方絶対値が2^192未満の加減算で結果が±Infになることはない
   336:     if (xd == 0.0) {  //±0
   337:       this.flg = QFP_P | QFP_Z;  //+0
   338:       return this;
   339:     }
   340:     this.flg = (int) (Double.doubleToLongBits (xd) >>> 36) & QFP_M;
   341:     //指数部の下駄を設定する
   342:     int e = Math.getExponent (xd) - QFP_GETA_BASE & -QFP_GETA_SIZE;
   343:     if (e != 0) {
   344:       xe += e;
   345:       xd = Math.scalb (xd, -e);
   346:       xc = Math.scalb (xc, -e);
   347:     }
   348:     //結果を格納する
   349:     this.epp = xe;
   350:     this.dvl = xd;
   351:     this.cvl = xc;
   352:     return this;
   353:   }  //qfp.add(QFP,QFP)
   354: 
   355:   //qfp = qfp.atanh (x)
   356:   private static final QFP[] QFP_ATH_T = {
   357:     //  1-2/(2^n*sqrt(2)+1) == (2^(2*n+1)+1-2^(n+1)*sqrt(2))/(2^(2*n+1)-1)
   358:     //  echo read("efp.gp");for(n=0,106,qfpmem(eval("1-2/(2^n*sqrt(2)+1)"),Str("1-2/(2^",n,"*sqrt(2)+1)"))) | gp-2.7 -q
   359:     new QFP (QFP_P, 0, 0x1.5f619980c4337p-3, -0x1.165f626cdd52bp-60),  //<1-2/(2^0*sqrt(2)+1)=0.17157287525380990239662255158060
   360:     new QFP (QFP_P, 0, 0x1.e90df15b89be3p-2, 0x1.1f99b9abc530dp-56),  //>1-2/(2^1*sqrt(2)+1)=0.47759225007251711497046358616589
   361:     new QFP (QFP_P, 0, 0x1.662c704e79f12p-1, 0x1.282af86097e1cp-55),  //>1-2/(2^2*sqrt(2)+1)=0.69955779035533030998666097439750
   362:     new QFP (QFP_P, 0, 0x1.acd734cfa2558p-1, -0x1.d43396def4576p-55),  //<1-2/(2^3*sqrt(2)+1)=0.83757939371677542692262189301295
   363:     new QFP (QFP_P, 0, 0x1.d4a917bee0f8ep-1, 0x1.a848f4c92e287p-56),  //>1-2/(2^4*sqrt(2)+1)=0.91535257535041283451731107793599
   364:     new QFP (QFP_P, 0, 0x1.e9dc9d2d2669p-1, 0x1.1f0d824ddcde5p-58),  //>1-2/(2^5*sqrt(2)+1)=0.95676127601764627106824226753814
   365:     new QFP (QFP_P, 0, 0x1.f4cf57477a9dfp-1, -0x1.58d05e5b8e13cp-57),  //>1-2/(2^6*sqrt(2)+1)=0.97814438579126405002483016033466
   366:     new QFP (QFP_P, 0, 0x1.fa5fcd25fa7dp-1, 0x1.d5e0dfd5720f4p-57),  //>1-2/(2^7*sqrt(2)+1)=0.98901215637783403019827166620693
   367:     new QFP (QFP_P, 0, 0x1.fd2deaca257dap-1, -0x1.b4d133e623a68p-55),  //>1-2/(2^8*sqrt(2)+1)=0.99449094503028873919489082537865
   368:     new QFP (QFP_P, 0, 0x1.fe9675ec66c9dp-1, 0x1.0314f0b167a78p-61),  //<1-2/(2^9*sqrt(2)+1)=0.99724167357216553275215115146172
   369:     new QFP (QFP_P, 0, 0x1.ff4b1b0724de6p-1, -0x1.4f6460da8573bp-55),  //>1-2/(2^10*sqrt(2)+1)=0.99861988508422135618277088368592
   370:     new QFP (QFP_P, 0, 0x1.ffa58585b10e2p-1, -0x1.e812fa2f88614p-56),  //<1-2/(2^11*sqrt(2)+1)=0.99930970437028696214760189421028
   371:     new QFP (QFP_P, 0, 0x1.ffd2c0c31c61fp-1, -0x1.82bc0859369f6p-55),  //>1-2/(2^12*sqrt(2)+1)=0.99965479261135554959522980932001
   372:     new QFP (QFP_P, 0, 0x1.ffe95fe196accp-1, 0x1.991433745dacbp-59),  //<1-2/(2^13*sqrt(2)+1)=0.99982738140837446316402328234885
   373:     new QFP (QFP_P, 0, 0x1.fff4afd0cc65ep-1, 0x1.4035d2da50783p-58),  //>1-2/(2^14*sqrt(2)+1)=0.99991368697937921695563957694829
   374:     new QFP (QFP_P, 0, 0x1.fffa57e06654ep-1, 0x1.9b510deb5e06ep-58),  //<1-2/(2^15*sqrt(2)+1)=0.99995684255842732237712901779679
   375:     new QFP (QFP_P, 0, 0x1.fffd2bee332ebp-1, -0x1.82787a0eb7631p-57),  //<1-2/(2^16*sqrt(2)+1)=0.99997842104639055378876556861730
   376:     new QFP (QFP_P, 0, 0x1.fffe95f69997ep-1, -0x1.fbebab6e35a44p-58),  //>1-2/(2^17*sqrt(2)+1)=0.99998921046498855802453226317115
   377:     new QFP (QFP_P, 0, 0x1.ffff4afb2cccp-1, 0x1.e5d3881a1907dp-59),  //>1-2/(2^18*sqrt(2)+1)=0.99999460521794248154009253151716
   378:     new QFP (QFP_P, 0, 0x1.ffffa57d8e66p-1, 0x1.2de46d4fa2f96p-56),  //>1-2/(2^19*sqrt(2)+1)=0.99999730260533327668252503098688
   379:     new QFP (QFP_P, 0, 0x1.ffffd2bec533p-1, 0x1.71c6418321567p-57),  //<1-2/(2^20*sqrt(2)+1)=0.99999865130175714547943879214761
   380:     new QFP (QFP_P, 0, 0x1.ffffe95f62198p-1, 0x1.82beb77000cdfp-58),  //>1-2/(2^21*sqrt(2)+1)=0.99999932565065119929426999557121
   381:     new QFP (QFP_P, 0, 0x1.fffff4afb0eccp-1, 0x1.86fcd50738a9dp-59),  //>1-2/(2^22*sqrt(2)+1)=0.99999966282526875625702343002444
   382:     new QFP (QFP_P, 0, 0x1.fffffa57d86e6p-1, 0x1.880c5c7086a0ap-60),  //<1-2/(2^23*sqrt(2)+1)=0.99999983141262016727739016874704
   383:     new QFP (QFP_P, 0, 0x1.fffffd2bec353p-1, 0x1.88503e4b4a1e5p-61),  //<1-2/(2^24*sqrt(2)+1)=0.99999991570630653092546549088408
   384:     new QFP (QFP_P, 0, 0x1.fffffe95f61a2p-1, -0x1.fcef3d927beep-55),  //>1-2/(2^25*sqrt(2)+1)=0.99999995785315237728436919619601
   385:     new QFP (QFP_P, 0, 0x1.ffffff4afb0cfp-1, -0x1.fcef3516408b1p-56),  //>1-2/(2^26*sqrt(2)+1)=0.99999997892657596659758669192677
   386:     new QFP (QFP_P, 0, 0x1.ffffffa57d867p-1, -0x1.fcef32f731b1fp-57),  //<1-2/(2^27*sqrt(2)+1)=0.99999998946328792778764299206308
   387:     new QFP (QFP_P, 0, 0x1.ffffffd2bec33p-1, 0x1.406219b212409p-55),  //>1-2/(2^28*sqrt(2)+1)=0.99999999473164395001603379788677
   388:     new QFP (QFP_P, 0, 0x1.ffffffe95f61ap-1, -0x1.7fcef324d7d0ep-55),  //<1-2/(2^29*sqrt(2)+1)=0.99999999736582197153856996069848
   389:     new QFP (QFP_P, 0, 0x1.fffffff4afb0dp-1, -0x1.8fcef324500d3p-56),  //<1-2/(2^30*sqrt(2)+1)=0.99999999868291098490192324407443
   390:     new QFP (QFP_P, 0, 0x1.fffffffa57d86p-1, 0x1.9a0c4336f478fp-55),  //<1-2/(2^31*sqrt(2)+1)=0.99999999934145549223412118775431
   391:     new QFP (QFP_P, 0, 0x1.fffffffd2bec3p-1, 0x1.990c4336f698p-56),  //>1-2/(2^32*sqrt(2)+1)=0.99999999967072774606285048527966
   392:     new QFP (QFP_P, 0, 0x1.fffffffe95f62p-1, -0x1.99dcef3242381p-55),  //<1-2/(2^33*sqrt(2)+1)=0.99999999983536387301787271548711
   393:     new QFP (QFP_P, 0, 0x1.ffffffff4afb1p-1, -0x1.99ecef32422f9p-56),  //>1-2/(2^34*sqrt(2)+1)=0.99999999991768193650554822595495
   394:     new QFP (QFP_P, 0, 0x1.ffffffffa57d8p-1, 0x1.9982c4336f74ap-55),  //<1-2/(2^35*sqrt(2)+1)=0.99999999995884096825192708003027
   395:     new QFP (QFP_P, 0, 0x1.ffffffffd2becp-1, 0x1.9981c4336f74cp-56),  //<1-2/(2^36*sqrt(2)+1)=0.99999999997942048412575178177833
   396:     new QFP (QFP_P, 0, 0x1.ffffffffe95f6p-1, 0x1.998144336f74dp-57),  //>1-2/(2^37*sqrt(2)+1)=0.99999999998971024206282295132996
   397:     new QFP (QFP_P, 0, 0x1.fffffffff4afbp-1, 0x1.998104336f74dp-58),  //>1-2/(2^38*sqrt(2)+1)=0.99999999999485512103139824077518
   398:     new QFP (QFP_P, 0, 0x1.fffffffffa57ep-1, -0x1.e667f1bcc908bp-55),  //>1-2/(2^39*sqrt(2)+1)=0.99999999999742756051569581166514
   399:     new QFP (QFP_P, 0, 0x1.fffffffffd2bfp-1, -0x1.e667f2bcc908bp-56),  //>1-2/(2^40*sqrt(2)+1)=0.99999999999871378025784707865196
   400:     new QFP (QFP_P, 0, 0x1.fffffffffe95fp-1, 0x1.86660330cdbddp-55),  //<1-2/(2^41*sqrt(2)+1)=0.99999999999935689012892333253082
   401:     new QFP (QFP_P, 0, 0x1.ffffffffff4bp-1, -0x1.3cccfe6f99211p-55),  //>1-2/(2^42*sqrt(2)+1)=0.99999999999967844506446161456663
   402:     new QFP (QFP_P, 0, 0x1.ffffffffffa58p-1, -0x1.3cccfe7399211p-56),  //>1-2/(2^43*sqrt(2)+1)=0.99999999999983922253223079435862
   403:     new QFP (QFP_P, 0, 0x1.ffffffffffd2cp-1, -0x1.3cccfe7599211p-57),  //>1-2/(2^44*sqrt(2)+1)=0.99999999999991961126611539394813
   404:     new QFP (QFP_P, 0, 0x1.ffffffffffe96p-1, -0x1.3cccfe7699211p-58),  //>1-2/(2^45*sqrt(2)+1)=0.99999999999995980563305769616627
   405:     new QFP (QFP_P, 0, 0x1.fffffffffff4bp-1, -0x1.3cccfe7719211p-59),  //>1-2/(2^46*sqrt(2)+1)=0.99999999999997990281652884788119
   406:     new QFP (QFP_P, 0, 0x1.fffffffffffa5p-1, 0x1.f619980c4536fp-55),  //<1-2/(2^47*sqrt(2)+1)=0.99999999999998995140826442389010
   407:     new QFP (QFP_P, 0, 0x1.fffffffffffd3p-1, -0x1.04f333f9dde48p-55),  //>1-2/(2^48*sqrt(2)+1)=0.99999999999999497570413221193243
   408:     new QFP (QFP_P, 0, 0x1.fffffffffffe9p-1, 0x1.7d86660310edcp-55),  //>1-2/(2^49*sqrt(2)+1)=0.99999999999999748785206610596306
   409:     new QFP (QFP_P, 0, 0x1.ffffffffffff5p-1, -0x1.413cccfe77912p-55),  //>1-2/(2^50*sqrt(2)+1)=0.99999999999999874392603305298074
   410:     new QFP (QFP_P, 0, 0x1.ffffffffffffap-1, 0x1.5f619980c4357p-55),  //>1-2/(2^51*sqrt(2)+1)=0.99999999999999937196301652649017
   411:     new QFP (QFP_P, 0, 0x1.ffffffffffffdp-1, 0x1.5f619980c4347p-56),  //>1-2/(2^52*sqrt(2)+1)=0.99999999999999968598150826324504
   412:     new QFP (QFP_P, 0, 0x1.fffffffffffffp-1, -0x1.a827999fcef3p-55),  //>1-2/(2^53*sqrt(2)+1)=0.99999999999999984299075413162251
   413:     new QFP (QFP_P, 0, 0x1.fffffffffffffp-1, 0x1.2bec333018867p-55),  //<1-2/(2^54*sqrt(2)+1)=0.99999999999999992149537706581125
   414:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bccp-55),  //>1-2/(2^55*sqrt(2)+1)=0.99999999999999996074768853290563
   415:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-56),  //<1-2/(2^56*sqrt(2)+1)=0.99999999999999998037384426645281
   416:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-57),  //<1-2/(2^57*sqrt(2)+1)=0.99999999999999999018692213322641
   417:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-58),  //<1-2/(2^58*sqrt(2)+1)=0.99999999999999999509346106661320
   418:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-59),  //<1-2/(2^59*sqrt(2)+1)=0.99999999999999999754673053330660
   419:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-60),  //<1-2/(2^60*sqrt(2)+1)=0.99999999999999999877336526665330
   420:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-61),  //<1-2/(2^61*sqrt(2)+1)=0.99999999999999999938668263332665
   421:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-62),  //<1-2/(2^62*sqrt(2)+1)=0.99999999999999999969334131666333
   422:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-63),  //<1-2/(2^63*sqrt(2)+1)=0.99999999999999999984667065833166
   423:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-64),  //<1-2/(2^64*sqrt(2)+1)=0.99999999999999999992333532916583
   424:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-65),  //<1-2/(2^65*sqrt(2)+1)=0.99999999999999999996166766458292
   425:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-66),  //<1-2/(2^66*sqrt(2)+1)=0.99999999999999999998083383229146
   426:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-67),  //<1-2/(2^67*sqrt(2)+1)=0.99999999999999999999041691614573
   427:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-68),  //<1-2/(2^68*sqrt(2)+1)=0.99999999999999999999520845807286
   428:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-69),  //<1-2/(2^69*sqrt(2)+1)=0.99999999999999999999760422903643
   429:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-70),  //<1-2/(2^70*sqrt(2)+1)=0.99999999999999999999880211451822
   430:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-71),  //<1-2/(2^71*sqrt(2)+1)=0.99999999999999999999940105725911
   431:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-72),  //<1-2/(2^72*sqrt(2)+1)=0.99999999999999999999970052862955
   432:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-73),  //<1-2/(2^73*sqrt(2)+1)=0.99999999999999999999985026431478
   433:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-74),  //<1-2/(2^74*sqrt(2)+1)=0.99999999999999999999992513215739
   434:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-75),  //<1-2/(2^75*sqrt(2)+1)=0.99999999999999999999996256607869
   435:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-76),  //<1-2/(2^76*sqrt(2)+1)=0.99999999999999999999998128303935
   436:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-77),  //<1-2/(2^77*sqrt(2)+1)=0.99999999999999999999999064151967
   437:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-78),  //<1-2/(2^78*sqrt(2)+1)=0.99999999999999999999999532075984
   438:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-79),  //<1-2/(2^79*sqrt(2)+1)=0.99999999999999999999999766037992
   439:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-80),  //<1-2/(2^80*sqrt(2)+1)=0.99999999999999999999999883018996
   440:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-81),  //<1-2/(2^81*sqrt(2)+1)=0.99999999999999999999999941509498
   441:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-82),  //<1-2/(2^82*sqrt(2)+1)=0.99999999999999999999999970754749
   442:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-83),  //<1-2/(2^83*sqrt(2)+1)=0.99999999999999999999999985377374
   443:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-84),  //<1-2/(2^84*sqrt(2)+1)=0.99999999999999999999999992688687
   444:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-85),  //<1-2/(2^85*sqrt(2)+1)=0.99999999999999999999999996344344
   445:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-86),  //<1-2/(2^86*sqrt(2)+1)=0.99999999999999999999999998172172
   446:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-87),  //<1-2/(2^87*sqrt(2)+1)=0.99999999999999999999999999086086
   447:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-88),  //<1-2/(2^88*sqrt(2)+1)=0.99999999999999999999999999543043
   448:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-89),  //<1-2/(2^89*sqrt(2)+1)=0.99999999999999999999999999771521
   449:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-90),  //<1-2/(2^90*sqrt(2)+1)=0.99999999999999999999999999885761
   450:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-91),  //<1-2/(2^91*sqrt(2)+1)=0.99999999999999999999999999942880
   451:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-92),  //<1-2/(2^92*sqrt(2)+1)=0.99999999999999999999999999971440
   452:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-93),  //<1-2/(2^93*sqrt(2)+1)=0.99999999999999999999999999985720
   453:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-94),  //<1-2/(2^94*sqrt(2)+1)=0.99999999999999999999999999992860
   454:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-95),  //<1-2/(2^95*sqrt(2)+1)=0.99999999999999999999999999996430
   455:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-96),  //<1-2/(2^96*sqrt(2)+1)=0.99999999999999999999999999998215
   456:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-97),  //<1-2/(2^97*sqrt(2)+1)=0.99999999999999999999999999999108
   457:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-98),  //<1-2/(2^98*sqrt(2)+1)=0.99999999999999999999999999999554
   458:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-99),  //<1-2/(2^99*sqrt(2)+1)=0.99999999999999999999999999999777
   459:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-100),  //<1-2/(2^100*sqrt(2)+1)=0.99999999999999999999999999999888
   460:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-101),  //<1-2/(2^101*sqrt(2)+1)=0.99999999999999999999999999999944
   461:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-102),  //<1-2/(2^102*sqrt(2)+1)=0.99999999999999999999999999999972
   462:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-103),  //<1-2/(2^103*sqrt(2)+1)=0.99999999999999999999999999999986
   463:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-104),  //<1-2/(2^104*sqrt(2)+1)=0.99999999999999999999999999999993
   464:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-105),  //<1-2/(2^105*sqrt(2)+1)=0.99999999999999999999999999999997
   465:     new QFP (QFP_P, 0, 0x1.0p0, -0x1.6a09e667f3bcdp-106),  //<1-2/(2^106*sqrt(2)+1)=0.99999999999999999999999999999998
   466:   };
   467:   //  echo read("efp.gp");eval("f(x)=atanh(x)");a=2*sqrt(2)-3;b=-a;n=31;qfpchebyshev("QFP_ATH_C",f,a,b,n) | gp-2.7 -q
   468:   private static final QFP QFP_ATH_C1 = new QFP (QFP_P, 0, 0x1.0p0, -0x1.f36f641e20d72p-114);  //>c1=1.00000000000000000000000000000
   469:   private static final QFP QFP_ATH_C3 = new QFP (QFP_P, 0, 0x1.5555555555555p-2, 0x1.5555555555611p-56);  //>c3=0.333333333333333333333333333334
   470:   private static final QFP QFP_ATH_C5 = new QFP (QFP_P, 0, 0x1.999999999999ap-3, -0x1.9999999a4175dp-57);  //>c5=0.199999999999999999999999998941
   471:   private static final QFP QFP_ATH_C7 = new QFP (QFP_P, 0, 0x1.2492492492492p-3, 0x1.24924b582faeep-57);  //>c7=0.142857142857142857142858053426
   472:   private static final QFP QFP_ATH_C9 = new QFP (QFP_P, 0, 0x1.c71c71c71c71cp-4, 0x1.c7140426f0f4ap-58);  //<c9=0.111111111111111111110664924582
   473:   private static final QFP QFP_ATH_C11 = new QFP (QFP_P, 0, 0x1.745d1745d1746p-4, -0x1.5fe638aeacf12p-59);  //<c11=0.0909090909090909092295808117708
   474:   private static final QFP QFP_ATH_C13 = new QFP (QFP_P, 0, 0x1.3b13b13b13b12p-4, -0x1.a3d69d8be81fp-58);  //>c13=0.0769230769230768939015654936276
   475:   private static final QFP QFP_ATH_C15 = new QFP (QFP_P, 0, 0x1.1111111111249p-4, 0x1.d1d90451b5d3fp-59);  //<c15=0.0666666666666709987679838034753
   476:   private static final QFP QFP_ATH_C17 = new QFP (QFP_P, 0, 0x1.e1e1e1e1d17b1p-5, -0x1.464b242b5513p-59);  //>c17=0.0588235294112985410236031906597
   477:   private static final QFP QFP_ATH_C19 = new QFP (QFP_P, 0, 0x1.af286bcf2dc59p-5, -0x1.87a381950827p-61);  //<c19=0.0526315789842832289529081789289
   478:   private static final QFP QFP_ATH_C21 = new QFP (QFP_P, 0, 0x1.8618605cae00ep-5, -0x1.338aeff0d76ebp-59);  //<c21=0.0476190454550647741781645777628
   479:   private static final QFP QFP_ATH_C23 = new QFP (QFP_P, 0, 0x1.642cb7cf14af4p-5, 0x1.e3d44dca87f29p-59);  //>c23=0.0434783544557343670267754616356
   480:   private static final QFP QFP_ATH_C25 = new QFP (QFP_P, 0, 0x1.47a7e890f5ccdp-5, 0x1.31e1771122396p-59);  //<c25=0.0399970571813184995379450688298
   481:   private static final QFP QFP_ATH_C27 = new QFP (QFP_P, 0, 0x1.2ff12aaefe92p-5, -0x1.6898bb0e3245dp-60);  //>c27=0.0371023019469751679228447470154
   482:   private static final QFP QFP_ATH_C29 = new QFP (QFP_P, 0, 0x1.129828d81e45cp-5, -0x1.e5836f6efc4bcp-62);  //<c29=0.0335198209513200835791642635031
   483:   private static final QFP QFP_ATH_C31 = new QFP (QFP_P, 0, 0x1.4cd3eb19356dep-5, -0x1.9b438b675ceacp-60);  //>c31=0.0406283942954084167863663627263
   484:   //  52738671016840858519654937824341/1298074214633706907132624082305024*x^31 + 174044861024189167268553365393105/5192296858534827628530496329220096*x^29 + 192646165843689783173132312173475/5192296858534827628530496329220096*x^27 + 51919148588299478172157465924043/1298074214633706907132624082305024*x^25 + 112876261627386639600856890638121/2596148429267413814265248164610048*x^23 + 123626110061379994152690835949845/2596148429267413814265248164610048*x^21 + 34159847777477682956032434436057/649037107316853453566312041152512*x^19 + 9544663342819264156052426631917/162259276829213363391578010288128*x^17 + 173076561951172167729046456327487/2596148429267413814265248164610048*x^15 + 6240741416508203917305250584545/81129638414606681695789005144064*x^13 + 118006746784882446282972010813559/1298074214633706907132624082305024*x^11 + 72115234146317050395967301846949/649037107316853453566312041152512*x^9 + 46359793379775246683308298435959/324518553658426726783156020576256*x^7 + 129807421463370690713262407542947/649037107316853453566312041152512*x^5 + 108172851219475575594385340192273/324518553658426726783156020576256*x^3 + 46768052394588893382517914646921052235912054765895/46768052394588893382517914646921056628989841375232*x
   485:   //  113.0359bit
   486:   public QFP atanh () {
   487:     return this.atanh (this);
   488:   }  //qfp.atanh()
   489:   public QFP atanh (QFP x) {
   490:     int xf = x.flg;
   491:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   492:       if ((xf & QFP_I) != 0) {  //±Inf
   493:         QFP.fpsr |= QFP_OE;
   494:         this.flg = QFP_N;  //atanh(±Inf)=NaN
   495:       } else {  //±0,NaN
   496:         this.flg = xf;  //atanh(±0)=±0,atanh(NaN)=NaN
   497:       }
   498:       return this;
   499:     }
   500:     //±0,±Inf,NaN以外
   501:     x = new QFP ().abs (x);  //|x|
   502:     {
   503:       int t = x.cmp1 ();
   504:       if (0 <= t) {  //1<=|x|
   505:         if (t == 0) {  //|x|==1
   506:           QFP.fpsr |= QFP_DZ;
   507:           this.flg = xf | QFP_I;
   508:         } else {  //1<|x|
   509:           QFP.fpsr |= QFP_OE;
   510:           this.flg = QFP_N;
   511:         }
   512:         return this;
   513:       }
   514:     }
   515:     //|x|<1
   516:     int n = 0;
   517:     if (QFP_ATH_T[0].lt (x)) {  //3-2*sqrt(2)<|x|
   518:       QFP t1 = new QFP ().negdec (x);  //t1=1-|x|
   519:       n = -Math.getExponent (t1.dvl);  //-floor(log2(1-|x|))。|x|<1なのでt1.eppは0
   520:       if (QFP_ATH_T[n].lt (x)) {  //1-2/(2^n*sqrt(2)+1)<|x|
   521:         n++;
   522:       }
   523:       t1.shl (n);  //t1=2^n*(1-|x|)
   524:       x.inc ();  //x=1+|x|
   525:       QFP t2 = new QFP ().add (x, t1);  //t2=1+|x|+2^n*(1-|x|)
   526:       x.sub (t1)  //x=1+|x|-2^n*(1-|x|)
   527:         .div (t2);  //x=(1+|x|-2^n*(1-|x|))/(1+|x|+2^n*(1-|x|))
   528:     }
   529:     QFP x2 = new QFP ().squ (x);  //x^2
   530:     this.mul (QFP_ATH_C31, x2)
   531:       .add (QFP_ATH_C29).mul (x2)
   532:         .add (QFP_ATH_C27).mul (x2)
   533:           .add (QFP_ATH_C25).mul (x2)
   534:             .add (QFP_ATH_C23).mul (x2)
   535:               .add (QFP_ATH_C21).mul (x2)
   536:                 .add (QFP_ATH_C19).mul (x2)
   537:                   .add (QFP_ATH_C17).mul (x2)
   538:                     .add (QFP_ATH_C15).mul (x2)
   539:                       .add (QFP_ATH_C13).mul (x2)
   540:                         .add (QFP_ATH_C11).mul (x2)
   541:                           .add (QFP_ATH_C9).mul (x2)
   542:                             .add (QFP_ATH_C7).mul (x2)
   543:                               .add (QFP_ATH_C5).mul (x2)
   544:                                 .add (QFP_ATH_C3).mul (x2)
   545:                                   .add (QFP_ATH_C1).mul (x);
   546:     if (n != 0) {
   547:       this.add (new QFP (n).mul (QFP_LN_2_2));  //+n*log(2)/2
   548:     }
   549:     return this.neg ((xf & QFP_M) != 0);
   550:   }  //qfp.atanh(QFP)
   551: 
   552:   //s = x.cmp (y)
   553:   public int cmp (QFP y) {
   554:     int xf = this.flg;
   555:     int yf = y.flg;
   556:     if (((xf | yf) & QFP_ZIN) != 0) {  //どちらかが±0,±Inf,NaN
   557:       return EFPBox.EFP_CMP_TABLE[xf >>> 24] << (yf >>> 24 - 1) >> 30;
   558:     }
   559:     //両方±0,±Inf,NaN以外
   560:     int s = (xf & QFP_M) != 0 ? -1 : 1;
   561:     return (xf != yf ? s :
   562:             this.epp != y.epp ? this.epp < y.epp ? -s : s :
   563:             this.dvl != y.dvl ? this.dvl < y.dvl ? -s : s :
   564:             this.cvl != y.cvl ? this.cvl < y.cvl ? -s : s :
   565:             0);
   566:   }  //qfp.cmp(QFP)
   567: 
   568:   //s = x.cmp1 ()
   569:   //  s=x<=>1
   570:   //  1との比較
   571:   public int cmp1 () {
   572:     return ((this.flg & QFP_MZIN) != 0 ?  //-x,±0,±Inf,NaN
   573:             (this.flg & QFP_N) != 0 ? 0 :  //NaN
   574:             (this.flg & (QFP_M | QFP_Z)) != 0 ? -1 :  //(-Inf,-x,-0,+0)<1
   575:             1  //1<+Inf
   576:             :  //+x
   577:             this.epp != 0 ? this.epp < 0 ? -1 : 1 :
   578:             this.dvl != 1.0 ? this.dvl < 1.0 ? -1 : 1 :
   579:             this.cvl < 0.0 ? -1 : 0.0 < this.cvl ? 1 : 0);
   580:   }  //qfp.cmp1()
   581: 
   582:   //x = x.dec ()
   583:   //y = y.dec (x)
   584:   //  y=x-1
   585:   //  デクリメント
   586:   public QFP dec () {
   587:     return this.sub (QFP_ONE);
   588:   }  //qfp.dec()
   589:   public QFP dec (QFP x) {
   590:     return this.sub (x, QFP_ONE);
   591:   }  //qfp.dec(QFP)
   592: 
   593:   //qfp = qfp.div (y)
   594:   //qfp = qfp.div (x, y)
   595:   //  除算
   596:   //  フラグ
   597:   //    (NaN)/(NaN,+Inf,0<y,+0,-0,y<0,-Inf)=NaN
   598:   //    (NaN,+Inf,0<x,+0,-0,x<0,-Inf)/(NaN)=NaN
   599:   //    (+0,-0)/(+0,-0)=NaN,OE
   600:   //    (+Inf,-Inf)/(+Inf,-Inf)=NaN,OE
   601:   //    (+Inf,0<x)/(+0)=+Inf,DZ
   602:   //    (x<0,-Inf)/(+0)=-Inf,DZ
   603:   //    (+Inf,0<x)/(-0)=-Inf,DZ
   604:   //    (x<0,-Inf)/(-0)=+Inf,DZ
   605:   //    (+Inf)/(0<y)=+Inf
   606:   //    (+Inf)/(y<0)=-Inf
   607:   //    (-Inf)/(0<y)=-Inf
   608:   //    (-Inf)/(y<0)=+Inf
   609:   //    (+0)/(0<y)=+0
   610:   //    (-0)/(0<y)=-0
   611:   //    (+0)/(y<0)=-0
   612:   //    (-0)/(y<0)=+0
   613:   //    (0<x,+0)/(+Inf)=+0
   614:   //    (-0,x<0)/(+Inf)=-0
   615:   //    (0<x,+0)/(-Inf)=-0
   616:   //    (-0,x<0)/(-Inf)=+0
   617:   //    (0<x)/(0<y)=+z
   618:   //    (0<x)/(y<0)=-z
   619:   //    (x<0)/(0<y)=-z
   620:   //    (x<0)/(y<0)=+z
   621:   //    +------+------+------+------+------+------+------+
   622:   //    |      | +Inf |  0<y |   +0 |   -0 |  y<0 | -Inf |
   623:   //    +------+------+------+------+------+------+------+
   624:   //    | +Inf |  NaN | +Inf | +Inf | -Inf | -Inf |  NaN |
   625:   //    |  0<x |   +0 |   +z | +Inf | -Inf |   -z |   -0 |
   626:   //    |   +0 |   +0 |   +0 |  NaN |  NaN |   -0 |   -0 |
   627:   //    |   -0 |   -0 |   -0 |  NaN |  NaN |   +0 |   +0 |
   628:   //    |  x<0 |   -0 |   -z | -Inf | +Inf |   +z |   +0 |
   629:   //    | -Inf |  NaN | -Inf | -Inf | +Inf | +Inf |  NaN |
   630:   //    +------+------+------+------+------+------+------+
   631:   private static final int[] QFP_DIV_FLG = {
   632:     //  perl -e "use strict;my($M,$Z,$I,$N,$OE,$DZ)=(0x08000000,0x04000000,0x02000000,0x01000000,0x00002000,0x00000400);for my$a(0..15){my$x=$a<<24;for my$b(0..15){my$y=$b<<24;($b&7)==0 and print'      ';printf('0x%08x,',0b11101010>>($x>>24&7)&1||0b11101010>>($y>>24&7)&1?$N:($x&$y)&($Z|$I)?$N|$OE:($x^$y)&$M|($x&$I&&$y&$Z?$I|$DZ:$x&$I?$I:$x&$Z||$y&$I?$Z:0));($b&7)==7 and print chr 10;}}"
   633:     0x00000000,0x01000000,0x04000000,0x01000000,0x00000000,0x01000000,0x01000000,0x01000000,
   634:     0x08000000,0x01000000,0x0c000000,0x01000000,0x08000000,0x01000000,0x01000000,0x01000000,
   635:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   636:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   637:     0x02000000,0x01000000,0x01002000,0x01000000,0x02000400,0x01000000,0x01000000,0x01000000,
   638:     0x0a000000,0x01000000,0x01002000,0x01000000,0x0a000400,0x01000000,0x01000000,0x01000000,
   639:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   640:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   641:     0x04000000,0x01000000,0x04000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   642:     0x0c000000,0x01000000,0x0c000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   643:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   644:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   645:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   646:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   647:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   648:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   649:     0x08000000,0x01000000,0x0c000000,0x01000000,0x08000000,0x01000000,0x01000000,0x01000000,
   650:     0x00000000,0x01000000,0x04000000,0x01000000,0x00000000,0x01000000,0x01000000,0x01000000,
   651:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   652:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   653:     0x0a000000,0x01000000,0x01002000,0x01000000,0x0a000400,0x01000000,0x01000000,0x01000000,
   654:     0x02000000,0x01000000,0x01002000,0x01000000,0x02000400,0x01000000,0x01000000,0x01000000,
   655:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   656:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   657:     0x0c000000,0x01000000,0x0c000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   658:     0x04000000,0x01000000,0x04000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   659:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   660:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   661:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   662:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   663:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   664:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   665:   };
   666:   public QFP div (QFP y) {
   667:     int zf = QFP_DIV_FLG[this.flg >>> 24 - 4 | y.flg >>> 24];
   668:     if ((zf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   669:       QFP.fpsr |= zf & QFP_EXC;
   670:       this.flg = zf & QFP_MZIN;
   671:       return this;
   672:     }
   673:     return this.mul (new QFP ().rcp (y));
   674:   }  //qfp.div(QFP)
   675:   public QFP div (QFP x, QFP y) {
   676:     int zf = QFP_DIV_FLG[x.flg >>> 24 - 4 | y.flg >>> 24];
   677:     if ((zf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   678:       QFP.fpsr |= zf & QFP_EXC;
   679:       this.flg = zf & QFP_MZIN;
   680:       return this;
   681:     }
   682:     return this.mul (x, new QFP ().rcp (y));
   683:   }  //qfp.div(QFP,QFP)
   684: 
   685:   //qfp = qfp.div2 ()
   686:   //qfp = qfp.div2 (x)
   687:   //  1/2倍
   688:   public QFP div2 () {
   689:     if ((this.flg & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   690:       return this;
   691:     }
   692:     //±0,±Inf,NaN以外
   693:     //1/2倍する
   694:     this.dvl *= 0.5;
   695:     this.cvl *= 0.5;
   696:     //指数部の下駄を設定する
   697:     int e = Math.getExponent (this.dvl) - QFP_GETA_BASE & -QFP_GETA_SIZE;
   698:     if (e != 0) {
   699:       this.epp += e;
   700:       this.dvl = Math.scalb (this.dvl, -e);
   701:       this.cvl = Math.scalb (this.cvl, -e);
   702:     }
   703:     return this;
   704:   }  //qfp.div2()
   705:   public QFP div2 (QFP x) {
   706:     this.flg = x.flg;
   707:     if ((this.flg & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   708:       return this;
   709:     }
   710:     //±0,±Inf,NaN以外
   711:     //1/2倍する
   712:     this.epp = x.epp;
   713:     this.dvl = x.dvl * 0.5;
   714:     this.cvl = x.cvl * 0.5;
   715:     //指数部の下駄を設定する
   716:     int e = Math.getExponent (this.dvl) - QFP_GETA_BASE & -QFP_GETA_SIZE;
   717:     if (e != 0) {
   718:       this.epp += e;
   719:       this.dvl = Math.scalb (this.dvl, -e);
   720:       this.cvl = Math.scalb (this.cvl, -e);
   721:     }
   722:     return this;
   723:   }  //qfp.div2(QFP)
   724: 
   725:   //qfp.dump ()
   726:   //qfp.dump (value)
   727:   //  ダンプ
   728:   public void dump () {
   729:     this.dump (null, false);
   730:   }  //qfp.dump()
   731:   public void dump (String name) {
   732:     this.dump (name, false);
   733:   }  //qfp.dump(String)
   734:   public void dump (boolean value) {
   735:     this.dump (null, value);
   736:   }  //qfp.dump(boolean)
   737:   public void dump (String name, boolean value) {
   738:     if (name != null) {
   739:       System.out.printf ("%s=", name);
   740:     }
   741:     if ((this.flg & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   742:       System.out.printf ("QFP{flg:0x%08x", this.flg);
   743:       if (value) {
   744:         System.out.printf (",val:%s", this.toString ());
   745:       }
   746:       System.out.println ("}");
   747:     } else {  //±0,±Inf,NaN以外
   748:       System.out.printf ("QFP{flg:0x%08x,epp:%d,dvl:%.16g=0x%016x,cvl:%.16g=0x%016x",
   749:                          this.flg,
   750:                          this.epp,
   751:                          this.dvl, Double.doubleToLongBits (this.dvl),
   752:                          this.cvl, Double.doubleToLongBits (this.cvl));
   753:       if (value) {
   754:         System.out.printf (",val:%s", this.toString ());
   755:       }
   756:       System.out.println ("}");
   757:     }
   758:   }  //qfp.dump()
   759: 
   760:   //i = qfp.getd ()
   761:   //  double取得
   762:   public double getd () {
   763:     int xf = this.flg;
   764:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   765:       return ((xf & QFP_Z) != 0 ? (xf & QFP_M) == 0 ? 0.0 : -0.0 :  //±0
   766:               (xf & QFP_I) != 0 ? (xf & QFP_M) == 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY :  //±Inf
   767:               Double.NaN);  //NaN
   768:     }
   769:     //±0,±Inf,NaN以外
   770:     double xd = this.dvl;
   771:     if (this.epp != 0) {
   772:       xd = Math.scalb (xd, this.epp);
   773:       if (Double.isInfinite (xd)) {  //オーバーフローした
   774:         QFP.fpsr |= QFP_OF;
   775:       } else if (xd == 0.0) {  //アンダーフローした
   776:         QFP.fpsr |= QFP_UF;
   777:       }
   778:     }
   779:     return xd;
   780:   }  //qfp.getd()
   781: 
   782:   //f = qfp.getf ()
   783:   //  float取得
   784:   public float getf () {
   785:     int xf = this.flg;
   786:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   787:       return ((xf & QFP_Z) != 0 ? (xf & QFP_M) == 0 ? 0.0F : -0.0F :  //±0
   788:               (xf & QFP_I) != 0 ? (xf & QFP_M) == 0 ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY :  //±Inf
   789:               Float.NaN);  //NaN
   790:     }
   791:     //±0,±Inf,NaN以外
   792:     return (float) this.dvl;
   793:   }  //qfp.getf()
   794: 
   795:   //i = qfp.geti ()
   796:   //  int取得
   797:   public int geti () {
   798:     int xf = this.flg;
   799:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   800:       if ((xf & QFP_Z) != 0) {  //±0
   801:         return 0;
   802:       }
   803:       QFP.fpsr |= QFP_OE;
   804:       return ((xf & QFP_I) != 0 ? (xf & QFP_M) == 0 ? 0x7fffffff : 0x80000000 :  //±Inf
   805:               0);  //NaN
   806:     }
   807:     //±0,±Inf,NaN以外
   808:     if (this.epp != 0) {
   809:       if (this.epp < 0) {
   810:         return 0;
   811:       }
   812:       QFP.fpsr |= QFP_OE;
   813:       return (xf & QFP_M) == 0 ? 0x7fffffff : 0x80000000;
   814:     }
   815:     double xd = this.dvl;
   816:     double xc = this.cvl;
   817:     //上位と下位の符号を合わせる
   818:     if (xd * xc < 0.0) {
   819:       double u = Math.ulp (xd) * Math.signum (xc);
   820:       xd += u;
   821:       xc -= u;
   822:     }
   823:     if (xd <= (double) Integer.MIN_VALUE - 1.0 || (double) Integer.MAX_VALUE + 1.0 <= xd) {
   824:       QFP.fpsr |= QFP_OE;
   825:       return 0.0 <= xd ? 0x7fffffff : 0x80000000;
   826:     }
   827:     return (int) xd;
   828:   }  //qfp.geti()
   829: 
   830:   //l = qfp.getl ()
   831:   //  long取得
   832:   public long getl () {
   833:     int xf = this.flg;
   834:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   835:       if ((xf & QFP_Z) != 0) {  //±0
   836:         return 0L;
   837:       }
   838:       QFP.fpsr |= QFP_OE;
   839:       return ((xf & QFP_I) != 0 ? (xf & QFP_M) == 0 ? 0x7fffffffffffffffL : 0x8000000000000000L :  //±Inf
   840:               0L);  //NaN
   841:     }
   842:     //±0,±Inf,NaN以外
   843:     if (this.epp != 0) {
   844:       if (this.epp < 0) {
   845:         return 0;
   846:       }
   847:       QFP.fpsr |= QFP_OE;
   848:       return (xf & QFP_M) == 0 ? 0x7fffffffffffffffL : 0x8000000000000000L;
   849:     }
   850:     double xd = this.dvl;
   851:     double xc = this.cvl;
   852:     //上位と下位の符号を合わせる
   853:     if (xd * xc < 0.0) {
   854:       double u = Math.ulp (xd) * Math.signum (xc);
   855:       xd += u;
   856:       xc -= u;
   857:     }
   858:     if (xd <= (double) Long.MIN_VALUE - 1.0 || (double) Long.MAX_VALUE + 1.0 <= xd) {
   859:       QFP.fpsr |= QFP_OE;
   860:       return 0.0 <= xd ? 0x7fffffffffffffffL : 0x8000000000000000L;
   861:     }
   862:     return (long) xd + (long) xc;
   863:   }  //qfp.getl()
   864: 
   865:   //x = x.inc ()
   866:   //  インクリメント
   867:   public QFP inc () {
   868:     return this.add (QFP_ONE);
   869:   }  //qfp.inc()
   870: 
   871:   //b = x.lt (y)
   872:   public boolean lt (QFP y) {
   873:     int xf = this.flg;
   874:     int yf = y.flg;
   875:     if (((xf | yf) & QFP_ZIN) != 0) {  //どちらかが±0,±Inf,NaN
   876:       return EFPBox.EFP_LT_TABLE[xf >>> 24] << (yf >>> 24 - 1) < 0;
   877:     }
   878:     //両方±0,±Inf,NaN以外
   879:     int s = (xf & QFP_M) != 0 ? -1 : 1;
   880:     return (xf != yf ? s :
   881:             this.epp != y.epp ? this.epp < y.epp ? -s : s :
   882:             this.dvl != y.dvl ? this.dvl < y.dvl ? -s : s :
   883:             this.cvl != y.cvl ? this.cvl < y.cvl ? -s : s :
   884:             0) < 0;
   885:   }  //qfp.lt(QFP)
   886: 
   887:   //qfp = qfp.mul (y)
   888:   //qfp = qfp.mul (x, y)
   889:   //  乗算
   890:   //  フラグ
   891:   //    (NaN)*(NaN,+Inf,0<y,+0,-0,y<0,-Inf)=NaN
   892:   //    (NaN,+Inf,0<x,+0,-0,x<0,-Inf)*(NaN)=NaN
   893:   //    (+0,-0)*(+Inf,-Inf)=NaN,OE
   894:   //    (+Inf,-Inf)*(+0,-0)=NaN,OE
   895:   //    (+Inf)*(+Inf,0<y)=+Inf
   896:   //    (+Inf)*(y<0,-Inf)=-Inf
   897:   //    (-Inf)*(+Inf,0<y)=-Inf
   898:   //    (-Inf)*(y<0,-Inf)=+Inf
   899:   //    (+Inf,0<x)*(+Inf)=+Inf
   900:   //    (x<0,-Inf)*(+Inf)=-Inf
   901:   //    (+Inf,0<x)*(-Inf)=-Inf
   902:   //    (x<0,-Inf)*(-Inf)=+Inf
   903:   //    (+0)*(0<y,+0)=+0
   904:   //    (+0)*(-0,y<0)=-0
   905:   //    (-0)*(0<y,+0)=-0
   906:   //    (-0)*(-0,y<0)=+0
   907:   //    (0<x,+0)*(+0)=+0
   908:   //    (-0,x<0)*(+0)=-0
   909:   //    (0<x,+0)*(-0)=-0
   910:   //    (-0,x<0)*(-0)=+0
   911:   //    (0<x)*(0<y)=+z
   912:   //    (0<x)*(y<0)=-z
   913:   //    (x<0)*(0<y)=-z
   914:   //    (x<0)*(y<0)=+z
   915:   //    +------+------+------+------+------+------+------+
   916:   //    |      | +Inf |  0<y |   +0 |   -0 |  y<0 | -Inf |
   917:   //    +------+------+------+------+------+------+------+
   918:   //    | +Inf | +Inf | +Inf |  NaN |  NaN | -Inf | -Inf |
   919:   //    |  0<x | +Inf |   +z |   +0 |   -0 |   -z | -Inf |
   920:   //    |   +0 |  NaN |   +0 |   +0 |   -0 |   -0 |  NaN |
   921:   //    |   -0 |  NaN |   -0 |   -0 |   +0 |   +0 |  NaN |
   922:   //    |  x<0 | -Inf |   -z |   -0 |   +0 |   +z | +Inf |
   923:   //    | -Inf | -Inf | -Inf |  NaN |  NaN | +Inf | +Inf |
   924:   //    +------+------+------+------+------+------+------+
   925:   private static final int[] QFP_MUL_FLG = {
   926:     //  perl -e "use strict;my($M,$Z,$I,$N,$OE,$DZ)=(0x08000000,0x04000000,0x02000000,0x01000000,0x00002000,0x00000400);for my$a(0..15){my$x=$a<<24;for my$b(0..15){my$y=$b<<24;($b&7)==0 and print'      ';printf('0x%08x,',0b11101010>>($x>>24&7)&1||0b11101010>>($y>>24&7)&1?$N:(($x|$y)&($Z|$I))==($Z|$I)?$N|$OE:($x^$y)&$M|(($x|$y)&$I?$I:($x|$y)&$Z?$Z:0));($b&7)==7 and print chr 10;}}"
   927:     0x00000000,0x01000000,0x02000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
   928:     0x08000000,0x01000000,0x0a000000,0x01000000,0x0c000000,0x01000000,0x01000000,0x01000000,
   929:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   930:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   931:     0x02000000,0x01000000,0x02000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   932:     0x0a000000,0x01000000,0x0a000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   933:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   934:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   935:     0x04000000,0x01000000,0x01002000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
   936:     0x0c000000,0x01000000,0x01002000,0x01000000,0x0c000000,0x01000000,0x01000000,0x01000000,
   937:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   938:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   939:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   940:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   941:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   942:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   943:     0x08000000,0x01000000,0x0a000000,0x01000000,0x0c000000,0x01000000,0x01000000,0x01000000,
   944:     0x00000000,0x01000000,0x02000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
   945:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   946:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   947:     0x0a000000,0x01000000,0x0a000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   948:     0x02000000,0x01000000,0x02000000,0x01000000,0x01002000,0x01000000,0x01000000,0x01000000,
   949:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   950:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   951:     0x0c000000,0x01000000,0x01002000,0x01000000,0x0c000000,0x01000000,0x01000000,0x01000000,
   952:     0x04000000,0x01000000,0x01002000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
   953:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   954:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   955:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   956:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   957:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   958:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
   959:   };
   960:   public QFP mul (QFP y) {
   961:     //フラグを設定する
   962:     int zf = QFP_MUL_FLG[this.flg >>> 24 - 4 | y.flg >>> 24];
   963:     if ((zf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
   964:       QFP.fpsr |= zf & QFP_EXC;
   965:       this.flg = zf & QFP_MZIN;
   966:       return this;
   967:     }
   968:     //乗算を行う
   969:     int ze = this.epp + y.epp;
   970:     double xd = this.dvl;
   971:     double xc = this.cvl;
   972:     double yd = y.dvl;
   973:     double yc = y.cvl;
   974:     {
   975:       double t1 = xd * yd;
   976:       double t2 = 0x1.0000002p27 * xd;
   977:       double t3 = 0x1.0000002p27 * yd;
   978:       double t4 = (xd - t2) + t2;
   979:       double t5 = (yd - t3) + t3;
   980:       t2 = xd - t4;
   981:       t3 = yd - t5;
   982:       t2 = ((xd + xc) * yc + xc * yd) + ((((t4 * t5 - t1) + t3 * t4) + t2 * t5) + t2 * t3);
   983:       xd = t1 + t2;
   984:       xc = (t1 - xd) + t2;
   985:     }
   986:     //指数部の下駄を設定する
   987:     int e = Math.getExponent (xd) - QFP_GETA_BASE & -QFP_GETA_SIZE;
   988:     if (e != 0) {
   989:       ze += e;
   990:       xd = Math.scalb (xd, -e);
   991:       xc = Math.scalb (xc, -e);
   992:     }
   993:     //結果を格納する
   994:     this.flg = zf;
   995:     this.epp = ze;
   996:     this.dvl = xd;
   997:     this.cvl = xc;
   998:     return this;
   999:   }  //qfp.mul(QFP)
  1000:   public QFP mul (QFP x, QFP y) {
  1001:     //フラグを設定する
  1002:     int zf = QFP_MUL_FLG[x.flg >>> 24 - 4 | y.flg >>> 24];
  1003:     if ((zf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1004:       QFP.fpsr |= zf & QFP_EXC;
  1005:       this.flg = zf & QFP_MZIN;
  1006:       return this;
  1007:     }
  1008:     //乗算を行う
  1009:     int ze = x.epp + y.epp;
  1010:     double xd = x.dvl;
  1011:     double xc = x.cvl;
  1012:     double yd = y.dvl;
  1013:     double yc = y.cvl;
  1014:     {
  1015:       double t1 = xd * yd;
  1016:       double t2 = 0x1.0000002p27 * xd;
  1017:       double t3 = 0x1.0000002p27 * yd;
  1018:       double t4 = (xd - t2) + t2;
  1019:       double t5 = (yd - t3) + t3;
  1020:       t2 = xd - t4;
  1021:       t3 = yd - t5;
  1022:       t2 = ((xd + xc) * yc + xc * yd) + ((((t4 * t5 - t1) + t3 * t4) + t2 * t5) + t2 * t3);
  1023:       xd = t1 + t2;
  1024:       xc = (t1 - xd) + t2;
  1025:     }
  1026:     //指数部の下駄を設定する
  1027:     int e = Math.getExponent (xd) - QFP_GETA_BASE & -QFP_GETA_SIZE;
  1028:     if (e != 0) {
  1029:       ze += e;
  1030:       xd = Math.scalb (xd, -e);
  1031:       xc = Math.scalb (xc, -e);
  1032:     }
  1033:     //結果を格納する
  1034:     this.flg = zf;
  1035:     this.epp = ze;
  1036:     this.dvl = xd;
  1037:     this.cvl = xc;
  1038:     return this;
  1039:   }  //qfp.mul(QFP,QFP)
  1040: 
  1041:   //qfp = qfp.mul2 ()
  1042:   //qfp = qfp.mul2 (x)
  1043:   //  2倍
  1044:   public QFP mul2 () {
  1045:     if ((this.flg & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1046:       return this;
  1047:     }
  1048:     //±0,±Inf,NaN以外
  1049:     //2倍する
  1050:     this.dvl *= 2.0;
  1051:     this.cvl *= 2.0;
  1052:     //指数部の下駄を設定する
  1053:     int e = Math.getExponent (this.dvl) - QFP_GETA_BASE & -QFP_GETA_SIZE;
  1054:     if (e != 0) {
  1055:       this.epp += e;
  1056:       this.dvl = Math.scalb (this.dvl, -e);
  1057:       this.cvl = Math.scalb (this.cvl, -e);
  1058:     }
  1059:     return this;
  1060:   }  //qfp.mul2()
  1061:   public QFP mul2 (QFP x) {
  1062:     this.flg = x.flg;
  1063:     if ((this.flg & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1064:       return this;
  1065:     }
  1066:     //±0,±Inf,NaN以外
  1067:     //2倍する
  1068:     this.epp = x.epp;
  1069:     this.dvl = x.dvl * 2.0;
  1070:     this.cvl = x.cvl * 2.0;
  1071:     //指数部の下駄を設定する
  1072:     int e = Math.getExponent (this.dvl) - QFP_GETA_BASE & -QFP_GETA_SIZE;
  1073:     if (e != 0) {
  1074:       this.epp += e;
  1075:       this.dvl = Math.scalb (this.dvl, -e);
  1076:       this.cvl = Math.scalb (this.cvl, -e);
  1077:     }
  1078:     return this;
  1079:   }  //qfp.mul2(QFP)
  1080: 
  1081:   //qfp = qfp.neg ()
  1082:   //qfp = qfp.neg (x)
  1083:   //qfp = qfp.neg (b)
  1084:   //qfp = qfp.neg (b, x)
  1085:   //  符号反転
  1086:   public QFP neg () {
  1087:     int xf = this.flg;
  1088:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1089:       if ((xf & (QFP_Z | QFP_I)) != 0) {  //±0,±Inf
  1090:         this.flg = xf ^ QFP_M;
  1091:       }
  1092:       return this;
  1093:     }
  1094:     //±0,±Inf,NaN以外
  1095:     this.flg = xf ^ QFP_M;
  1096:     this.dvl = -this.dvl;
  1097:     this.cvl = -this.cvl;
  1098:     return this;
  1099:   }  //qfp.neg()
  1100:   public QFP neg (QFP x) {
  1101:     int xf = x.flg;
  1102:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1103:       if ((xf & (QFP_Z | QFP_I)) != 0) {  //±0,±Inf
  1104:         xf ^= QFP_M;
  1105:       }
  1106:       this.flg = xf;
  1107:       return this;
  1108:     }
  1109:     //±0,±Inf,NaN以外
  1110:     this.flg = xf ^ QFP_M;
  1111:     this.dvl = -this.dvl;
  1112:     this.cvl = -this.cvl;
  1113:     return this;
  1114:   }  //qfp.neg(QFP)
  1115:   public QFP neg (boolean b) {
  1116:     return b ? this.neg () : this;
  1117:   }  //qfp.neg(boolean)
  1118:   public QFP neg (boolean b, QFP x) {
  1119:     return b ? this.neg (x) : this.setq (x);
  1120:   }  //qfp.neg(boolean,QFP)
  1121: 
  1122:   //x = x.negdec ()
  1123:   //  x=1-x
  1124:   //  逆デクリメント
  1125:   public QFP negdec () {
  1126:     return this.sub (QFP_ONE, this);
  1127:   }  //qfp.dec()
  1128:   public QFP negdec (QFP x) {
  1129:     return this.sub (QFP_ONE, x);
  1130:   }  //qfp.dec(QFP)
  1131: 
  1132:   //qfp = qfp.rcp ()
  1133:   //qfp = qfp.rcp (x)
  1134:   //  逆数
  1135:   //  z[n+1]=2*z[n]-x*z[n]^2
  1136:   public QFP rcp () {
  1137:     int xf = this.flg;
  1138:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1139:       if ((xf & QFP_Z) != 0) {  //±0
  1140:         QFP.fpsr |= QFP_DZ;
  1141:         xf ^= QFP_Z | QFP_I;
  1142:       } else if ((xf & QFP_I) != 0) {  //±Inf
  1143:         xf ^= QFP_Z | QFP_I;
  1144:       }
  1145:       this.flg = xf;
  1146:       return this;
  1147:     }
  1148:     QFP t1 = new QFP (xf, -this.epp, 1.0 / this.dvl, 0.0);
  1149:     QFP t2 = new QFP ().squ (t1).mul (this);
  1150:     return this.mul2 (t1).sub (t2);
  1151:   }  //qfp.rcp()
  1152:   public QFP rcp (QFP x) {
  1153:     int xf = x.flg;
  1154:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1155:       if ((xf & QFP_Z) != 0) {  //±0
  1156:         QFP.fpsr |= QFP_DZ;
  1157:         xf ^= QFP_Z | QFP_I;
  1158:       } else if ((xf & QFP_I) != 0) {  //±Inf
  1159:         xf ^= QFP_Z | QFP_I;
  1160:       }
  1161:       this.flg = xf;
  1162:       return this;
  1163:     }
  1164:     QFP t1 = new QFP (xf, -x.epp, 1.0 / x.dvl, 0.0);
  1165:     QFP t2 = new QFP ().squ (t1).mul (x);
  1166:     return this.mul2 (t1).sub (t2);
  1167:   }  //qfp.rcp(QFP)
  1168: 
  1169:   //qfp = qfp.set0 ()
  1170:   //  0代入
  1171:   public final QFP set0 () {
  1172:     this.flg = QFP_P | QFP_Z;
  1173:     return this;
  1174:   }  //qfp.set0()
  1175: 
  1176:   //qfp = qfp.setd (d)
  1177:   //  double代入
  1178:   public final QFP setd (double d) {
  1179:     int xf = (int) (Double.doubleToLongBits (d) >>> 36) & QFP_M;
  1180:     if (d == 0.0) {  //±0
  1181:       this.flg = xf | QFP_Z;
  1182:       return this;
  1183:     }
  1184:     if (Double.isInfinite (d)) {  //±Inf
  1185:       this.flg = xf | QFP_I;
  1186:       return this;
  1187:     }
  1188:     if (Double.isNaN (d)) {  //NaN
  1189:       this.flg = QFP_N;
  1190:       return this;
  1191:     }
  1192:     //±0,±Inf,NaN以外
  1193:     this.flg = xf;
  1194:     this.epp = Math.getExponent (d) - QFP_GETA_BASE & -QFP_GETA_SIZE;
  1195:     this.dvl = Math.scalb (d, -this.epp);
  1196:     this.cvl = 0.0;
  1197:     return this;
  1198:   }  //qfp.setd(double)
  1199: 
  1200:   //qfp = qfp.setf (f)
  1201:   //  float代入
  1202:   public final QFP setf (float f) {
  1203:     int xf = Float.floatToIntBits (f) >>> 4 & QFP_M;
  1204:     if (f == 0.0) {  //±0
  1205:       this.flg = xf | QFP_Z;
  1206:       return this;
  1207:     }
  1208:     if (Float.isInfinite (f)) {  //±Inf
  1209:       this.flg = xf | QFP_I;
  1210:       return this;
  1211:     }
  1212:     if (Float.isNaN (f)) {  //NaN
  1213:       this.flg = QFP_N;
  1214:       return this;
  1215:     }
  1216:     //±0,±Inf,NaN以外
  1217:     this.flg = xf;
  1218:     this.epp = 0;  //512<=QFP_GETA_SIZEのときfloatの指数部はepp==0に入り切る
  1219:     this.dvl = (double) f;
  1220:     this.cvl = 0.0;
  1221:     return this;
  1222:   }  //qfp.setf(float)
  1223: 
  1224:   //qfp = qfp.seti (i)
  1225:   //  int代入
  1226:   public final QFP seti (int i) {
  1227:     if (i == 0) {
  1228:       this.flg = QFP_P | QFP_Z;
  1229:       return this;
  1230:     }
  1231:     this.flg = i >>> 4 & QFP_M;
  1232:     this.epp = 0;
  1233:     this.dvl = (double) i;
  1234:     this.cvl = 0.0;
  1235:     return this;
  1236:   }  //qfp.seti(int)
  1237: 
  1238:   //qfp = qfp.setl (l)
  1239:   //  long代入
  1240:   public final QFP setl (long l) {
  1241:     if (l == 0L) {  //+0
  1242:       this.flg = QFP_P | QFP_Z;
  1243:       return this;
  1244:     }
  1245:     this.flg = (int) (l >>> 36) & QFP_M;
  1246:     this.epp = 0;
  1247:     //  Javaのlongからdoubleへのキャストはto nearest evenなので誤差の絶対値はulp/2以下
  1248:     //  http://docs.oracle.com/javase/specs/jls/se8/html/jls-5.html#jls-5.1.2
  1249:     //  doubleに変換しても整数であることに変わりないのでlongにキャストすればほとんど元に戻るが、
  1250:     //  longからdoubleへのキャストで2^63になったときだけdoubleからlongへのキャストが飽和変換で2^63-1になってしまう
  1251:     this.dvl = (double) l;  //2^63になる場合がある。longにキャストすると
  1252:     this.cvl = (double) (l - (this.dvl == 0x1p63 ? 0x8000000000000000L : (long) this.dvl));
  1253:     return this;
  1254:   }  //qfp.setl(long)
  1255: 
  1256:   //qfp = qfp.setq (x)
  1257:   //  QFP代入
  1258:   public final QFP setq (QFP x) {
  1259:     this.flg = x.flg;
  1260:     this.epp = x.epp;
  1261:     this.dvl = x.dvl;
  1262:     this.cvl = x.cvl;
  1263:     return this;
  1264:   }  //qfp.setq(QFP)
  1265: 
  1266:   //qfp = qfp.shl (n)
  1267:   //  左シフト
  1268:   public QFP shl (int n) {
  1269:     if ((this.flg & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1270:       return this;
  1271:     }
  1272:     //±0,±Inf,NaN以外
  1273:     //2^n倍する
  1274:     int e = Math.getExponent (this.dvl);  //元の指数部の下位
  1275:     n += this.epp + e;  //新しい指数部の全体
  1276:     this.epp = n - QFP_GETA_BASE & -QFP_GETA_SIZE;  //新しい指数部の上位
  1277:     n -= this.epp + e;  //新しい指数部の下位-元の指数部の下位
  1278:     this.dvl = Math.scalb (this.dvl, n);
  1279:     this.cvl = Math.scalb (this.cvl, n);
  1280:     return this;
  1281:   }  //qfp.shl(int)
  1282: 
  1283:   //qfp = qfp.shr (n)
  1284:   //  右シフト
  1285:   public QFP shr (int n) {
  1286:     if ((this.flg & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1287:       return this;
  1288:     }
  1289:     //±0,±Inf,NaN以外
  1290:     //2^-n倍する
  1291:     n = -n;
  1292:     int e = Math.getExponent (this.dvl);  //元の指数部の下位
  1293:     n += this.epp + e;  //新しい指数部の全体
  1294:     this.epp = n - QFP_GETA_BASE & -QFP_GETA_SIZE;  //新しい指数部の上位
  1295:     n -= this.epp + e;  //新しい指数部の下位-元の指数部の下位
  1296:     this.dvl = Math.scalb (this.dvl, n);
  1297:     this.cvl = Math.scalb (this.cvl, n);
  1298:     return this;
  1299:   }  //qfp.shr(int)
  1300: 
  1301:   //qfp = qfp.sqrt ()
  1302:   //qfp = qfp.sqrt (x)
  1303:   //  平方根
  1304:   //  z[n+1]=(z[n]+x/z[n])/2
  1305:   public QFP sqrt () {
  1306:     int xf = this.flg;
  1307:     if ((xf & QFP_MZIN) != 0) {  //-x,±0,±Inf,NaN
  1308:       if (xf == QFP_M || xf == (QFP_M | QFP_I)) {  //sqrt(-x)=NaN, sqrt(-Inf)=NaN
  1309:         QFP.fpsr |= QFP_OE;
  1310:         this.flg = QFP_N;
  1311:       } else {
  1312:         this.flg = xf;  //sqrt(±0)=±0, sqrt(+Inf)=+Inf, sqrt(NaN)=NaN
  1313:       }
  1314:       return this;
  1315:     }
  1316:     //-x,±0,±Inf,NaN以外
  1317:     QFP t = new QFP (Math.sqrt (this.dvl)).shl (this.epp >> 1);
  1318:     return this.div (this, t).add (t).div2 ();
  1319:   }  //qfp.sqrt()
  1320:   public QFP sqrt (QFP x) {
  1321:     int xf = x.flg;
  1322:     if ((xf & QFP_MZIN) != 0) {  //-x,±0,±Inf,NaN
  1323:       if (xf == QFP_M || xf == (QFP_M | QFP_I)) {  //sqrt(-x)=NaN, sqrt(-Inf)=NaN
  1324:         QFP.fpsr |= QFP_OE;
  1325:         this.flg = QFP_N;
  1326:       } else {
  1327:         this.flg = xf;  //sqrt(±0)=±0, sqrt(+Inf)=+Inf, sqrt(NaN)=NaN
  1328:       }
  1329:       return this;
  1330:     }
  1331:     //-x,±0,±Inf,NaN以外
  1332:     QFP t = new QFP (Math.sqrt (x.dvl)).shl (x.epp >> 1);
  1333:     return this.div (x, t).add (t).div2 ();
  1334:   }  //qfp.sqrt(QFP)
  1335: 
  1336:   //qfp = qfp.squ ()
  1337:   //qfp = qfp.squ (x)
  1338:   //  2乗
  1339:   public QFP squ () {
  1340:     return this.mul (this, this);
  1341:   }  //qfp.squ()
  1342:   public QFP squ (QFP x) {
  1343:     return this.mul (x, x);
  1344:   }  //qfp.squ(QFP)
  1345: 
  1346:   //qfp = qfp.sub (y)
  1347:   //qfp = qfp.sub (x, y)
  1348:   //  減算
  1349:   //  フラグ
  1350:   //    (NaN)-(NaN,+Inf,0<y,+0,-0,y<0,-Inf)=NaN
  1351:   //    (NaN,+Inf,0<x,+0,-0,x<0,-Inf)-(NaN)=NaN
  1352:   //    (+Inf)-(+Inf)=NaN,OE
  1353:   //    (-Inf)-(-Inf)=NaN,OE
  1354:   //    (+Inf)-(0<y,+0,-0,y<0,-Inf)=+Inf
  1355:   //    (-Inf)-(+Inf,0<y,+0,-0,y<0)=-Inf
  1356:   //    (0<x,+0,-0,x<0,-Inf)-(+Inf)=-Inf
  1357:   //    (+Inf,0<x,+0,-0,x<0)-(-Inf)=+Inf
  1358:   //    (+0)-(+0)=+0
  1359:   //    (+0)-(-0)=+0
  1360:   //    (-0)-(+0)=-0
  1361:   //    (-0)-(-0)=+0
  1362:   //    (+0,-0)-(0<y,y<0)=-y
  1363:   //    (0<x,x<0)-(+0,-0)=x
  1364:   //    (0<x,x<0)-(0<y,y<0)=z
  1365:   //    +------+------+------+------+------+------+------+
  1366:   //    |      | +Inf |  0<y |   +0 |   -0 |  y<0 | -Inf |
  1367:   //    +------+------+------+------+------+------+------+
  1368:   //    | +Inf |  NaN | +Inf | +Inf | +Inf | +Inf | +Inf |
  1369:   //    |  0<x | -Inf |    z |    x |    x |    z | +Inf |
  1370:   //    |   +0 | -Inf |   -y |   +0 |   +0 |   -y | +Inf |
  1371:   //    |   -0 | -Inf |   -y |   -0 |   +0 |   -y | +Inf |
  1372:   //    |  x<0 | -Inf |    z |    x |    x |    z | +Inf |
  1373:   //    | -Inf | -Inf | -Inf | -Inf | -Inf | -Inf |  NaN |
  1374:   //    +------+------+------+------+------+------+------+
  1375:   private static final int[] QFP_SUB_FLG = {
  1376:     //  perl -e "use strict;my($M,$Z,$I,$N,$OE,$DZ)=(0x08000000,0x04000000,0x02000000,0x01000000,0x00002000,0x00000400);for my$a(0..15){my$x=$a<<24;for my$b(0..15){my$y=$b<<24;($b&7)==0 and print'      ';printf('0x%08x,',0b11101010>>($x>>24&7)&1||0b11101010>>($y>>24&7)&1?$N:($x&$y)&$I&&$x==$y?$N|$OE:($x|$y)&$I?$x&$I?$x:$y^$M:($x&$y)&$Z?($x&~$y)&$M|$Z:$x&$Z?0xc0000000:$y&$Z?0x80000000:0);($b&7)==7 and print chr 10;}}"
  1377:     0x00000000,0x01000000,0x0a000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
  1378:     0x00000000,0x01000000,0x02000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
  1379:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1380:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1381:     0x02000000,0x01000000,0x01002000,0x01000000,0x02000000,0x01000000,0x01000000,0x01000000,
  1382:     0x02000000,0x01000000,0x02000000,0x01000000,0x02000000,0x01000000,0x01000000,0x01000000,
  1383:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1384:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1385:     0xc0000000,0x01000000,0x0a000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
  1386:     0xc0000000,0x01000000,0x02000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
  1387:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1388:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1389:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1390:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1391:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1392:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1393:     0x00000000,0x01000000,0x0a000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
  1394:     0x00000000,0x01000000,0x02000000,0x01000000,0x80000000,0x01000000,0x01000000,0x01000000,
  1395:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1396:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1397:     0x0a000000,0x01000000,0x0a000000,0x01000000,0x0a000000,0x01000000,0x01000000,0x01000000,
  1398:     0x0a000000,0x01000000,0x01002000,0x01000000,0x0a000000,0x01000000,0x01000000,0x01000000,
  1399:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1400:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1401:     0xc0000000,0x01000000,0x0a000000,0x01000000,0x0c000000,0x01000000,0x01000000,0x01000000,
  1402:     0xc0000000,0x01000000,0x02000000,0x01000000,0x04000000,0x01000000,0x01000000,0x01000000,
  1403:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1404:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1405:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1406:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1407:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1408:     0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,0x01000000,
  1409:   };
  1410:   public QFP sub (QFP y) {
  1411:     int xf = this.flg;
  1412:     int xe = this.epp;
  1413:     double xd = this.dvl;
  1414:     double xc = this.cvl;
  1415:     int yf = y.flg;
  1416:     int ye = y.epp;
  1417:     double yd = y.dvl;
  1418:     double yc = y.cvl;
  1419:     //フラグを設定する
  1420:     int zf = QFP_SUB_FLG[xf >>> 24 - 4 | yf >>> 24];
  1421:     if ((zf & (0xc0000000 | QFP_ZIN)) != 0) {  //x,y,±0,±Inf,NaN
  1422:       if (0 <= zf) {  //±0,±Inf,NaN
  1423:         QFP.fpsr |= zf & QFP_EXC;
  1424:         this.flg = zf & QFP_MZIN;
  1425:       } else if (zf << 1 < 0) {  //-y
  1426:         this.flg = yf ^ QFP_M;
  1427:         this.epp = ye;
  1428:         this.dvl = yd;
  1429:         this.cvl = yc;
  1430:         //} else {  //x
  1431:         //  this.flg = xf;
  1432:         //  this.epp = xe;
  1433:         //  this.dvl = xd;
  1434:         //  this.cvl = xc;
  1435:       }
  1436:       return this;
  1437:     }
  1438:     //両方±0,±Inf,NaN以外
  1439:     //指数部の下駄を合わせる
  1440:     if (xe != ye) {
  1441:       if (xe + 128 == ye) {
  1442:         yd = Math.scalb (yd, 128);
  1443:         yc = Math.scalb (yc, 128);
  1444:       } else if (xe - 128 == ye) {
  1445:         yd = Math.scalb (yd, -128);
  1446:         yc = Math.scalb (yc, -128);
  1447:       } else {
  1448:         if (xe < ye) {  //yが大きすぎるとき-y
  1449:           this.flg = yf ^ QFP_M;
  1450:           this.epp = ye;
  1451:           this.dvl = yd;
  1452:           this.cvl = yc;
  1453:           //} else {  //xが大きすぎるときx
  1454:           //  this.flg = xf;
  1455:           //  this.epp = xe;
  1456:           //  this.dvl = xd;
  1457:           //  this.cvl = xc;
  1458:         }
  1459:         return this;
  1460:       }
  1461:     }
  1462:     //加減算を行う
  1463:     {
  1464:       double t1 = xd - yd;
  1465:       double t2 = xd - t1;
  1466:       t2 = (((xd - (t1 + t2)) + (t2 - yd)) + xc) - yc;
  1467:       xd = t1 + t2;
  1468:       xc = (t1 - xd) + t2;
  1469:     }
  1470:     //フラグを設定する
  1471:     //  両方±0,±Inf,NaN以外の加減算で結果がNaNになることはない
  1472:     //  両方絶対値が2^192未満の加減算で結果が±Infになることはない
  1473:     if (xd == 0.0) {  //±0
  1474:       this.flg = QFP_P | QFP_Z;  //+0
  1475:       return this;
  1476:     }
  1477:     this.flg = (int) (Double.doubleToLongBits (xd) >>> 36) & QFP_M;
  1478:     //指数部の下駄を設定する
  1479:     int e = Math.getExponent (xd) - QFP_GETA_BASE & -QFP_GETA_SIZE;
  1480:     if (e != 0) {
  1481:       xe += e;
  1482:       xd = Math.scalb (xd, -e);
  1483:       xc = Math.scalb (xc, -e);
  1484:     }
  1485:     //結果を格納する
  1486:     this.epp = xe;
  1487:     this.dvl = xd;
  1488:     this.cvl = xc;
  1489:     return this;
  1490:   }  //qfp.sub(QFP)
  1491:   public QFP sub (QFP x, QFP y) {
  1492:     int xf = x.flg;
  1493:     int xe = x.epp;
  1494:     double xd = x.dvl;
  1495:     double xc = x.cvl;
  1496:     int yf = y.flg;
  1497:     int ye = y.epp;
  1498:     double yd = y.dvl;
  1499:     double yc = y.cvl;
  1500:     //フラグを設定する
  1501:     int zf = QFP_SUB_FLG[xf >>> 24 - 4 | yf >>> 24];
  1502:     if ((zf & (0xc0000000 | QFP_ZIN)) != 0) {  //x,y,±0,±Inf,NaN
  1503:       if (0 <= zf) {  //±0,±Inf,NaN
  1504:         QFP.fpsr |= zf & QFP_EXC;
  1505:         this.flg = zf & QFP_MZIN;
  1506:       } else if (zf << 1 < 0) {  //-y
  1507:         this.flg = yf ^ QFP_M;
  1508:         this.epp = ye;
  1509:         this.dvl = yd;
  1510:         this.cvl = yc;
  1511:       } else {  //x
  1512:         this.flg = xf;
  1513:         this.epp = xe;
  1514:         this.dvl = xd;
  1515:         this.cvl = xc;
  1516:       }
  1517:       return this;
  1518:     }
  1519:     //両方±0,±Inf,NaN以外
  1520:     //指数部の下駄を合わせる
  1521:     if (xe != ye) {
  1522:       if (xe + 128 == ye) {
  1523:         yd = Math.scalb (yd, 128);
  1524:         yc = Math.scalb (yc, 128);
  1525:       } else if (xe - 128 == ye) {
  1526:         yd = Math.scalb (yd, -128);
  1527:         yc = Math.scalb (yc, -128);
  1528:       } else {
  1529:         if (xe < ye) {  //yが大きすぎるとき-y
  1530:           this.flg = yf ^ QFP_M;
  1531:           this.epp = ye;
  1532:           this.dvl = yd;
  1533:           this.cvl = yc;
  1534:         } else {  //xが大きすぎるときx
  1535:           this.flg = xf;
  1536:           this.epp = xe;
  1537:           this.dvl = xd;
  1538:           this.cvl = xc;
  1539:         }
  1540:         return this;
  1541:       }
  1542:     }
  1543:     //加減算を行う
  1544:     {
  1545:       double t1 = xd - yd;
  1546:       double t2 = xd - t1;
  1547:       t2 = (((xd - (t1 + t2)) + (t2 - yd)) + xc) - yc;
  1548:       xd = t1 + t2;
  1549:       xc = (t1 - xd) + t2;
  1550:     }
  1551:     //フラグを設定する
  1552:     //  両方±0,±Inf,NaN以外の加減算で結果がNaNになることはない
  1553:     //  両方絶対値が2^192未満の加減算で結果が±Infになることはない
  1554:     if (xd == 0.0) {  //±0
  1555:       this.flg = QFP_P | QFP_Z;  //+0
  1556:       return this;
  1557:     }
  1558:     this.flg = (int) (Double.doubleToLongBits (xd) >>> 36) & QFP_M;
  1559:     //指数部の下駄を設定する
  1560:     int e = Math.getExponent (xd) - QFP_GETA_BASE & -QFP_GETA_SIZE;
  1561:     if (e != 0) {
  1562:       xe += e;
  1563:       xd = Math.scalb (xd, -e);
  1564:       xc = Math.scalb (xc, -e);
  1565:     }
  1566:     //結果を格納する
  1567:     this.epp = xe;
  1568:     this.dvl = xd;
  1569:     this.cvl = xc;
  1570:     return this;
  1571:   }  //qfp.sub(QFP,QFP)
  1572: 
  1573:   //s = qfp.toString ()
  1574:   //  10進数文字列化
  1575:   public static final QFP[] QFP_TEN_P16QR = {
  1576:     //  echo read("efp.gp");for(i=0,50,q=floor(i/16);r=i-16*q;s=Str("10^(16^",q,"*",r,")");qfpmem(eval(s),s)) | gp-2.7 -q
  1577:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(16^0*0)=1
  1578:     new QFP (QFP_P, 0, 0x1.4p3, 0.0),  //10^(16^0*1)=10
  1579:     new QFP (QFP_P, 0, 0x1.9p6, 0.0),  //10^(16^0*2)=100
  1580:     new QFP (QFP_P, 0, 0x1.f4p9, 0.0),  //10^(16^0*3)=1000
  1581:     new QFP (QFP_P, 0, 0x1.388p13, 0.0),  //10^(16^0*4)=10000
  1582:     new QFP (QFP_P, 0, 0x1.86ap16, 0.0),  //10^(16^0*5)=100000
  1583:     new QFP (QFP_P, 0, 0x1.e848p19, 0.0),  //10^(16^0*6)=1000000
  1584:     new QFP (QFP_P, 0, 0x1.312dp23, 0.0),  //10^(16^0*7)=10000000
  1585:     new QFP (QFP_P, 0, 0x1.7d784p26, 0.0),  //10^(16^0*8)=100000000
  1586:     new QFP (QFP_P, 0, 0x1.dcd65p29, 0.0),  //10^(16^0*9)=1000000000
  1587:     new QFP (QFP_P, 0, 0x1.2a05f2p33, 0.0),  //10^(16^0*10)=10000000000
  1588:     new QFP (QFP_P, 0, 0x1.74876e8p36, 0.0),  //10^(16^0*11)=100000000000
  1589:     new QFP (QFP_P, 0, 0x1.d1a94a2p39, 0.0),  //10^(16^0*12)=1000000000000
  1590:     new QFP (QFP_P, 0, 0x1.2309ce54p43, 0.0),  //10^(16^0*13)=10000000000000
  1591:     new QFP (QFP_P, 0, 0x1.6bcc41e9p46, 0.0),  //10^(16^0*14)=100000000000000
  1592:     new QFP (QFP_P, 0, 0x1.c6bf52634p49, 0.0),  //10^(16^0*15)=1000000000000000
  1593:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(16^1*0)=1
  1594:     new QFP (QFP_P, 0, 0x1.1c37937e08p53, 0.0),  //10^(16^1*1)=10000000000000000
  1595:     new QFP (QFP_P, 0, 0x1.3b8b5b5056e17p106, -0x1.3107fp52),  //10^(16^1*2)=1.00000000000000000000000000000 e32
  1596:     new QFP (QFP_P, 0, 0x1.5e531a0a1c873p159, -0x1.14b4c7a76a406p105),  //10^(16^1*3)=1.00000000000000000000000000000 e48
  1597:     new QFP (QFP_P, 0, 0x1.84f03e93ff9f5p212, -0x1.2ac340948e389p157),  //10^(16^1*4)=1.00000000000000000000000000000 e64
  1598:     new QFP (QFP_P, 512, 0x1.afcef51f0fb5fp-247, -0x1.08f322e84da1p-308),  //10^(16^1*5)=1.00000000000000000000000000000 e80
  1599:     new QFP (QFP_P, 512, 0x1.df67562d8b363p-194, -0x1.ae9d180b58861p-248),  //10^(16^1*6)=1.00000000000000000000000000000 e96
  1600:     new QFP (QFP_P, 512, 0x1.0a1f5b8132466p-140, 0x1.4f01f167b5e3p-194),  //10^(16^1*7)=1.00000000000000000000000000000 e112
  1601:     new QFP (QFP_P, 512, 0x1.27748f9301d32p-87, -0x1.901cc86649e4ap-141),  //10^(16^1*8)=1.00000000000000000000000000000 e128
  1602:     new QFP (QFP_P, 512, 0x1.4805738b51a75p-34, -0x1.18a0e9df9363ap-89),  //10^(16^1*9)=1.00000000000000000000000000000 e144
  1603:     new QFP (QFP_P, 512, 0x1.6c2d4256ffcc3p19, -0x1.56a2119e533adp-38),  //10^(16^1*10)=1.00000000000000000000000000000 e160
  1604:     new QFP (QFP_P, 512, 0x1.945145230b378p72, -0x1.b20a11c22bf0cp15),  //10^(16^1*11)=1.00000000000000000000000000000 e176
  1605:     new QFP (QFP_P, 512, 0x1.c0e1ef1a724ebp125, -0x1.4abd220ed605cp71),  //10^(16^1*12)=1.00000000000000000000000000000 e192
  1606:     new QFP (QFP_P, 512, 0x1.f25c186a6f04cp178, 0x1.45a7709a56ccep123),  //10^(16^1*13)=1.00000000000000000000000000000 e208
  1607:     new QFP (QFP_P, 512, 0x1.14a52dffc6799p232, 0x1.2f82bd6b70d9ap177),  //10^(16^1*14)=1.00000000000000000000000000000 e224
  1608:     new QFP (QFP_P, 1024, 0x1.33234de7ad7e3p-227, -0x1.34a66b24bc3ebp-283),  //10^(16^1*15)=1.00000000000000000000000000000 e240
  1609:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(16^2*0)=1
  1610:     new QFP (QFP_P, 1024, 0x1.54fdd7f73bf3cp-174, -0x1.7222446fe467p-229),  //10^(16^2*1)=1.00000000000000000000000000000 e256
  1611:     new QFP (QFP_P, 1536, 0x1.c633415d4c1d2p164, 0x1.c6cc655c54bc5p109),  //10^(16^2*2)=1.00000000000000000000000000000 e512
  1612:     new QFP (QFP_P, 2560, 0x1.2e7f39519a015p-9, -0x1.739372c035fd3p-64),  //10^(16^2*3)=1.00000000000000000000000000000 e768
  1613:     new QFP (QFP_P, 3584, 0x1.92eceb0d02ea2p-183, -0x1.f44d79616b874p-237),  //10^(16^2*4)=1.00000000000000000000000000000 e1024
  1614:     new QFP (QFP_P, 4096, 0x1.0c59181dd70aep156, -0x1.34f7a4332a3fap101),  //10^(16^2*5)=1.00000000000000000000000000000 e1280
  1615:     new QFP (QFP_P, 5120, 0x1.65706a7673275p-18, -0x1.bbda853c4b3e3p-77),  //10^(16^2*6)=1.00000000000000000000000000000 e1536
  1616:     new QFP (QFP_P, 6144, 0x1.dc1bbb0924957p-192, 0x1.185a80e1e6764p-248),  //10^(16^2*7)=1.00000000000000000000000000000 e1792
  1617:     new QFP (QFP_P, 6656, 0x1.3d1676bb8a7acp147, -0x1.0dac596b98e6bp93),  //10^(16^2*8)=1.00000000000000000000000000000 e2048
  1618:     new QFP (QFP_P, 7680, 0x1.a65c406483e9p-27, 0x1.bcfd432008104p-84),  //10^(16^2*9)=1.00000000000000000000000000000 e2304
  1619:     new QFP (QFP_P, 8704, 0x1.194aa9804143ep-200, 0x1.4cbbfda13244dp-256),  //10^(16^2*10)=1.00000000000000000000000000000 e2560
  1620:     new QFP (QFP_P, 9216, 0x1.76ae153537b2fp138, -0x1.9d9bfbc56226fp80),  //10^(16^2*11)=1.00000000000000000000000000000 e2816
  1621:     new QFP (QFP_P, 10240, 0x1.f312ba4bb116bp-36, 0x1.457ee84626d0ap-90),  //10^(16^2*12)=1.00000000000000000000000000000 e3072
  1622:     new QFP (QFP_P, 11264, 0x1.4c61defaad34p-209, -0x1.baaf1c99fdcacp-264),  //10^(16^2*13)=1.00000000000000000000000000000 e3328
  1623:     new QFP (QFP_P, 11776, 0x1.babb91457e4fep129, 0x1.fbcad508c760fp75),  //10^(16^2*14)=1.00000000000000000000000000000 e3584
  1624:     new QFP (QFP_P, 12800, 0x1.26dc0ee6fb8cap-44, -0x1.24e7716893ec1p-100),  //10^(16^2*15)=1.00000000000000000000000000000 e3840
  1625:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(16^3*0)=1
  1626:     new QFP (QFP_P, 13824, 0x1.88c0a40514413p-218, -0x1.94dacfab01fedp-275),  //10^(16^3*1)=1.00000000000000000000000000000 e4096
  1627:     new QFP (QFP_P, 27136, 0x1.2d4743a2ff5e4p77, 0x1.1a0c7c2892306p22),  //10^(16^3*2)=1.00000000000000000000000000000 e8192
  1628:   };
  1629:   public static final QFP[] QFP_TEN_M16QR = {
  1630:     //  echo read("efp.gp");for(i=0,50,q=floor(i/16);r=i-16*q;s=Str("10^(-16^",q,"*",r,")");qfpmem(eval(s),s)) | gp-2.7 -q
  1631:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(-16^0*0)=1
  1632:     new QFP (QFP_P, 0, 0x1.999999999999ap-4, -0x1.999999999999ap-58),  //10^(-16^0*1)=0.100000000000000000000000000000
  1633:     new QFP (QFP_P, 0, 0x1.47ae147ae147bp-7, -0x1.eb851eb851eb8p-63),  //10^(-16^0*2)=0.0100000000000000000000000000000
  1634:     new QFP (QFP_P, 0, 0x1.0624dd2f1a9fcp-10, -0x1.89374bc6a7efap-66),  //10^(-16^0*3)=0.00100000000000000000000000000000
  1635:     new QFP (QFP_P, 0, 0x1.a36e2eb1c432dp-14, -0x1.6a161e4f765fep-68),  //10^(-16^0*4)=0.000100000000000000000000000000000
  1636:     new QFP (QFP_P, 0, 0x1.4f8b588e368f1p-17, -0x1.ee78183f91e64p-71),  //10^(-16^0*5)=1.00000000000000000000000000000 e-5
  1637:     new QFP (QFP_P, 0, 0x1.0c6f7a0b5ed8dp-20, 0x1.b5a63f9a49c2cp-75),  //10^(-16^0*6)=1.00000000000000000000000000000 e-6
  1638:     new QFP (QFP_P, 0, 0x1.ad7f29abcaf48p-24, 0x1.5e1e99483b023p-78),  //10^(-16^0*7)=1.00000000000000000000000000000 e-7
  1639:     new QFP (QFP_P, 0, 0x1.5798ee2308c3ap-27, -0x1.03023df2d4c94p-82),  //10^(-16^0*8)=1.00000000000000000000000000000 e-8
  1640:     new QFP (QFP_P, 0, 0x1.12e0be826d695p-30, -0x1.34674bfabb83bp-84),  //10^(-16^0*9)=1.00000000000000000000000000000 e-9
  1641:     new QFP (QFP_P, 0, 0x1.b7cdfd9d7bdbbp-34, -0x1.20a5465df8d2cp-88),  //10^(-16^0*10)=1.00000000000000000000000000000 e-10
  1642:     new QFP (QFP_P, 0, 0x1.5fd7fe1796495p-37, 0x1.7f7bc7b4d28aap-91),  //10^(-16^0*11)=1.00000000000000000000000000000 e-11
  1643:     new QFP (QFP_P, 0, 0x1.19799812dea11p-40, 0x1.97f27f0f6e886p-96),  //10^(-16^0*12)=1.00000000000000000000000000000 e-12
  1644:     new QFP (QFP_P, 0, 0x1.c25c268497682p-44, -0x1.ecd79a5a0df95p-99),  //10^(-16^0*13)=1.00000000000000000000000000000 e-13
  1645:     new QFP (QFP_P, 0, 0x1.6849b86a12b9bp-47, 0x1.ea70909833de7p-107),  //10^(-16^0*14)=1.00000000000000000000000000000 e-14
  1646:     new QFP (QFP_P, 0, 0x1.203af9ee75616p-50, -0x1.937831647f5ap-104),  //10^(-16^0*15)=1.00000000000000000000000000000 e-15
  1647:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(-16^1*0)=1
  1648:     new QFP (QFP_P, 0, 0x1.cd2b297d889bcp-54, 0x1.5b4c2ebe68799p-109),  //10^(-16^1*1)=1.00000000000000000000000000000 e-16
  1649:     new QFP (QFP_P, 0, 0x1.9f623d5a8a733p-107, -0x1.a2cc10f3892d4p-161),  //10^(-16^1*2)=1.00000000000000000000000000000 e-32
  1650:     new QFP (QFP_P, 0, 0x1.7624f8a762fd8p-160, 0x1.595560c018581p-215),  //10^(-16^1*3)=1.00000000000000000000000000000 e-48
  1651:     new QFP (QFP_P, 0, 0x1.50ffd44f4a73dp-213, 0x1.a53f2398d747bp-268),  //10^(-16^1*4)=1.00000000000000000000000000000 e-64
  1652:     new QFP (QFP_P, -512, 0x1.2f8ac174d6123p246, 0x1.a5dccd879fc96p191),  //10^(-16^1*5)=1.00000000000000000000000000000 e-80
  1653:     new QFP (QFP_P, -512, 0x1.116805effaeaap193, 0x1.cd88ede5810c7p139),  //10^(-16^1*6)=1.00000000000000000000000000000 e-96
  1654:     new QFP (QFP_P, -512, 0x1.ec866b79e0cbap139, 0x1.bea6a30bdaffap85),  //10^(-16^1*7)=1.00000000000000000000000000000 e-112
  1655:     new QFP (QFP_P, -512, 0x1.bba08cf8c979dp86, -0x1.afa9c1a60497dp32),  //10^(-16^1*8)=1.00000000000000000000000000000 e-128
  1656:     new QFP (QFP_P, -512, 0x1.8f9574dcf8a7p33, 0x1.647f32529774bp-21),  //10^(-16^1*9)=1.00000000000000000000000000000 e-144
  1657:     new QFP (QFP_P, -512, 0x1.67e9c127b6e74p-20, 0x1.26b3da42cecadp-76),  //10^(-16^1*10)=1.00000000000000000000000000000 e-160
  1658:     new QFP (QFP_P, -512, 0x1.442e4fb67196p-73, 0x1.7dc56d0072d28p-131),  //10^(-16^1*11)=1.00000000000000000000000000000 e-176
  1659:     new QFP (QFP_P, -512, 0x1.23ff06eea847ap-126, -0x1.fcc24e7cae5cfp-180),  //10^(-16^1*12)=1.00000000000000000000000000000 e-192
  1660:     new QFP (QFP_P, -512, 0x1.0701bd527b498p-179, -0x1.cfdedc1a30d3p-233),  //10^(-16^1*13)=1.00000000000000000000000000000 e-208
  1661:     new QFP (QFP_P, -512, 0x1.d9ca79d89462ap-233, -0x1.425b0740a9caep-288),  //10^(-16^1*14)=1.00000000000000000000000000000 e-224
  1662:     new QFP (QFP_P, -1024, 0x1.aac0bf9b9e65cp226, 0x1.d6fb1e4a9a909p171),  //10^(-16^1*15)=1.00000000000000000000000000000 e-240
  1663:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(-16^2*0)=1
  1664:     new QFP (QFP_P, -1024, 0x1.8062864ac6f43p173, 0x1.39fa911155ffp118),  //10^(-16^2*1)=1.00000000000000000000000000000 e-256
  1665:     new QFP (QFP_P, -1536, 0x1.2093dc65b647ap-165, 0x1.0e3899699972p-219),  //10^(-16^2*2)=1.00000000000000000000000000000 e-512
  1666:     new QFP (QFP_P, -2560, 0x1.b14cda94a0bbdp8, 0x1.6b5ae1b259461p-47),  //10^(-16^2*3)=1.00000000000000000000000000000 e-768
  1667:     new QFP (QFP_P, -3584, 0x1.454d054bb4af8p182, 0x1.7b0f4c02b0d7ap126),  //10^(-16^2*4)=1.00000000000000000000000000000 e-1024
  1668:     new QFP (QFP_P, -4096, 0x1.e870ba12ebdb7p-157, 0x1.5f0fa5fb329e7p-211),  //10^(-16^2*5)=1.00000000000000000000000000000 e-1280
  1669:     new QFP (QFP_P, -5120, 0x1.6eb2893ea54e2p17, 0x1.b268e1eb7534p-38),  //10^(-16^2*6)=1.00000000000000000000000000000 e-1536
  1670:     new QFP (QFP_P, -6144, 0x1.134c7749892f7p191, -0x1.78c9be2976cf8p137),  //10^(-16^2*7)=1.00000000000000000000000000000 e-1792
  1671:     new QFP (QFP_P, -6656, 0x1.9d5ca69e686c6p-148, -0x1.0ddb6d7695869p-202),  //10^(-16^2*8)=1.00000000000000000000000000000 e-2048
  1672:     new QFP (QFP_P, -7680, 0x1.3655081e5142cp26, 0x1.c7f1c9d4e0198p-28),  //10^(-16^2*9)=1.00000000000000000000000000000 e-2304
  1673:     new QFP (QFP_P, -8704, 0x1.d1f6fb85bd815p199, -0x1.fdd33889c15f4p145),  //10^(-16^2*10)=1.00000000000000000000000000000 e-2560
  1674:     new QFP (QFP_P, -9216, 0x1.5dd2e72224515p-139, 0x1.e5718c3e1a28p-193),  //10^(-16^2*11)=1.00000000000000000000000000000 e-2816
  1675:     new QFP (QFP_P, -10240, 0x1.06a17e7922aebp35, 0x1.43f3cf11b5fc6p-19),  //10^(-16^2*12)=1.00000000000000000000000000000 e-3072
  1676:     new QFP (QFP_P, -11264, 0x1.8a57514d5d62cp208, -0x1.36b0b39a33d88p154),  //10^(-16^2*13)=1.00000000000000000000000000000 e-3328
  1677:     new QFP (QFP_P, -11776, 0x1.280d5f1f07facp-130, 0x1.32a5a6f1a7076p-184),  //10^(-16^2*14)=1.00000000000000000000000000000 e-3584
  1678:     new QFP (QFP_P, -12800, 0x1.bc85ff1a6f95bp43, 0x1.fc518fbd22355p-14),  //10^(-16^2*15)=1.00000000000000000000000000000 e-3840
  1679:     new QFP (QFP_P, 0, 0x1.0p0, 0.0),  //10^(-16^3*0)=1
  1680:     new QFP (QFP_P, -13824, 0x1.4dba0991a59d4p217, -0x1.0e90e3f6e2f1ep159),  //10^(-16^3*1)=1.00000000000000000000000000000 e-4096
  1681:     new QFP (QFP_P, -27136, 0x1.b30d8416d0db5p-78, 0x1.a574753f616c1p-135),  //10^(-16^3*2)=1.00000000000000000000000000000 e-8192
  1682:   };
  1683:   public String toString () {
  1684:     int xf = this.flg;
  1685:     if ((xf & QFP_ZIN) != 0) {  //±0,±Inf,NaN
  1686:       return (xf == (QFP_P | QFP_Z) ? "0" :
  1687:               xf == (QFP_M | QFP_Z) ? "-0" :
  1688:               xf == (QFP_P | QFP_I) ? "Infinity" :
  1689:               xf == (QFP_M | QFP_I) ? "-Infinity" :
  1690:               "NaN");
  1691:     }
  1692:     StringBuilder sb = new StringBuilder ();
  1693:     QFP x = new QFP (this);
  1694:     if ((xf & QFP_M) != 0) {  //-x
  1695:       sb.append ('-');
  1696:       x.neg ();
  1697:     }
  1698:     QFP t = new QFP ();
  1699:     //10進数で表現したときの指数部を求める
  1700:     //  10^e<=x<10^(e+1)となるeを求める
  1701:     int e = (int) Math.floor ((double) (x.epp + Math.getExponent (x.dvl)) * 0.30102999566398119521373889472);  //log10(2)
  1702:     //10^-eを掛けて1<=x<10にする
  1703:     if (0 < e) {
  1704:       x.mul (QFP_TEN_M16QR[e & 15]);
  1705:       if (16 <= e) {
  1706:         x.mul (QFP_TEN_M16QR[16 + (e >> 4 & 15)]);
  1707:         if (256 <= e) {
  1708:           x.mul (QFP_TEN_M16QR[32 + (e >> 8 & 15)]);
  1709:           if (4096 <= e) {
  1710:             x.mul (QFP_TEN_M16QR[48 + (e >> 12)]);
  1711:           }
  1712:         }
  1713:       }
  1714:     } else if (e < 0) {
  1715:       x.mul (QFP_TEN_P16QR[-e & 15]);
  1716:       if (e <= -16) {
  1717:         x.mul (QFP_TEN_P16QR[16 + (-e >> 4 & 15)]);
  1718:         if (e <= -256) {
  1719:           x.mul (QFP_TEN_P16QR[32 + (-e >> 8 & 15)]);
  1720:           if (e <= -4096) {
  1721:             x.mul (QFP_TEN_P16QR[48 + (-e >> 12)]);
  1722:           }
  1723:         }
  1724:       }
  1725:     }
  1726:     //整数部2桁、小数部32桁の10進数に変換する
  1727:     //  1<=x<10なのでw[1]が先頭になるはずだが誤差で前後にずれる可能性がある
  1728:     char[] w = new char[34];
  1729:     {
  1730:       int num = x.geti ();
  1731:       int bcd = XEiJ.FMT_BCD4[num];
  1732:       w[0] = (char) ('0' | bcd >> 4     );
  1733:       w[1] = (char) ('0' | bcd      & 15);
  1734:       for (int i = 2; i < 34; i += 4) {
  1735:         x.sub (t.seti (num)).mul (QFP_TENTO4);
  1736:         num = x.geti ();
  1737:         bcd = XEiJ.FMT_BCD4[num];
  1738:         w[i    ] = (char) ('0' | bcd >> 12     );
  1739:         w[i + 1] = (char) ('0' | bcd >>  8 & 15);
  1740:         w[i + 2] = (char) ('0' | bcd >>  4 & 15);
  1741:         w[i + 3] = (char) ('0' | bcd       & 15);
  1742:       }
  1743:     }
  1744:     //先頭の位置を確認する
  1745:     //  w[h]が先頭(0でない最初の数字)の位置
  1746:     int h = w[0] != '0' ? 0 : w[1] != '0' ? 1 : 2;
  1747:     //30+1桁目を四捨五入する
  1748:     int o = h + 30;  //w[o]は四捨五入する桁の位置。w[]の範囲内
  1749:     if ('5' <= w[o]) {
  1750:       int i = o;
  1751:       while ('9' < ++w[--i]) {
  1752:         w[i] = '0';
  1753:       }
  1754:       if (i < h) {  //先頭から繰り上がった。このとき新しい先頭は1でそれ以外はすべて0
  1755:         h--;  //先頭を左にずらす
  1756:         o--;  //末尾を左にずらす
  1757:       }
  1758:     }
  1759:     //先頭の位置に応じて指数部を更新する
  1760:     //  w[h]が整数部、w[h+1..o-1]が小数部。10^eの小数点はw[h]の右側。整数部の桁数はe+1桁
  1761:     e -= h - 1;
  1762:     //末尾の位置を確認する
  1763:     //  w[o-1]が末尾(0でない最後の数字)の位置
  1764:     while (w[o - 1] == '0') {  //全体は0ではないので必ず止まる。小数点よりも左側で止まる場合があることに注意
  1765:       o--;
  1766:     }
  1767:     //指数形式にするかどうか選択して文字列に変換する
  1768:     if (0 <= e && e < 30) {  //1<=x<10^30。指数形式にしない
  1769:       sb.append (w, h, e + 1);  //整数部。末尾の位置に関係なく1の位まで書く
  1770:       h += e + 1;
  1771:       if (h < o) {  //小数部がある
  1772:         sb.append ('.')  //小数部があるときだけ小数点を書く
  1773:           .append (w, h, o - h);  //小数部
  1774:       }
  1775:     } else if (-4 <= e && e < 0) {  //10^-4<=x<1。指数形式にしない
  1776:       sb.append ('0')  //整数部の0
  1777:         .append ('.');  //小数点
  1778:       while (++e < 0) {
  1779:         sb.append ('0');  //小数部の先頭の0の並び
  1780:       }
  1781:       sb.append (w, h, o - h);  //小数部。全体は0ではないので必ず小数部がある
  1782:     } else {  //x<10^-4または10^30<=x。指数形式にする
  1783:       sb.append (w[h++]);  //整数部
  1784:       if (h < o) {  //小数部がある
  1785:         sb.append ('.')  //小数部があるときだけ小数点を書く
  1786:           .append (w, h, o - h);  //小数部
  1787:       }
  1788:       sb.append ('e')  //指数部の始まり
  1789:         .append (e);  //指数部。正符号は省略する
  1790:     }
  1791:     return sb.toString ();
  1792:   }  //qfp.toString()
  1793: 
  1794: }  //class QFP
  1795: 
  1796: 
  1797: