Suppose there are numbers a b c d on the stack, then
locals| d c b a |
pop the values on the stack and store them in the locals. Normally there is no real need of locals in Forth, when factoring optimally, but when the stack is used for counted number series
n1 n2 ... nk k
locals is handy. And of course, locals could be used to uncomplicate algorithms.
The Pollard rho factoring routine for single cell numbers is fast, even in 64 bit systems, and can be used to define number theoretical functions. In BigZ the word
pollard# \ n -- p1 ... pk k
factorize the number and present it in a form that can be sorted by the word
sort \ n1 ... nk k -- m1 ... mk k
: maxprimefactor \ n -- p
pollard# sort
over >r drops r> ;
drops \ n1 ... nk k --
The radical of a number can be defined easily by factoring, dropping all copies of prime factors and multiply the rest of the factors:
: radical \ n -- r
1 swap \ just a value to be dropped at the end
pollard# sort \ p1 ... pk k sorted with largest on top of stack
1 swap 0 \ p1 ... pk 1 k 0
do undequ \ is two primes eual?
if nip \ drop the first of them (second on the stack)
else * \ multiply single prime
then
loop nip ; \ drop the number "1" used by undequ
The word
undequ \ a b c -- a b c flag
compares the second and the third values on the stack and flag is true if a=b else false.
{ 1000000 1001000 | all } ok
utime function radical transform-set utime d- d. -118123 ok
That is, transforming the set of the numbers {1000000,1000001,...,1000999} to the set of their radicals takes about a tenth of a second.
Now zdup cardinality . cr zet. gives
1000
{10,42,1034,1158,1954,3910,4119, ... ,1000995,1000997,1000999} ok
(Non of the numbers appears to have the same radical).
Also, it is easy to define a test for square free numbers
: sqrfree \ n -- flag
dup radical = ;
that's fairly fast
utime { 1 10000 | sqrfree } utime d- d. -2750161 ok
zdup cardinality . 6083 ok
cr zet.
{1,2,3,5,6,7,10,11,13,14,15, ... ,9993,9994,9995,9997,9998} ok
A nice word to analyse a sorted counted bundle of numbers of the stack is
: hist \ a1 ... ak k -- a1 ... ai i ak nk
2dup 0 locals| n k1 a k |
begin dup a = k1 and
while n 1+ to n
k1 1- to k1 drop
repeat k n - a n ;
that counts the uppermost copies of the same number, leaving the remaining counted bundle under the histogram value ak nk on the top of the stack, indicating nk copies of the number ak.
For example, define a function theta that gives the greatest factor of n that is a sum of two squares.
Facts:
any prime of the form 4n+1 can be written as a sum of two squares;
the product of two squaresums is a squaresum;
for primes p of the form 4n+3, p^2i is of the form 0²+b².
The word squsumfac gives the contribution from the prime factor pk.
: squsumfac \ pk nk -- fac fac=a²+b²
over 3mod4
if dup odd if 1- then
then ** ;
: theta \ n -- m
dup 1 = if exit then
1 locals| m |
oddpart dup 1 = \ r s flag, where n=s*2^r, s odd.
if swap lshift exit then
pollard# sort
begin hist squsumfac
m * to m dup 0=
until drop m swap lshift ;
The sets in BigZ can't have big number members (yet) but it might be interesting to create sets of single number factors of big numbers.
\ testing for small (fast) single number divisors
\ of big number w in the intervall n,m-1
: sfac \ w -- w ?v | m n -- f
beven if 2drop 2 bdup b2/ exit then
0 locals| flag |
0 locals| flag |
do bdup i pnr@ bs/mod 0=
if i pnr@ to flag leave then bdrop
loop flag ;
: sfacset \ b -- set
0 \ count of the number of elements
begin pi_plim 1 sfac ?dup
while >zst 2 - bnip
repeat bdrop >zst reduce ;
Testing a conjecture about divisibility of Fibonacci numbers:
: bsfib \ n -- F(n) single input and big output
dup 2 < if s>b exit then
bzero bone 1
do btuck b+ loop bnip ;
Conjecture:
dup 2 < if s>b exit then
bzero bone 1
do btuck b+ loop bnip ;
Conjecture:
Any prime number p<n divide some Fibonacci number F(m), 0<m≤n.
: fibtest \ m n -- flag
false locals| flag |
do i pnr@ 1+ 1
do j pnr@ i bsfib sfacset smember
if true to flag leave then
loop flag if leave then
loop flag ;
pi_plim . 1077871 ok
pi_plim . 1077871 ok
utime pi_plim 1 fibtest . utime d- d. -1 -12836944 ok (t<13 sec).
That is, the conjecture is true for all primes Pn where n<1077871.
That is, the conjecture is true for all primes Pn where n<1077871.
__________
Due to Wikipedia there is a formula:
p|F(p-i), where i=(5/p), the Legendre symbol.
: test \ p -- flag
dup
dup 5 over legendre -
bsfib
bs/mod bdrop 0= ;
100000 random nextprime test . -1 ok
Due to Wikipedia there is a formula:
p|F(p-i), where i=(5/p), the Legendre symbol.
: test \ p -- flag
dup
dup 5 over legendre -
bsfib
bs/mod bdrop 0= ;
100000 random nextprime test . -1 ok
