include mpfr.e
atom t0 = time()
mpz z = mpz_init()
randstate state = gmp_randinit_mt()
for i=18 to 25 do
    integer pi = get_prime(i)
    mpz_ui_pow_ui(z,2,pi)
    mpz_sub_ui(z,z,1)
--  sequence f = mpz_prime_factors(z,20000000)
    sequence f = mpz_prime_factors(z)
    string zs = mpz_get_str(z),
           fs = mpz_factorstring(f)
    if fs!=zs then zs &= " = "&fs end if
    if length(f[$])=1 then -- (as per docs)
        mpz_set_str(z,f[$][1])
        if not mpz_probable_prime_p(z, state) then
            if length(f)=1 then
                zs &= " (not prime)"
            else
                zs &= " (last factor is not prime)"
            end if
        end if
    end if
    printf(1,"2^%d-1 = %s\n",{pi,zs})
end for
?elapsed(time()-t0)