What words should be included in a minimal number theory package?

[lehs]

While developing my BigZ (https://forthmath.blogspot.se/ https://github.com/Lehs/BigZ) I often create integer sequences which I try to find on https://oeis.org/. Now and then I propose a new sequence and have to design a program that generate the sequence.

However, I always use Forth which probably is very unknown by OEIS users and lacks the most math words that has to be used, as SQRTF ( u -- floor ) and PRIME ( u -- flag ).

(ANS-Forth)
\ prime \ n -- flag
: Gauss_twins \ n -- a(n)
dup * locals| square | 0
begin 1+ dup dup * square + prime
over 2 + dup * square + prime and
until ;

So far I've solved the problem by state the stack diagram for words not included in the ANS standard, but it would be nice to refer to a minimal standard package and what about:

sqrtf \ u -- floor
umod \ u q -- u(mod q)
u*mod \ u v q -- uv(mod q)
u**mod \ u v q -- u**v(mod q)
ugcd \ u v -- w
prime \ u -- flag
nextprime \ u -- p

[krishna...@ccreweb.org]

On Sunday, March 26, 2017 at 5:13:43 AM UTC-5, lehs wrote:
> While developing my BigZ (https://forthmath.blogspot.se/ https://github.com/Lehs/BigZ) I often create integer sequences which I try to find on https://oeis.org/. Now and then I propose a new sequence and have to design a program that generate the sequence.
>
> However, I always use Forth which probably is very unknown by OEIS users and lacks the most math words that has to be used, as SQRTF ( u -- floor ) and PRIME ( u -- flag ).

...

Why not use the name PRIME? or ?PRIME since it's a query?

Krishna

[lehs]

PRIME? and ISPRIME are perhaps better than PRIME, I'm just used to the latter. I think that ISPRIME is the most common name.

Forth is not oriented for math, but with a few extensions it works perfect. The main problem is that words as * may or may not result in overflow. There is a lot of nice words in ANS-Forth to use to override that problem: UM* UM/MOD M*/ etc, but it's not always obvious how to use those words to solve the problem. One has to make a lot of reinvention all the time.

With words like UMOD, U*MOD U**MOD a lot of routines are within reach for defining more sophisticated primarily tests and factoring (Pollard rho) which make it easy to define a lot of usual number theoretical functions as totients, big omega and Möbius.

[NN]

Hi Lehs,

Most number theory is beyond me but I wondered if you had looked at

http://mathworld.wolfram.com/NumberTheory.html.

for inspiration ?

NN

[lehs]

Thanks for the link! Since there is a lot of number theory beyond my head too, there are functions I will never implement. Most interesting are integer valued functions which also are possible to realize with varied efficiency. Standard Forth isn't perfect for integer functions but can be extended to be good enough for my purpose. Single cell arithmetic is fast but limited because of overflow problems. Words must be carefully defined and should be collected in a libraries ready to be used. If there was an interest those libraries would grow worldwide, but Forth isn't the first choice for number theorists.

My minimal package will be a minimalistic extension of Forth, with words that probably never will be standard Forth, extentions that makes it more easy to define other words.

[Elizabeth D. Rather]

On 3/26/17 2:03 PM, lehs wrote:
> Thanks for the link! Since there is a lot of number theory beyond my head too, there are functions I will never implement. Most interesting are integer valued functions which also are possible to realize with varied efficiency. Standard Forth isn't perfect for integer functions but can be extended to be good enough for my purpose. Single cell arithmetic is fast but limited because of overflow problems. Words must be carefully defined and should be collected in a libraries ready to be used. If there was an interest those libraries would grow worldwide, but Forth isn't the first choice for number theorists.
>
> My minimal package will be a minimalistic extension of Forth, with words that probably never will be standard Forth, extentions that makes it more easy to define other words.
>

Chuck intended Forth to be wonderful for integer math; it was designed
that way in days before hardware FP was common. He wrote extensive math
functions with integer arithmetic, including 16- and 32-bit fixed-point
fraction arithmetic, trig, least-squares fitting algorithms, and FFTs.

Cheers,
Elizabeth

--
Elizabeth R. Conklin
FORTH, Inc.
6080 Center Drive, Suite 600
Los Angeles, CA 90045
USA

[hughag...@gmail.com]

On Sunday, March 26, 2017 at 5:16:35 PM UTC-7, Elizabeth D. Rather wrote:
> Chuck intended Forth to be wonderful for integer math; it was designed
> that way in days before hardware FP was common. He wrote extensive math
> functions with integer arithmetic, including 16- and 32-bit fixed-point
> fraction arithmetic, trig, least-squares fitting algorithms, and FFTs.

You are a professional name-dropper --- you don't have any accomplishments of your own --- you have made an entire career out of bragging about how you knew Charles Moore in the 1970s.

I'm glad that Charles Moore wrote extensive math functions with integer arithmetic before you kicked him out of Forth Inc. in 1982 --- good for him! --- but I don't think this gives you the right to say that I'm not a Forth programmer.

> Cheers,
> Elizabeth

Well, it was after 5:00 pm when you posted this --- so you have already begun your drinking --- or did you get a head-start on that at 5:00 am?

Cheers!

[lehs]

Forth has all what is needed except libraries for functions like:

: isprime \ u -- flag
dup 5 u< if 2 4 within exit then
dup 1 and 0= if 0= exit then
true over sqrtf 3 umax 1+ 3
do over i umod 0=
if 0= leave then 2
+loop nip ;

: ugcd \ u1 u2 -- u \ Greatest common divisor
2dup u< if swap then \ smallest of a b on top of stack (tos)
?dup 0= if exit then \ return second value if tos is zero
begin tuck \ y x y first time b a b
0 swap um/mod \ y x 0 y --> y r q
drop dup 0= \ y r [r=0]
until drop ; \ y

: u*mod \ u1 u2 u3 -- u1*u2(mod u3)
>r r@ umod swap r@ umod um*
r> um/mod drop ;

: u**mod \ b a m -- x
locals| m a b |
1 \ preparation for repeated multiplication
a log~ 0 \ n 0 n is the number of binary digits in a
?do a 1 i lshift and \ flag the i-th digit of a is 0/1
if b m u*mod \ multiply if flag
then b dup m \ square b for each i: b b^2 b^4 b^8 ...
u*mod to b \ new squared number to b
loop ; \ result of the repeated multiplication

With such libraries Forth would be perfect for integer math. But as it is every user of standard Forth has to reinvent all those words that is necessary building stones for all integer math that might interest mathematicians.

[Elizabeth D. Rather]

I'm not a mathematician myself (databases and UI were my thing when I
was programming), but I don't recall Chuck or the math guys I worked
with ever using any of these. What do you do with them?

[Paul Rubin]

"Elizabeth D. Rather" <era...@forth.com> writes:
> I don't recall Chuck or the math guys I worked with ever using any of
> these. What do you do with them?

They show up in number-theoretic and combinatorial math, and these days
in cryptography. They wouldn't show up in the traditional physics or
engineering type math that I think was more Chuck's area.

[lehs]

I guess the single cell arithmetic is mostly for efficient testing of conjectures and similar, but the big integer variants

bgcd \ b1 b2 -- b3
b*mod \ b1 b2 b3 -- b1*b2(mod b3)
b**mod \ b1 b2 b3 -- b1**b2(mod b3)

seems to be useful in certain applications - if they are fast enough.

Thanks to all the double and mixed precision words in ANS-Forth it's possible to develop fast number theoretical words in Forth, but it's inconvenient to do from scratch and words like this ought to be collected in certain libraries, that doesn't have to contain hundreds of words but just the most basic words in the same style as ANS-Forth extensions.

[Elizabeth D. Rather]

There is an existing math library that you're welcome to contribute to.
I'm not sure who's administering it these days, though...

[minf...@arcor.de]

Abd there is the Forth Scientific Library #47
https://www.taygeta.com/fsl/library/big.fth

Depending on your Forth you might better link with the GNU BigNum Library.

Multiplication and particularly division/modulus are tricky and there are many non-trivial optimized algorithms. Karatsuba's one of the more simple ones.

But feel free to reinvent the wheel as long as it is a fun project. ;-)

[lehs]

I already had that fun: https://github.com/Lehs/BIG-INTEGER-ANS-FORTH and I do think that b/mod and b**mod works ok. The advantage of my system is that it use big integer stacks which make programming easy. The disadvantage of my system is that big integer stacks may not be optimal efficient. And that my system only is for unsigned integers.

Also, I have difficulties loading most of the code at https://www.taygeta.com/fsl/library/big.fth

What I'm doing right now is an experimental minimalist library for single cell arithmetic extensions, perhaps including the most well known functions as RADICAL, TOTIENTS, MÖBIUS, BIGOMEGA, SMALLOMEGA, LEGENDRE etc.

I might organize an experimental library on github, with the ambition being possible to load an all ANS-Forth systems (32/64 bits).

[krishna...@ccreweb.org]

On Sunday, March 26, 2017 at 8:13:17 PM UTC-5, lehs wrote:
> Forth has all what is needed except libraries for functions like:
>
> : isprime \ u -- flag

...

A good example to demonstrate your library is to write a sample program to illustrate factoring a number into its composite primes. An efficient factoring word would actually be a useful addition to a scientific/numeric library. BTW, you didn't include NEXT-PRIME in your library?

Krishna

[lehs]

I've put it in
https://github.com/Lehs/ANS-Forth-libraries/blob/master/basicsinglext.4th
and it contain NEXTPRIME and PRIMEFACTORS plus some number theoretical functions as RADICAL TOTIENTS MOBIUS BIGOMEGA SMALLOMEGA LEGENDRE TAU and SIGMA.

I think it will work for all true ANS Forth systems. I use LOCALS| since that is the standard.

From this point any barefoot mathematician forther could start experimenting in number theory and extend on their own without having to do the cumbersome part direct from the standard.

[Albert van der Horst]

In article <162e9326-e7e5-4f7b...@googlegroups.com>,
lehs <skydda...@gmail.com> wrote:
>While developing my BigZ (https://forthmath.blogspot.se/ https://github.com=
>/Lehs/BigZ) I often create integer sequences which I try to find on https:/=
>/oeis.org/. Now and then I propose a new sequence and have to design a prog=
>ram that generate the sequence.=20
>
>However, I always use Forth which probably is very unknown by OEIS users an=
>d lacks the most math words that has to be used, as SQRTF ( u -- floor ) an=

>d PRIME ( u -- flag ).
>
>(ANS-Forth)
>\ prime \ n -- flag
>: Gauss_twins \ n -- a(n)
> dup * locals| square | 0
> begin 1+ dup dup * square + prime
> over 2 + dup * square + prime and
> until ;
>

>So far I've solved the problem by state the stack diagram for words not inc=
>luded in the ANS standard, but it would be nice to refer to a minimal stand=

>ard package and what about:
>

>sqrtf \ u -- floor=20

>umod \ u q -- u(mod q)
>u*mod \ u v q -- uv(mod q)
>u**mod \ u v q -- u**v(mod q)
>ugcd \ u v -- w
>prime \ u -- flag
>nextprime \ u -- p

Get rid of the u. They are numbers.
When was the last time that you had a number between
18446744073709551615 and half of that?

I have
FACTOR \ For n and p return smallest factor larger than p.
PRIMES? \ For p return flag whether it is prime
GCD \ Your GCD
SQRT \ Your sqrtf
XGCD \ Extended gcd, for A B return C and GCD such that A*C + B*X is GCD
CHS \ m n -- m "over" n (choose function)
PYR \ pyramidal function
TRI \ triangular function
MU \ moebius function
PHI \ totient function
FAC \ faculty (!)

Also a modulo package ( all arithmetic operators modulo a fixed
number, a global that you can fill in.)
Code can be stolen form the lab file of one of the ciforth's.

Groetjes Albert
--
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

[Albert van der Horst]

In article <M6GdnfHLpdzF60XF...@supernews.com>,

This screen has modulo arithmetic for all the usuals operations.
It is one of the most heavily used. E.g. the latest euler problem
ask for the solution to be calculated modulo 1,000,000,007
(which is a prime).
So a solution starts with
1,000,000,007 _m !
and then I can use `` +m '' instead of `` + '' etc.

0 ( +m -m *m /m %:m **m x^x ) \ AvdH B6mar23
1 \ Can be INLINE d. Suffix n : normalized number.
2 VARIABLE _m ( Modulo number) WANT XGCD
3 : _norm_-m DUP 0< _m @ AND + ; ( x -- xn ) \ -m<xn<+m
4 : +m + _m @ - _norm_-m ; ( an bn -- sumn )
5 : -m - _norm_-m ; ( an bn -- diffn)
6 : *m M* _m @ SM/REM DROP ; ( an bn -- prodn)
7 : /m _m @ XGCD DROP _norm_-m *m ; ( a b -- quotn)
8 : %:m S>D _m @ SM/REM DROP _norm_-m ; ( a -- an)
9
10 \ Both steps: For A B and C: return A B en C. Invariant A*B^C.
11 : _reduce_1- 1- >R >R R@ *m R> R> ;
12 : _reduce_2/ 2/ >R DUP *m R> ;
13 : **m 1 ROT ROT BEGIN DUP 1 AND IF _reduce_1- THEN
14 _reduce_2/ DUP 0= UNTIL 2DROP ; ( a b -- apowbn )
15 : x^x _m @ >R _m ! **m R> _m ! ; ( a b m -- a^b mod m ) OK


P.S. Division only works if _m contains a prime.

>
>Cheers,
>Elizabeth

[lehs]

No, I focus on unsigned integers.

\ Unsigned integer number theory basicsinglext.4th
\
\ isprime \ u -- flag
\ nextprime \ u -- p
\ primefactors \ u -- q1 q2 ... qn n
\ radical \ n -- r
\ totients \ n -- t
\ mobius \ n -- m
\ bigomega \ n -- b
\ smallomega \ n -- s
\ legendre \ a p -- i
\ tau \ u -- n The number of divisors of u
\ sigma \ u -- s The sum of all divisors of u
\ choose \ n k -- nCk
\ sqrtf \ u -- floor
\ log~ \ n -- 1+log n
\ umin \ u1 u2 -- u
\ umax \ u1 u2 -- u
\ umod \ u q -- u(mod q)
\ ugcd \ u1 u2 -- u
\ u*mod \ u1 u2 u3 -- u1*u2(mod u3)
\ u**mod \ b a m -- b^a(mod m)

https://github.com/Lehs/ANS-Forth-libraries

[hughag...@gmail.com]

On Monday, March 27, 2017 at 11:33:50 PM UTC-7, lehs wrote:
> Den tisdag 28 mars 2017 kl. 01:44:40 UTC+2 skrev Albert van der Horst:
> > Get rid of the u. They are numbers.
> > When was the last time that you had a number between
> > 18446744073709551615 and half of that?
>

> No, I focus on unsigned integers.
>
> \ Unsigned integer number theory basicsinglext.4th
> \
> \ isprime \ u -- flag
> \ nextprime \ u -- p
> \ primefactors \ u -- q1 q2 ... qn n
> \ radical \ n -- r
> \ totients \ n -- t
> \ mobius \ n -- m
> \ bigomega \ n -- b
> \ smallomega \ n -- s
> \ legendre \ a p -- i
> \ tau \ u -- n The number of divisors of u
> \ sigma \ u -- s The sum of all divisors of u
> \ choose \ n k -- nCk
> \ sqrtf \ u -- floor
> \ log~ \ n -- 1+log n
> \ umin \ u1 u2 -- u
> \ umax \ u1 u2 -- u
> \ umod \ u q -- u(mod q)
> \ ugcd \ u1 u2 -- u
> \ u*mod \ u1 u2 u3 -- u1*u2(mod u3)
> \ u**mod \ b a m -- b^a(mod m)
>
> https://github.com/Lehs/ANS-Forth-libraries

I have double-precision GCD etc. in the novice-package.

I have CF.4TH in the novice-package that calculates continued fractions. This was written by Nathaniel Grosman in a "Forth Dimensions" article for 16-bit Forth-83 --- I upgraded it slightly to work under 32-bit ANS-Forth.

Porting this to ANS-Forth was one of my very first ANS-Forth efforts. Elizabeth Rather told me that I can never be a Forth programmer because I complained about the lack of double-precision division in ANS-Forth.

Here is the first half of the CF.4TH code:
-------------------------------------------------------------------------

\ ******
\ ****** This is the first section, which provides UD/MOD.
\ ******

: LPswap ( d1 d2 d3 d4 -- d3 d4 d1 d2 )
>r >r 2swap >r >r 2swap r> r> r> r> 2swap
>r >r 2swap r> r> ;

: LPdup ( d1 d2 -- d1 d2 d1 d2 )
2dup >r >r 2over r> r> ;

: LPover ( d1 d2 d3 d4 -- d1 d2 d3 d4 d1 d2 )
>r >r >r >r LPdup r> r> r> r> LPswap ;

: LProt ( d1 d2 d3 d4 d5 d6 -- d3 d4 d5 d6 d1 d2 )
>r >r >r >r LPswap r> r> r> r> LPswap ;

: LPdrop ( d1 d2 -- )
drop drop drop drop ;

: um/ ( ud un -- un ) \ divide UD by UN and drop remainder
um/mod swap drop ;

: u*/ ( Umultiplier Udividend Udivisor -- Uquotient )
>r um* r> um/mod nip ;

: t* ( ud un -- ut )
dup rot um* >r >r
um*
0 r> r> d+ ;

: t/ ( ut un -- ud )
>r r@ um/mod swap
rot 0 r@ um/mod swap
rot r> um/mod swap drop
0 2swap swap d+ ;

: ut*/ ( ud un un -- ud ) \ this was called U*/ in Grossman's article, but I'm already using that name
>r t* r> t/ ;

: d- ( d1 d2 -- d1-d2 )
dnegate d+ ;

: narrow-ud/mod ( UDividend UNdivisor -- UDremainder UDquotient )
drop >r 2dup r@ \ shuck high cell of divisor
1 swap ut*/ \ -- UDquotient
2swap 2over r> 1 ut*/ d- 2swap ; \ -- UDrem UDquot

variable numH
variable denH
variable denL
variable denscale
2variable num
2variable den
0 1 2constant superbase \ this was 65536. in Grossman's article

: us>d \ convert unsigned single-cell to double-cell
0 ;

: ud/mod-tuck ( ud ud -- ) \ save parts of num and den
2dup den 2! denH ! denL !
2dup num 2! numH ! drop ;

: ud/mod-denscale ( -- un ) \ for scaling-up den
superbase denH @ 1+ um/
denscale ! ;

: scale-den ( -- ) \ multiply denominator by scale factor
den 2@ denscale @ 1 ut*/
denH ! denL ! ;

: wide-quot ( -- ud ) \ if divisor needs more than one cell
num 2@ numH @ us>d
denL @ denH @ ut*/ d-
denscale @ denH @ ut*/ swap drop ;

: ?narrow-divisor ( d -- ? ) \ is divisor < superbase?
dup 0= ;

: wide-rem ( -- ud ) \ remainder in wide division
dup num 2@ rot den 2@
rot 1 ut*/ d- rot us>d ;

: wide-ud/mod ( UDdividend UDdivisor -- UDremainder UDquotient )
ud/mod-tuck
ud/mod-denscale
scale-den
wide-quot
wide-rem ;

: ud/mod ( UDdividend UDdivisor -- UDremainder UDquotient )
?narrow-divisor if narrow-ud/mod
else wide-ud/mod then ;

: udmod ( UDdividend UDdivisor -- UDremainder )
ud/mod 2drop ;

: ud/ ( UDdividend UDdivisor -- UDquotient )
ud/mod 2swap 2drop ;

[lehs]

Impressive with double precision division words, but it doesn't solve the problems with single precision number theoretical words.

If someone loaded https://github.com/Lehs/ANS-Forth-libraries there was an error with TOTIENTS TAU SIGMA SIGMA** which is corrected now.

It really doesn't matter where libraries are stored as long as they are described appropriate enough to be possible to google when there is a certain need.

[krishna...@ccreweb.org]

Thanks for the nice code. I'll try it out.

Krishna

[lehs]

Thank YOU! I've changed the address: https://github.com/Lehs/ANS-Forth-libraries

[krishna...@ccreweb.org]

On Saturday, April 1, 2017 at 9:41:51 AM UTC-5, lehs wrote:
...

>
> Thank YOU! I've changed the address: https://github.com/Lehs/ANS-Forth-libraries

I've tried it now. Although I don't know what all of the words do, here are a few comments:

1) The word NEXTPRIME returns the same number if the argument is prime, e.g.

19 nextprime .
19 ok

2) Consider using the common usage word(s), [DEFINED] and/or [UNDEFINED] , instead of devising your own.

3) Maybe an implementation of a prime counting function would be useful to include, for small upper bound, e.g.

: pcf ( u -- ucount )
0 swap 0 do i isprime if 1+ then loop ;

That's quick and dirty, and not efficient at all once "u" is on the order of 1 million.

4) Some test code, perhaps in a separate file, would help the user understand if the words are working as intended.

Thanks very much for this!

: show-primes ( u -- )
0 do i isprime if i 6 .r then i 1+ 10 mod 0= if cr then loop ;

100 show-primes
2 3 5 7
11 13 17 19
23 29
31 37
41 43 47
53 59
61 67
71 73 79
83 89
97
ok

Krishna

[lehs]

I decided that nextprime ( n -- p ) should give n if n was prime. I don't know if it's unconventional?

If I use [defined] etc I still have to include the code since those are not ANS standard.

There is a fast prime counting function PI in primecounting.4th which use squeezed tables described in http://forthmath.blogspot.se/2015/12/prime-tables-and-prime-counting-function.html and even a function PNR@ ( n -- Pn ) that gives the nth prime.

I shall include examples.

Lars-Erik Svahn
https://forthmath.blogspot.se
https://github.com/Lehs/ANS-Forth-libraries

[Julian Fondren]

On Monday, April 17, 2017 at 6:34:51 PM UTC-5, lehs wrote:
> If I use [defined] etc I still have to include the code since those are not ANS standard.

It's not Forth-94, but it's Forth-2012: http://forth-standard.org/standard/tools/

[lehs]

Thanks! This seems to be the first place to look after names of non ANS words. Nice that even the code is given as definition since formal human language isn't quit as unambiguos as code.

[lehs]

Den tisdag 18 april 2017 kl. 01:08:27 UTC+2 skrev krishna...@ccreweb.org:

Okay, I've looked around and it seems to be a convention that

nextprime ( n -- p )

should give the next prime greater than n.

[m...@iae.nl]

That is a sweeping statement. Good Forths have all these simple
words and many more, some (i.e. iForth) even with arbitrary
precision. And there *are* (Forth) libraries and publications
of this stuff.

So, in the future, please write:

"But as it is, I THINK THAT every user of standard Forth has to

reinvent all those words that is necessary building stones for
all integer math that might interest mathematicians.

I MIGHT BE WRONG."

-marcel

[lehs]

Yes, as it seems there are fantastic libraries hidden on the Internet which some of them seems to be difficult to use.

After all, 'easy to use' is an advantages among with 'extremely efficient'.

It's a dull thing to do, looking around for libraries that, if you found them, are hard to use.

[krishna...@ccreweb.org]

On Monday, April 17, 2017 at 6:34:51 PM UTC-5, lehs wrote:

> Den tisdag 18 april 2017 kl. 01:08:27 UTC+2 skrev

> ...

> I decided that nextprime ( n -- p ) should give n if n was prime. I don't know if it's unconventional?
>

Ok. It was unexpected but the behavior is fine if it was by design -- it's possibly more useful the way you have define it.


> ... There is a fast prime counting function PI in primecounting.4th which use squeezed tables described in http://forthmath.blogspot.se/2015/12/prime-tables-and-prime-counting-function.html and even a function PNR@ ( n -- Pn ) that gives the nth prime.

Thanks for the link.

Krishna