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

Saturday, May 14, 2016

Simple graphs

The idea of using stacks for sets is not a bad idea - in my opinion. Besides from the simplified garbage collecton there is the benefit that the data can be accessed in arrays (in the stacks), which at least can be used to define some fast primitive routines. The set routines are not slow. When enumerating all the cards in a deck from 0 to 51, Zet can calculate the set of all 1326 possible hold cards in Texas hold´em poker in a few hundredths of a second. On my Android:

{ 0 52 | all } utime 2 power# utime 2swap d- cr d. cardinality cr .

29099 
1326

Here utime counts in μs.


Formally a simple graph is a set of vertices V and a subset of all unordered pairs of vertices E. Visually, a simple graph is a collection of verticies joined by zero or one edges and which consist no loops.


A subgraph (V,E) of a simple graph (V',E') is a graph such that V is a subset of V' and E of E'.


: subgraph \ -- flag | (V,E) (V',E') -- 

  unfence zrot unfence 
  zrot subset
  zswap subset and ;

There is a maximal subgraph for each subset V generated by E'.


\ E = intersection of E' and power#(2,V)

: edges~ \ E' V -- E
  2 power# intersection ;

But this straightforward implementation is inefficient. About 20 times faster is:


: edges \ E' V -- E

  0 >xst 
  zst yst setmove 
  foreach \ {u,v}∈E'
  ?do zdup unfence yzcopy1 member 
     if yzcopy1 member
        if zfence xzmerge
        else zdrop
        then
     else zst> drop zdrop
     then
  loop yst setdrop xst zst setmove ;

To make a random simple graph with v vertices and with an edge between two vertices in m cases of n:


: randgraph \ m n v -- | -- (V,E)

  loc{ m n v } 
  0 >xst
  { v 0 do i 1+ loop } 
  zdup 2 power# foreach \ {u,v}
  do n random m < 
     if zfence xzmerge
     else zdrop
     then
  loop xst zst setmove pair ;

The word extend creates a superset to the graph created by edges where all edges connected to V is submitted plus all edges connected to the submitted points.


\ V={x∈V'|y∈V" & {x,y}∈E'}

\ E={{x,y}∈E'|x∈V & y∈V}
: extend \ E' V" -- (V,E)
  zswap zst yst setmove 
  zst xst setcopy 
  foreach \ v∈V"
  do zst> yzcopy1
     begin zst@ 
     while zsplit zdup dup smember
        if xzmerge
        else zdrop
        then
     repeat zet> 2drop
  loop xst zst setmove 
  set-sort reduce
  yst zst setmove 
  zover edges pair ;

Counts all isolated points in a graph:

  
: isolated-vertices# \ -- n | (V,E) --
  unfence 0 dup loc{ flag }
  zst yst setmove
  foreach
  do zst> yzcopy1 true to flag
     begin zst@ flag and
     while zsplit dup smember 0= to flag
     repeat zdrop drop flag -
  loop yst zst setmove zdrop ;

4 5 9 randgraph zdup cr zet.

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

isolated-vertices# . 0  ok

Counts all isolated components in a graph:
  
: components# \ -- n | (V,E) --
  zdup 0 >xst
  unfence 
  znip zst yst setcopy
  foreach
  do begin yzcopy1 zover
        extend unfence zdrop ztuck zet=
     until zfence xzmerge
  loop yst setdrop 
  xst zst setmove reduce cardinality
  isolated-vertices# + ;

Due to the formula for vertices, edges and components for a forest, that is, a graph without circuits, v=e+c:


: forest? \ -- flag | (V,E) -- 

  zdup unfence 
  cardinality \ e
  cardinality \ v
  components# \ c
  rot + = ;

4 5 9 randgraph zdup cr zet.
({1,2,3,4,5,6,7,8,9},{{7,9},{6,9},{5,9},{4,9},{3,9},{2,9},{6,8},{5,8},{4,8},{3,8},{2,8},{1,8},{5,7},{4,7},{3,7},{1,7},{5,6},{4,6},{3,6},{2,6},{1,6},{4,5},{3,5},{2,5},{1,5},{3,4},{2,4},{1,4},{2,3},{1,2}}) ok

forest? . 0  ok

\ Using set-sort to sort a vector

: vector-sort \ s -- s'
  set-sort zst> 1- >zst ;


\ check if E is a cycle 

: cycle \ -- flag | E --
  zdup multiunion
  zdup cardinality true loc{ v flag }
  zover zdup cardinality v = 0=
  if triplet zdrop false exit
  then pair components# 1 >
  if zdrop false exit
  then 0 >xst foreach
  do xzmerge
  loop xst zst setmove
  zet> cs sort 2 - 0
  do over = flag and to flag
     over > flag and to flag
  +loop = flag and ;


: clear-table \ s --
  pad 0 foreach
  do zst> max
  loop cells erase ;
 
: cyc!check \ n -- flag
  cells pad + 1 over +! @ 2 > ;


\ Test if (V,E) is 2-regular 

: 2-regular \ -- flag | (V,E) --
  unfence zswap clear-table
  begin zst@
  while zsplit unfence
     zst> cyc!check if zst> drop zdrop false exit then
     zst> cyc!check if zdrop false exit then
  repeat zdrop true ;

4 5 9 randgraph zdup cr zet.
({1,2,3,4,5,6,7,8,9},{{8,9},{7,9},{6,9},{4,9},{3,9},{2,9},{7,8},{6,8},{5,8},{4,8},{2,8},{1,8},{6,7},{4,7},{2,7},{1,7},{4,6},{2,6},{1,6},{4,5},{3,5},{1,5},{3,4},{1,4},{1,3},{1,2}}) ok

2-regular . 0  ok

Tuesday, May 3, 2016

Fast generation of the symmetric and alternating groups

From Wikipedia I got this algorithm, how to generate all permutation in alphabetical order:

1. Find the largest index k such that a[k]<a[k+1]. If no such index
   exists, the permutation is the last.
2. Find the largest index l greater than k such that a[k]<a[l].
3. Swap the value of a[k] with that of a[l].
4. Reverse the sequence from a[k+1] up to and including the final
   element a[n].

First a word that reverse the order of all n characters starting at address ad.

: reverse-string \ ad n --
  2dup + 1- loc{ ad1 n ad2 } n 2/ 0
  ?do ad1 i + c@ ad2 i - c@ 
     ad1 i + c! ad2 i - c!
  loop ; 

Then the 1'st part of the algorithm, returning the address corresponding to the index k if it exists or else return 0.

: lex-perm1 \ ad n -- a1
  0 loc{ a1 } 2 - over + 
  do i c@ i 1+ c@ <
     if i to a1 leave then -1
  +loop a1 ;

Find the largest address a2 greater than a1 such that [a1]<[a2].

: lex-perm2 \ ad n a1 -- a2
  0 loc{ a1 a2 } 1- over +
  do a1 c@ i c@ <
     if i to a2 leave then -1
  +loop a2 ;

Swap the values at addresses a1 and a2.

: lex-perm3 \ a1 a2 --
  over c@ over c@
  swap rot c!
  swap c! ;

Reverse the order of the last characters, from address a1 to the end.

: lex-perm4 \ ad n a1 -- 
  reverse from a1+1 to ad+n-1 
  1+ -rot            \ a1+1 ad n
  + over -           \ a1+1 ad+n-(a1+1) 
  reverse-string ; 

Calculate the next permutation:
  
: nextp \ ad n -- 
  2dup 2dup          \ ad n ad n ad n
  lex-perm1 dup 0=
  if 2drop 2drop drop exit 
  then dup >r        \ ad n ad n a1
  lex-perm2 r>       \ ad n a2 a1
  tuck swap          \ ad n a1 a1 a2
  lex-perm3          \ ad n a1
  lex-perm4 ;

Create the string 123...n:

: n>str \ n -- ad n
  dup 0 do i 49 + pad i + c! loop pad swap ;

Create a vector on the z-stack from the string.

: str>vect \ ad n -- | -- s
  loc{ ad n } n dup 0
  do ad i + c@ 15 and >zst loop 2* 1+ negate >zst ;

Fast calculation of the symmetry group of n! permutations.

: sym \ n -- | -- s
  n>str loc{ ad n }
  n dup ufaculty dup 0
  do ad n str>vect 
     ad n nextp
  loop swap 1+ * 2* negate >zst ;

utime 7 sym cardinality . utime d- d. 5040 -3931  ok

What would take hours with straight forward generation is now done in 4 milliseconds.

Next word calculates how many components in the vector s that is greater than the number m:

: perm> \ m -- n | s --
  loc{ m } 0
  foreach do zst> m > + loop negate ;

This is used to calculate the number of pairs of components in the vector s that is unsorted:
  
: #perm \ -- n | s -- 
  0
  begin zst@ -3 <
  while zsplit zst> zdup perm> +
  repeat zdrop ;

Which determine if the vector correspond to an odd permutation:

: oddperm \ -- flag | s --
  #perm 1 and ; 

: alt \ n -- | -- s
  n>str loc{ ad n }
  n dup ufaculty dup 0
  do ad n str>vect zdup oddperm
     if zdrop then ad n nextp
  loop swap 1+ * negate >zst ;

utime 7 alt cardinality . utime d- d. 2520 -35424  ok

To filter out the odd permutations takes some time, so the alternating group of n!/2 even permutations runs in 35 ms.

What is left is to figure out how to generate general groups fast. And to write a manual!