## Saturday, January 18, 2020

The Pollard rho algorithm factorize a composite number n in a time proportional to √p, where p is the smallest prime that divides n. In worst case p≈√n so the algorithm factorize proportional to ∜n. Compare with trial-and-error that factorize proportional to √n.

The algorithm is easy to implement and is built upon an the greatest common divisor gcd (which can be calculated by the algorithm of Euclide) and a polynomial P(x) of degree >1. Most common polynomial to use is P(x)=x²+1. Let Pn(x)=P(x) mod n, that is the rest when P(x) is divided by n. Pn is a function Pn:ℕ→ℤn, where ℤn={0,1,2,...,n-1}. The function Pn will act like a simple pseudo random generator in the algorithm.

Define Xi+1=Pn(Xi). Since ℤn is a finite set there are smallest i,j such that Xi=Xi+j, when the sequence start to repeat itself. Mostly the sequence Xi isn't cyclic from the start Xo, but later on at some Xi. There is at trick to find out if Xi is in the loop or not. It's like sending away a turtle and a hare on the same track at the same time. When the hare again is comming besides the turtle they must both be in the loop.

Therefore, define a second sequence Yi+1=Pn(Pn(Yi)), where Xo=Yo. When Yi=Xi, i>0, Xi and Yi are in the cycle. But for Pollard rho the real important sequences are the uncalculated sequences Vi+1=Pm(Vi) and Wi+1=Pm(Pm(Wi)) where Vo=Wo=Xo and m|n (which aren't calculated since we don't know any m|n yet). When Wi=Vi, i>0, then m|Xi-Yi and gcd(Xi-Yi,n)=m. If m=1 the result is neglected and the process go on with Xi+1 and Yi+1. If m=n the start values Xo=Yo is said to fail and may be increamented by 1 for an other try. Now, gcd(Xi-Yi,n)=n when Xi and Yi become equal before any of the uncalculated sequences Vi and Wi become equal for some non trivial divisor m>1 of n. There could still be non trivial divisors to obtain for gcd(Xj-Yj,n) for j>i when Xi=Yi, but there is no guarranty and Xi=Yi is a natural terminal case.

An example with P(x)=x²+1 that is factorized while continuing after failure for Xo=2 is n=4294952621. An example that can't be factorized after failure with Xo=2 but with Xo=3 is n=4294939069.

Increasing Xo after termination when gcd(Xi-Yi,n)=n will always find a factor of n, because if m is a non trivial factor of n, then Xo=m immediately gives m as a factor. That is, for
P(x)=x²+1.

## Friday, September 1, 2017

### BigZ - new top level instructions

I've reorganized and extended some top level words:

create-set \ m n xt -- set
filter-set \ set1 xt -- set2
build-set \ m n xt -- set
transform-set \ set1 xt -- set2

This words should be used with the word :| to define and form sets.

1 10 :| b a | a ^2 b ^2 + isprime ; create-set zdup cr zet.
{(1,1),(2,1),(4,1),(6,1),(1,2),(3,2),(5,2),(7,2),(2,3),(8,3),(1,4),(5,4),(9,4),(2,5),(4,5),(6,5),(8,5),(1,6),(5,6),(2,7),(8,7),(3,8),(5,8),(7,8),(4,9)} ok

:| b a | a ^2 b ^2 + ; transform-set zdup cr zet. {2,5,13,17,29,37,41,53,61,73,89,97,113} ok

:| n | n 4 mod 3 = ; filter-set zet. 0 ok

The word :| define a nameless word and count the number of parameters.

: bl# \ ad n -- m    count the number of spaces in the string
over + swap 0 -rot
do i c@ bl = -
loop ;

variable loc#    \ the number of parameters
variable sta#    \ the number of outputs on the stack

: locals# \ --
>in @ >r
[char] | parse bl# loc# !
r> >in !
1 sta# ! ; immediate

: :| \ --
:noname
postpone locals#
postpone locals| ;

The nameless words, represented by xt on stack, could be of two types:

1. Taking parameters and leaving a flag.
2. Taking parameters and leaving a non negative integer.

In the first case a set with the "dimension" stated by the parameters is the result. This works with CREATE-SET and FILTER-SET.

In the second case a set of values is the result. This works with BUILD-SET and TRANSFORM-SET.

1 10 :| x | x ^2 ; build-set cr zet.
{1,4,9,16,25,36,49,64,81} ok

Normally these nameless words leave one integer on the stack, but for TRANSFORM-SET there is an option (using 2; instead of ;) when the word leaves two non negative integers on the stack

{ 0 10 | all }

creates the set {0,1,2,3,4,5,6,7,8,9} which can be transformed to a two dimensional set

:| n | n n ^2 2; transform-set cr zet.
{(0,0),(1,1),(2,4),(3,9),(4,16),(5,25),(6,36),(7,49),(8,64),(9,81)} ok

The purpose with this words is to be able to quickly inspect conjectures and Diophantine equations.

a and b are coprime if there exists values x and y such that

ax+by=1

Since BigZ just deals with sets of non negative numbers we can search x and y in the equation

ax-by=1

1 10 :| b a y x | a x * b y * - 1 = ; create-set cr zet.
{(2,1,1,1),(3,2,1,1),(4,3,1,1),(5,4,1,1),(6,5,1,1),(7,6,1,1),(8,7,1,1),(9,8,1,1),(1,1,2,1),(2,3,2,1),(3,5,2,1),(4,7,2,1),(5,9,2,1),(1,2,3,1),(2,5,3,1),(3,8,3,1),(1,3,4,1),(2,7,4,1),(1,4,5,1),(2,9,5,1),(1,5,6,1),(1,6,7,1),(1,7,8,1),(1,8,9,1),(3,1,1,2),(5,2,1,2),(7,3,1,2),(9,4,1,2),(1,1,3,2),(3,4,3,2),(5,7,3,2),(1,2,5,2),(3,7,5,2),(1,3,7,2),(1,4,9,2),(4,1,1,3),(7,2,1,3),(2,1,2,3),(5,3,2,3),(8,5,2,3),(1,1,4,3),(4,5,4,3),(7,9,4,3),(2,3,5,3),(5,8,5,3),(1,2,7,3),(4,9,7,3),(2,5,8,3),(5,1,1,4),(9,2,1,4),(3,2,3,4),(7,5,3,4),(1,1,5,4),(5,6,5,4),(3,5,7,4),(1,2,9,4),(6,1,1,5),(3,1,2,5),(8,3,2,5),(2,1,3,5),(7,4,3,5),(4,3,4,5),(9,7,4,5),(1,1,6,5),(6,7,6,5),(3,4,7,5),(2,3,8,5),(4,7,9,5),(7,1,1,6),(5,4,5,6),(1,1,7,6),(7,8,7,6),(8,1,1,7),(4,1,2,7),(5,2,3,7),(2,1,4,7),(9,5,4,7),(3,2,5,7),(6,5,6,7),(1,1,8,7),(8,9,8,7),(4,5,9,7),(9,1,1,8),(3,1,3,8),(5,3,5,8),(7,6,7,8),(1,1,9,8),(5,1,2,9),(7,3,4,9),(2,1,5,9),(4,3,7,9),(8,7,8,9)} ok

The vectors in the set are of the form (x,y,a,b), that is the opposite of its appearance as parameters. (Which is logical in Forth with postfix notations and values on a stack, but still a bit awkward).

The idea is a one row code and for that purpose there is a need of short words:

: isp  isprime ;         \ n -- flag
: isq  sqr ;             \ is perfect square: n -- flag
: isqf  sqrfree ;        \ is square free: n -- flag
: isem  bigomega 2 = ;   \ is semi prime: n -- flag
: ispp  smallomega 1 = ; \ is prime power: n -- flag

: 2sqs  2sqsum ;         \ square sum: a b -- sum
: 3sqs  3sqsum ;         \ square sum: a b c -- sum
: 4sqs  4sqsum ;         \ square sum: a b c d -- sum

: cop  coprime ;         \ are coprime: m n -- flag
: div  swap mod 0= ;     \ divides: m n -- flag

: <<  \ i j k -- flag  i<j<k
over > -rot < and ;

: <<=  \ i j k -- flag  i<=j<=k
over >= -rot <= and ;

: z. zet. ;

: fi postpone else postpone false postpone then ; immediate
\ a short version of ELSE FALSE THEN.

When inspecting a Diophantine equation there might be symmetries and trivial cases to weed out.

1 100 :| c b a | a b c << a b cop and if a b 2sqs c ^2 = else false then ; create-set

can be shortened to

1 100 :| c b a | a b c << a b cop and if a b 2sqs c ^2 = fi ; create-set

cr zet.
{(3,4,5),(5,12,13),(8,15,17),(7,24,25),(20,21,29),(12,35,37),(9,40,41),(28,45,53),(11,60,61),(33,56,65),(16,63,65),(48,55,73),(36,77,85),(13,84,85),(39,80,89),(65,72,97)} ok

Experimenting with sets gives a great opportunity to find and test conjectures:

1 50 :| b a | a b 2sqsum isprime ; create-set  ok
:| b a | a b + ; transform-set cr z.
{2,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79,83,85,87,89,91,93,95} ok

This suggests that all odd numbers larger than 1 can be written as a sum a+b where a²+b² is prime.

Define

: sqeq \ m n a b -- ma²+nb²
dup * rot * -rot
dup * * + ;

and the set {a+b<50|a,b>0 & ma²+nb² is prime} for some different values of m and n

5 value m
7 value n

1 50 :| b a | m n a b sqeq isprime ; create-set
:| b a | a b + ; transform-set cr zet.
{5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,89,91} ok

3 to m
4 to n

1 50 :| b a | m n a b sqeq isprime ; create-set
:| b a | a b + ; transform-set cr zet.
{2,3,4,5,6,8,9,10,11,12,13,15,16,17,18,19,20,22,23,24,25,26,27,29,30,31,32,33,34,36,37,38,39,40,41,43,44,45,46,47,48,50,51,52,53,54,55,57,58,59,60,61,62,64,65,66,67,68,69,71,72,73,74,75,76,78,79,80,81,82,83,86,87,88,89,93,94,95,96} ok

This will show a pattern which suggest that
{a+b|a,b>0 & ma²+nb² is prime}={k>1|gcd(k,m+n)=1}

However, further tests will show that some changes are needed:

Routines to check the conjecture and print the table:

: prime_partition \ m n k -- flag
false locals| flag k |
k 1
?do 2dup k i - i sqeq isprime
if true to flag leave then
loop 2drop flag ;

\ 100000 value nx     Gives the same table as
1000 value nx

:
maximal_exception \ m n -- mx
2dup + 1 locals| mx m+n n m |
nx 3
do m+n i coprime
if m n i
prime_partition 0=
if i to mx then
then
loop mx ;

: table \ m2 n2 --
locals| n2 | cr
2 spaces dup 1 do i 3 .r loop 1
do n2 1 cr i 2 .r
do j i coprime
if j i
maximal_exception
else 0
then 3 .r
loop
loop ;

The source code for BigZ can be loaded here:

https://github.com/Lehs/BigZ/blob/master/bigzet.txt

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

(x0*y1+x1*y0)=(x0+x1)(y0+y1)-x0*y0-x1*y1

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 |
mb bvp @ +! ;
\ C is the number of digits in a cell

: bcells/ \ big m -- big/C^m
cells top\$ locals| n ad mb |
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 ;

ok

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

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

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

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

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

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)