is_prime for INTI

[K.Kodama]

class INIT_EXT is

        --- Use is_prime as a function.

        is_spsp(x,a:INTI):BOOL is
                -- is prime or strong pseudo prime of base "a"?
                s:INT:=0; d:INTI:=x-1.inti;
                loop while!(d.is_even); s:=s+1; d:=d/2.inti; end;
                ad:INTI:=1.inti; --  ad=(a^d)%x
                loop while!(d.is_pos);
                        if d.is_odd then ad:=(ad*a)%x; end;
                        d:=d/2.inti; a:=(a*a)%x;
                end;
                if ad=1.inti then return true; end;
                loop r::=0.upto!(s-1);
                        if ad=(x-1.inti) then return true; end;
                        ad:=(ad*ad)%x; -- ad=a^(d*2^r)
                end;
                return false;
        end;

        is_fpsp(x,a:INTI):BOOL is
                -- is prime or Fermat's pseudo prime of base "a"?
                -- return (a^(x-1.inti)%x=1.inti);
                d:INTI:=x-1.inti; ad:INTI:=1.inti; an:INTI:=a;
                loop while!(d.is_pos);
                        if d.is_odd then ad:=(ad*an)%x; end;
                        d:=d/2.inti; an:=(an*an)%x;
                end;
                return ad=1.inti;
        end;

        is_prime(x:INTI):BOOL is
                -- Is __x__  prime?
                if x.is_even then
                        if x=2.inti then return true; else return false; end;
                elsif 3.inti.evenly_divides(x) then
                        if x=3.inti then return true; else return false; end;
                elsif 5.inti.evenly_divides(x) then
                        if x=5.inti then return true; else return false; end;
                elsif 7.inti.evenly_divides(x) then 
                        if x=7.inti then return true;else return false; end;
                elsif 11.inti.evenly_divides(x) then
                        if x=11.inti then return true; else return false; end;
                elsif 13.inti.evenly_divides(x) then 
                        if x=13.inti then return true;else return false; end;
                elsif x<2.inti then return false;
                end;
                -- Check if spsp or fpsp.
                if (x>1000000000.inti)and(~is_fpsp(x,13.inti)) then return 
false; end;
                d:INTI:=x.sqrt;
                -- r1:=5; r2:=7; -- \pm 1 (mod 6)
                -- r1:=11; r2:=13;
                r1:INTI:=17.inti; r2:INTI:=19.inti;
                loop while!(r1<=d);
                        if r1.evenly_divides(x) then return false;
                        elsif r2.evenly_divides(x) then return false;
                        end;
                        r1:=r1+6.inti; r2:=r2+6.inti;
                end;
                return true;
        end;

end;

-- 
K.Kodama(address@hidden)