Best Euler #13 solution?

I don't know that this is any good, but it's something I saved a while
back. It has some unused hacks like finding the largest intermediate
result, which is around 2**34. Because of that, it requires 64 bit
cells. It runs in about 0.8 sec on my laptop with gforth-fast.

: cs clearstack ;

: cnext ( n -- n ) dup 1 and if 3 * 1+ else 1 rshift then ;

100000 constant msize

CREATE memo msize cells allot
DOES> ( n a ) swap cells + ;

: mclear ( -- ) \ zero out the table
0 memo msize cells erase ;
mclear

: mlook ( n -- n ) dup msize >= if drop 0 else memo @ then ;
: mset ( v k -- ) dup msize >= if 2drop else memo ! then ;

: mlook1 { n -- n } n msize >= if 0 else n memo @ then ;
: mset1 { v k -- } k msize < if v k memo ! then ;

: mlook2 ( n -- n flag ) dup msize >= IF false ELSE memo @ true THEN ;

1 1 mset

: clen3 ( n -- n ) >r
r@ 1 = IF
r> drop 0
ELSE
r@ mlook dup
0= IF
drop r@ cnext RECURSE 1+
dup r@ mset
THEN
r> drop
THEN ;

: clen2 { n -- n } n mlook2 0= ~~ IF
drop
n cnext RECURSE 1+ { cn } ~~
cn n mset
cn
THEN ;

: clen1 { n -- n } n mlook { cn } cn 0= IF
n cnext RECURSE 1+ to cn
cn n mset
THEN
cn ;

: clen1a ( n -- n ) dup ( n n ) mlook ( n cn ) dup 0= IF
drop ( n ) dup cnext RECURSE 1+ ( n cn )
dup -rot ( cn cn n ) mset
ELSE ( n cn ) nip
THEN ;

: euler13 ( -- n )
0
1e6 f>s 1 DO i clen1 max LOOP ;

euler13 . cr

Ron AARON

2024-05-02 11:41:50 UTC

Post by minforth
Inspired by the Euler problem #13, I dug out one of my old
Forth systems after a long time. See below for the solution.
Can you do the same in your Forth system?

Similar idea using 8th:

: doit
'# parse "\n" s:/
( >n n:+ ) 0 a:reduce

Post by minforth
s 10 s:lsub . cr ;

doit
37107287533902102798797998220837590246510135740250
46376937677490009712648124896970078050417018260538
74324986199524741059474233309513058123726617309629
91942213363574161572522430563301811072406154908250
23067588207539346171171980310421047513778063246676
89261670696623633820136378418383684178734361726757
28112879812849979408065481931592621691275889832738
44274228917432520321923589422876796487670272189318
47451445736001306439091167216856844588711603153276
70386486105843025439939619828917593665686757934951
62176457141856560629502157223196586755079324193331
64906352462741904929101432445813822663347944758178
92575867718337217661963751590579239728245598838407
58203565325359399008402633568948830189458628227828
80181199384826282014278194139940567587151170094390
35398664372827112653829987240784473053190104293586
86515506006295864861532075273371959191420517255829
71693888707715466499115593487603532921714970056938
54370070576826684624621495650076471787294438377604
53282654108756828443191190634694037855217779295145
36123272525000296071075082563815656710885258350721
45876576172410976447339110607218265236877223636045
17423706905851860660448207621209813287860733969412
81142660418086830619328460811191061556940512689692
51934325451728388641918047049293215058642563049483
62467221648435076201727918039944693004732956340691
15732444386908125794514089057706229429197107928209
55037687525678773091862540744969844508330393682126
18336384825330154686196124348767681297534375946515
80386287592878490201521685554828717201219257766954
78182833757993103614740356856449095527097864797581
16726320100436897842553539920931837441497806860984
48403098129077791799088218795327364475675590848030
87086987551392711854517078544161852424320693150332
59959406895756536782107074926966537676326235447210
69793950679652694742597709739166693763042633987085
41052684708299085211399427365734116182760315001271
65378607361501080857009149939512557028198746004375
35829035317434717326932123578154982629742552737307
94953759765105305946966067683156574377167401875275
88902802571733229619176668713819931811048770190271
25267680276078003013678680992525463401061632866526
36270218540497705585629946580636237993140746255962
24074486908231174977792365466257246923322810917141
91430288197103288597806669760892938638285025333403
34413065578016127815921815005561868836468420090470
23053081172816430487623791969842487255036638784583
11487696932154902810424020138335124462181441773470
63783299490636259666498587618221225225512486764533
67720186971698544312419572409913959008952310058822
95548255300263520781532296796249481641953868218774
76085327132285723110424803456124867697064507995236
37774242535411291684276865538926205024910326572967
23701913275725675285653248258265463092207058596522
29798860272258331913126375147341994889534765745501
18495701454879288984856827726077713721403798879715
38298203783031473527721580348144513491373226651381
34829543829199918180278916522431027392251122869539
40957953066405232632538044100059654939159879593635
29746152185502371307642255121183693803580388584903
41698116222072977186158236678424689157993532961922
62467957194401269043877107275048102390895523597457
23189706772547915061505504953922979530901129967519
86188088225875314529584099251203829009407770775672
11306739708304724483816533873502340845647058077308
82959174767140363198008187129011875491310547126581
97623331044818386269515456334926366572897563400500
42846280183517070527831839425882145521227251250327
55121603546981200581762165212827652751691296897789
32238195734329339946437501907836945765883352399886
75506164965184775180738168837861091527357929701337
62177842752192623401942399639168044983993173312731
32924185707147349566916674687634660915035914677504
99518671430235219628894890102423325116913619626622
73267460800591547471830798392868535206946944540724
76841822524674417161514036427982273348055556214818
97142617910342598647204516893989422179826088076852
87783646182799346313767754307809363333018982642090
10848802521674670883215120185883543223812876952786
71329612474782464538636993009049310363619763878039
62184073572399794223406235393808339651327408011116
66627891981488087797941876876144230030984490851411
60661826293682836764744779239180335110989069790714
85786944089552990653640447425576083659976645795096
66024396409905389607120198219976047599490197230297
64913982680032973156037120041377903785566085089252
16730939319872750275468906903707539413042652315011
94809377245048795150954100921645863754710598436791
78639167021187492431995700641917969777599028300699
15368713711936614952811305876380278410754449733078
40789923115535562561142322423255033685442488917353
44889911501440648020369068063960672322193204149535
41503128880339536053299340368006977710650566631954
81234880673210146739058568557934581403627822703280
82616570773948327592232845941706525094512325230608
22918802058777319719839450180888072429661980811197
77158542502016545090413245809786882778948721859617
72107838435069186155435662884062257473692284509516
20849603980134001723930671666823555245252804609722
53503534226472524250874054075591789781264330331690#

Gerry Jackson

2024-05-02 17:08:06 UTC

Post by minforth
Inspired by the Euler problem #13, I dug out one of my old
Forth systems after a long time. See below for the solution.
Can you do the same in your Forth system?
MinForth 3.6 (32 bit)
# fload euler13.mf
5537376230. solution: 5537376230. ok
This was made possible by a) high fp stack depth > 100, and
b) overflow when recognising integers (eg. within >NUMBER)
is treated as non-recognition, so the parsed string is passed
on to fp literal recognition.
Strictly speaking, this is obviously not a solution but an
approximation. Do you know a real solution, perhaps even
without loading a fat bignum library?
\ Work out the first ten digits of the sum of the
\ following one-hundred 50-digit numbers.

I'm no mathematician and don't understand Paul's or Ron's solutions.
Assuming the phrase 'the first 10 digits" means the 10 most significant
digits of the sum of the 100 numbers I think I would tackle it like this:

Use 64 bit GForth

1. Save the 100 numbers as strings with a 100 element array of addresses
of the strings, we know they are each 50 digits long so we don't need to
store the length.
2. For each integer string in turn
- select the least significant 10 digits (ca+40 10) where ca is a
string address in the array
- convert the 10 digit strings to integers
- accumulate the sum of the 100 integers
- divide the sum by 10^10 the quotient of which is the carry to the
next iteration
- repeat this using the next 10 digits i.e. (ca+30 10) adding in the
carry from above
- then repeat with (ca+20 10), (ca+10 10) (ca 10)
3. After the last iteration convert the sum of the 10 m.s. digits &
carry to a numeric string and select the 10 m.s. digits.

The maximum possible carry is when the digits are all 9s when the sum is
999999999900 and the carry to the second iteration is 99. So overflow
shouldn't be a problem.

Sorry I haven't the time to code a solution

--
Gerry

minforth

2024-05-02 21:14:13 UTC

My first similar idea was:
50 decimal digits are equivalent to 167 binary digits. So splitting
the numbers in half would allow the use double number arithmetic
of a 64-bit Forth system for each half column.

Krishna Myneni

2024-05-03 03:07:05 UTC

Post by minforth
50 decimal digits are equivalent to 167 binary digits. So splitting
the numbers in half would allow the use double number arithmetic
of a 64-bit Forth system for each half column.

Here's the first part of this to read in the numbers into an array, as
double length integers on a 64-bit Forth system (kforth64).

-- Krishna

=== begin code ===
\ euler13.4th
\
\ Work out the first ten digits of the sum of the
\ following one-hundred 50-digit numbers.

1 CELLS 8 < ABORT" Needs 64-bit Forth!"

10 constant LF
100 constant N
50 constant Ndig

create $nums N Ndig * allot
create Dnums[ N 2* 16 * allot
: ]D@ ( a idx -- ud ) 16 * + 2@ ;
: ]D! ( ud a idx -- ) 16 * + 2! ;
: $>ud ( a u -- ud ) 0 s>d 2swap >number 2drop ;
: $>2ud ( a u -- ud_low ud_high )
drop Ndig 2/ 2dup + over $>ud 2swap $>ud ;

\ Read N big numbers as strings and then parse into
\ double length array
: read-numbers ( -- )
N 0 DO
refill IF
LF parse dup Ndig <> ABORT" String Length Error!"
$nums Ndig I * + swap move
ELSE ABORT
THEN
LOOP
N 0 DO
$nums Ndig I * + Ndig
$>2ud Dnums[ I 2* 1+ ]D! Dnums[ I 2* ]D!
LOOP ;

\ read into array of 2*N doubles
read-numbers
37107287533902102798797998220837590246510135740250
46376937677490009712648124896970078050417018260538
74324986199524741059474233309513058123726617309629
...
=== end code ===

And, here's the check to ensure the numbers loaded into the array Dnums[

=== test above code ===
include euler13
ok
Dnums[ 0 ]D@ ud.
8220837590246510135740250 ok
Dnums[ 1 ]D@ ud.
3710728753390210279879799 ok
Dnums[ 2 ]D@ ud.
4896970078050417018260538 ok
Dnums[ 3 ]D@ ud.
4637693767749000971264812 ok
\ ...
Dnums[ 198 ]D@ ud.
4075591789781264330331690 ok
Dnums[ 199 ]D@ ud.
5350353422647252425087405 ok
=== end test ===

minforth

2024-05-03 03:33:00 UTC

My, this is cute! READ-NUMBERS using REFILL while including the file
(which uses REFILL itself)!

minforth

2024-05-03 03:41:05 UTC

Post by minforth
My, this is cute! READ-NUMBERS using REFILL while including the file
(which uses REFILL itself)!

To give credit where credit is due, Ron's solution follows a similar
path (although I have to admit that the ultra-compact 8th code looks
rather cryptic to me).

Ron AARON

2024-05-03 05:56:21 UTC

Post by minforth
My, this is cute! READ-NUMBERS using REFILL while including the file
(which uses REFILL itself)!

To give credit where credit is due, Ron's solution follows a similar
path (although I have to admit that the ultra-compact 8th code looks
rather cryptic to me).

Here's a reworking of it being much less "clever", and with commentary:

\ The numbers are in an array:
[37107287533902102798797998220837590246510135740250,
46376937677490009712648124896970078050417018260538,
74324986199524741059474233309513058123726617309629,
91942213363574161572522430563301811072406154908250,
23067588207539346171171980310421047513778063246676,
89261670696623633820136378418383684178734361726757,
28112879812849979408065481931592621691275889832738,
44274228917432520321923589422876796487670272189318,
47451445736001306439091167216856844588711603153276,
70386486105843025439939619828917593665686757934951,
62176457141856560629502157223196586755079324193331,
64906352462741904929101432445813822663347944758178,
92575867718337217661963751590579239728245598838407,
58203565325359399008402633568948830189458628227828,
80181199384826282014278194139940567587151170094390,
35398664372827112653829987240784473053190104293586,
86515506006295864861532075273371959191420517255829,
71693888707715466499115593487603532921714970056938,
54370070576826684624621495650076471787294438377604,
53282654108756828443191190634694037855217779295145,
36123272525000296071075082563815656710885258350721,
45876576172410976447339110607218265236877223636045,
17423706905851860660448207621209813287860733969412,
81142660418086830619328460811191061556940512689692,
51934325451728388641918047049293215058642563049483,
62467221648435076201727918039944693004732956340691,
15732444386908125794514089057706229429197107928209,
55037687525678773091862540744969844508330393682126,
18336384825330154686196124348767681297534375946515,
80386287592878490201521685554828717201219257766954,
78182833757993103614740356856449095527097864797581,
16726320100436897842553539920931837441497806860984,
48403098129077791799088218795327364475675590848030,
87086987551392711854517078544161852424320693150332,
59959406895756536782107074926966537676326235447210,
69793950679652694742597709739166693763042633987085,
41052684708299085211399427365734116182760315001271,
65378607361501080857009149939512557028198746004375,
35829035317434717326932123578154982629742552737307,
94953759765105305946966067683156574377167401875275,
88902802571733229619176668713819931811048770190271,
25267680276078003013678680992525463401061632866526,
36270218540497705585629946580636237993140746255962,
24074486908231174977792365466257246923322810917141,
91430288197103288597806669760892938638285025333403,
34413065578016127815921815005561868836468420090470,
23053081172816430487623791969842487255036638784583,
11487696932154902810424020138335124462181441773470,
63783299490636259666498587618221225225512486764533,
67720186971698544312419572409913959008952310058822,
95548255300263520781532296796249481641953868218774,
76085327132285723110424803456124867697064507995236,
37774242535411291684276865538926205024910326572967,
23701913275725675285653248258265463092207058596522,
29798860272258331913126375147341994889534765745501,
18495701454879288984856827726077713721403798879715,
38298203783031473527721580348144513491373226651381,
34829543829199918180278916522431027392251122869539,
40957953066405232632538044100059654939159879593635,
29746152185502371307642255121183693803580388584903,
41698116222072977186158236678424689157993532961922,
62467957194401269043877107275048102390895523597457,
23189706772547915061505504953922979530901129967519,
86188088225875314529584099251203829009407770775672,
11306739708304724483816533873502340845647058077308,
82959174767140363198008187129011875491310547126581,
97623331044818386269515456334926366572897563400500,
42846280183517070527831839425882145521227251250327,
55121603546981200581762165212827652751691296897789,
32238195734329339946437501907836945765883352399886,
75506164965184775180738168837861091527357929701337,
62177842752192623401942399639168044983993173312731,
32924185707147349566916674687634660915035914677504,
99518671430235219628894890102423325116913619626622,
73267460800591547471830798392868535206946944540724,
76841822524674417161514036427982273348055556214818,
97142617910342598647204516893989422179826088076852,
87783646182799346313767754307809363333018982642090,
10848802521674670883215120185883543223812876952786,
71329612474782464538636993009049310363619763878039,
62184073572399794223406235393808339651327408011116,
66627891981488087797941876876144230030984490851411,
60661826293682836764744779239180335110989069790714,
85786944089552990653640447425576083659976645795096,
66024396409905389607120198219976047599490197230297,
64913982680032973156037120041377903785566085089252,
16730939319872750275468906903707539413042652315011,
94809377245048795150954100921645863754710598436791,
78639167021187492431995700641917969777599028300699,
15368713711936614952811305876380278410754449733078,
40789923115535562561142322423255033685442488917353,
44889911501440648020369068063960672322193204149535,
41503128880339536053299340368006977710650566631954,
81234880673210146739058568557934581403627822703280,
82616570773948327592232845941706525094512325230608,
22918802058777319719839450180888072429661980811197,
77158542502016545090413245809786882778948721859617,
72107838435069186155435662884062257473692284509516,
20849603980134001723930671666823555245252804609722,
53503534226472524250874054075591789781264330331690]
\ "reduce" the array, adding each value
\ to the accumulator (0 to begin with)
' n:+ 0 a:reduce
\ convert the number to a string, lop off the 10
\ left-most characters and print them

Post by minforth
s 10 s:lsub . cr

bye

Krishna Myneni

2024-05-03 21:11:51 UTC

Post by minforth
My, this is cute! READ-NUMBERS using REFILL while including the file
(which uses REFILL itself)!

Nothing new. This is a commonly used technique in Forth to parse data
which needs special handling, without doing file i/o.

--
KM

Krishna Myneni

2024-05-03 22:02:55 UTC

Post by Krishna Myneni

Post by minforth
50 decimal digits are equivalent to 167 binary digits. So splitting
the numbers in half would allow the use double number arithmetic
of a 64-bit Forth system for each half column.

Here's the first part of this to read in the numbers into an array, as
double length integers on a 64-bit Forth system (kforth64).
-- Krishna
=== begin code ===
\ euler13.4th
\
\ Work out the first ten digits of the sum of the
\ following one-hundred 50-digit numbers.
1 CELLS 8 < ABORT" Needs 64-bit Forth!"
10 constant LF
100 constant N
50 constant Ndig
create $nums N Ndig * allot
create Dnums[ N 2* 16 * allot
: ]D! ( ud a idx -- ) 16 * + 2! ;
: $>ud ( a u -- ud ) 0 s>d 2swap >number 2drop ;
: $>2ud ( a u -- ud_low ud_high )
   drop Ndig 2/ 2dup + over $>ud 2swap $>ud ;
\ Read N big numbers as strings and then parse into
\ double length array
: read-numbers ( -- )
    N 0 DO
        refill IF
          LF parse dup Ndig <> ABORT" String Length Error!"
          $nums Ndig I * + swap move
        ELSE ABORT
        THEN
    LOOP
    N 0 DO
      $nums Ndig I * + Ndig
      $>2ud Dnums[ I 2* 1+ ]D! Dnums[ I 2* ]D!
    LOOP ;
\ read into array of 2*N doubles
read-numbers
37107287533902102798797998220837590246510135740250
46376937677490009712648124896970078050417018260538
74324986199524741059474233309513058123726617309629
...
=== end code ===
And, here's the check to ensure the numbers loaded into the array Dnums[
=== test above code ===
include euler13
ok
8220837590246510135740250 ok
3710728753390210279879799 ok
4896970078050417018260538 ok
4637693767749000971264812 ok
\ ...
4075591789781264330331690 ok
5350353422647252425087405 ok
=== end test ===

Here's the full program to solve Euler13 using addition of double
numbers on a 64-bit Forth system using only standard words and no
floating point. With a little work, it can be used to print the full sum
as well.

-- KM

=== begin code ===
\ euler13.4th
\
\ Work out the first ten digits of the sum of the
\ following one-hundred 50-digit numbers.

1 CELLS 8 < ABORT" Needs 64-bit Forth!"

10 constant LF
100 constant N
50 constant Ndig

create $nums N Ndig * allot
create Dnums[ N 2* 16 * allot
: ]D@ ( a idx -- ud ) 16 * + 2@ ;
: ]D! ( ud a idx -- ) 16 * + 2! ;
: $>ud ( a u -- ud ) 0 s>d 2swap >number 2drop ;
: $>2ud ( a u -- ud_low ud_high )
drop Ndig 2/ 2dup + over $>ud 2swap $>ud ;

\ Read N big numbers as strings and then parse into
\ double length array
: read-numbers ( -- )
N 0 DO
refill IF
LF parse dup Ndig <> ABORT" String Length Error!"
$nums Ndig I * + swap move
ELSE ABORT
THEN
LOOP
N 0 DO
$nums Ndig I * + Ndig
$>2ud Dnums[ I 2* 1+ ]D! Dnums[ I 2* ]D!
LOOP ;

read-numbers
37107287533902102798797998220837590246510135740250
46376937677490009712648124896970078050417018260538
74324986199524741059474233309513058123726617309629
91942213363574161572522430563301811072406154908250
23067588207539346171171980310421047513778063246676
89261670696623633820136378418383684178734361726757
28112879812849979408065481931592621691275889832738
44274228917432520321923589422876796487670272189318
47451445736001306439091167216856844588711603153276
70386486105843025439939619828917593665686757934951
62176457141856560629502157223196586755079324193331
64906352462741904929101432445813822663347944758178
92575867718337217661963751590579239728245598838407
58203565325359399008402633568948830189458628227828
80181199384826282014278194139940567587151170094390
35398664372827112653829987240784473053190104293586
86515506006295864861532075273371959191420517255829
71693888707715466499115593487603532921714970056938
54370070576826684624621495650076471787294438377604
53282654108756828443191190634694037855217779295145
36123272525000296071075082563815656710885258350721
45876576172410976447339110607218265236877223636045
17423706905851860660448207621209813287860733969412
81142660418086830619328460811191061556940512689692
51934325451728388641918047049293215058642563049483
62467221648435076201727918039944693004732956340691
15732444386908125794514089057706229429197107928209
55037687525678773091862540744969844508330393682126
18336384825330154686196124348767681297534375946515
80386287592878490201521685554828717201219257766954
78182833757993103614740356856449095527097864797581
16726320100436897842553539920931837441497806860984
48403098129077791799088218795327364475675590848030
87086987551392711854517078544161852424320693150332
59959406895756536782107074926966537676326235447210
69793950679652694742597709739166693763042633987085
41052684708299085211399427365734116182760315001271
65378607361501080857009149939512557028198746004375
35829035317434717326932123578154982629742552737307
94953759765105305946966067683156574377167401875275
88902802571733229619176668713819931811048770190271
25267680276078003013678680992525463401061632866526
36270218540497705585629946580636237993140746255962
24074486908231174977792365466257246923322810917141
91430288197103288597806669760892938638285025333403
34413065578016127815921815005561868836468420090470
23053081172816430487623791969842487255036638784583
11487696932154902810424020138335124462181441773470
63783299490636259666498587618221225225512486764533
67720186971698544312419572409913959008952310058822
95548255300263520781532296796249481641953868218774
76085327132285723110424803456124867697064507995236
37774242535411291684276865538926205024910326572967
23701913275725675285653248258265463092207058596522
29798860272258331913126375147341994889534765745501
18495701454879288984856827726077713721403798879715
38298203783031473527721580348144513491373226651381
34829543829199918180278916522431027392251122869539
40957953066405232632538044100059654939159879593635
29746152185502371307642255121183693803580388584903
41698116222072977186158236678424689157993532961922
62467957194401269043877107275048102390895523597457
23189706772547915061505504953922979530901129967519
86188088225875314529584099251203829009407770775672
11306739708304724483816533873502340845647058077308
82959174767140363198008187129011875491310547126581
97623331044818386269515456334926366572897563400500
42846280183517070527831839425882145521227251250327
55121603546981200581762165212827652751691296897789
32238195734329339946437501907836945765883352399886
75506164965184775180738168837861091527357929701337
62177842752192623401942399639168044983993173312731
32924185707147349566916674687634660915035914677504
99518671430235219628894890102423325116913619626622
73267460800591547471830798392868535206946944540724
76841822524674417161514036427982273348055556214818
97142617910342598647204516893989422179826088076852
87783646182799346313767754307809363333018982642090
10848802521674670883215120185883543223812876952786
71329612474782464538636993009049310363619763878039
62184073572399794223406235393808339651327408011116
66627891981488087797941876876144230030984490851411
60661826293682836764744779239180335110989069790714
85786944089552990653640447425576083659976645795096
66024396409905389607120198219976047599490197230297
64913982680032973156037120041377903785566085089252
16730939319872750275468906903707539413042652315011
94809377245048795150954100921645863754710598436791
78639167021187492431995700641917969777599028300699
15368713711936614952811305876380278410754449733078
40789923115535562561142322423255033685442488917353
44889911501440648020369068063960672322193204149535
41503128880339536053299340368006977710650566631954
81234880673210146739058568557934581403627822703280
82616570773948327592232845941706525094512325230608
22918802058777319719839450180888072429661980811197
77158542502016545090413245809786882778948721859617
72107838435069186155435662884062257473692284509516
20849603980134001723930671666823555245252804609722
53503534226472524250874054075591789781264330331690

2variable dsum_high
2variable dsum_low

: sumDnums ( -- )
0 s>d 2dup dsum_high 2! dsum_low 2!
N 2* 0 DO
I 1 and IF
dsum_high 2@ Dnums[ I ]D@ D+ dsum_high 2!
ELSE
dsum_low 2@ Dnums[ I ]D@ D+ dsum_low 2!
THEN
LOOP ;

sumDnums
\ dsum_low 2@ ud. cr
\ dsum_high 2@ ud. cr

\ add carry portion of low to high for unsigned 25 digit
\ numbers to find the new high portion.

: udsum_25 ( udhigh udlow -- udhigh' )
10000000000000000000 um/mod nip 1000000 / s>d d+ ;

\ print the leading 10 digits of the sum
: print-leading10 ( -- )
dsum_high 2@ dsum_low 2@ udsum_25
100000000000000000 um/mod nip u. ;

print-leading10 cr

=== end code ===

=== run output ===
include euler13
5537376230
ok
=== end output ===

Ahmed

2024-05-04 14:21:39 UTC

Hi,

What about this program. It is written in gforth.
the numbers are in the file "euler13.input" (Anton's file).
I reversed the strings holding the numbers and used a decimal-based adder
with carry (digit by digit) (as in primary school) but the summation goes
from left to right.
I didn't use foats or double numbers (just strings).
I think it is generalizable for any size (digits or numbers)

\ here, the program begins

s" euler13.input" slurp-file 2constant input

50 constant ns
ns 2 + constant ns+2
100 constant N

20 value ndigits

create input_rev N ns * allot
input_rev N ns * erase

: >input_rev
N 0 do
ns 0 do
input drop ns 1- i - j ns+2 * + + c@ input_rev i + ns j
* + c!
loop
loop ;

input_rev

create sum ns+2 allot sum ns+2 erase
create sv ns+2 allot sv ns+2 erase
create cry ns+2 allot cry ns+2 erase
create num ns+2 allot num ns+2 erase

: init_sum sum ns+2 erase ;
: init_cry cry ns+2 erase ;
: init_num num ns+2 erase ;

: d>c 48 + ;
: c>d 48 - ;

: >num ( n --) ns * input_rev + ns 0 do dup i + c@ c>d num i + c! loop
drop ;

: dda ( r --)
dup cry + c@ over num + c@ + over sum + c@ + 10 /mod rot swap over 1+ cry
+ c! sum + c! ;

: (nsad) ns+2 0 do i dda loop ;
: nsad init_sum N 0 do init_num init_cry i >num (nsad) loop ;

: >sv ns+2 0 do i sum + c@ d>c sv ns+2 1 - i - + c! loop ;

: .sum >sv sv ns+2 type ;

: 1st_nz_digit 0 ns+2 0 do sv i + c@ 0 d>c = if 1+ else unloop exit then
loop ;

: .sum_ndigits >sv sv 1st_nz_digit + ndigits type ;

nsad
cr .sum
cr 10 to ndigits .sum_ndigits
cr 20 to ndigits .sum_ndigits

\ here, the end of the program.

the results:
5537376230390876637302048746832985971773659831892672 for the whole sum
5537376230 for 10 digits
55373762303908766373 for 20 digits

Ahmed

Krishna Myneni

2024-05-04 17:12:00 UTC

Post by Ahmed
Hi,
What about this program. It is written in gforth.
the numbers are in the file "euler13.input" (Anton's file).
I reversed the strings holding the numbers and used a decimal-based adder
with carry (digit by digit) (as in primary school) but the summation goes
from left to right.
I didn't use foats or double numbers (just strings).
I think it is generalizable for any size (digits or numbers)
\ here, the program begins
s" euler13.input" slurp-file 2constant input
50 constant ns
ns 2 + constant ns+2
100 constant N
20 value ndigits
create input_rev N ns * allot
input_rev N ns * erase
: >input_rev    N 0 do            ns 0 do                   input drop
* + c!
           loop
   loop ;

Post by Ahmed
input_rev

create sum ns+2 allot sum ns+2 erase
create sv ns+2 allot sv ns+2 erase
create cry ns+2 allot cry ns+2 erase
create num ns+2 allot num ns+2 erase
: init_sum sum ns+2 erase ;
: init_cry cry ns+2 erase ;
: init_num num ns+2 erase ;
: d>c 48 + ;
: c>d 48 - ;
drop ;
: dda ( r --)
+ c! sum + c! ;
: (nsad) ns+2 0 do i dda loop ;
: nsad init_sum N 0 do init_num init_cry i >num (nsad) loop ;
: .sum >sv sv ns+2 type ;
loop ;
: .sum_ndigits >sv sv 1st_nz_digit + ndigits type ;
nsad
cr .sum
cr 10 to ndigits .sum_ndigits
cr 20 to ndigits .sum_ndigits
\ here, the end of the program.
5537376230390876637302048746832985971773659831892672 for the whole sum
5537376230 for 10 digits
55373762303908766373 for 20 digits

Nice. However, it's probably not a practical way to do arithmetic of big
numbers, compared to adding binary representations of the number in
larger chunks than single decimal digits. The machine addition
instructions are suited for 32/64 bit additions. Even my approach using
25 decimal digits at a time is slow because of the need to compute the
carry bits for adding to the high 25 digits. Better to have a contiguous
binary representation of the decimal number across the required memory
size needed to hold the sum.

--
Krishna

Ahmed

2024-05-04 17:56:42 UTC

You are right about number representations in memory (binary
representation).
But I find that the execution of nsad is fast (nsad calculates the sum).

for the speed:
tested with gforth:
utime nsad utime d>f d>f f- 1e-6 f* f. ." seconds"
gives about 0.334 ms

tested with vfxforth:
: timing_1000 timer-reset 1000 0 do nsad loop .elapsed ;
gives 125 ms for 1000 times nasd execution
so 0.125 ms for one execution of nsad

Ahmed

Anton Ertl

2024-05-04 17:42:37 UTC

Post by Krishna Myneni
Nice. However, it's probably not a practical way to do arithmetic of big
numbers, compared to adding binary representations of the number in
larger chunks than single decimal digits. The machine addition
instructions are suited for 32/64 bit additions. Even my approach using
25 decimal digits at a time is slow because of the need to compute the
carry bits for adding to the high 25 digits. Better to have a contiguous
binary representation of the decimal number across the required memory
size needed to hold the sum.

In the Euler 13 case, each input number is used only once, and as I
explained in <***@mips.complang.tuwien.ac.at>, converting
an input number to a binary representation costs more than adding it
digit-by-digit (i.e., column by column), which my program from that
posting does.

Concerning Ahmed's program, the description sounds like it does the
same, but there is the somewhat expensive "10 /mod" in DDA, which is
called by the inner loop. Does this program try to perform the
addition rowwise rather than columnwise?

- anton

--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
New standard: https://forth-standard.org/
EuroForth 2023: https://euro.theforth.net/2023

Krishna Myneni

2024-05-04 18:02:45 UTC

In the Euler 13 case, each input number is used only once, and as I
an input number to a binary representation costs more than adding it
digit-by-digit (i.e., column by column), which my program from that
posting does.

Thanks, I'll try it.

--
KM

Ahmed

2024-05-04 18:26:14 UTC

Post by Anton Ertl
Does this program try to perform the
addition rowwise rather than columnwise?

The sum is done as we do it manually with pencil.
the array: sum works as an accumulator.
we add a 50-digit number at once.
the array: cry holds the carries.

it is a decimal copy of the binary adder based on the binary (1-bit)
full-adder
____
c --->| |---> s
a --->| |
b --->| |---> c new carry
----

c
a
+ b
-------
c s
where the sum s = a XOR b and the carry c = a AND b
So based on this, a decimal full-adder (if I can call it like this) (one
decimal digit at once)

and the trick is that I reversed the strings holding the numbers to add
together and the addition is applied from left to right digit after digit.

Ahmed

Anton Ertl

2024-05-05 14:32:54 UTC

Post by Ahmed

Post by Anton Ertl
Does this program try to perform the
addition rowwise rather than columnwise?

The sum is done as we do it manually with pencil.
the array: sum works as an accumulator.
we add a 50-digit number at once.

But when I sum multiple numbers by hand, I don't use an accumulator
for whole numbers. I add up all the digits of one column and the
carry from the previous digit colum (which is in the range 0..99 in
the Euler #13 case). In the Euler #13 case this produces a number in
the range 0..999, and then I do the 10 /mod (via # in my program),
with the remainder being the result digit for the column, and the
quotient (in the range 0..99) being the carry for the next column.

- anton

Ahmed

2024-05-05 17:36:49 UTC

Post by Anton Ertl
But when I sum multiple numbers by hand, I don't use an accumulator
for whole numbers. I add up all the digits of one column and the
carry from the previous digit colum (which is in the range 0..99 in
the Euler #13 case). In the Euler #13 case this produces a number in
the range 0..999, and then I do the 10 /mod (via # in my program),
with the remainder being the result digit for the column, and the
quotient (in the range 0..99) being the carry for the next column.

Yes, you are right, we can do it that way.
In my program, you can consider the array sum as the result where we add up the 50-digit numbers one by one.

Here is a pseudo-code:

let n[0...99] given numbers (50-digit numbers)
let s = 0 initlalize the result which is the sum itself

for i from 0 to 99 do
s = s + n[i] here add up the numbers (accumulation)(but here, I do the sum digit by digit (and carries)
end for

display s

so you can consider what I called accumulator as the result holder

Ahmed

a***@spenarnc.xs4all.nl

2024-05-05 18:42:47 UTC

Post by Ahmed

Post by Anton Ertl
Does this program try to perform the
addition rowwise rather than columnwise?

The sum is done as we do it manually with pencil.
the array: sum works as an accumulator.
we add a 50-digit number at once.

The problem should include some numbers less that 50 digits,
then use your method, running in a 16 bit Forth.
That would at least more interesting than resorting to infinite
precision numbers. At the time it was intended that the numbers
would overflow precision.

Post by Anton Ertl
- anton

Groetjes Albert

--
Don't praise the day before the evening. One swallow doesn't make spring.
You must not say "hey" before you have crossed the bridge. Don't sell the
hide of the bear until you shot it. Better one bird in the hand than ten in
the air. First gain is a cat purring. - the Wise from Antrim -

Ahmed

2024-05-04 18:51:14 UTC

Post by Anton Ertl
but there is the somewhat expensive "10 /mod" in DDA, which is
called by the inner loop.

I replaced the 10 /mod by: cs defined as:
: cs ( n -- s c) dup 10 < if 0 exit then 10 - 1 ;

The code works fine but with gforth the execution time of nsad changes to
405ms (instead of 334ms with 10 /mod). So with 10 /mod the code is faster!
With vfxforth, no change in the execution time between the two versions.

Ahmed

Ahmed

2024-05-04 18:57:24 UTC

Post by Ahmed
The code works fine but with gforth the execution time of nsad changes to
405ms (instead of 334ms with 10 /mod). So with 10 /mod the code is faster!
With vfxforth, no change in the execution time between the two versions.

I forgot to mention that these timings are for 1000 executions of nsad.
So for just one execution of nsad: 0.406 ms and 0.125 ms for gforth and
vfxforth respectively

Ahmed

minforth

2024-05-05 05:23:17 UTC

Post by Ahmed
What about this program. It is written in gforth.
the numbers are in the file "euler13.input" (Anton's file).
I reversed the strings holding the numbers and used a decimal-based adder
with carry (digit by digit) (as in primary school) but the summation goes
from left to right.
I didn't use floats or double numbers (just strings).
I think it is generalizable for any size (digits or numbers)

This reminds me of BCD coded arithmetic. Sadly, BCD is hardly used
today or is only directly supported by a few CPUs, despite its great
advantages, especially for financial calculations. Some programming
languages
therefore support the numeric type BigDecimal (in addition to BigNum and
sometimes BigFloat, usually via the GMP library), such as Java, JS or
Python.

The approach of adding decimal digits individually in sufficiently large
buffers (for solving Euler problem #13) can be extended with little effort
to
multiplication, with the basic operation "multiplication of a large
decdigit
buffer by a single decimal digit".

This would also have been my approach to solving Euler problem #20 in
Forth:
How many decimal digits are in the number factorial(100)? (648)

Paul Rubin

2024-05-05 06:38:49 UTC

Post by minforth
This would also have been my approach to solving Euler problem #20 in
How many decimal digits are in the number factorial(100)? (648)

Brute force multiplication with 64 bit floats works fine:

: foo 101 1 ?do i s>f f* loop ; ok
1e foo f.s <1> 9.3326215444E157 ok

Ahmed

2024-05-05 06:56:30 UTC

Post by minforth
How many decimal digits are in the number factorial(100)? (648)

The question is:
Find the sum of the digits in the number factorial(100)
and not the number of digits in factorial(100)

With python:

def fact(n):
p = 1
for i in range(2,n+1):
p = p*i
return p

def sdf(n):
return sum([int(i) for i in list(str(fact(n)))])

then
sdf(100) gives 684
and sdf(1000) gives 10539

I didn't use recurson fin the definition of fact because of the depth of
recursion is not sufficent for 1000

Ahmed

Paul Rubin

2024-05-05 07:53:24 UTC

Post by Ahmed
Find the sum of the digits in the number factorial(100)
and not the number of digits in factorial(100)

Ah ok, I missed that part. Yes, doing it with bignums is easy. Maybe
there is some clever trick for doing it with machine integers.

minforth

2024-05-05 06:59:18 UTC

Post by Paul Rubin

Post by minforth
This would also have been my approach to solving Euler problem #20 in
How many decimal digits are in the number factorial(100)? (648)

: foo 101 1 ?do i s>f f* loop ; ok
1e foo f.s <1> 9.3326215444E157 ok

.. with the minor nuisance that the exponent should be 647

Paul Rubin

2024-05-05 08:00:59 UTC

Post by Paul Rubin
1e foo f.s <1> 9.3326215444E157 ok

.. with the minor nuisance that the exponent should be 647

No, that is the right value of 100!. You left out the part about
summing the digits:

digits n | n == 0 = [0]
| otherwise = b : digits a
where (a,b) = n `quotRem` 10

main = print . sum . digits . product $ [1..10000]

prints 648.

mhx

2024-05-05 08:41:28 UTC

(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : Collection at
http://projecteuler.net/index.php?section=problems
* CATEGORY : Contests
* AUTHOR : Marcel Hendrix
* LAST CHANGE : April 26, 2008, Marcel Hendrix
*)

NEEDS -miscutil
NEEDS -bignum

REVISION -euler20 "--- Sum of digits in n! Version 1.00 ---"

PRIVATES

DOC
(*
n! means n * (n - 1) * ... * 3 * 2 * 1

Find the sum of the digits in the number 100!
*)
ENDDOC

MAX.DIGITS BIGNUM n! PRIVATE

: Euler20 ( -- )
TIMER-RESET
1 n! V! #101 2 DO n! I VS* LOOP
CR ." 100! = " n! .Vnum
CR ." The sum of the digits of the number 100! = "
0 BEGIN
n! #10 VS/MOD +
n! VS0=
UNTIL
U>D (n,3)
CR .ELAPSED ;

:ABOUT CR ." Euler20 -- Finds the sum of the digits in the number 100!" ;

.ABOUT -euler20 CR
DEPRIVE

(* End of Source *)

FORTH> in euler20
Euler20 -- Finds the sum of the digits in the number 100!
ok
FORTH> Euler20
100! = 93,326,215,443,944,152,681,699,238,856,266,700,490,715,968,
264,381,621,468,592,963,895,217,599,993,229,915,608,941,463,976,
156,518,286,253,697,920,827,223,758,251,185,210,916,864,000,000,
000,000,000,000,000,000

The sum of the digits of the number 100! = 648

Slightly interesting that the last 24 digits are '0'.

-marcel

Anton Ertl

2024-05-05 11:24:44 UTC

Post by mhx
The sum of the digits of the number 100! = 648
Slightly interesting that the last 24 digits are '0'.

10 from 10 20 ... 100.
10 from 5 15 ... 95; there are enough even numbers to produce a 0 from that.
4 additional ones from 25 50 75 100

An additional one would come from multiples of 125, but that's not
included in 100!.

- anton

a***@spenarnc.xs4all.nl

2024-05-05 08:43:50 UTC

Post by Ahmed
What about this program. It is written in gforth.
the numbers are in the file "euler13.input" (Anton's file).
I reversed the strings holding the numbers and used a decimal-based

adder

Post by Ahmed
with carry (digit by digit) (as in primary school) but the summation

goes

Post by Ahmed
from left to right.
I didn't use floats or double numbers (just strings).
I think it is generalizable for any size (digits or numbers)

This reminds me of BCD coded arithmetic. Sadly, BCD is hardly used
today or is only directly supported by a few CPUs, despite its great
advantages, especially for financial calculations.

You mean to imply that DAA DAS AAA AAS AAM AAD instructions are not
available to 32 and/or 64 bits intel/AMD CPU's? To spare a few hundred
of the millions of transistors on the die? At the cost of complaining
by MSDOS emulator builders?
That has to be corroborated, but I think it unlikely.

Groetjes Albert

mhx

2024-05-05 10:24:10 UTC

Post by a***@spenarnc.xs4all.nl

[..]

Post by minforth
This reminds me of BCD coded arithmetic. Sadly, BCD is hardly used
today or is only directly supported by a few CPUs, despite its great
advantages, especially for financial calculations.

ISTR that AAM is just a renamed opcode for "mul immediate by 10" and
something simular for AAD.
This led me to believe that these instructions are just for
convenience, and receive no (hardware) optimization (anymore).
There is a nice (F.)-like opcode also, but it does not support
80-bit numbers, only 56 DP (?) or even less (I did not check
if somebody hacked around that)?

-marcel

minforth

2024-05-05 15:15:24 UTC

Post by a***@spenarnc.xs4all.nl
You mean to imply that DAA DAS AAA AAS AAM AAD instructions are not
available to 32 and/or 64 bits intel/AMD CPU's? To spare a few hundred
of the millions of transistors on the die? At the cost of complaining
by MSDOS emulator builders?
That has to be corroborated, but I think it unlikely.

There are indeed a few BCD artefacts in Intel CPU instruction sets.
But they are as good as useless.
https://en.wikipedia.org/wiki/Intel_BCD_opcodes

Anton Ertl

2024-05-05 11:31:32 UTC

What are these advantages supposed to be? The reason why this stuff
has vanished is that it has no significant advantages. The only
advantage is in cases like Euler 13 where there is hardly any
arithmetic between input and output, but such problems use so little
CPU time that this advantage is not significant. For financial
calculations use binary fixed point with an appropriate (decimal)
scale factor.

Post by minforth
Some programming
languages
therefore support the numeric type BigDecimal (in addition to BigNum and
sometimes BigFloat, usually via the GMP library), such as Java, JS or
Python.

Looking at
<https://docs.oracle.com/javase/8/docs/api/java/math/BigDecimal.html>,
this looks like binary fixed point with a decimal scale factor.

- anton

minforth

2024-05-05 13:28:01 UTC

Post by Anton Ertl
Looking at
<https://docs.oracle.com/javase/8/docs/api/java/math/BigDecimal.html>,
this looks like binary fixed point with a decimal scale factor.

It is not quite clear to me from the documentation. I think
bignum = arbitrary precision math in base 2 (for the mantissa)
bigdecimal = arbitrary precision math in base 10 (for the mantissa)
The latter to avoid rounding errors like
0.1 + 0.2 -> 0.30000000000000004

Of course at the bottom, bigdecimal mantissas are not necesssarily
ASCII digit fields, but can be encoded into packed binary fields.

I don't know whether above assumed 'binary fixed point with scaling'
could eliminate the rounding errors.

Anton Ertl

2024-05-05 17:12:11 UTC

Post by Anton Ertl
Looking at
<https://docs.oracle.com/javase/8/docs/api/java/math/BigDecimal.html>,
this looks like binary fixed point with a decimal scale factor.

...

Post by minforth
bigdecimal = arbitrary precision math in base 10 (for the mantissa)

What makes you think so?

Post by minforth
The latter to avoid rounding errors like
0.1 + 0.2 -> 0.30000000000000004

There is no such rounding error for

1/10 + 2/10 -> 3/10

i.e., fixed point with a scale factor 1/10.

That's irrespective of whether the mantissa is represented in binary
or decimal. And given that

1) binary is faster

2) and they already have a binary arbitrary-length implementation

it would be stupid of them to use something else for the mantissa.

Post by minforth
I don't know whether above assumed 'binary fixed point with scaling'
could eliminate the rounding errors.

Of course not. There is no way to represent 1/3 as a decimal
fraction. The point of decimal arithmetic stuff for financial
computations is not to eliminate rounding errors, but to make the same
rounding errors as the decimal calculating machines of yesteryear, the
limitations of which found their way into various regulations.
Decimally scaled fixed point is exactly what you need for that.

Incidentially, on Monday I have been in the Arithmeum in Bonn, which
displays various calculating machines (and you can even play with some
of them). I also saw a Zuse Z25 in operation there (but only used as
a calculator).

- anton

minforth

2024-05-05 19:31:08 UTC

Post by Anton Ertl
Looking at
<https://docs.oracle.com/javase/8/docs/api/java/math/BigDecimal.html>,
this looks like binary fixed point with a decimal scale factor.

....

Post by minforth
bigdecimal = arbitrary precision math in base 10 (for the mantissa)

What makes you think so?

Thank you for the clarification. I must have been on the wrong path for a
few years. BigFloat seems to be a term only in Julia (Julia uses the GMP
and MPFR libraries). Regarding BigDecimal, I seem to be confused by a very
old perhaps naive implementation I had seen many years ago. In this
encoding of mantissas, the input/output of numbers did not require
conversion of the number base, similar to BCD arithmetic.

I did some research, and there does indeed seem to be no significant
difference (except in rounding behaviour) between BigFloat and BigDecimal.
It's a bit confusing because everyone can use the names however they want,
similar to BigInt and BigNum.

Translated with DeepL.com (free version)

Post by minforth
The latter to avoid rounding errors like
0.1 + 0.2 -> 0.30000000000000004

There is no such rounding error for
1/10 + 2/10 -> 3/10
i.e., fixed point with a scale factor 1/10.

This is actually from Javascript documentation
https://www.npmjs.com/package/bigdecimal

minforth

2024-05-05 19:48:25 UTC

I forgot to mention the actual Javascript proposal for base-10 encoding
of BigDecimals
https://github.com/tc39/proposal-decimal

It is already available in the libbf library by Fabrice Bellard.

minforth

2024-05-09 06:02:49 UTC

There is an optional libdecnumber library for gcc, with DPD packed decimal
digits in 128 bit mantissa acc. to IEEE754-2008.

Could also be implemented in Forth, but who uses softfloats anyway ...

Krishna Myneni

2024-05-04 17:57:33 UTC

On 5/3/24 17:02, Krishna Myneni wrote:
...

Post by Krishna Myneni
Here's the full program to solve Euler13 using addition of double
numbers on a 64-bit Forth system using only standard words and no
floating point. With a little work, it can be used to print the full sum
as well.
-- KM
=== begin code ===
\ euler13.4th
\
\ Work out the first ten digits of the sum of the
\ following one-hundred 50-digit numbers.
1 CELLS 8 < ABORT" Needs 64-bit Forth!"
10 constant LF
100 constant N
50 constant Ndig
create $nums N Ndig * allot
create Dnums[ N 2* 16 * allot
: ]D! ( ud a idx -- ) 16 * + 2! ;
: $>ud ( a u -- ud ) 0 s>d 2swap >number 2drop ;
: $>2ud ( a u -- ud_low ud_high )
   drop Ndig 2/ 2dup + over $>ud 2swap $>ud ;
\ Read N big numbers as strings and then parse into
\ double length array
: read-numbers ( -- )
    N 0 DO
        refill IF
          LF parse dup Ndig <> ABORT" String Length Error!"
          $nums Ndig I * + swap move
        ELSE ABORT
        THEN
    LOOP
    N 0 DO
      $nums Ndig I * + Ndig
      $>2ud Dnums[ I 2* 1+ ]D! Dnums[ I 2* ]D!
    LOOP ;
read-numbers
37107287533902102798797998220837590246510135740250
...
2variable dsum_high
2variable dsum_low
: sumDnums ( -- )
    0 s>d 2dup dsum_high 2! dsum_low 2!
    N 2* 0 DO
      I 1 and IF
      ELSE
      THEN
    LOOP ;
sumDnums
\ add carry portion of low to high for unsigned 25 digit
\ numbers to find the new high portion.
: udsum_25 ( udhigh udlow -- udhigh' )
   10000000000000000000 um/mod nip 1000000 / s>d d+ ;
\ print the leading 10 digits of the sum
: print-leading10 ( -- )
    100000000000000000 um/mod nip u. ;
print-leading10 cr
=== end code ===
=== run output ===
include euler13
5537376230
ok
=== end output ===

Now that we have seen the code works, it's time to dispense with the
string storage and the array storage, which were useful for intermediate
diagnostics, and compact the code.

=== begin code ===
\ euler13-2.4th
\
\ Print the first ten digits of the sum of a
\ sequence of N 50-digit numbers.

1 CELLS 8 < ABORT" Needs 64-bit Forth!"

10 constant LF
50 constant Ndig

: $>ud ( a u -- ud ) 0 s>d 2swap >number 2drop ;
: $>2ud ( a u -- ud_low ud_high )
drop Ndig 2/ 2dup + over $>ud 2swap $>ud ;
: d+! ( d a -- ) dup >r 2@ d+ r> 2! ;

\ add carry portion of udlow to udhigh for 25 decimal
\ digit numbers; return new high portion.
: udsum_25 ( udhigh udlow -- udhigh' )
10000000000000000000 um/mod nip 1000000 / s>d d+ ;

\ print the leading 10 digits
: print-leading10 ( udhigh udlow -- )
udsum_25 100000000000000000 um/mod nip u. ;

\ Sum n big numbers as two doubles, holding the
\ high and low 25 decimal digits
2variable dsum_high
2variable dsum_low
: sum-numbers ( n -- udhigh udlow )

Post by Krishna Myneni
r 0 s>d 2dup dsum_high 2! dsum_low 2!

dsum_high dsum_low
r> 0 DO
refill IF
LF parse dup Ndig <> ABORT" String Length Error!"
$>2ud 5 pick D+! 2 pick D+!
ELSE ABORT
THEN
LOOP 2drop dsum_high 2@ dsum_low 2@ ;

100 sum-numbers
37107287533902102798797998220837590246510135740250
\ ...
print-leading10 cr

=== end code ===

=== run output ===
include euler13-2
5537376230
ok
.s
<empty>
ok
=== end output ===

Paul Rubin

2024-05-03 03:54:10 UTC

Post by Gerry Jackson
I'm no mathematician and don't understand Paul's or Ron's
solutions. Assuming the phrase 'the first 10 digits" means the 10 most
significant digits of the sum of the 100 numbers I think I would

Wait what, I think I have #13 as the wrong problem? Yes, I had a saved
euler13.fs file that I posted, but it was actually for problem 14. Not
sure how that happened, sorry.

Yeah for #13, I would just add up the leftmost 14 or so digits of each
number as 64-bit ints, then check that the result was not anywhere near
causing an overflow in the top 10 digits. In the unlikely case where
that is an issue, use multi-precision. Or in a language with native
bignums, the whole thing becomes trivial.

minforth

2024-05-03 06:53:48 UTC

Post by Paul Rubin
Yeah for #13, I would just add up the leftmost 14 or so digits of each
number as 64-bit ints, then check that the result was not anywhere near
causing an overflow in the top 10 digits. In the unlikely case where
that is an issue, use multi-precision.

I did check that with quad precision floats:
MinForth 3.6 (64 bit) (fp 128 bit)
# fload euler13.mf
5537376230. answer: 5537376230. ok
#

However even 112 mantissa bits are not sufficient for lossless conversion
of integers with 50 decimal digits.

Paul Rubin

2024-05-03 07:46:23 UTC

Post by minforth
However even 112 mantissa bits are not sufficient for lossless conversion
of integers with 50 decimal digits.

You don't have to add up all 50. You only have to be sure that the
lower 40 don't overflow into the top 10 without your noticing it.

mhx

2024-05-03 08:06:49 UTC

MinForth 3.6 (64 bit) (fp 128 bit)
# fload euler13.mf
5537376230. answer: 5537376230. ok
#
However even 112 mantissa bits are not sufficient for lossless

conversion

Post by minforth
of integers with 50 decimal digits.

Hmm, the solution must be somewhere in your idea: when summing
N numbers with many digits and then inspecting the K highest bits
of the sum, the MSB's won't be influenced by the M LSB's of *one* of
those N numbers. Now, what happens when you have N numbers, all with
M of the worst possible sets of LSBs? There has to be some arithmetic
rule that explains this. (There are people who have worked out all
of the Euler problems with pencil and paper).
The worst I can imagine is that problem is defined to be just on the
boundary what is possible, or that some other trick is needed to
prove that the approach works in the given case although it violates
the main rule.

It is nice that your double-double floats signal overflow, than you
don't need to find the rule :--)

I think this gives me the idea how to *experimentally* check how
many bits are needed.

-marcel

a***@spenarnc.xs4all.nl

2024-05-03 09:01:07 UTC

In article <***@www.novabbs.com>,
mhx <***@iae.nl> wrote:
<SNIP>

Post by mhx
M of the worst possible sets of LSBs? There has to be some arithmetic
rule that explains this. (There are people who have worked out all
of the Euler problems with pencil and paper).

Make that "some". There are certainly problems that cannot be handled
this way.

Post by mhx
-marcel

mhx

2024-05-03 09:39:30 UTC

By text editor (replace '\n' by 'e f+\n' and insert a
14th column of full-stops):

0e
3710728753390.2102798797998220837590246510135740250e f+
4637693767749.0009712648124896970078050417018260538e f+
7432498619952.4741059474233309513058123726617309629e f+
..
2084960398013.4001723930671666823555245252804609722e f+
5350353422647.2524250874054075591789781264330331690e f+
cr +e.

\ 5.5373762303908766364e+0014

For 100 numbers as shown, their fractions can sum to
at most 100. As the numbers have 13 digit significands,
summing 100 of them needs only 12 digits. So the first
10 digits of 5.5373762303908766364e+0014 should be correct:
"5537376230". (As already shown by minforth.)

-marcel

a***@spenarnc.xs4all.nl

2024-05-03 08:57:33 UTC

Post by Paul Rubin

Wait what, I think I have #13 as the wrong problem? Yes, I had a saved
euler13.fs file that I posted, but it was actually for problem 14. Not
sure how that happened, sorry.
Yeah for #13, I would just add up the leftmost 14 or so digits of each
number as 64-bit ints, then check that the result was not anywhere near
causing an overflow in the top 10 digits. In the unlikely case where
that is an issue, use multi-precision. Or in a language with native
bignums, the whole thing becomes trivial.

That is sensible. First inspect that the file contain numbers of the
same length. Then trim it at length 30.

I even used a totally superfluous 2VARIABLE of sorts, thus
preventing stack juggling.

Using scripting:
------------------
#!/usr/bin/lina64 -s
CREATE total 2 CELLS ALLOT 0. total 2!
1 ARG[] GET-FILE
BEGIN OVER WHILE
^J $/ DROP 30 0. 2SWAP >NUMBER 2DROP total 2@ D+ total 2!
REPEAT 2DROP
total 2@ D. CR
------------------

Trimmed:
------------------
#!/usr/bin/lina -s
0.0 1 ARG[] GET-FILE
BEGIN OVER WHILE
^J $/ DROP 30 0. 2SWAP >NUMBER 2DROP 2SWAP D+ 2SWAP
REPEAT 2DROP
D. CR
------------------
Then round the resulting number by hand.

P.S. It is time that the usefulness of $/ finally sinks in.

Groetjes Albert

mhx

2024-05-03 18:55:34 UTC

Paul Rubin wrote:
[..]

Post by Paul Rubin
Wait what, I think I have #13 as the wrong problem? Yes, I had a saved
euler13.fs file that I posted, but it was actually for problem 14. Not
sure how that happened, sorry.

Look Mum, no recursion.

(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : Collection at
http://projecteuler.net/index.php?section=problems
* CATEGORY : Contests
* AUTHOR : Marcel Hendrix
* LAST CHANGE : April 21, 2008, Marcel Hendrix
*)

NEEDS -miscutil

REVISION -euler14 "--- Longest sequence Version 1.00 ---"

PRIVATES

DOC
(*
The following iterative sequence is defined for the set of positive
integers:

n -> n/2 (n is even)
n -> 3n + 1 (n is odd)

Using the rule above and starting with 13, we generate the following
sequence:

13 -> 40 -> 20 -> 10 -> 5 -> 16 -> 8 -> 4 -> 2 -> 1

It can be seen that this sequence (starting at 13 and finishing at 1)
contains
10 terms. Although it has not been proved yet (Collatz Problem), it is
thought
that all starting numbers finish at 1.

Which starting number, under one million, produces the longest chain?

NOTE: Once the chain starts the terms are allowed to go above one
million.
*)
ENDDOC

0. DVALUE greatest PRIVATE

: iteration ( dn - dm )
OVER 1 AND
IF ( odd) 2DUP D2* D+ 1. D+
2DUP greatest DMAX TO greatest
ELSE ( even) D2/
ENDIF ; PRIVATE

: sequence ( n -- dk ) \ dk = #terms
1. DLOCAL dk
S>D
BEGIN ( dn -- )
iteration ( dm -- )
1. +TO dk ( dm -- )
2DUP 1. D=
UNTIL 2DROP dk ; PRIVATE

0. DVALUE maxchain

: Euler14
MS? LOCAL btime
CLEAR maxchain CLEAR greatest
CR ." iter# length elapsed"
CR ." --------------------------"
#1000000
1 DO
I sequence ( dk -- )
2DUP maxchain D>
IF ( dk -- )
2DUP TO maxchain
CR I 6 .R 3 SPACES 2DUP D. #16 HTAB MS? btime - 6 .R ." ms "
ENDIF 2DROP
LOOP
CR ." greatest number encountered = " greatest (n,3) ;

:ABOUT CR ." Try: Euler14 -- Find longest sequence"
CR
CR ." \ ...37min 53sec till ...39min 30sec = 01min30sec"
CR ." \ PowerBook G4 (867 MHz PowerPC G4)" ;

.ABOUT -euler14 CR
DEPRIVE

(* End of Source *)

Anton Ertl

2024-05-03 08:29:27 UTC

***@gmx.net (minforth) writes:
[print the first 10 digits of the sum of 100 50-digit numbers]

Post by minforth
Strictly speaking, this is obviously not a solution but an
approximation. Do you know a real solution, perhaps even
without loading a fat bignum library?

Gerry Jackson's answer inspired my solution: sum up each column, the
way you would do it by hand. This not only avoids the need for
bignums and the uncertainty about the correctness that some of the
other solutions have, it is also the most efficient: In this digitwise
solution, every digit is only added, while in the conversion
solutions, each digit costs at least an add, a subtract, and a
multiply (and more in the case of FP); my solution does the subtracts
for the whole column once instead of 100 times:

s" euler13.input" slurp-file 2constant input

: sumcolumn ( ucarry1 c-addr -- ucarry2 )
5100 bounds ?do
i c@ +
51 +loop
'0' 100 * - 0 # drop ;

: sum ( -- c-addr2 u2 )
0 <<#
input drop 50 1 mem-do
i sumcolumn
loop
0 #s #> ;

sum drop 10 type cr

This outputs "5537376230".

I still use # for convenience; that could be optimized (e.g., divide
by 10000000000 every 10 columns for the last 40 columns, then use
um/mod and add '0' explicitly to the remainder).

- anton

minforth

2024-05-03 09:52:24 UTC

Your solution is good.
I had toyed with the idea of building a tiny bignum library myself.
Binary operations, and with intrinsics, addition/subtraction are also
trivial. Multiplication wouldn't be difficult either, for division
the shift-subtract algorithm would be sufficient for me. Let's see,
maybe something for the next rainy season... ;-)

a***@spenarnc.xs4all.nl

2024-05-03 11:00:00 UTC

The modulo operation is a real hassle. (We did it for the transputer.)

A more fruitful idea is to open a pipe to dc.

The heart of it is:

: cmd1 dc-EMIT ;
: command dc-TYPE ^J cmd1 ;
...
: +h &+ cmd1 ;
: -h &- cmd1 ;
: *h &* cmd1 ;
: /h &/ cmd1 ;
: MODh &% cmd1 ;
: /MODh &~ cmd1 &r cmd1 ;
: **h &^ cmd1 ;
: */h "s_ * l_ /" command ;
: |h &| cmd1 ;

...

I started dc using a shell instead of a straight fork,
it worked, but for this reason it was not good enough to publish it.

In gforth it is probably easier to fork off dc, and you can
get a powerful easy wordlist.

Groetjes Albert

minforth

2024-05-03 11:39:37 UTC

Post by a***@spenarnc.xs4all.nl

Post by minforth
I had toyed with the idea of building a tiny bignum library myself.
Binary operations, and with intrinsics, addition/subtraction are also
trivial. Multiplication wouldn't be difficult either, for division
the shift-subtract algorithm would be sufficient for me. Let's see,
maybe something for the next rainy season... ;-)

The modulo operation is a real hassle. (We did it for the transputer.)

Out of the blue: IIRC the shift-subtract algo generates quotient AND
modulus as stop criterion.

a***@spenarnc.xs4all.nl

2024-05-05 09:20:25 UTC

Post by a***@spenarnc.xs4all.nl

The modulo operation is a real hassle. (We did it for the transputer.)

Out of the blue: IIRC the shift-subtract algo generates quotient AND
modulus as stop criterion.

This is one of the algorithms of Knuth. It is in fact a shift-subtract
algorithm but using 32 bits chunks instead of 1 bit.
I would be interested to see the result of a simpler shift-subtract
to compare speed.

Groetjes Albert

Anton Ertl

2024-05-03 12:22:51 UTC

Post by Anton Ertl
Gerry Jackson's answer inspired my solution: sum up each column, the
way you would do it by hand. This not only avoids the need for
bignums and the uncertainty about the correctness that some of the
other solutions have, it is also the most efficient: In this digitwise
solution, every digit is only added, while in the conversion
solutions, each digit costs at least an add, a subtract, and a
multiply (and more in the case of FP); my solution does the subtracts
s" euler13.input" slurp-file 2constant input
: sumcolumn ( ucarry1 c-addr -- ucarry2 )
5100 bounds ?do
51 +loop
'0' 100 * - 0 # drop ;
: sum ( -- c-addr2 u2 )
0 <<#
input drop 50 1 mem-do
i sumcolumn
loop
0 #s #> ;
sum drop 10 type cr
This outputs "5537376230".

One additional observation inspires an optimization: You can add up to
25 digits without overflowing a byte. So you can do the sum of up to
25 sequences of digits in SIMD style with the ordinary +; with 8-byte
cells, you can add 8 digits with one +. With AVX-512, you can add all
the digits of one line at once, but unfortunately, Gforth does not
support this properly (and my vector package is designed for long
vectors and has too much overhead for short vectors). So I
implemented this for working with 8-byte cells;

here 5114 allot 14 + constant input

s" euler13.input" r/o open-file throw
input 5100 rot read-file throw 5100 <> throw

create sum25buf 256 allot

: sum25 ( -- )
sum25buf 8 + input 6 - 4 0 do { d s }
s 0 25 0 do ( addr w )
over @ $ffff000000000000 and + swap 51 + swap loop
$3030000000000000 25 * - d ! drop
d cell+ s 56 + s cell+ do
i 0 25 0 do ( daddr saddr w )
over @ + swap 51 + swap loop
nip $3030303030303030 25 * - over !
cell+
cell +loop
drop
d 64 + s 25 51 * +
loop
2drop ;

: sumcolumn ( ucarry1 c-addr -- ucarry2 )
256 bounds ?do
i c@ +
64 +loop
0 # drop ;

: sum ( -- c-addr2 u2 )
sum25
0 <<#
sum25buf 14 + 50 1 mem-do
i sumcolumn
loop
0 #s #> ;

sum drop 10 type cr

This is quite messy, especially dealing with the two extra bytes per
line. Does it at least pay off in speed? Calling it with

perf stat ~/gforth/gforth-fast <euler13>.4th -e "#>> : foo 1000000 0 do sum 2drop #>> loop ; foo bye"

i.e., 1M calls to SUM costs (on a 4GHz Skylake):

euler13 euler13-faster
30_647_182_330 9_653_276_959 cycles
75_201_467_616 24_413_536_115 instructions

7.859516000 2.411344000 seconds user

So more than a factor of 3. 9653 cycles for a call to SUM (in
euler13-faster) is not bad. A quick measurement indicates that about
1400 cycles of that are spent in # and #S.

- anton

PMF

2024-05-03 16:40:07 UTC

One additional observation inspires an optimization: You can add up to
25 digits without overflowing a byte. So you can do the sum of up to
25 sequences of digits in SIMD style with the ordinary +; with 8-byte
cells, you can add 8 digits with one +. With AVX-512, you can add all
the digits of one line at once, but unfortunately, Gforth does not
support this properly (and my vector package is designed for long
vectors and has too much overhead for short vectors). So I
implemented this for working with 8-byte cells;
here 5114 allot 14 + constant input
s" euler13.input" r/o open-file throw
input 5100 rot read-file throw 5100 <> throw
create sum25buf 256 allot
: sum25 ( -- )
sum25buf 8 + input 6 - 4 0 do { d s }
s 0 25 0 do ( addr w )
$3030000000000000 25 * - d ! drop
d cell+ s 56 + s cell+ do
i 0 25 0 do ( daddr saddr w )
nip $3030303030303030 25 * - over !
cell+
cell +loop
drop
d 64 + s 25 51 * +
loop
2drop ;
: sumcolumn ( ucarry1 c-addr -- ucarry2 )
256 bounds ?do
64 +loop
0 # drop ;
: sum ( -- c-addr2 u2 )
sum25
0 <<#
sum25buf 14 + 50 1 mem-do
i sumcolumn
loop
0 #s #> ;
sum drop 10 type cr
This is quite messy, especially dealing with the two extra bytes per
line. Does it at least pay off in speed? Calling it with
perf stat ~/gforth/gforth-fast <euler13>.4th -e "#>> : foo 1000000 0 do

sum 2drop #>> loop ; foo bye"

Post by Anton Ertl
euler13 euler13-faster
30_647_182_330 9_653_276_959 cycles
75_201_467_616 24_413_536_115 instructions
7.859516000 2.411344000 seconds user
So more than a factor of 3. 9653 cycles for a call to SUM (in
euler13-faster) is not bad. A quick measurement indicates that about
1400 cycles of that are spent in # and #S.
- anton

The last few month I have been learning Python. My daughter is taking an
university course in it.
To help her I have also followed her course.

Seeing this problem looked like a school example for Python
''' Solve Euler 13 '''

fh = open('euler13.dat', 'r')

sum = 0

for line in fh :
sum = sum + int(line)

print(sum)

Executing this gives

D:\dev\Python>python3 euler13.py
5537376230390876637302048746832985971773659831892672

BR
Peter

mhx

2024-05-03 18:44:52 UTC