labunix's blog

labunixのラボUnix

gawk4.1.4+MPFRでフィボナッチ数列をF[0]からF[998]まで出力してみる。

■gawk4.1.4+MPFRでフィボナッチ数列をF[0]からF[998]まで出力してみる。

$ awk -V | grep MPFR
GNU Awk 4.1.4, API: 1.1 (GNU MPFR 3.1.5, GNU MP 6.1.2)

$ ldd /usr/bin/gawk 
	linux-vdso.so.1 (0x00007ffcb7b73000)
	libsigsegv.so.2 => /usr/lib/x86_64-linux-gnu/libsigsegv.so.2 (0x00007f32e7839000)
	libreadline.so.7 => /lib/x86_64-linux-gnu/libreadline.so.7 (0x00007f32e75ec000)
	libmpfr.so.4 => /usr/lib/x86_64-linux-gnu/libmpfr.so.4 (0x00007f32e7385000)
	libgmp.so.10 => /usr/lib/x86_64-linux-gnu/libgmp.so.10 (0x00007f32e7102000)
	libdl.so.2 => /lib/x86_64-linux-gnu/libdl.so.2 (0x00007f32e6efe000)
	libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007f32e6bfa000)
	libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f32e685b000)
	libtinfo.so.5 => /lib/x86_64-linux-gnu/libtinfo.so.5 (0x00007f32e6631000)
	/lib64/ld-linux-x86-64.so.2 (0x00007f32e7ce5000)

■gawkの53bit制限とMPFR拡張

$ echo -e "53\n54" | awk '{print 2^$1-1}'
9007199254740991
180143985094819840からはじめて1000までを試してみる。
 制御文字の「^I」を消して最初に差異が出るのは79番目。

$ echo "0 1000" | \
    awk '{for(i=$1;i<($2-2);i++) \
         {f[i]+=1;f[(i+2)]=f[i]+f[(i+1)]}} \
      END{for(n in f){print "F("n+1") =",f[n]}}' | sort -uV > fib3.txt

$ echo "0 1000" | \
    awk -M '{for(i=$1;i<($2-2);i++) \
            {f[i]+=1;f[(i+2)]=f[i]+f[(i+1)]}} \
         END{for(n in f){print "F("n+1") =",f[n]}}' | sort -uV > fib4.txt 

$ sdiff -w 1000 -s fib[34].txt | sed -e 's/\x9//g' | head -2
F(79) = 14472334024676220   |F(79) = 14472334024676221
F(80) = 23416728348467684   |F(80) = 23416728348467685

■MPFR拡張を使用して、clispの実行結果との最初の差異が出るのは999番目。

$ sdiff -s -w 2000 fib[24].txt | sed -e 's/\x9//g' | head -2
F(999) = 26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626      |F(999) = 26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174625
F(1000) = 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875      |F(1000) = 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228874

■元の定義からF(n)=F(n-1)+F(n-2)を計算してみる。

$ grep -B 2 'F(999)' fib4.txt 
F(997) = 10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377
F(998) = 16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249
F(999) = 26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174625

$ echo '10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377+16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249' | bc
26863810024485359386146727202142923967616609318986952340123175997617\
98170024788168933836965448335656419182785616144335631297667364221035\
03246348504103776803673341511728991697231970827639856157644500784741\
74626

■当然だけど下一桁を足せば6になるので、5になる方が誤り。
 同じ209桁なのに、どこかで捨てられているようだ。

$ echo $(((7+9)%10))
6

$ echo '10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377+16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249' | awk -F'+' '{print (substr($1,length($1))+substr($2,length($2)))%10}'
6

$ echo '10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377+16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249' | awk -F'+' '{print length($1),length($2),length($1+$2)}'
209 209 209

■両方の最上位の桁を外せば確かに正しい。

$ echo '10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377+16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249' | awk -F'+' -M '{print (substr($1,2)+substr($2,2))}'
6863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626

$ echo '10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377+16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249' | awk -F'+' -M '{print ($1+substr($2,2))}'
168638100244853593861467272021429239676166093189869523401231759976179817002478816893383696544833565641918278561614433563129766736422103503246348504103776803673341511728991697231970827639856157644500784741746262^10000の実行結果の3012桁でも問題無いので加算に問題があるように思える。

$ echo "2" | awk -M '{print $1^10000}' > 1.txt
$ echo "2^10000" | bc | tr -d '\n\|\\' | xargs echo  > 2.txt
$ sdiff -s [12].txt;wc -c [12].txt
3012 1.txt
3012 2.txt
6024 合計

■F(1000)-F(999)=F(998)についての減算は正しく求められるので、
 MPFRではなくgawk自体の加算時限定のバグな気がする。
 gawk-4.2.1やgawk-5.0.1が出ているけど、MPFRを組み込むのが大変そうななので後日試してみることにする。

$ echo '43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626' | awk -M '{print $1-$2}'
16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249