お題というか、協力してほしい感じなんですが、素因数分解関数をHaskellで書いて色んな数を素因数分解して遊んでいたら確認したい事実に出くわしたので。

31 <- 素数
331 <- 素数
3331 <- 素数
33331 <- 素数
333331 <- 素数

と、３が５個並んで末尾が１の数字までは素数という事が分かりましたが、いかんせん、ノートだと力不足。
それにCとかで書き直したらもっと先まで行けるかも？という事で、この先、どこまで33...31が素数なのかを調べて欲しいのです。
協力お願いします<(_ _)>

一応、Haskellではこんなコードです。

factorization n = f primes n
　where　primes = 2:(sieve [3,5..])
　　　　　　where sieve (p:xs) = p:(sieve [x | x <- xs, x `mod` p /= 0])

　　　　　f (p:ps) n | n <= p = [n]
　　　　　f (p:ps) n | n `mod` p == 0 = p:f (p:ps) (n `div` p)
　　　　　f (p:ps) n = f ps n