多倍長計算

 投稿者：しばっち  投稿日：2020年 1月 5日(日)19時06分46秒
  GMP互換のmpirライブラリー、mpfrライブラリーを使用して多倍長計算を行います。

https://gmplib.org
http://mpir.org
https://www.mpfr.org

mpirには整数型、浮動小数型、有理数型がありますが四則計算しかサポートされていません。
そこでmpirを使って関数を使えるようにしたのがmpfrです。
初等関数等がサポートされています。
mpfrを使って複素数を使えるようにしたのがmpcです。
自前で定義するより高速に計算できます。
mpfr及びmpcはC++インターフェイスがありません(Cインターフェイスのみ)のでラッパーにboostライブラリーを使用しています。
(※C++インターフェイスは関数定義でCインターフェイスは副プログラム定義のようなもの)

なお、これらのDLLは自前でライブラリーをビルドしたものではなく、ネットよりダウンロードして入手したものです。
予めご了承ください。

https://www.dll4free.com


1000桁モードで使用できるように関数として定義する場合と
2進モードや10進モードで使用する副プログラムとして定義する2通りがあります。

1000桁モードで使用する場合は下記のように関数として定義します。

OPTION ARITHMETIC DECIMAL_HIGH
INPUT X
PRINT LEXP(X)
END

EXTERNAL  FUNCTION LEXP(X)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION CHARACTER BYTE
LET KETA=1000
LET X$=STR$(X)
LET B$=REPEAT$(CHR$(32),KETA+100)
CALL EXP1000(KETA,X$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET LEXP=VAL(B$(1:I))

SUB EXP1000(KETA,X$,Y$)
   IF POS(X$,"/")>0 OR POS(X$,"(")>0 OR X$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\exp1000.dll","exp1000"
END SUB
END FUNCTION

1000桁モードではほぼ数式通りの記述ができます。
LET N=LEXP(-X)*X

副プログラムとして定義する場合は

OPTION CHARACTER BYTE
LET X$="1.04719755119659774615421446109316762806572313312503527365831486410260546876206966620934494178070568932738269550442743554903128153651686074390845313604282703915009470090064617370185321487431631831012732147627032522197781537615854941126226105509040063638188285564115344953681810888273779786908674971375790819566886877186272496050697365427641803057178812263086345333711017684960682217379471565064717053647768575678858653065103072870579397753726436837284935815412665424985578396191757496374264606100398304327789112081355221436200713164879840824573023405995364790092351307239209772558412822493948922313504400018937571508785360926192378091926320305787905957382281363374165114338218319512368359742656308630784733998537070967398695467813938660454325825710332017290240378333333279099268331701991057760536543953167481981844896943421417410275111489501175397706272367000104594625096219584440279380687239255638243453275116347625182291038652095462745126253125065259395259351072374226887100064262553706530307214006633333333333333"
CALL LSIN(1000,X$,RESULT$)
PRINT RESULT$
END

EXTERNAL  SUB LSIN(KETA,X$,RESULT$)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),KETA+100)
CALL SIN1000(KETA,X$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB SIN1000(KETA,X$,Y$)
   IF POS(X$,"/")>0 OR POS(X$,"(")>0 OR X$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\sin1000.dll","sin1000"
END SUB
END SUB

副プログラムとして定義する場合は数式通りには記述できません。
LET N=LEXP(-X)*X
PRINT N

これを記述し直すと

LET KETA=4000
!!! CALL LMUL(KETA,X$,STR$(-1),XX$)
LET XX$="-" & X$
CALL LEXP(KETA,XX$,Y$) ! Y=EXP(-X)
CALL LMUL(KETA,Y$,X$,N$) ! N=Y*X
PRINT N$
! CALL DISPLAY(N$)

記述通りに書けない代わりに1000桁である必要もないので
例では4000桁を計算します。

あまり意味はありませんが、桁数を16桁とすると2進モード、10進モードでも
使用できます。

PRINT LSIN(PI/6)
END

EXTERNAL  FUNCTION LCOS(X)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),32)
CALL COS1000(16,STR$(X),B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET LCOS=VAL(B$(1:I))

SUB COS1000(KETA,X$,Y$)
   IF POS(X$,"/")>0 OR POS(X$,"(")>0 OR X$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\cos1000.dll","cos1000"
END SUB
END FUNCTION

複素数モードで使用できる関数も定義してみました。

OPTION ARITHMETIC COMPLEX
PRINT CSIN(COMPLEX(PI/4,0))
END

EXTERNAL  FUNCTION CSIN(X)
OPTION CHARACTER BYTE
OPTION ARITHMETIC COMPLEX
LET X$=REPEAT$(CHR$(0),32)
LET Y$=REPEAT$(CHR$(0),32)
CALL CSIN1000(16,STR$(RE(X)),STR$(IM(X)),X$,Y$)
FOR I=LEN(X$) TO 1 STEP -1
   IF X$(I:I)<="9" AND X$(I:I)>="0" THEN EXIT FOR
NEXT I
FOR I=LEN(Y$) TO 1 STEP -1
   IF Y$(I:I)<="9" AND Y$(I:I)>="0" THEN EXIT FOR
NEXT I
LET CSIN=COMPLEX(RE(VAL(X$(1:I))),IM(VAL(Y$(1:I))))

SUB CSIN1000(KETA,A$,B$,X$,Y$)
   IF POS(A$,"/")>0 OR POS(B$,"/")>0 OR POS(A$,"(")>0 OR POS(B$,"(")>0 OR A$="" OR B$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\complex1000.dll","csin1000"
END SUB
END FUNCTION

副プログラムとして定義すれば、複素数多倍長で計算できます。

OPTION CHARACTER BYTE
LET A=1
LET B=2
CALL CEXP(1000,STR$(A),STR$(B),RE$,IM$)
PRINT "(";RE$;" , ";IM$;")"
END

EXTERNAL  SUB CEXP(KETA,A$,B$,RE$,IM$)
OPTION CHARACTER BYTE
LET X$=REPEAT$(CHR$(32),KETA+100)
LET Y$=REPEAT$(CHR$(32),KETA+100)
CALL CEXP1000(KETA,A$,B$,X$,Y$)
FOR I=LEN(X$) TO 1 STEP -1
   IF X$(I:I)<="9" AND X$(I:I)>="0" THEN EXIT FOR
NEXT I
LET RE$=X$(1:I)
FOR I=LEN(Y$) TO 1 STEP -1
   IF Y$(I:I)<="9" AND Y$(I:I)>="0" THEN EXIT FOR
NEXT I
LET IM$=Y$(1:I)

SUB CEXP1000(KETA,A$,B$,X$,Y$)
   IF POS(A$,"/")>0 OR POS(B$,"/")>0 OR POS(A$,"(")>0 OR POS(B$,"(")>0 OR A$="" OR B$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\cexp1000.dll","cexp1000"
END SUB
END SUB


三角関数(SIN,COS等)、逆三角関数(ASIN,ACOS等)、双曲線関数(COSH,TANH等)、逆双曲線関数(ACOSH,ASINH等)、指数関数(EXP)、対数関数(LOG)等
定数計算として副プログラムで円周率(π)、ネイピア数(e)、log(2)等を定義しています。

また、ガウス・ルジャンドル則(有限区間:-1～1)、ガウス・ラゲール則(半無限区間:0～∞)、ガウス・エルミート則(無限区間:-∞～∞)等を使用し
高精度数値積分をサポートしました。

#214
#215
#216

これは、1000桁モード等で使用するためにその零点及び重み係数を1000桁等の精度で算出し次数も1000次(1000点)等という超高精度で数値積分を行おうというものです。

但し、DATA文として記述するにはサイズが大きすぎるため、データはファイルからの読み出しになります。
(1000桁で零点と重み係数がそれぞれ1000点では2MB近いファイルサイズになります)

ガウス・ルジャンドル則で被積分関数 f(x)=1/(1+x*x) を積分区間A～Bで数値積分します。

OPTION ARITHMETIC DECIMAL_HIGH
LET N=1000 !次数
LET A=0 !下限
LET B=1 !上限
!'INPUT B
LET U=(B+A)/2
LET V=(B-A)/2
OPEN #1:NAME "..\data\legendre1000_"&STR$(N)&".txt"
FOR I=1 TO N
   LINE INPUT #1:X$
   LINE INPUT #1:W$
   LET X=VAL(X$)
   LET WEIGHT=VAL(W$)
   LET S=S+F(U+V*X)*V*WEIGHT
NEXT I
PRINT S*4 !ATN(B)
PRINT PI
END

EXTERNAL  FUNCTION F(X)
OPTION ARITHMETIC DECIMAL_HIGH
LET F=1/(1+X*X)
END FUNCTION

                              実行結果

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095324
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199

注意点として1000桁精度で数値積分を行っても、計算結果も1000桁の精度があるわけではありません。

計算結果が既知である場合は、その値と比較して計算精度を確認してください。
計算結果が未知である場合は、数値積分を精度を変えて2度(800点と1000点等)行い、どこまで一致しているかで
計算精度を確認してください。

下記は副プログラムで定義して計算しています。
被積分関数 f(x)=exp(-x*x)として無限区間積分をガウス・エルミート則で行っています。
例えばファイルに2000桁で係数を求めていれば2000桁までの数値積分ができます。


OPTION CHARACTER BYTE
LET KETA=1000
LET N=1000 !次数
OPEN #1:NAME "..\data\hermite1000_"&STR$(N)&".txt"
FOR I=1 TO N
   LINE INPUT #1:X$
   LINE INPUT #1:WEIGHT$
   CALL LMUL(KETA,X$,X$,W$) !'W=X*X
   CALL LEXP(KETA,"-"&W$,F$) !'F=EXP(-W)
   CALL LMUL(KETA,F$,WEIGHT$,S$) !' S=F*WEIGHT
   CALL LADD(KETA,SS$,S$,TOTAL$) !'TOTAL=SS+S
   LET SS$=TOTAL$
NEXT I
CALL DISPLAY(SS$) ! SQR(PI)
END

EXTERNAL  SUB LEXP(KETA,X$,RESULT$)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),KETA+100)
CALL EXP1000(KETA,X$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB EXP1000(KETA,X$,Y$)
   IF POS(X$,"/")>0 OR POS(X$,"(")>0 OR X$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\exp1000.dll","exp1000"
END SUB
END SUB

EXTERNAL  SUB LADD(KETA,X$,Y$,RESULT$)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),KETA+100)
IF X$="" THEN LET X$="0"
IF Y$="" THEN LET Y$="0"
CALL LADD1000(KETA,X$,Y$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB LADD1000(KETA,X$,Y$,RESULT$)
   IF POS(X$,"/")>0 OR POS(Y$,"/")>0 OR POS(X$,"(")>0 OR POS(Y$,"(")>0 OR X$="" OR Y$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\calc1000.dll","add1000"
END SUB
END SUB

EXTERNAL  SUB LSUB(KETA,X$,Y$,RESULT$)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),KETA+100)
IF X$="" THEN LET X$="0"
IF Y$="" THEN LET Y$="0"
CALL LSUB1000(KETA,X$,Y$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB LSUB1000(KETA,X$,Y$,RESULT$)
   IF POS(X$,"/")>0 OR POS(Y$,"/")>0 OR POS(X$,"(")>0 OR POS(Y$,"(")>0 OR X$="" OR Y$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\calc1000.dll","sub1000"
END SUB
END SUB

EXTERNAL  SUB LMUL(KETA,X$,Y$,RESULT$)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),KETA+100)
IF X$="" THEN LET X$="1"
IF Y$="" THEN LET Y$="1"
CALL LMUL1000(KETA,X$,Y$,B$)
FOR I=LEN(B$) TO 1 STEP -1
   IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
      EXIT FOR
   END IF
NEXT I
LET RESULT$=B$(1:I)

SUB LMUL1000(KETA,X$,Y$,RESULT$)
   IF POS(X$,"/")>0 OR POS(Y$,"/")>0 OR POS(X$,"(")>0 OR POS(Y$,"(")>0 OR X$="" OR Y$="" THEN
      PRINT "ERROR"
      STOP
   END IF
   ASSIGN ".\DLL\calc1000.dll","mul1000"
END SUB
END SUB

EXTERNAL  SUB DISPLAY(X$)
OPTION CHARACTER BYTE
LET N=POS(X$,".")
IF N>0 THEN
   FOR I=1 TO N
      PRINT X$(I:I);
      IF MOD(I,5)=0 THEN PRINT " ";
      IF MOD(I,50)=0 THEN PRINT "   ";
      IF MOD(I,100)=0 THEN PRINT
   NEXT I
   PRINT
END IF
FOR I=0 TO LEN(X$)-N STEP 5
   PRINT X$(I+N+1:I+N+5);" ";
   IF MOD(I+5,100)<>0 AND MOD(I+5,50)=0 THEN PRINT "   ";
   IF MOD(I+5,100)=0 THEN PRINT ":";I+5
   IF MOD(I+5,1000)=0 THEN PRINT
   IF MOD(I+5,10000)=0 THEN PRINT
NEXT I
END SUB

パーサーも作ってみました。sin,cos,tan,sqrt,%,^,abs,exp,log,max,acos,sinhなどの関数が使えます。
EXPRESSION1$にはf(x) (1変数)
EXPRESSION2$にはg(x,y) (2変数)を
EXPRESSION3$にはh(x,y,z) (3変数)の関数定義します。
FUNC$に呼び出す関数を入れます。
パラメーターにu,v,w,x,y,zが使用できます。
EXPRESSION1$に被積分関数 1/(1+x*x)を定義し
EXPRESSION2$で和を求めています。


OPTION CHARACTER BYTE
LET KETA=1000
LET A=0 !下限
LET B=1 !上限
LET U=(B+A)/2
LET V=(B-A)/2
LET EXPRESSION1$="1/(1+x*x)"
LET EXPRESSION2$="y+x"
LET N=800
OPEN #1:NAME "..\data\legendre1000_"+STRT$(N)+".txt"
FOR I=1 TO N
   LINE INPUT #1:X$
   LINE INPUT #1:WEIGHT$
   LET FUNC$="f(u+v*x)*v*w"
   CALL PARSER(KETA,FUNC$,STR$(U),STR$(V),WEIGHT$,X$,"","",EXPRESSION1$,EXPRESSION2$,EXPRESSION3$,OUTPUT$)
   LET FUNC$="g(x,y)"
   CALL PARSER(KETA,FUNC$,"","","",OUTPUT$,Y$,"",EXPRESSION1$,EXPRESSION2$,EXPRESSION3$,S$)
   LET Y$=S$
NEXT  I
CLOSE #1
CALL PARSER(KETA,"f(x)",U$,V$,W$,S$,Y$,Z$,"4*x",EXP2$,EXP3$,S$)
PRINT S$
END

EXTERNAL  SUB PARSER(KETA,INPUT$,U$,V$,W$,X$,Y$,Z$,EXPRESSION1$,EXPRESSION2$,EXPRESSION3$,OUTPUT$)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),KETA+100)
IF EXPRESSION1$="" THEN LET EXPRESSION1$="x" !' F(X)
IF EXPRESSION2$="" THEN LET EXPRESSION2$="x" !' G(X,Y)
IF EXPRESSION3$="" THEN LET EXPRESSION3$="x" !' H(X,Y,Z)
IF INPUT$="" THEN
   PRINT "ERROR"
   STOP
END IF
IF U$="" THEN LET U$="0"
IF V$="" THEN LET V$="0"
IF W$="" THEN LET W$="0"
IF X$="" THEN LET X$="0"
IF Y$="" THEN LET Y$="0"
IF Z$="" THEN LET Z$="0"
CALL PARSER1000(KETA,LCASE$(INPUT$),U$,V$,W$,X$,Y$,Z$,LCASE$(EXPRESSION1$),LCASE$(EXPRESSION2$),LCASE$(EXPRESSION3$),B$)
IF B$(1:5)="error" THEN
   PRINT "ERROR!!"
   STOP
ELSE
   FOR I=LEN(B$) TO 1 STEP -1
      IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
         EXIT FOR
      END IF
   NEXT I
   LET OUTPUT$=B$(1:I)
END IF

SUB PARSER1000(KETA,INPUT$,U$,V$,W$,X$,Y$,Z$,EXP1$,EXP2$,EXP3$,OUTPUT$)
   ASSIGN ".\DLL\parser1000_3.dll","parser1000"
END SUB
END SUB
 		 

Re: 多倍長計算

 投稿者：しばっち  投稿日：2020年 1月 5日(日)19時08分52秒
  > No.4760[元記事へ]

なお、係数算出にはコンソールアプリですが作成しました。
win32版、win64版を用意しました。
---------------------------------------------------------------------
                         legendre.cpp


#include <iostream>
#include <iomanip>
#include <string>
#include <cstdio>
#include <boost/multiprecision/mpfr.hpp>

#pragma comment(lib, "mpir.lib")
#pragma comment(lib, "mpfr.lib")

using namespace std;
using namespace boost::multiprecision;

int flg=0;
mpfr_float eps;

void derivative(int n,mpfr_float *a,mpfr_float *b)
{
    int i;
    for(i=0; i<=n; i++) b[i]=0.0;
    for(i=n; i>=1; i--) b[i-1]=(mpfr_float)i*a[i];
}

void div(int n,mpfr_float *a,mpfr_float p)
{
    mpfr_float *c;
    c=new mpfr_float [n+1];
    int i;
    for(i=0; i<=n; i++) c[i]=0.0;
    for(i=n; i>=1; i--) c[i-1]=a[i]+c[i]*p;
    for(i=0; i<=n; i++) a[i]=c[i];
    delete [] c;
}

void legendrepoly(int kk,mpfr_float *newp)
{
    mpfr_float *p,*oldp;
    p=new mpfr_float [kk+1];
    oldp=new mpfr_float [kk+1];
    int k,j,i;
    for(i=0; i<=kk; i++) {
        newp[i]=0.0;
        p[i]=0.0;
        oldp[i]=0.0;
    }
    oldp[0]=1.0;
    p[1]=1.0;
    for(k=2; k<=kk; k++) {
        for(j=1; j<=k; j++) {
            newp[j]=newp[j]+((mpfr_float)2.0*(mpfr_float)k-1.0)/(mpfr_float)k*p[j-1];
            newp[j-1]=newp[j-1]-((mpfr_float)k-1.0)/(mpfr_float)k*oldp[j-1];
        }
        if(k<kk)
            for(i=0; i<=k; i++) {
                oldp[i]=p[i];
                p[i]=newp[i];
                newp[i]=0.0;
            }
    }
    delete [] p;
    delete [] oldp;
}

mpfr_float horner(int n,mpfr_float *a,mpfr_float xx)
{
    int i;
    mpfr_float y;
    y=a[n];
    for(i=n-1; i>=0; i--) y=y*xx+a[i];
    return y;
}

mpfr_float legendre(int k,mpfr_float x)
{
    mpfr_float oldp,newp,p;
    int n;
    oldp=1.0;
    p=x;
    for(n=1; n<=k-1; n++) {
        newp=((2.0*(mpfr_float)n+1.0)*x*p-(mpfr_float)n*oldp)/((mpfr_float)n+1.0);
        oldp=p;
        p=newp;
    }
    return p;
}

mpfr_float legendrediff(int n,mpfr_float x)
{
    return ((mpfr_float)n*x*legendre(n,x)-(mpfr_float)n*legendre(n-1,x))/(x*x-(mpfr_float)1.0);
}

mpfr_float weight(int n,mpfr_float x)
{
    mpfr_float s;
    s=(mpfr_float)n*legendre(n-1,x);
    return (mpfr_float)2.0*(1.0-x*x)/s/s;
}

void legendrepara(int n,mpfr_float *a,mpfr_float *w)
{
    mpfr_float *p,*d,x,xx;
    p=new mpfr_float [n+1];
    d=new mpfr_float [n+1];
    int i,j;
    for(i=0; i<=n; i++) {
        p[i]=0.0;
        d[i]=0.0;
    }
    legendrepoly(n,p);
    for(i=1; i<=(n+1)/2; i++)
    {
        if (flg==0) cerr << "あと " << (n+1)/2-i << endl;
        if (n % 2==1 && i==(n+1)/2) {
            a[(n+1)/2]=0.0;
            w[(n+1)/2]=weight(n,0.0);
        } else {
            for(j=n; j>=0; j--) if (p[j]!=0.0) break;
            for(int k=0; k<j; k++) p[k]/=p[j];
            p[j]=1.0;
            derivative(n,p,d);
            xx=-1.0;
            for(j=0; j<1000; j++) {
                x=xx;
                xx=x-horner(n,p,x)/horner(n,d,x);
                if (fabs(x-xx)<eps) break;
            }
            a[i]=xx;
            a[n-i+1]=-xx;
            w[i]=weight(n,xx);
            w[n-i+1]=w[i];
            div(n,p,xx);
        }
    }
    delete [] p;
    delete [] d;
}

int main(int argc, char *argv[])
{
    int keta=0,prec=0,maxlevel=0;
    if (argc==4 || argc==5)
    {
        keta=atoi(argv[1]);
        prec=atoi(argv[2]);
        maxlevel=atoi(argv[3]);
    }
    else {
        cerr << "計算精度(桁) ";
        cin >> keta;
        cerr << "出力精度(桁) ";
        cin >> prec;
        cerr << "求める次数 ";
        cin >> maxlevel;
    }
    if (argc==5) flg=1;
    if (keta<=0) keta=100;
    if (prec<=0) prec=16;
    if (maxlevel<=0) maxlevel=20;
    mpfr_float::default_precision(keta);
    mpfr_float *x,*w;
    eps=pow(10.0,-prec);
    x=new mpfr_float [maxlevel+1];
    w=new mpfr_float [maxlevel+1];
    legendrepara(maxlevel,x,w);
    for(int i=1; i<=maxlevel; i++) {
        if (prec<=32)
            cout << std::scientific << setprecision(prec+1) << x[i] << " , " << w[i] << endl;
        else {
            cout << std::scientific << setprecision(prec+1) << x[i] << endl;
            cout << std::scientific << setprecision(prec+1) << w[i] << endl;
        }
    }
    cout << endl;
    delete [] x;
    delete [] w;
    return 0;
}

バッチファイル(.bat)を作成すると簡単です。
標準出力に書き出すのでリダイレクトします。

                makedata.bat


legendre 2000 1008 500 > legendre1000_500.txt
legendre 2000 1008 800 > legendre1000_800.txt
legendre 2000 1008 1000 > legendre1000_1000.txt
legendre 2000 1008 1200 > legendre1000_1200.txt

計算精度や出力精度、次数を指定します。
ニュートン法を使用し、出力精度分収束したらループを打ち切るようにしています。
誤差が累積していくので計算精度は出力精度以上(出力精度の2倍近く)で計算させます。

laguerre 128 16 8 > laguerre16_8.txt
laguerre 128 16 15 > laguerre16_15.txt
laguerre 128 16 20 > laguerre16_20.txt

のように出力精度16桁にすると2進モード、10進モードでも利用できます。
※計算精度を高めにしています。

legendre 4000 2000 1500 * > legendre2000_1500.txt
のように次数指定のあとに何か書くと残数表示しません。

時間がかかりすぎる等途中で止めたい時は、CTRL+C で実行中断します。


次にライブラリーを呼び出しているだけで計算方法は全く不明なのですが...

OPTION CHARACTER BYTE
!!LET KETA=100000000 !1億桁
LET KETA=1000000 !100万桁
LET T=TIME
CALL PI1000(KETA,RESULT$)
PRINT TIME-T
OPEN #1:NAME "pi.txt"
ERASE #1
PRINT #1:RESULT$
CLOSE #1
END

EXTERNAL  SUB PI1000(KETA,RESULT$)
OPTION CHARACTER BYTE
LET B$=REPEAT$(CHR$(32),KETA+100)
CALL PI_CALC(KETA,B$)
IF B$(1:5)="error" THEN
   PRINT "ERROR"
   STOP
ELSE
   FOR I=LEN(B$) TO 1 STEP -1
      IF B$(I:I)<="9" AND B$(I:I)>="0" THEN
         EXIT FOR
      END IF
   NEXT I
   LET RESULT$=B$(1:I)
END IF

SUB PI_CALC(KETA,X$)
   ASSIGN ".\DLL\pi1000_2.dll","pi1000"
END SUB
END SUB
------------------------------------------------------------------
                  pi1000_2.cpp

#include <string>
#include <mpreal.h>
#include <sstream>

#pragma comment(lib, "mpir.lib")
#pragma comment(lib, "mpfr.lib")

using namespace std;
using mpfr::mpreal;

extern "C"  __declspec(dllexport)  void pi1000(int keta,char *b)
{
    mpreal::set_default_prec(mpfr::digits2bits(keta));
    string s;
    ostringstream oss;
    try {
        const mpreal pi=mpfr::const_pi();
        oss.precision(keta);
        oss << pi;
        s=oss.str();
    } catch (...) {
        s="error";
    }
    strcpy(b,s.c_str());
}


桁数に1億桁の指定ができます。(mpreal型)
(※mpfr_float型ではどうやら400万桁程で打ち切りのようです。未確認)

BASICAcc(x64版)で 712.827秒
実行するとおよそ100MB程のファイル(pi.txt)ができました。

なお、BASIC(x86版)だと 1687.12秒
かかりました。

今回は64ビット版(BASICAcc.exe)がリリースされましたので
BASICAcc(x64)用にx64版DLL 及び 32ビット版BASIC用にx86版DLLを公開します。
x64版はVC++2019(x64)でx86版はVC++2012(x86)でコンパイルしています。


BASICAccには1000桁モードはありませんので副プログラム定義のみです。
32ビット版BASICには1000モード用の関数定義と副プログラム定義を入れています。

下記からダウンロードしてください。
●BASICAcc(win64) x64版 (multi precision(x64).zip)

https://5.gigafile.nu/0305-cf927d401fe60d0851521225f84509659

ダウンロード期限：2020年3月5日(木)
ダウンロードパスワード：未設定


●BASIC(win32) x86版 (multi precision(x86).zip)

https://5.gigafile.nu/0305-c438fcda253ac5ecae0510162529420fc

ダウンロード期限：2020年3月5日(木)
ダウンロードパスワード：未設定


解凍してできたdataフォルダやdllフォルダ、mpir.dll,mpfr.dll,libgmp-10.dll,libmpfr-4.dll,libmpc-3.dll等は
BASIC.exe(x86) 又はBASICAcc.exe(x64)と同じフォルダに入れてください。
実行時にmpir.dll,mpfr.dll 又はlibgmp-10.dll,libmpfr-4.dll,libmpc-3.dllが必要です(x86版とx64版があります)

BASICAcc1202(x64)ではコンパイルするとoutputフォルダに実行ファイルが出力されるためエラーが出ます。
そのまま終了させてoutputフォルダから取り出してから実行し直してください。
 		 

戻る