\ extended gcd algorithm
: extended-gcd ( b a -- y x d )
    {: | _b _a :}
    over            ( b a -- b a b )
    0 = if          ( b a b -- b a )
        swap drop   ( b a -- a )
        0 swap      ( a -- 0 a )
        1 swap      ( 0 a -- 0 1 a )
    else            ( b a b -- b a )
        dup to _a   \ store a for later
        swap        ( b a -- a b )
        dup to _b   \ store b for later
        swap        ( a b -- b a )
        over        ( b a -- b a b )
        mod         ( b a b -- b a%b )
        swap        ( b a%b -- a%b b )
        recurse     ( a%b b -- y x d )
        rot rot     ( y x d -- d y x )
        over        ( d y x -- d y x y )
        _a _b /     ( d y x y -- d y x y a/b )
        *           ( d y x y a/b -- d y x [y * [a / b]] = e )
        -           ( d y x e -- d y [x - e] = f )
        swap        ( d y f -- d f y )
        rot         ( d f y -- f y d )
    then ;

\ returns -1 if both are equal to 0, 0 otherwise
: check-2-top-0 ( a b -- result )
    0 = if 0 = else drop 0 then ;

\ computes the solution of a linear diophantine equation
\ ax + by = c
\ if a solution exists this returns y x -1
\ otherwise 0 0 0
: diophantine ( c b a -- y x result )
    {: | _d _c _b _a :}
    over over                   ( c b a -- c b a b a )
    check-2-top-0 if            ( c b a b a -- c b a )
        drop drop drop          ( c b a -- )
        0 0 0                   ( -- 0 0 0 )
    else                        ( c b a b a -- c b a )
        rot to _c               ( c b a -- b a )
        over over to _a to _b   ( b a -- b a )
        extended-gcd            ( b a -- y x d ) \ b a extended-gcd 
        dup to _d               ( y x d -- y x d )
        _c swap mod             ( y x d -- y x c%d )
        0 = invert if           ( y x c%d -- y x )
            drop drop           ( y x -- )
            0 0 0               ( -- 0 0 0 )
        else                    ( y x c%d -- y x )
            _c _d /             ( y x -- y x [c/d] = e )
            dup rot             ( y x e -- y e e x )
            *                   ( y e e x -- y e [e*x] = f )
            rot rot             ( y e f -- f y e )
            *                   ( f y e -- f [y*e] = g )
            _b 0 < if -1 * then ( x0 y0 -- x0 -y0 )
            swap                ( x0 y0 -- y0 x0 )
            _a 0 < if -1 * then ( y0 x0 -- y0 -x0 )
            -1                  ( y0 x0 -- y0 x0 -1 ) \ x = x0; y = y0; result = -1 - true
        then
    then ;

\ prints and consumes the result of a diophantine equation
\ examples:
\ 0 result = → there's no solution
\ -1 result = → x = 0 ∧ y = 1 
: print-diophantine-result ( y x result -- )
    invert if ." there's no solution" drop drop else
        ." x = " . ." ∧ y = " .
    then ;

\ prints and consumes the number, adds a space if positive
\ examples:
\    3: ` 3`
\   -5: `-5`
: print-number ( a -- )
    dup 0 < invert if 32 emit else 45 emit -1 * then 48 + emit
    ;

\ prints and consumes the input of the diophantine equation
\ ax + by = c
: print-diophantine-input ( c b a -- )
    print-number ." x + "
    print-number ." y = "
    print-number
    ;

\ 3 nested loops in the range [-i, i]
\ diophantine equation results of all those
\ formatted nicely :-)
: diophantine-loop ( i -- )
    dup 1 + swap -1 * ( i -- i+1 -i )
    cr
    over over do
        over over do
            over over do
                i j k print-diophantine-input ."  → " i j k diophantine print-diophantine-result cr
                ." ———————————————————————————————————————" cr
            loop
        loop
    loop
    drop drop ( i+1 -i -- )
    ;

\ vim: ft=forth