d*mod

: d*mod { D: a D: b D: m -- d }     \ multiplication modulo m
    a m dmod to a            \ normalize a
    b m dmod to b            \ normalize b
    a 1. d= if b exit then
    b 1. d= if a exit then
    a m a d- dumin { D: aA }
    b m b d- dumin { D: bB }
    aA d0= bB d0= or if 0. exit then
    aA 1. d= if m b d- exit then
    bB 1. d= if m a d- exit then
    aA a d= bB b d= and aA a d<> bB b d<> and or  { pos }  \ pos is True if positive, False otherwise
    aA UMAX bB d/ du<= if
        aA bB d* m dmod pos 0= if m 2swap d- m dmod then exit
    then
    aA d2/          { D: a0 }
    aA a0 d-        { D: a1 }
    bB d2/          { D: b0 }
    bB b0 d-        { D: b1 }
    a1 b1 m recurse { D: p4 }
    0. 0. 0.        { D: p1 D: p2 D: p3 }
    a0 a1 d= b0 b1 d= and if
        p4 to p1
        p4 to p2
        p4 to p3
    else
        a0 a1 d= if
            p4 m a1 d- m dmod m d+mod 2dup   to p3
                                             to p1
            p4                               to p2
        else
            p4 m b1 d- m dmod m d+mod        to p2
            b0 b1 d= if
                p2                           to p1
                p4                           to p3
            else
                p4 m a1 d- m dmod m b1 d- m dmod 1. m d+mod m d+mod m d+mod    to p1
                p4 m a1 d- m dmod m d+mod                                      to p3
            then
        then
    then
    p1 p2 p3 p4 m d+mod m d+mod m d+mod pos 0= if m 2swap d- m dmod then
;