Saturday, March 18, 2017

How to use BigZ - part 2

Gaps in increasing sequences of natural numbers

Strictly increasing sequences as

(1,2,3,...), (2,3,5,7,11,13,...) and (1,2,4,8,16,...)


can partly be represented as sorted sets in BigZ. Differences of sequences is an analogy to differentiation of functions. When defining a function gapz, to be the sorted set of all gaps between consecutive numbers of a set, also gapz become an analogy of the derivative of functions.


Apply gapz on an arithmetic sequence gives a set with a single element. If n times apply gapz on an increasing polynomial series of degree n gives a singleton set. On a sequence of exponential function with base two does nothing to the infinite set.


: gapz \ s -- s'
  0 locals| n |        \ counts the number of gaps in s'
  foreach 1+           \ prepare elements of s for the do-loop
  do zst> zst@ - >xst  \ the gap between the largest consecutive's
     n 1+ to n
  loop zst> drop       \ drop the smallest element of s
  n 2* negate >xst     \ calculate the set-count for s'
  xst zst setmove      \ move the set to zst
  set-sort reduce ;    \ sort and eliminate copies

{ 1 1000 | prime } gapz cr zet.

{1,2,4,6,8,10,12,14,18,20} ok



Partitions of a number n into distinct primes

: collincl \ s n -- s'
  0 >xst
  begin zst@
  while zsplit 
     dup >zst zfence zmerge
     set-sort reduce zfence
     xst zst setmove zmerge
     zst xst setmove
  repeat zdrop drop
  xst zst setmove
  reduce ;
\ include n in all sets in s

: xunion \ set --
  xst zst setmove union 
  zst xst setmove ;
\ Union of the top sets on the xst- and zst-stacks
\ is put on the xst-stack

: primeset \ m -- set
  pi dup 1+ 1
  ?do i pnr@ >zst 
  loop 2* negate >zst ;
\ Create the set of all primes < m+1

: memb \ s n -- flag
  false swap
  adn1 over + swap
  ?do dup i @ =
     if -1 under+ leave then cell
  +loop drop ;
\ Faster test if n is a member in the sorted number set s

For T being the set of primes:



The algorithm can be used with corrections for n=2p.

: termcase \ n -- flag
  case 2 of true endof
       3 of true endof
      11 of true endof
     dup of false endof
  endcase ; 
\ terminal cases: prime numbers without additional partitions

I have no proof that there are additional partitions for all primes greater than 11, but as far as the algorithm will go the terminal cases above are correct.

: z2@ \ set -- set n
  zst> zst@ swap >zst ;
\ read the largest element in the set

: lowlim \ set n -- set p
  0 swap adn1 over + swap
  ?do i @ under+ 2dup < 0= 
     if 2drop i @ leave then cell
  +loop ;
\ p is the smallest prime such that 2+3+5+...+p > n
  
: setsum \ set -- sum
  0 foreach ?do zst> + loop ;
\ The sum of all elements in set

: sumcorr \ s n -- s'
  locals| n |
  0 >xst
  begin zst@ 
  while zsplit zdup setsum n =
     if zfence xunion
     else zdrop
     then
  repeat zst> drop 
  xst zst setmove ;
\ Removes all partitions from s such that the sum < n

: dps \ n -- set
  dup 2 < if drop 0 >zst exit then 
  dup termcase if >zst -2 >zst -4 >zst exit then
  0 >xst
  dup primeset
  dup lowlim locals| low n |
  begin zst@ 
     if z2@ low <
        if false else true then
     else false 
     then
  while zsplit n zst> dup >r - ?dup 
     if recurse
        zst@ 
        if r> collincl n sumcorr xunion 
        else zst> drop r> drop 
        then 
     else { { r> } } xunion  
     then 
  repeat zdrop
  xst zst setmove 

  set-sort reduce ;
\ The set of partitions of n>0 into distinct primes


20 dps cr zet.

{{2,7,11},{2,5,13},{7,13},{3,17}} ok



50 dps cr zet.
{{2,7,11,13,17},{2,5,11,13,19},{7,11,13,19},{2,5,7,17,19},{3,11,17,19},{2,5,7,13,23},{3,11,13,23},{2,3,5,17,23},{3,7,17,23},{3,5,19,23},{2,3,5,11,29},{3,7,11,29},{3,5,13,29},{2,19,29},{3,5,11,31},{2,17,31},{19,31},{2,11,37},{13,37},{2,7,41},{2,5,43},{7,43},{3,47}} ok

: A000586 \ n -- 
  ." 1," 1+ 1 
  ?do i dps cardinality 0
     <# [char] , hold #s #> type 
  loop ;
\ List A000586 

100 a000586 cr
1,0,1,1,0,2,0,2,1,1,2,1,2,2,2,2,3,2,4,3,4,4,4,5,5,5,6,5,6,7,6,9,7,9,9,9,11,11,11,13,12,14,15,15,17,16,18,19,20,21,23,22,25,26,27,30,29,32,32,35,37,39,40,42,44,45,50,50,53,55,57,61,64,67,70,71,76,78,83,87,89,93,96,102,106,111,114,119,122,130,136,140,147,150,156,164,170,178,183,188,198, ok


Partitions of a number n into distinct non composites

A variant of the above.

: termcase1 \ n -- flag
  case 1 of true endof
       2 of true endof
     dup of false endof
  endcase ; 

: dps1 \ n -- set
  dup 0= if >zst exit then 
  dup termcase1 if >zst -2 >zst -4 >zst exit then
  0 >xst
  dup { 1 } primeset zmerge
  dup lowlim locals| low n |
  begin zst@ 
     if z2@ low <
        if false else true then
     else false 
     then
  while zsplit n zst> dup >r - ?dup 
     if recurse
        zst@ 
        if r> collincl n sumcorr xunion 
        else zst> drop r> drop 
        then 
     else { { r> } } xunion  
     then 
  repeat zdrop
  xst zst setmove 
  set-sort reduce ;

50 dps1 cr zet. 
{{1,3,5,11,13,17},{2,7,11,13,17},{1,2,3,5,7,13,19},{2,5,11,13,19},{7,11,13,19},{2,5,7,17,19},{1,2,11,17,19},{3,11,17,19},{1,13,17,19},{1,3,5,7,11,23},{2,5,7,13,23},{1,2,11,13,23},{3,11,13,23},{2,3,5,17,23},{1,2,7,17,23},{3,7,17,23},{1,2,5,19,23},{3,5,19,23},{1,7,19,23},{2,3,5,11,29},{1,2,7,11,29},{3,7,11,29},{1,2,5,13,29},{3,5,13,29},{1,7,13,29},{1,3,17,29},{2,19,29},{1,2,5,11,31},{3,5,11,31},{1,7,11,31},{1,2,3,13,31},{1,5,13,31},{2,17,31},{19,31},{1,2,3,7,37},{1,5,7,37},{2,11,37},{13,37},{1,3,5,41},{2,7,41},{2,5,43},{7,43},{1,2,47},{3,47}} ok

: test \ n --    n>0
  1+ 1 
  ?do i dps1 cardinality 0
     <# [char] , hold #s #> type 
  loop ;

100 cr test
1,1,2,1,2,2,2,3,2,3,3,3,4,4,4,5,5,6,7,7,8,8,9,10,10,11,11,11,13,13,15,16,16,18,18,20,22,22,24,25,26,29,30,32,33,34,37,39,41,44,45,47,51,53,57,59,61,64,67,72,76,79,82,86,89,95,100,103,108,112,118,125,131,137,141,147,154,161,170,176,182,189,198,208,217,225,233,241,252,266,276,287,297,306,320,334,348,361,371,386, ok

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

Friday, August 26, 2016

The abc-conjecture 1

Suppose a, b and c are natural numbers such that a,b,c are mutual co-prime and a+b=c. Those triples are called abc-triplets. The abc-conjecture concern the unusual possibility that

(1) a+b>rad(ab(a+b))

where rad is the radical, the product of all unique prime factors of a number. I.e. rad(4)=2, rad(6)=6, rad(60)=30, rad(81)=3, rad(101)=101, ...

There is an infinite number of pairs (a,b) as (1), but given a real number epsilon>0 there seems to be only a finite number of abc-triplets (a,b,a+b) such that

(2) a+b>rad(ab(a+b))^(1+epsilon) 

and that's one version of the famous abc-conjecture.

: abcpair \ a b -- flag
  locals| b a |
  a b ugcd 1 =
  a b + 
  dup a ugcd 1 =
  swap b ugcd 1 =
  and and ;

test if (a,b,a+b) is a abc-triplet, and

: unusual \ a b -- flag
  locals| b a |
  a b 2dup + * * radical
  a b + < ;

test if a+b>rad(ab(a+b)).

1 1000 condition non create-set zdup cardinality . 999  ok

creates the set {1,...,999}.

zdup cartprod zdup cardinality . 998001  ok

create the Cartesian product {1,...,999}x{1,...,999} and 

2dim < filter-set zdup cardinality . 498501  ok

filter the set so that the first component is less than the second.

2dim abcpair filter-set zdup cardinality . 303791  ok

This skip all (a,b) but those where (a,b,a+b) is a abc-triplet.

2dim unusual filter-set zdup cardinality . 32  ok

This is the remaining set of pairs such that (1):

zdup cr zet.
{(1,8),(1,48),(1,63),(1,80),(1,224),(1,242),(1,288),(1,512),(1,624),(1,675),(1,728),(1,960),(2,243),(3,125),(4,121),(5,27),(5,507),(7,243),(13,243),(25,704),(27,512),(32,49),(32,343),(49,576),(81,175),(81,544),(100,243),(104,625),(169,343),(200,529),(343,625),(640,729)} ok

Considering the pairs as Gaussian integers and transform the set of unusual pairs to their Gaussian norms give:

zdup 2dim gnorm transform-set zdup cardinality . cr zet. 32
{65,754,2305,3425,3970,6401,14657,15634,37186,50177,58565,59053,59098,59218,69049,82945,118673,146210,257074,262145,262873,302497,319841,334177,389377,401441,455626,496241,508274,529985,921601,941041} ok

Since the both sets have 32 elements I hasten to raise the conjecture:

(3) All (ordered) unusual pairs has unique Gaussian norms. 

To test the conjecture for different limits without stack overflow, conj3 works:

: abcunusual \ ab -- flag
  2dup abcpair 0= 
  if 2drop false
  else unusual
  then ;

: conj3 \ n -- set flag
  true locals| flag |
  0 >zst \ empty set on zst stack
  2 
  ?do i 1 
     ?do i j abcunusual 
        if i j gnorm dup zdup smember
           if false to flag drop i j pad 2! leave
           else >zst zfence union
           then
        then
     loop flag 0= if leave then
  loop flag ;

100 conj3 . -1  ok

zet. {65,754,2305,3425,3970,6401} ok

5000 conj3 . -1  ok
zdup cardinality . 87  ok
cr zet.
{65,754,2305,3425,3970,6401,14657,15634,37186,50177,58565,59053,59098,59218,69049,82945,118673,146210,257074,262145,262873,302497,319841,334177,389377,401441,455626,496241,508274,529985,921601,941041,1048577,1048601,1242793,1476226,1569061,1750393,1837097,2566561,2944705,3067769,3317074,4093154,4101154,4194385,4213625,4484017,4584929,4735097,4783069,4798594,4916545,5303810,5592434,5646001,5760001,5774602,5831545,8977273,8998393,9144577,9439993,9765746,9976306,10185529,11944561,12379505,13302409,13986466,14548594,15108770,15745025,15784466,15818497,15944098,16769026,16777337,16778441,17155426,18548777,19131877,23070401,23660897,23819585,24153953,33667138} ok

True so far. This take some time but I try 10000:

10000 conj3 . -1  ok
zdup cardinality . 129  ok
cr zet.
{65,754,2305,3425,3970,6401,14657,15634,37186,50177,58565,59053,59098,59218,69049,82945,118673,146210,257074,262145,262873,302497,319841,334177,389377,401441,455626,496241,508274,529985,921601,941041,1048577,1048601,1242793,1476226,1569061,1750393,1837097,2566561,2944705,3067769,3317074,4093154,4101154,4194385,4213625,4484017,4584929,4735097,4783069,4798594,4916545,5303810,5592434,5646001,5760001,5774602,5831545,8977273,8998393,9144577,9439993,9765746,9976306,10185529,11944561,12379505,13302409,13986466,14548594,15108770,15745025,15784466,15818497,15944098,16769026,16777337,16778441,17155426,18548777,19131877,23070401,23660897,23819585,24153953,26040898,31640674,31706945,32283521,33667138,34000562,37515986,37520281,38950162,39052481,39421505,40947202,40985921,43033601,43050817,43445377,44289026,44446210,47045882,47046137,47048690,56811506,64000361,64235537,65713618,66928882,66961570,68374489,68508353,70761674,73260281,73530626,85470281,87890626,88510465,93655426,94008377,96040001,96060226,118771553,123820633,140639489,141533305} ok

(To be continued)

Saturday, August 13, 2016

Instructions to create and manipulate sets

The basic words for sets are 

member \ element set -- flag
set= \ set1 set2 -- flag
subset \ set1 set2 -- flag

The parameters on the left sides above are located on the zst-stack. The 'element' could be a single (on the zst-stack) or a set. The flag is true or false on the ordinary stack. To check if a non negative single number on the parameter stack is a member in a set (located on the zst-stack), use the word

smember \ n set -- flag

3 { 1 2 3 4 } smember . -1  ok

Using member 

3 >zst { 1 2 3 4 } member . -1  ok

Except from set operations like

union \ s1 s2 -- s3
intersection \ s1 s2 -- s3
diff \ s1 s2 -- s3
powerset \ s1 -- s2     (s2 = set of all subsets of s1)
cartprod \ s1 s2 -- s3  (Cartesian product)
cardinality \ s1 -- n   (n = number of elements in s1)

and the possibility to create small sets

{ 10 10000 | prime } cardinality . 1225  ok

there are a lot of cryptic words to manipulate even rather big sets of non negative integers

intcond \ low hi xt -- | -- s   "intervall condition"
setcond \ xt -- | s -- s'       "set condition"
intimage \ low hi xt -- | -- s  "intervall image"
setimage \ xt -- | s -- s'      "set image"
paircond \ xt -- | s -- s'
pairimage \ xt -- | s -- s'
int2cond \ low hi n xt -- | -- s   "intervall two-condition"
int2image \ low hi n xt -- | -- s  "intervall image"
set2cond \ n xt -- | s -- s'       "set condition"
set2image \ n xt -- | s -- s'      "set image"

Inspired by the idea of orthogonal instructions I created a simple syntax just for set manipulations, to summarize and simplify the use of the words above:


variable zp

variable cf2

: condition  ' 0 cf2 ! sp@ zp ! ;
: function  ' 2 cf2 ! sp@ zp ! ;
: 2dim  ' -1 cf2 ! sp@ zp ! ;
: syntax  sp@ zp @ - 0= if 0 0 else 1 then cf2 @ ;

The word syntax check the number of input parameters and put a dummy parameter on the stack when needed.

(I have started to use the ANS-Forth notation for locals and will reform the code and the blog. It's awkward, but "forthish" awkward since it loads the parameters in the reversed direction: 

k l m n locals| n m l k |

It has some advantages even if it looks strange.)

The ten cryptic words are now squeezed into the three words

create-set
filter-set
transform-set


\ e.g. 1 20 condition < 7 create-set

: create-set \ m n xt nr -- set
  syntax locals| cf k nr xt n m |
  k cf or
  case 0 of m n xt intcond endof
       1 of m n nr xt int2cond endof
       2 of m n xt intimage endof
       3 of m n nr xt int2image endof
  endcase ;

\ e.g. condition > 5 filter-set
: filter-set \ set xt nr -- set'
  syntax locals| cf k nr xt |
  cf 0< if xt paircond exit then k
  case 0 of xt setcond endof
       1 of nr xt set2cond endof
  endcase ;

\ e.g. 2dim + transform-set
: transform-set \ set xt nr -- set'
  syntax locals| cf k nr xt |
  cf 0< if xt pairimage exit then k
  case 0 of xt setimage endof
       1 of nr xt set2image endof
  endcase ;

Creating sets:


1 25 condition prime create-set cr zet.

{2,3,5,7,11,13,17,19,23} ok

10 29 condition coprime 6 create-set cr zet.
{11,13,17,19,23,25} ok


1 10 function 2* create-set cr zet.

{2,4,6,8,10,12,14,16,18} ok

1 10 function * 3 create-set cr zet.
{3,6,9,12,15,18,21,24,27} ok

The word after condition/function must be a defined condition or function for one or two parameters.

Filter sets:

{ 1 25 | all } condition odd filter-set cr zet.
{1,3,5,7,9,11,13,15,17,19,21,23} ok

{ 1 25 | all } condition < 10 filter-set cr zet.
{1,2,3,4,5,6,7,8,9} ok

{ 1 2 3 } zdup cartprod zdup cr zet.
{(3,3),(3,2),(3,1),(2,3),(2,2),(2,1),(1,3),(1,2),(1,1)} ok
2dim < filter-set cr zet.
{(1,2),(1,3),(2,3)} ok

Transform sets:

{ 1 100 | prime } function 2* transform-set  ok
function 1- transform-set  ok
condition prime filter-set cr zet.
{3,5,13,37,61,73,157,193} ok

1 1000000 condition prime create-set  ok
condition < 50 filter-set  ok
function + 2 transform-set cr zet.
{4,5,7,9,13,15,19,21,25,31,33,39,43,45,49} ok

{ 1 10 | odd } zdup cartprod  ok
2dim + transform-set cr zet.
{2,4,6,8,10,12,14,16,18} ok

The word 2dim works like condition and function but acts on a set of pairs and can be used both for 2-dim conditions and 2-dim functions.

So far those instructions don't works in definitions and they just works for numbers and pairs of numbers.

Saturday, June 18, 2016

String operations and bioinformatics

Strings makes it possible to generalize the concept of sets. In BigZ a set is a nested set of nested sets and lists like

{ 123 ( 234 { 345 456 { 567 678 } } ) } cr zet. 
{123,(234,{345,456,{567,678}})} ok

and the only lack of generality concern the atomic elements, which must be non negative single numbers. But virtually anything can be denoted as a string which can be interpreted as a list of characters:

s" {Hello world!,How are you?}" >str stringset>zet cr zet.
{(72,101,108,108,111,32,119,111,114,108,100,33),(72,111,119,32,97,114,101,32,121,111,117,63)} ok
In this way also sets of big integers, Gaussian integers etc can be elements of sets.

A nice way to handle strings in Forth is using a string stack, which in this implementation consists of two stacks, one for the arrays of ASCII signs and one for addresses to the arrays of signs.

>str \ ad n -- string    Push a string on the stack
str> \ string -- ad n    Pop a string from the stack
str@ \ string -- string | -- ad n

sempty \ string -- string | -- flag
.str \ --  Prints the stack without changing it
str. \ str --  Print and drop the topmost element

sdup sdrop sover snip sswap srot stuck spick does the normal operations.
soover \ str1 str2 str3 -- str1 str2 str3 str1
A shorter way to enter strings from commando line is
s Hello world"
However, in definitions one must use 
s" Hello world" >str
Some words for string manipulations:
s& \ s1 s2 -- s1&s2   Concatenation
sleft \ s1 -- s2 | n --  Skip all but the n leftmost characters
sright \ s1 -- s2 | n -- The samr for the n rightmost chars
ssplit \ s -- s' s" | n --  split string after the nth letter
sanalyze \ s1 s2 -- s1 s3 s1 s4 / s2 | -- flag 
split s2 if s1 is a part of s2 and if true flag then s2=s3&s1&s4.
substring \ s1 s2 -- s1 s2 | -- flag
sreplace \ s1 s2 s3 -- s4    Replace s2 with s1 in s3
scomp \ s1 s2 -- | -- n    -1:s1>s2, +1:s1<s2, 0:s1=s2
snull \ -- emptystring
schr& \ s -- s' | ch --   Concatenate ch to top string
slen= \ s1 s2 -- | -- flag   Test if same length
strail \ s -- s'  Remove trailing spaces
>capital \ ch -- ch'  Change common to capital
>common \ ch -- ch'  The oposite 
capital \ ch --flag  Test if capital letter
common \ ch -- flag  Test if common letter
slower \ s -- s'  Change to lower in string
supper \ s -- s'  Opposite as above
str>ud \ s -- s' | -- ud flag   Unsigned double from string
str>d \ s -- s' | -- d flag     Double from string
snobl \ s -- s'      Remove all blanks
sjustabc \ s -- s'   Remove all signs but eng. letters
alphabet \ s -- s'   Gives the alphabet of string
zet>stringset \ set -- string
stringset>zet \ string -- set
sunion \ str1 str2 -- str3
sintersection \ str1 str2 -- str3
sdiff \ str1 str2 -- str3
s {brown,red,orange,yellow,green}"  ok
s {blue,violet,brown,black}"  ok
sunion str. {black,brown,violet,blue,green,yellow,orange,red} ok
hamming \ s1 s2 -- s1 s2 | n   The Hamming distance
editdistance \ s1 s2 -- s1 s2 | n   The Levenshtein distance

This code is now included in the BigZ code.

Friday, May 27, 2016

BigZ, Zet with big integers included

Forth is a very special computer language - a kind of smart macro assembler. Even if the programming is on high level you always program directly upon an addressable part of the memory. The visible stack orientation is the simple way to handle data. Postfix notation makes brackets and operator priorities unnecessary. There is no black box parser limiting what is allowed and possible. The programmer decide virtually everything about the procedures, input, output and memory. And this makes it very easy to extend the system with new data types with postfix algebra.

In BigZ (see the top bar) a system for big integers is included to Zet. It's simple, efficient and rather complete, in spite of my limited programming ability. As for integers, floating point numbers and sets, there are stacks for numbers of dynamical length. And as for sets, there is no other limit of the size than the allocated memory. When needed the memory is reallocated. 

The stack for big numbers is really two stacks in the same area of memory: one growing towards high memory (the numbers) and one growing towards low memory (the addresses to the numbers). Writing

b 702486742867487684278678476028746724601 

and pressing enter, reads the number string and convert it to a multi-decimal number A_0*B^0+...+A_r*B^r where B is 2^b and b is the maximal number of bits of a single cell integer. The single cell numbers A_i are stored in the stack with i growing towards hi memory.

The operators have the usual names but with the prefix b.

b 2000 bfaculty cr b.

fill the screen with a 5736 figure number within some tenth of a second.

The word 

b**mod \ -- | a n m -- b

counts aⁿ(mod m) and 

b 26359783991551070871965201979080333254038743646746158379582192038055842791146833506745978666309678710238746262325665407448047112858614221184120023774728850927701745782077979943165434355776993447809155163506304287949484786229043007193369097865681445643720004387345872800008950502312482268122160708155160328564   ok
b 39327375191467048647521018841730348998598651522372708158670691892420060531890655533991461398550696324722351925617448324344083141484661951392820800479947685042791549748743564268081958080246132723174581232969062661990556972176861792341905425252382562697686127413259201904144867482279552760624394742040590855602   ok
b 16306675630939784626502807161425612212340131109932460322238576590610976549235432503764166590338557086651302692446204937053886057954003032814801648230668894753863150108029570966337696939560565101312815273803547555619029325583815565767168841836143511237512006023630352590540140638620838094344583215840893265550   ok
b**mod cr b.
13589152355661418800480923196880892958242540305591911682228289531245831992616543540153046075650497423341388067834762342417455873016510076223412593080539625696104324076486611754715629613440368922879919737699812253207474005960378943204378352266794545844684513362891786202126506784931492440650088987888029494246  ok

executes in about 0.5 seconds, which might be fast enough for some encryption experiments.

Some of the most important words are:

bdrop, bdup, bover, brot, bnip, btuck, bswap, b+, b-, b*, b/, b=, b<, b>, b0=, bmod, b/mod, bsqrtf, bgcd and b**mod.

with obvious operations. The word .b prints the b-stack. There are also mixed words:

bs* \ n -- | b -- n*b
bs/mod \ n -- r | b -- q, where b=nq+r

which is much faster than the corresponding b* and b/mod.

Ahead I will try so submit adequate prime tests and factoring functions.

Sunday, May 15, 2016

Words for testing conjectures

Making code for testing conjectures can be cumbersome even in the case of easy programming. So far in Zet there is the possibility of make and calculate with sets of numbers interactively.

{ 1 1000 | prime } ok
{ 1 1000 | pairprime } { 1 1000 | notpairprime }
union ok

zet= . -1 ok

Conditions so far are
: all dup = ;
: odd 1 and ;
: 1mod4 4 mod 1 = ;
: 3mod4 4 mod 3 = ;
: sqr dup sqrtf dup * = ;
: sqrfree dup radical = ;
: pairprime dup prime over 2 + prime rot 2 - prime or and ; 
: notpairprime dup prime swap pairprime 0= and ;
: semiprime bigomega 2 = ;  \ A product of two primes?
: uniprime smallomega 1 = ; \ Only divisional by one prime?
: biprime smallomega 2 = ;  \ Exact two different primes?


The construction { 1 10000 | pairprime } is fancy but slow and risk overflow in data stack. All the pairprimes in the intervall will first be created on the stack and then be moved to the zst-stack. It's better to check number for number and create the set directly on the zst-stack.

: intcond \ low hi xt -- | -- s   "intervall condition"

  loc{ xt } 
  swap 0 -rot
  do i xt execute 
     if i >zst 1+ then
  loop 2* negate >zst ;

utime 1 100000 ' pairprime intcond utime cr d- d. cardinality .


-35954 2447  ok

A set of 2447 primes is created in about 0.04 seconds. This construction is also possible to use in definitions, then using ['] instead of '.

To filtrate a set on the zst-stack:

: setcond \ xt -- | s -- s'       "set condition"
  loc{ xt } 0
  foreach
  do zst> dup xt execute
     if >xst 1+ else drop then
  loop dup 0
  do xst> >zst 
  loop 2* negate >zst ;

{ 1 100 | prime } ' 1mod4 setcond cr zet.


{5,13,17,29,37,41,53,61,73,89,97} ok

It's also nice to be able to create the image of a function:

: intimage \ low hi xt -- | -- s  "intervall image"
  loc{ xt } 
  swap 2dup
  do i xt execute >zst
  loop - 2* negate >zst
  set-sort reduce ;

: setimage \ xt -- | s -- s'      "set image"

  loc{ xt } 0
  foreach 
  do zst> xt execute >xst 1+
  loop dup 0
  do xst> >zst
  loop 2* negate >zst
  set-sort reduce ;


Functions so far are: 

log~ ( n -- nr ) where nr=1+²log n
random ( u1 -- u2 ) where 0≤u2<u1
nextprime ( numb -- prime )
prevprime ( numb -- prime )
sqrtf ( m -- n ) "floor"
sqrtc ( m -- n ) "ceiling"
radical ( n -- r )
totients ( n -- t )
bigomega ( n -- b )
smallomega ( n -- s )
ufaculty ( u -- u! )
pnr@ ( n -- p ) prime number n
pi ( x -- n ) number of primes ≤ x

Functions and conditions both must have the stackdiagram ( m -- n ), but the concept will be generalized.

1 20 ' radical intimage zet. {1,2,3,5,6,7,10,11,13,14,15,17,19} ok

Some test functions:

: square dup * ;                \ x → x²
: sqr>prime square nextprime ;  \ x → nextprime(x²)
: sqr<prime square prevprime ;  \ x → prevprime(x²)
: foo dup totients mod ;        \ x → x(mod ϕ(x)) Euler's totient.

{ 1 100 | all } ' foo setimage cr zet.
{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,21,22,23,25,27,31,33,35,39} ok

1 100 ' square intimage ' foo setimage cr zet.
{0,3,5,7,11,13,17,19,20,23,27,28,29,31,37,41,43,44,47,52,53,59,61,67,68,71,73,76,79,80,83,89,92,97,105,112,116,124,125,148,164,172,176,180,188,189,208,243,252,272,304,320,343,368,385,396,429,448,468,500,585,704,720,825,945,969,1008,1105,1197,1280,1309,1372,1540,1620,1701,1725,1729,1785,2185,2187,2625,2697,3069,3861} ok

Hmm, it seems like all odd primes less than 100 belongs to the image...

1 10000 ' square intimage ' foo setimage ok
1 10000 ' prime intcond ok
zswap diff zet. {2} ok

So I asked Mathematics stack exchange about it. (: 






Well, it might be sound to expect non dramatic explanations to conjectures, especially conjectures concerning primes.

To check relations R m there is a need for testing subsets of Cartesian products, sets of pairs of integers.

: paircond \ xt -- | s -- s'
  loc{ xt } 0
  foreach
  do zdup zet> drop xt execute
     if zst xst setmove 1+ else zdrop then
  loop 6 * negate >xst
  xst zst setmove ;

{ 1 10 | all } zdup cartprod ' = paircond cr zet.

{(1,1),(2,2),(3,3),(4,4),(5,5),(6,6),(7,7),(8,8),(9,9)} ok

: pairimage \ xt -- | s -- s'
  loc{ xt } 0
  foreach
  do 1+ zet> drop xt execute >xst
  loop dup 0 
  do xst> >zst
  loop 2* negate >zst
  set-sort reduce ;

{ 2 10 | all } zdup cartprod ' * pairimage cr zet.
{4,6,8,9,10,12,14,15,16,18,20,21,24,25,27,28,30,32,35,36,40,42,45,48,49,54,56,63,64,72,81} ok

Some conditions and functions N²→N are =, <>, <, >, >=, <=, +, *, /, mod, **, ugcd, -, invmod, legendre, jacobi, kronecker, gnorm, choose, where m ≥ n for m n - and m n must be coprime for m n invmod. the gnorm of two integers m n is the norm of the gaussian integer m+in, that is the number m²+n².

: coprime ugcd 1 = ;
: divide swap mod 0= ; 

{ 1 10 | all } zdup cartprod ' coprime paircond cr zet.
{(1,1),(1,2),(1,3),(1,4),(1,5),(1,6),(1,7),(1,8),(1,9),(2,1),(2,3),(2,5),(2,7),(2,9),(3,1),(3,2),(3,4),(3,5),(3,7),(3,8),(4,1),(4,3),(4,5),(4,7),(4,9),(5,1),(5,2),(5,3),(5,4),(5,6),(5,7),(5,8),(5,9),(6,1),(6,5),(6,7),(7,1),(7,2),(7,3),(7,4),(7,5),(7,6),(7,8),(7,9),(8,1),(8,3),(8,5),(8,7),(8,9),(9,1),(9,2),(9,4),(9,5),(9,7),(9,8)} ok

{ 1 10 | all } zdup cartprod ' divide paircond cr zet.
{(1,1),(1,2),(1,3),(1,4),(1,5),(1,6),(1,7),(1,8),(1,9),(2,2),(2,4),(2,6),(2,8),(3,3),(3,6),(3,9),(4,4),(4,8),(5,5),(6,6),(7,7),(8,8),(9,9)} ok

{ 1 10 | all } zdup cartprod ' coprime paircond ' gnorm pairimage cr zet.
{2,5,10,13,17,25,26,29,34,37,41,50,53,58,61,65,73,74,82,85,89,97,106,113,130,145} ok