Tuesday, May 30, 2017

Karatsuba multiplication

The time for direct multiplication is proportional to n² where n is the number of figures in the multiplicands x,y. When choosing a big base B=2ⁿ and and writing

xy=(x0+B*x1)(y0+B*y1)=x0y0 + (x0*y1+x1*y0)B + x1*y1*B²

there are four multiplications of smaller numbers, also here the calculation time is proportional to n². However


why it's enough to calculate three smaller multiplications 

x0*y0, x1*y1 and (x0+x1)(y0+y1). 

The multiplication with B are fast left shifting and if the shifting and the addition where cost free, the recursive Karatsuba multiplication would be very efficient

but unfortunately the extra math (and in Forth also some stack juggling) takes a lot of time and the method is efficient only for rather big numbers. For very big numbers, however, it's very efficient.

Here is the way I implemented it in ANS Forth:

: bcells* \ big m -- big*C^m
  cells top$ locals| n ad mb |
  ad ad mb + n move
  ad mb erase
  mb bvp @ +! ;
\ C is the number of digits in a cell

: bcells/ \ big m -- big/C^m
  cells top$ locals| n ad mb |
  ad mb + ad n move
  mb negate bvp @ +! ;

: bsplit \ w ad -- u v 
  dup nextfree < 
  if bvp @ dup @ vp+ bvp @ ! ! 
  else drop bzero
  then ;
\ A big number is split on the big stack at address ad

: btransmul \ x y -- x0 x1 y0 y1 m     B=2^bits 
  len1 len2 max cell + lcell 1+ rshift     \ m
  dup >r cells 
  >bx first over + bsplit 
  bx> first + bsplit r> ; 
\ x=x0+x1*B^m  y=y0+y1*B^m 

0x84 value karalim \ break point byte length for termination.

: b* \ x y -- xy
  len1 len2 max karalim < 
  if b* exit then
  btransmul >r                   \ x0 x1 y0 y1
  3 bpick 2 bpick recurse >bx    \ bx: x0*y0
  2 bpick 1 bpick recurse >bx    \ bx: x0*y0 x1*y1
  b+ >bx b+ bx>   recurse        \ (x0+x1)(y0+y1)
  bx b- by b- r@ bcells*         \ z1*C^m
  bx> r> 2* bcells* bx> b+ b+ <top ;
\ Karatsuba multiplication

Monday, May 29, 2017

How to use BigZ - part 3

The binomial coefficients for big integers

The number of possibilities to choose k objects from n objects soon get to big for a single cell number. The word bschoose gives a big integer result for single cell inputs.

: bschoose \ n k -- b

  bone 0
  ?do dup i - bs*
     i 1+ bs/mod drop
  loop drop ;


2000 500 bschoose cr b.
5648284895675941420424412140748481039502890353942825357221051675360331984776743417002364625179991976070068866284527555107208940603781511988000970130381311935878493235111594076219803768997324618773852975824828528735285833615310777764160933348372329757027402537319600321600269195597902747298520883357267710485334098751949232380773741897267988881873218260056305793069941805234442045890109611836653468404129012879905442075185208447514284775689056520318572740750419026192611832748925888424320  ok

This word produce big integers with single cell factors that can be analysed by the word 

sfacset \ b -- b' set

2000 1000 bschoose sfacset bdrop cr zet.

{2,5,7,11,13,17,19,23,37,41,43,53,59,67,73,79,101,103,113,127,131,149,151,167,173,179,181,211,251,257,263,269,271,277,281,283,337,347,349,353,359,367,373,379,383,389,397,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,1153,1163,1171,1181,1187,1193,1201,1213,1217,1223,1229,1231,1237,1249,1259,1277,1279,1283,1289,1291,1297,1301,1303,1307,1319,1321,1327,1361,1367,1373,1381,1399,1409,1423,1427,1429,1433,1439,1447,1451,1453,1459,1471,1481,1483,1487,1489,1493,1499,1511,1523,1531,1543,1549,1553,1559,1567,1571,1579,1583,1597,1601,1607,1609,1613,1619,1621,1627,1637,1657,1663,1667,1669,1693,1697,1699,1709,1721,1723,1733,1741,1747,1753,1759,1777,1783,1787,1789,1801,1811,1823,1831,1847,1861,1867,1871,1873,1877,1879,1889,1901,1907,1913,1931,1933,1949,1951,1973,1979,1987,1993,1997,1999} ok

The word bsetprod calculates the big product of the singles in set

: bsetprod \ set -- b
  bone                   \ big one
  foreach                \ make ready for do-loop
  ?do zst> bs* loop ;

and can be used to calculate the radical for big integers with single cell factors:

: bsradical \ b -- b'

  sfacset bdrop 
  bsetprod ;

50 25 bschoose bsradical cr b.

1504888171878  ok

Erdős squarefree conjecture (proved 1996) states that the central binomial coefficient (2n)Cn is not squarefree if n>4. The word sqrfacset calculates the set of all factors that occurs more than once: 

: sqrfacset \ b -- set

  bdup bsradical b/
  sfacset bdrop ;

20000 10000 bschoose sqrfacset cr zet.

{2,3,7,11,23,29,41,47,53,61,71,73,79,109,127,137,139} ok

The word

: maxel \ set -- n   non e
  zst> zst@ swap >zst zdrop ;

gives the maximal element in a set of integers.

: erdprime \ n -- p
  dup 2* swap bschoose
  sqrfacset maxel ;

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
  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 
  while zsplit n zst> dup >r - ?dup 
     if recurse
        if r> collincl n sumcorr xunion 
        else zst> drop r> drop 
     else { { r> } } xunion  
  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 
  while zsplit n zst> dup >r - ?dup 
     if recurse
        if r> collincl n sumcorr xunion 
        else zst> drop r> drop 
     else { { r> } } xunion  
  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 
  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