*> Combination

*> 全体でn個の選択から縦(または横)を
*> r個選ぶ組み合わせ

*> nCr
*> = nPr / r!
*> = n! / {(n - r)! * r!}

*> (W + H)C(W > H ? W : H)
*> = (W + H)! / {(W + H - biggerWH)! * biggerWH!}

PROGRAM-ID.                      ABC_999_D.
DATA                             DIVISION.
WORKING-STORAGE                  SECTION.
    01 INF     EXTERNAL PIC 9(10).

*> Fact
    01 F1 EXTERNAL.
        03 F PIC 9(10) OCCURS 300000 INDEXED BY FIX.

    01 INF2    PIC 9(10).

    01 INP     PIC X(34).

    01 N       PIC 9(6).
    01 AX      PIC 9(6).
    01 BX      PIC 9(6).
    01 K       PIC 9(12).

    01 a       PIC 9(6) BINARY.
    01 b       PIC 9(6) BINARY.

    01 tmp     PIC 9(12) BINARY.

    01 rm      PIC 9(12).
    01 dv      PIC 9(12).

    01 dvx     PIC 9(19).

    01 b_txt   PIC X(32).
    01 d_num   PIC 9(10).
    01 d_ret   PIC 9(10).

    01 i       PIC 9(18) BINARY.
    01 j       PIC 9(18) BINARY.

*> N!
    01 x       PIC 9(10).

*> smalla!
    01 smalla  PIC 9(10).
*> biga!
    01 biga    PIC 9(10).

*> smallb!
    01 smallb  PIC 9(10).
*> bigb!
    01 bigb    PIC 9(10).

*> smalla! * biga!
    01 y       PIC 9(10).

*> smallb! * bigb!
    01 z       PIC 9(10).

*> inverse a
    01 inva    PIC 9(10).

*> inverse b
    01 invb    PIC 9(10).

*> combination a
    01 comba   PIC 9(10).

*> combination b
    01 combb   PIC 9(10).

    01 m       PIC 9(10).
    01 sm      PIC 9(18).

    01 ans     PIC X(11).

    01 ZS      PIC Z(18)9.
    01 DUMMY   PIC X(1).

    01 XXX     PIC 9(30).

PROCEDURE                        DIVISION.
    ACCEPT INP.

    UNSTRING INP
    DELIMITED BY SPACE
    INTO N AX BX K.

*> 外部データ
*> 実行単位中の複数のプログラムで使用する定数を設定
    CALL "SET_EXTERNAL".

    SUBTRACT 2 FROM INF GIVING INF2.

    MOVE N TO d_num.

*> 階乗計算
    CALL "FACT_CALC"    USING BY CONTENT d_num
                              BY REFERENCE d_ret.

*> (a + b)!
    MOVE d_ret TO x.

    MOVE ZERO TO sm.

    PERFORM VARYING i FROM ZERO BY 1 UNTIL N < i

        IF ZERO = K
            ADD 1 TO sm
            EXIT PERFORM
        END-IF

        MOVE i TO a
        MULTIPLY AX BY a GIVING tmp

        IF K < tmp
            EXIT PERFORM
        END-IF

        SUBTRACT tmp FROM K GIVING rm

        DIVIDE BX INTO rm GIVING dv REMAINDER rm

        IF ZERO = rm AND dv <= N
            MOVE dv TO b
*> a!
            MOVE 1 TO smalla
            IF ZERO NOT = a
                SET FIX TO a
                MOVE F(FIX) TO smalla
            END-IF
*> (N - a)!
            MOVE 1 TO biga
            IF a NOT = N
                SUBTRACT a FROM N GIVING tmp

                SET FIX TO tmp
                MOVE F(FIX) TO biga
            END-IF
*> b!
            MOVE 1 TO smallb
            IF ZERO NOT = b
                SET FIX TO b
                MOVE F(FIX) TO smallb
            END-IF
*> (N - b)!
            MOVE 1 TO bigb
            IF b NOT = N
                SUBTRACT b FROM N GIVING tmp

                SET FIX TO tmp
                MOVE F(FIX) TO bigb
            END-IF
*> a! * (N - a)!
            MULTIPLY smalla BY biga GIVING dvx
            DIVIDE INF INTO dvx GIVING dvx REMAINDER y
*> 逆元
*> 1 / a! * (N - a)!
            CALL "MOD_CALC"     USING BY CONTENT y
                                      BY CONTENT INF2
                                      BY REFERENCE inva
*> b! * (N - b)!
            MULTIPLY smallb BY bigb GIVING dvx
            DIVIDE INF INTO dvx GIVING dvx REMAINDER z
*> 逆元
*> 1 / b! * (N - b)!
            CALL "MOD_CALC"     USING BY CONTENT z
                                      BY CONTENT INF2
                                      BY REFERENCE invb
*> combination a
            MULTIPLY x BY inva GIVING dvx
            DIVIDE INF INTO dvx GIVING dvx REMAINDER comba
*> combination b
            MULTIPLY x BY invb GIVING dvx
            DIVIDE INF INTO dvx GIVING dvx REMAINDER combb
*> combination a * combination b
            MULTIPLY comba BY combb GIVING dvx
            DIVIDE INF INTO dvx GIVING dvx REMAINDER m

            ADD m TO sm
            IF INF < sm
                SUBTRACT INF FROM sm
            END-IF

        END-IF

    END-PERFORM.

    MOVE sm TO ZS.
    PERFORM UNANS.

    DISPLAY ans(1:FUNCTION STORED-CHAR-LENGTH(ans)).
    STOP RUN.

UNANS                            SECTION.
    UNSTRING
        ZS DELIMITED BY ALL SPACE
        INTO DUMMY ans
    END-UNSTRING.

END PROGRAM ABC_999_D.

PROGRAM-ID.                      SET_EXTERNAL.
*> 外部データ
*> 実行単位中の複数のプログラムで使用する定数を設定
*> external data
*> Set constants used for multiple programs in running unit
DATA                             DIVISION.
WORKING-STORAGE                  SECTION.
    01 INF     EXTERNAL PIC 9(10).

    01 BNY     EXTERNAL PIC 9(1).
    01 BIT4    EXTERNAL PIC 9(1).
    01 BYT4    EXTERNAL PIC 9(1).
    01 BYT32   EXTERNAL PIC 9(2).
    01 BYT256  EXTERNAL PIC 9(3).
    01 HXD     EXTERNAL PIC 9(2).
    01 NUM_OFF EXTERNAL PIC 9(1).
    01 NUM_ON  EXTERNAL PIC 9(1).
    01 NUM_X   EXTERNAL PIC 9(1).
    01 OCT     EXTERNAL PIC 9(1).

PROCEDURE                        DIVISION.
    MOVE 998244353 TO INF.

    MOVE 2         TO BNY.
    MOVE 4         TO BIT4.
    MOVE 4         TO BYT4.
    MOVE 32  TO BYT32.
    MOVE 256 TO BYT256.
    MOVE 16  TO HXD.
    MOVE 0   TO NUM_OFF.
    MOVE 1   TO NUM_ON.
    MOVE 9   TO NUM_X.
    MOVE 8   TO OCT.

END PROGRAM SET_EXTERNAL.

PROGRAM-ID.                      FACT_CALC.
*> 階乗を計算する (a!)
DATA                             DIVISION.
WORKING-STORAGE                  SECTION.
    01 INF     EXTERNAL PIC 9(10).

    01 F1 EXTERNAL.
        03 F PIC 9(10) OCCURS 300000 INDEXED BY FIX.

*> 2
    01 BNY     EXTERNAL PIC 9(1).

    01 i       PIC 9(10) BINARY.

    01 tmp     PIC 9(19).
    01 d       PIC 9(19).
    01 m       PIC 9(10).

LINKAGE                          SECTION.
*> (IN)  d_num
*> (OUT) d_ret
    01 d_num   PIC 9(10).
    01 d_ret   PIC 9(10).

PROCEDURE                        DIVISION USING d_num d_ret.
    MOVE 1 TO i.
    SET FIX TO i.
    MOVE 1 TO F(FIX).

*> 2
    PERFORM VARYING i FROM BNY BY 1 UNTIL d_num < i

        MULTIPLY i BY F(FIX) GIVING tmp
        DIVIDE INF INTO tmp GIVING d REMAINDER m

        SET FIX UP BY 1
        MOVE m TO F(FIX)

    END-PERFORM.

    MOVE F(FIX) TO d_ret.

END PROGRAM FACT_CALC.

PROGRAM-ID.                      MOD_MULTI.
*> aとbを掛けた値をmodする (a * b mod p)
DATA                             DIVISION.
WORKING-STORAGE                  SECTION.
    01 INF     EXTERNAL PIC 9(10).

*> 4
    01 BIT4    EXTERNAL PIC 9(1).
*> 2
    01 BNY     EXTERNAL PIC 9(1).
*> 4
    01 BYT4    EXTERNAL PIC 9(1).
*> 32
    01 BYT32   EXTERNAL PIC 9(2).
*> 256
    01 BYT256  EXTERNAL PIC 9(3).
*> 16
    01 HXD     EXTERNAL PIC 9(2).
*> 8
    01 OCT     EXTERNAL PIC 9(1).
*> 1
    01 NUM_ON  EXTERNAL PIC 9(1).

    01 d       PIC 9(10).
    01 md      PIC 9(10).
    01 m       PIC 9(1).

LINKAGE                          SECTION.
*> (IN)  d_a
*>       d_b
*> (OUT) d_ret
    01 d_a     PIC 9(10).
    01 d_b     PIC 9(10).
    01 d_ret   PIC 9(10).

PROCEDURE                        DIVISION USING d_a d_b d_ret.
    MOVE ZERO TO d_ret.
    DIVIDE INF INTO d_a GIVING d REMAINDER md.

    PERFORM UNTIL d_b <= ZERO

        DIVIDE BNY INTO d_b GIVING d REMAINDER m

        IF NUM_ON = m
            ADD md TO d_ret
            IF INF < d_ret
                SUBTRACT INF FROM d_ret
            END-IF
        END-IF

*> 2
*> left shift
        MULTIPLY BNY BY md

        IF INF < md
            SUBTRACT INF FROM md
        END-IF

*> 2
*> right shift
        DIVIDE BNY INTO d_b
    END-PERFORM.
END PROGRAM MOD_MULTI.

PROGRAM-ID.                      MOD_CALC.
*> aのb乗をpで割った余りを計算する (a ^ b mod p)
DATA                             DIVISION.
WORKING-STORAGE                  SECTION.
    01 INF     EXTERNAL PIC 9(10).

    01 tmp     PIC 9(10).

    01 d       PIC 9(10).
    01 m       PIC 9(1).

*> 4
    01 BIT4    EXTERNAL PIC 9(1).
*> 2
    01 BNY     EXTERNAL PIC 9(1).
*> 4
    01 BYT4    EXTERNAL PIC 9(1).
*> 32
    01 BYT32   EXTERNAL PIC 9(2).
*> 256
    01 BYT256  EXTERNAL PIC 9(3).
*> 16
    01 HXD     EXTERNAL PIC 9(2).
*> 8
    01 OCT     EXTERNAL PIC 9(1).
*> 0
    01 NUM_OFF EXTERNAL PIC 9(1).

LINKAGE                          SECTION.
*> (IN)  d_a
*>       d_b
*> (OUT) d_ret
    01 d_a     PIC 9(10).
    01 d_b     PIC 9(10).
    01 d_ret   PIC 9(10).

PROCEDURE                        DIVISION USING d_a d_b d_ret.
*> O(log b)
    MOVE ZERO TO tmp.

    IF ZERO = d_b
        MOVE 1 TO d_ret
    ELSE
        DIVIDE BNY INTO d_b GIVING d REMAINDER m

        IF NUM_OFF = m
*> 偶数 bを半分にできる
            DIVIDE BNY INTO d_b

            CALL "MOD_CALC"     USING BY CONTENT   d_a
                                      BY CONTENT   d_b
                                      BY REFERENCE tmp

            CALL "MOD_MULTI"    USING BY CONTENT   tmp
                                      BY CONTENT   tmp
                                      BY REFERENCE d_ret
        ELSE
*> 奇数 -1乗から1乗をかける
            SUBTRACT 1 FROM d_b

            CALL "MOD_CALC"     USING BY CONTENT   d_a
                                      BY CONTENT   d_b
                                      BY REFERENCE tmp

            CALL "MOD_MULTI"    USING BY CONTENT   d_a
                                      BY CONTENT   tmp
                                      BY REFERENCE d_ret
        END-IF
    END-IF.
END PROGRAM MOD_CALC.
