Petitverse

https://qiita.com:443/qnighy/items/ba7e6cb6e321972bebf2#%E6%95%B4%E6%95%B0%E3%81%AB%E3%82%88%E3%82%8B%E8%BF%91%E4%BC%BC%E8%A8%88%E7%AE%97
AKSに使うのはこれですね。
後で読むとして、SQRのROUNDは多倍長には使えないので、これを使うしかなさそう。

近似とあるので、無理数を誤差なしで表現できるとかではないんでしょうか。(無理数を誤差なしで表現できると嬉しい。Sinc()は無理数なので)

いくつかの無理数自体を有限の桁数で表現できるのかな…と思ったんですけど、そうではない感じなんですかね？
いや、3^(1/3)を2進数で表してる画像が流れてきて無理数を表現できないかと思ったんですけど、よく見たらこれも無限に続きそうですね。
https://twitter.com:443/fmathsecond/status/1886069803510903153

p進数、剰余はおろか除算もよく分かってない(除算が定義されているかどうかすら知らない)ので、剰余が分かりやすいってのはピンと来てないですね。乗算も怪レい。

そんなわけで、調べてみます。

https://blog.miraikan.jst.go.jp/articles/20190315post-339.html
例えば、2の平方根は7進数で表現できるっぽくて、それは無限桁が必要？
実数の
1.41421356…
と何が違うんだろう？
循環小数…と言って良いのかな…で表せるってこと？

-1も2進数で無限桁必要だっけ…
…11111
奇しくも2の補数と同じ形。いや、導出方法を考えれば当然か。

https://ja.wikipedia.org/wiki/P%E9%80%B2%E6%95%B0#%E7%95%A5%E5%BC%8F%E3%81%AE%E8%A7%A3%E8%AA%AC
p進数で四則演算も存在するっぽくて、任意のpについて負数の概念も存在するらしい。
2進数の実家感よ…
有限桁の小数範囲で必ず逆数が存在する…らしい(整数範囲が有限桁とは言ってない)
というより、逆数が必ず存在するためにp(素数)なのか…(これはmod pと似てる)

逆数を実際に計算してみる。
&b11(3)の逆数は
…010101011
らしい
3倍するには元の数に1桁右シフトした数を足せば良い
…010101011+
…101010110=
…000000001
本当だ。

プログラムで計算してみる。
プログラムの場合、上位のビットを切り捨てれば良い。
OPTION DEFINT
VAR A#=&b010101010101011
VAR B#=3
PRINT FORMAT$("%d*%d=%d",A#,B#,A#*B#AND 32767)

やってることはmod 32768上でのモジュラ逆数と同じ。
https://keisan.casio.jp/exec/system/1589337561

mod 32768だと偶数の逆数は存在しないんだけど、その場合は小数点を使えば良いのか。
例えば、&b10(2)の逆数は普通に0.1
つまり、&b110(6)の逆数は2の逆数と3の逆数を掛けたものなので
…101010101010101.1


OPTION DEFINT
VAR A#=21845.5
VAR B#=6
PRINT FORMAT$("%f*%d=%d",A#,B#,ROUND(A#*B#)AND 32767)

なかなか面白い。
無理数を誤差なしで表現できる気がしてきた。
正直、(整数)/(2^n)を誤差なしで表現できて離散フーリエ変換が定義できる時点で十分期待できる。(ガウシアンフィルタは(整数)/(2^n)の形なので)

モジュラ逆数は拡張ユークリッドの互除法というアルゴリズムで求まるので、頑張れば2進数上の逆数を求めるプログラムを書けそう。

3^(1/3)を計算してみる。(もちろん実数では無理数)
https://twitter.com:443/fmathsecond/status/1886069803510903153

OPTION DEFINT
VAR A#=379
PRINT FORMAT$("%d*%d*%d=%d",A#,A#,A#,A#*A#*A#AND 1023);

言われて見ればmod 1024上での3の3乗根なんだけど、合成数mod上でn乗根を計算したことないので新鮮(実はフィボナッチ数を計算するためにmod p上で5の平方根を計算しようとしたことはあるけど、計算方法が分からなくて挫折したことがある)


とりあえず、2進数についてはある程度理解できたし、当初の疑問も解決しそうです。

印象としては、剰余環(mod n)と有限小数のハイブリットって感じですね。なんで今まで知らなかったんだろう…

まぁ名前ややこしいので無理もないですね。
ボクも最初AKS実装しようとしたとき、p-aryと思って気にしてなかったので…

しかし、p進数、普通に生きていれば触れることすらない概念なのでどういう書籍で勉強すればいいのか分からないですね。

とりあえず、当初の予定通りp進ニュートン法とp進FFTを理解したい。
https://qiita.com/qnighy/items/ba7e6cb6e321972bebf2
記事のコードをSmileBasicで書き直してみる。

'立っているビットの数を数える
DEF BITCNT(N)
　N=(N AND &H55555555)+(N AND &HAAAAAAAA)>>1
　N=(N AND &H33333333)+(N AND &HCCCCCCCC)>>2
　N=(N AND &H0F0F0F0F)+(N AND &HF0F0F0F0)>>4
　N=(N AND &H00FF00FF)+(N AND &HFF00FF00)>>8
　N=(N AND &H00FF00FF)+(N AND &HFF00FF00)>>8
END

'奇数Xに対しての平方数判定
DEF IS_SQR_ODD(X)
　VAR I%,Y%=1
　'逆平方根を求める
　FOR I%=1 TO 3
　　Y%=((3-Y%*Y%*X)*Y%)>>1
　　Y%=Y% AND 32767
　NEXT I%
　Y%=X*Y% AND 32767
　IF Y% AND 16384 THEN
　　Y%=32768-Y%
　ENDIF
　RETURN Y%*Y%==X
END

'平方数判定(奇偶数判定部)
DEF IS_SQR(X)
　VAR K%=BITCNT((X&-X)-1)
　IF 0==X THEN
　　RETURN 1
　ELSE IF K% AND 1 THEN
　　RETURN 0
　ELSE
　　RETURN IS_SQR_ODD(X>>K%)
　ENDIF
END