105 lines
3.7 KiB
Forth
105 lines
3.7 KiB
Forth
\ 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
|