' Copyrighted by YMOjisan
' パスカルの三角形 Part2

Dim M As Long
Dim P As Long
Dim N As Long


Sub Gの計算()
 P = InputBox("素数P=", "PASCAL part2 のデータ入力　1/3")
 If 素数ではない(P) Then End
 
 M = InputBox("素数べきM=", "PASCAL part2 のデータ入力　2/3")
 X = P ^ M
 If (X > 10000) Then MsgBox "大きすぎ": End
 
 N = InputBox("N=", "PASCAL part2 のデータ入力　3/3")
 Y = Int(Log(N) / Log(P)) + 1
 
 ReDim UP_K_M(Y, M)
 ReDim DP_K_M(Y, M)
 ReDim UP_K_M_Rk(Y, M)
 ReDim DP_K_M_Rk(Y, M)
 ReDim Rk(Y)
 
'------------------------
   PP1 = P * (P + 1) / 2
   PM1 = P * (P - 1) / 2
   
  K = 1
   UP_K_M(K, 0) = PP1
   DP_K_M(K, 0) = 0

   UP_K_M(K, 1) = 0
   DP_K_M(K, 1) = PM1
   
   For sm = 2 To M
    UP_K_M(K, sm) = 0
    DP_K_M(K, sm) = 0
   Next sm
   
  For K = 1 To Y - 1
    For sm = 0 To M
      UP_K_M(K + 1, sm) = PP1 * UP_K_M(K, sm) + PM1 * DP_K_M(K, sm)
      If sm > 0 Then
        DP_K_M(K + 1, sm) = PM1 * UP_K_M(K, sm - 1) + PP1 * DP_K_M(K, sm - 1)
       Else
        DP_K_M(K + 1, sm) = 0
      End If
    Next sm
    DP_K_M(K + 1, M) = DP_K_M(K + 1, M) + PM1 * UP_K_M(K, M) + PP1 * DP_K_M(K, M)
  Next K
'--------------------
 
 NXX = N
 For i = 0 To Y - 1
  Rk(i) = NXX Mod P
  NXX = Int(NXX / P)
 Next i
 
 K = 1
   Dk = Rk(0)
   UP_K_M_Rk(K, 0) = Dk * (Dk + 1) / 2
   DP_K_M_Rk(K, 0) = 0
   
   UP_K_M_Rk(K, 1) = 0
   DP_K_M_Rk(K, 1) = Dk * (P - 1 + P - Dk) / 2
   
   For sm = 2 To M
    UP_K_M_Rk(K, sm) = 0
    DP_K_M_Rk(K, sm) = 0
   Next sm

   For K = 1 To Y - 1
    Dk = Rk(K)
    DP1 = Dk * (Dk + 1) / 2
    DM1 = Dk * (Dk - 1) / 2
    DP2 = Dk * (2 * P + 1 - Dk) / 2
    DM2 = Dk * (2 * P - 1 - Dk) / 2
    
    For sm = 0 To M
      UP_K_M_Rk(K + 1, sm) = DP1 * UP_K_M(K, sm) + DM1 * DP_K_M(K, sm) _
                            + (Dk + 1) * UP_K_M_Rk(K, sm) + Dk * DP_K_M_Rk(K, sm)
      If sm > 0 Then
        DP_K_M_Rk(K + 1, sm) = DM2 * UP_K_M(K, sm - 1) + DP2 * DP_K_M(K, sm - 1) _
                            + (P - Dk) * DP_K_M_Rk(K, sm - 1) + (P - 1 - Dk) * UP_K_M_Rk(K, sm - 1)
       Else
        DP_K_M_Rk(K + 1, sm) = 0
      End If
    Next sm
    
    DP_K_M_Rk(K + 1, M) = DP_K_M_Rk(K + 1, M) + _
                          DM2 * UP_K_M(K, M) + DP2 * DP_K_M(K, M) _
                          + (P - Dk) * DP_K_M_Rk(K, M) + (P - 1 - Dk) * UP_K_M_Rk(K, M)
   Next K
'---------------------------------------------------
答え = UP_K_M_Rk(Y, M)

総数 = 0
For sm = 0 To M
 総数 = 総数 + UP_K_M_Rk(Y, sm)
Next sm

正しい総数 = N * (N + 1) / 2

If 正しい総数 <> 総数 Then MsgBox ("桁あふれか?")

MsgBox (N & "段のパスカル△の" & X & "の倍数の数は " & 答え)
正解 = 単純計算(X, N)
If 正解 = 答え Then MsgBox ("単純な方法でも " & 正解) _
Else MsgBox ("??? 単純な方法では " & 正解)
End Sub

Function 単純計算(X, N)
 Xの倍数の数 = 0
 ReDim 組合わせの数(N)
 
 組合わせの数(1) = 1
 組合わせの数(2) = 1
 
 For i = 3 To N
  For J = i - 1 To 2 Step -1
    組合わせの数(J) = (組合わせの数(J) + 組合わせの数(J - 1)) Mod X
    If 組合わせの数(J) = 0 Then Xの倍数の数 = Xの倍数の数 + 1
  Next J
  組合わせの数(i) = 1
 Next i
 
 単純計算 = Xの倍数の数
End Function

Function 素数ではない(P候補)
 素数ではない = True
 探索上限 = Sqr(P候補)
 
 For i = 2 To 探索上限
     If (P候補 Mod i) = 0 Then Exit Function
 Next i
 
 素数ではない = False
End Function