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 
  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

{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. 

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 ;

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 5 over legendre -
  bs/mod bdrop 0= ;

100000 random nextprime test . -1  ok