## Monday, February 13, 2017

### How to use BigZ - part 1

I've started to use the standard ANS-Forth notation for locals. It's a bit awkward but awkward in a forthish way. When I started this blog I wasn't aware of this notation.

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

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 |
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:
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
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.
__________

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