オリジナル
整数の性質
(1) 5**4*x-2**4*y-1
(3) 5**5*x-2**5*y-1
(4)11**5*x-2**5*y-1
解法1:手計算とwolframのPowerModで
手計算なのに？逆元の計算にwolframalphaのPowerModを使いました。
(例)2元1次不定方程式
X=ai(mod mi) i=1,2
X=0 (mod 5^4)
X=1 (mod 2^4)
x=X/5^4
(1)
a m (M/m) (M/m)^-1 a*(M/m)*(M/m)^-1
0 ＜以下省略＞
1 2^4 5^4 1 5^4*1
x= 5^4 /5^4 = 1
y=(5^4-1)/2^4 = 39
x=2^4*k+ 1
y=5^4*k+39
(2)省略
(3)
1 2^5 5^5 29 5^5*29
x= 5^5*29 /5^5 = 29
y=(5^5*29-1)/2^5 =2832
x=2^5*k+ 29
y=5^5*k+2832
(4)
1 2^5 11^5 19 11^5*19
x= 11^5*19 /11^5= 19
y=(11^5*19-1)/2^5=95624
x= 2^5*k+ 19
y=11^5*k+95624
逆元の計算
結果
{1, 29, 19}
解法2: 解法1(手計算とwolframalphaのPowerModで)をvbaで
Private Function mul_inv(a As Long, n As Long) As Variant
If n < 0 Then n = -n
If a < 0 Then a = n - ((-a) Mod n)
Dim t As Long: t = 0
Dim nt As Long: nt = 1
Dim r As Long: r = n
Dim nr As Long: nr = a
Dim q As Long
Do While nr <> 0
q = r \ nr
tmp = t
t = nt
nt = tmp - q * nt
tmp = r
r = nr
nr = tmp - q * nr
Loop
If r > 1 Then
mul_inv = "a is not invertible"
Else
If t < 0 Then t = t + n
mul_inv = t
End If
End Function
Sub myLine(ir, aa As Long, mm As Long, bb As Long)
Cells(ir, 4) = mul_inv(bb, mm)
Cells(ir, 5) = aa * bb * Cells(ir, 4)
Cells(ir, 6) = (Cells(ir, 5) - 1) / mm
End Sub
Sub aaa_mi()
ActiveSheet.Cells.Clear
Cells(1, 1) = 1
Cells(1, 2) = 2 ^ 4
Cells(1, 3) = 5 ^ 4
Cells(2, 1) = 1
Cells(2, 2) = 2 ^ 5
Cells(2, 3) = 5 ^ 5
Cells(3, 1) = 1
Cells(3, 2) = 2 ^ 5
Cells(3, 3) = 11 ^ 5
For i = 1 To 3
Call myLine(i, Cells(i, 1), Cells(i, 2), Cells(i, 3))
Next
End Sub
'1 16 625 1 625 39
'1 32 3125 29 90625 2832
'1 32 161051 19 3059969 95624
解法3:WolframAlphaのChineseRemainderで
(1)ア～ウ
y = 39
(1)エ～ク
解
y = 664
(3)サ～ツ
解
y = 12207
(4)テ～ノ
y = 95624
解法4:Sageのcrtで
var('x,y,k,kk')
def mySage_Solve2(myA,myB,myC,myXmin):
mymy=mySage_Solve(myA,myB,myC)
kk=solve(abs(myB)*k+mymy[0]>=myXmin,k)[0][0].rhs()
kk=math.ceil(float(kk))
xo=abs(myB)*kk+mymy[0]
yo=solve(myA*xo+myB*y+myC,y)[0].rhs()
return xo,yo
def mySage_Solve(myA,myB,myC):
xo=crt([0,-myC],[myA,-myB])/myA
yo=solve(myA*xo+myB*y+myC,y)[0].rhs()
return xo,yo
print("#(1)ア～ウ",mySage_Solve ( 5**4,-2**4,-1) )
print("#(1)エ～ク",mySage_Solve2( 5**4,-2**4,-1, 10))
print("#(3)サ～ツ",mySage_Solve2( 5**5,-2**5,-1,100))
print("#(4)テ～ノ",mySage_Solve (11**5,-2**5,-1) )
m =mySage_Solve (5**4,-2**4,-1)[1]
KE=solve(625**2-5**x,x)[0].rhs()
KO=solve(625**2-(2**KE*m**2+2**x*m+1),x)[0].rhs()
print("#(2)ケ～コ",KE,KO)
#(1)ア～ウ (1, 39)
#(1)エ～ク (17, 664)
#(3)サ～ツ (125, 12207)
#(4)テ～ノ (19, 95624)
#(2)ケ～コ 8 5
解法5:sympyのsolve_congruenceで
from sympy import *
from sympy.ntheory.modular import solve_congruence
import math
var('x,y,k,kk')
def mySym_Solve2(myA,myB,myC,myXmin):
mymy=mySym_Solve(myA,myB,myC)
kk=solve(abs(myB)*k+mymy[0]-myXmin,k)[0]
kk=math.ceil(float(kk))
xo=abs(myB)*kk+mymy[0]
yo=solve(myA*xo+myB*y+myC,y)[0]
return xo,yo
def mySym_Solve(myA,myB,myC):
xo=solve_congruence((0, myA),(-myC,-myB))[0]//myA
yo=solve(myA*xo+myB*y+myC,y)[0]
return xo,yo
print("#(1)ア～ウ",mySym_Solve ( 5**4,-2**4,-1) )
print("#(1)エ～ク",mySym_Solve2( 5**4,-2**4,-1, 10))
print("#(3)サ～ツ",mySym_Solve2( 5**5,-2**5,-1,100))
print("#(4)テ～ノ",mySym_Solve (11**5,-2**5,-1) )
m =mySym_Solve (5**4,-2**4,-1)[1]
KE=solve(625**2-5**x,x)[0]
KO=solve(625**2-(2**KE*m**2+2**x*m+1),x)[0]
print("#(2)ケ～コ",KE,KO)
#(1)ア～ウ (1, 39)
#(1)エ～ク (17, 664)
#(3)サ～ツ (125, 12207)
#(4)テ～ノ (19, 95624)
#(2)ケ～コ 8 5
参考
(20220511)
ChineseRemainder「整数合同式...を満たすx>=dの最小のxを与える。」ありました。
原点の移動ができました。「原点の移動」は、世界基準？かもしれません。
(20220512)
solve_congruence