Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #10531 > unrolled thread

Sum of squares

Started bymhx@iae.nl (Marcel Hendrix)
First post2012-03-26 20:50 +0200
Last post2012-03-29 09:05 +0200
Articles 20 on this page of 24 — 10 participants

Back to article view | Back to comp.lang.forth

Contents

  Sum of squares mhx@iae.nl (Marcel Hendrix) - 2012-03-26 20:50 +0200
    Re: Sum of squares "A. K." <akk@nospam.org> - 2012-03-26 21:43 +0200
    Re: Sum of squares Paul Rubin <no.email@nospam.invalid> - 2012-03-26 16:52 -0700
      Re: Sum of squares Helmar Wodtke <helmwo@gmail.com> - 2012-03-27 03:15 -0700
        Re: Sum of squares Paul Rubin <no.email@nospam.invalid> - 2012-03-28 08:16 -0700
      Re: Sum of squares BruceMcF <agila61@netscape.net> - 2012-03-28 06:06 -0700
    Re: Sum of squares Zbiggy <zbigniew2011REMOVE@gmail.REMOVE.com> - 2012-03-27 12:13 +0100
      Re: Sum of squares Helmar Wodtke <helmwo@gmail.com> - 2012-03-27 04:03 -0700
        Re: Sum of squares Zbiggy <zbigniew2011REMOVE@gmail.REMOVE.com> - 2012-03-27 13:55 +0100
    Re: Sum of squares Albert van der Horst <albert@spenarnc.xs4all.nl> - 2012-03-27 11:05 +0000
      Re: Sum of squares Albert van der Horst <albert@spenarnc.xs4all.nl> - 2012-03-28 11:41 +0000
        Re: Sum of squares Paul Rubin <no.email@nospam.invalid> - 2012-03-28 07:41 -0700
        Re: Sum of squares mhx@iae.nl (Marcel Hendrix) - 2012-03-29 22:06 +0200
        Re: Sum of squares Albert van der Horst <albert@spenarnc.xs4all.nl> - 2012-03-31 11:48 +0000
      Re: Sum of squares Alexander Adler <alexander.adler@stud.uni-frankfurt.de> - 2012-03-29 19:33 +0200
        Re: Sum of squares Helmar Wodtke <helmwo@gmail.com> - 2012-03-29 11:16 -0700
          Re: Sum of squares Albert van der Horst <albert@spenarnc.xs4all.nl> - 2012-04-03 09:47 +0000
      Re: Sum of squares Helmar Wodtke <helmwo@gmail.com> - 2012-03-28 09:42 -0700
      Re: Sum of squares Helmar Wodtke <helmwo@gmail.com> - 2012-03-28 11:08 -0700
    Re: Sum of squares Ecki <ecki@intershop.de> - 2012-03-27 16:15 +0200
      Re: Sum of squares Helmar Wodtke <helmwo@gmail.com> - 2012-03-27 07:35 -0700
        Re: Sum of squares hwfwguy@gmail.com - 2012-03-27 08:00 -0700
          Re: Sum of squares Helmar Wodtke <helmwo@gmail.com> - 2012-03-27 08:21 -0700
          Re: Sum of squares Ecki <ecki@intershop.de> - 2012-03-29 09:05 +0200

Page 1 of 2  [1] 2  Next page →

#10531 — Sum of squares

Frommhx@iae.nl (Marcel Hendrix)
Date2012-03-26 20:50 +0200
SubjectSum of squares
Message-ID<10851509008435@frunobulax.edu>
Admit it (when trying it yourself), this exercise is 
surprisingly tricky. 

-marcel

-- -------------------------------------------
(*
 * LANGUAGE    : ANS Forth with extensions
 * PROJECT     : Forth Environments
 * DESCRIPTION : Sum of squares of two largest of three values
 * CATEGORY    : Programming exercise 
 * AUTHOR      : Marcel Hendrix 
 * LAST CHANGE : March 25, 2012, Marcel Hendrix 
 *)



	NEEDS -miscutil

	REVISION -sumsquares "--- Sum of Squares      Version 0.01 ---"

	PRIVATES

DOC
(*

  Sum of squares of two largest of three values

  This exercise comes to us from the book "Structure and Interpretation of 
  Computer Programs" by Abelson and Sussman (exercise 1.3):
  
  Define a procedure that takes three numbers as arguments and returns 
  the sum of the squares of the two larger numbers.
*)
ENDDOC

: DSQR ( u -- du ) DUP M* ; PRIVATE
: .ANSWER ( ud -- ) CR ." The sum of the squares of the two larger numbers is " (n,3) ; PRIVATE

: SOS ( n1 n2 n3 -- )
	2DUP > IF  SWAP  ENDIF ( x1) DSQR 2>R
	MAX ( x2) DSQR 2R> D+ .ANSWER ;


: SQR ( u1 -- u2 ) DUP * ; PRIVATE
: SUM-OF-SQUARES ( n1 n2 n3 -- ) 2DUP > IF  SWAP  ENDIF  SQR -ROT MAX SQR + . ;

:ABOUT	CR ." Print the sum of the squares of the two larger numbers from {n1, n2, n3}."
	CR ." The numbers should have less than 18 digits."
     	CR ." Try: 1 2 3 SOS            -- print 13, double precision." 
	CR ."      1 2 3 SUM-OF-SQUARES -- single precision." ;

		.ABOUT -sumsquares CR
		DEPRIVE

                              (* End of Source *)

[toc] | [next] | [standalone]

#10534

From"A. K." <akk@nospam.org>
Date2012-03-26 21:43 +0200
Message-ID<4f70c6e1$0$6635$9b4e6d93@newsspool2.arcor-online.net>
In reply to#10531
On 26.03.2012 20:50, Marcel Hendrix wrote:
> Admit it (when trying it yourself), this exercise is
> surprisingly tricky.

Not when you do it with only one argument.
:-)

[toc] | [prev] | [next] | [standalone]

#10541

FromPaul Rubin <no.email@nospam.invalid>
Date2012-03-26 16:52 -0700
Message-ID<7x62dqapjr.fsf@ruckus.brouhaha.com>
In reply to#10531
mhx@iae.nl (Marcel Hendrix) writes:
>   Define a procedure that takes three numbers as arguments and returns 
>   the sum of the squares of the two larger numbers.

: 2sort ( x y -- largest smallest ) 2dup < if swap endif ;
: 3sort ( x y z -- largest middle smallest )
    2sort -rot 2sort rot 2sort ;
: square ( n -- n^2 ) dup * ;
: sum2 ( x y z -- largest^2 + middle^2 )
    3sort drop square swap square + ;

[toc] | [prev] | [next] | [standalone]

#10560

FromHelmar Wodtke <helmwo@gmail.com>
Date2012-03-27 03:15 -0700
Message-ID<23403720.182.1332843312154.JavaMail.geo-discussion-forums@ynlw24>
In reply to#10541
Am Dienstag, 27. März 2012 01:52:08 UTC+2 schrieb Paul Rubin:
> mhx@iae.nl (Marcel Hendrix) writes:
> >   Define a procedure that takes three numbers as arguments and returns 
> >   the sum of the squares of the two larger numbers.
> 
> : 2sort ( x y -- largest smallest ) 2dup < if swap endif ;
> : 3sort ( x y z -- largest middle smallest )
>     2sort -rot 2sort rot 2sort ;
> : square ( n -- n^2 ) dup * ;
> : sum2 ( x y z -- largest^2 + middle^2 )
>     3sort drop square swap square + ;

Mhm, why not (using your idea):

: keep    2dup > if swap then ;
: ^2      dup * ;
: result  keep -rot keep nip
          ^2 swap ^2 + ;

Regards,
-Helmar

[toc] | [prev] | [next] | [standalone]

#10628

FromPaul Rubin <no.email@nospam.invalid>
Date2012-03-28 08:16 -0700
Message-ID<7xd37wwybx.fsf@ruckus.brouhaha.com>
In reply to#10560
Helmar Wodtke <helmwo@gmail.com> writes:
> Mhm, why not (using your idea):
>
> : keep    2dup > if swap then ;
> : ^2      dup * ;
> : result  keep -rot keep nip
>           ^2 swap ^2 + ;

I like this.  I wrote 3sort as a factor that was very easy to specify
and test, but that did more work than the problem actually required.
Your version has better Zen.

[toc] | [prev] | [next] | [standalone]

#10723

FromBruceMcF <agila61@netscape.net>
Date2012-03-28 06:06 -0700
Message-ID<1e74021b-3100-4043-bcfc-98f8cf931336@z5g2000yqj.googlegroups.com>
In reply to#10541
On Mar 26, 7:52 pm, Paul Rubin <no.em...@nospam.invalid> wrote:
> m...@iae.nl (Marcel Hendrix) writes:
> >   Define a procedure that takes three numbers as arguments and returns
> >   the sum of the squares of the two larger numbers.
>
> : 2sort ( x y -- largest smallest ) 2dup < if swap endif ;
> : 3sort ( x y z -- largest middle smallest )
>     2sort -rot 2sort rot 2sort ;
> : square ( n -- n^2 ) dup * ;
> : sum2 ( x y z -- largest^2 + middle^2 )
>     3sort drop square swap square + ;

You only need max2 ~ the two inputs to the sum of squares don't have
to be in sorted order

: 2sort ( x y -- largest smallest ) 2dup < if swap then ;

: max2 ( n1 n2 n3 -- n4 n5 )
   \ n4,n5 largest two of n1,n2,n3
   sort2 rot max ;

: sum2squares ( n1 n2 -- n3 ) dup * swap dup * + ;

: sum2 ( n1 n2 n3 -- n4 ) max2 sum2squares ;

[toc] | [prev] | [next] | [standalone]

#10559

FromZbiggy <zbigniew2011REMOVE@gmail.REMOVE.com>
Date2012-03-27 12:13 +0100
Message-ID<slrnjn3bv2.2uh.zbigniew2011REMOVE@Tichy.myhome.org>
In reply to#10531
In comp.lang.forth, Marcel Hendrix wrote:

>   Define a procedure that takes three numbers as arguments and returns 
>   the sum of the squares of the two larger numbers.

: 3dup ( n1 n2 n3 -- n1 n2 n3 n1 n2 n3 )   dup 2over rot ;
: square ( n -- n^2 )   dup * ;
: solution
  3dup min min swap 2swap
  square swap square rot square + + swap square -
;

-- 
Forth is a preserver of health (Hippocrates)

[toc] | [prev] | [next] | [standalone]

#10565

FromHelmar Wodtke <helmwo@gmail.com>
Date2012-03-27 04:03 -0700
Message-ID<27162337.0.1332846184146.JavaMail.geo-discussion-forums@yngr3>
In reply to#10559
Am Dienstag, 27. März 2012 13:13:15 UTC+2 schrieb Zbiggy:
> In comp.lang.forth, Marcel Hendrix wrote:
> 
> >   Define a procedure that takes three numbers as arguments and returns 
> >   the sum of the squares of the two larger numbers.
> 
> : 3dup ( n1 n2 n3 -- n1 n2 n3 n1 n2 n3 )   dup 2over rot ;
> : square ( n -- n^2 )   dup * ;
> : solution
>   3dup min min swap 2swap
>   square swap square rot square + + swap square -
> ;
> 
> -- 
> Forth is a preserver of health (Hippocrates)

I've just for fun also a SWAP- and ROT-free version ;)

: ^2      dup * ;
: result  2dup max ^2 >r min max ^2 r> + ;

Regards,
-Helmar

[toc] | [prev] | [next] | [standalone]

#10567

FromZbiggy <zbigniew2011REMOVE@gmail.REMOVE.com>
Date2012-03-27 13:55 +0100
Message-ID<slrnjn3hve.3uv.zbigniew2011REMOVE@Tichy.myhome.org>
In reply to#10565
In comp.lang.forth, Helmar Wodtke wrote:

> I've just for fun also a SWAP- and ROT-free version ;)
>
>: ^2      dup * ;
>: result  2dup max ^2 >r min max ^2 r> + ;

You're right: smart use of min/max does the trick.
-- 
Forth is a preserver of health (Hippocrates)

[toc] | [prev] | [next] | [standalone]

#10562

FromAlbert van der Horst <albert@spenarnc.xs4all.nl>
Date2012-03-27 11:05 +0000
Message-ID<m1jito.mmz@spenarnc.xs4all.nl>
In reply to#10531
In article <10851509008435@frunobulax.edu>, Marcel Hendrix <mhx@iae.nl> wrote:
>Admit it (when trying it yourself), this exercise is
>surprisingly tricky.
>
>-marcel
>
>-- -------------------------------------------
>(*
> * LANGUAGE    : ANS Forth with extensions
> * PROJECT     : Forth Environments
> * DESCRIPTION : Sum of squares of two largest of three values
> * CATEGORY    : Programming exercise
> * AUTHOR      : Marcel Hendrix
> * LAST CHANGE : March 25, 2012, Marcel Hendrix
> *)
>
>
>
>       NEEDS -miscutil
>
>       REVISION -sumsquares "--- Sum of Squares      Version 0.01 ---"
>
>       PRIVATES
>
>DOC
>(*
>
>  Sum of squares of two largest of three values
>
>  This exercise comes to us from the book "Structure and Interpretation of
>  Computer Programs" by Abelson and Sussman (exercise 1.3):
>
>  Define a procedure that takes three numbers as arguments and returns
>  the sum of the squares of the two larger numbers.
>*)
>ENDDOC
>
>: DSQR ( u -- du ) DUP M* ; PRIVATE
>: .ANSWER ( ud -- ) CR  ." The sum of the squares of the two larger numbers is " (n,3) ; PRIVATE
>
>: SOS ( n1 n2 n3 -- )
>       2DUP > IF  SWAP  ENDIF ( x1) DSQR 2>R
>       MAX ( x2) DSQR 2R> D+ .ANSWER ;
>
>
>: SQR ( u1 -- u2 ) DUP * ; PRIVATE
>: SUM-OF-SQUARES ( n1 n2 n3 -- ) 2DUP > IF  SWAP  ENDIF  SQR -ROT MAX SQR + . ;
>
>:ABOUT CR ." Print the sum of the squares of the two larger numbers from {n1, n2, n3}."
>       CR ." The numbers should have less than 18 digits."
>       CR ." Try: 1 2 3 SOS            -- print 13, double precision."
>       CR ."      1 2 3 SUM-OF-SQUARES -- single precision." ;
>
>               .ABOUT -sumsquares CR
>               DEPRIVE
>
>                              (* End of Source *)
>

Why not start with
( n1 n2 n3 -- n4 n5 )   2MAX leave the two largest numbers of 3.
( a b -- csquare )   HYP2 leave the sum of the squares of a and b.

: SUM-OF-SQUARES 2MAX HYP2 ;

A much more interesting problem is:
   It is well-known that a prime is the sum of 2 squares in exactly
   one way 1]. Find them.

1] Unless of course it leaves 3 by division by 4. Those numbers have
   no solutions.

   Example:
   101 = 10^2 + 1^2
    97 = 9^2  + 4^2

This sometimes comes up as a sub problem in Euler, where it has to
be executed millions of times if not more.  So it must be fast.

First all , let us find some primes.
"
  AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
  ...
   OK
  INCLUDE rabbin.frt WANT NEW-IF
  ...
  AGAIN : ISN'T UNIQUE
   OK
  DO-DEBUG

  S[ ] OK 1,000,000,000,000,000,001 BEGIN 4 + DUP rprime? UNTIL
  S[ 1000000000000000009 ] OK BEGIN 4 + DUP rprime? UNTIL
  S[ 1000000000000000177 ] OK BEGIN 4 + DUP rprime? UNTIL
  S[ 1000000000000000201 ] OK BYE
"
Now finding solutions:

"
  AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
  ...
   OK
  INCLUDE findsq.frt WANT MARK-TIME
  ...
   OK
  MARK-TIME 1,000,000,000,000,000,009 solve-p . . ELAPSED .uS
  3 1000000000  517 uS  OK
  MARK-TIME 1,000,000,000,000,000,177 solve-p . . ELAPSED .uS
  633708561 773571884  535 uS  OK
  633708561 SQ DUP . 773571884 SQ DUP . + .
  401586540284690721 598413459715309456 1000000000000000177  OK
  BYE
"
The first solution is easy to find, but the second?
I won't spoil the fun by showing findsq.frt already.
It is non-trivial indeed and surprisingly short, once you
have enough modulo arithmetic in place.

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

[toc] | [prev] | [next] | [standalone]

#10623

FromAlbert van der Horst <albert@spenarnc.xs4all.nl>
Date2012-03-28 11:41 +0000
Message-ID<m1lf4n.9d@spenarnc.xs4all.nl>
In reply to#10562
In article <10851509008435@frunobulax.edu>, Marcel Hendrix <mhx@iae.nl> wrote:
>Admit it (when trying it yourself), this exercise is
>surprisingly tricky.
>
>-marcel
<SNIP>
>
>  Define a procedure that takes three numbers as arguments and returns
>  the sum of the squares of the two larger numbers.

A much more interesting problem is the other way around:
Define a procedure that takes a prime 1 mod 4 and returns two squares
whose sum is that prime.

(It is known that there is exactly one solution.)

EXAMPLE

1,000,000,000,000,000,177 is a prime

These are the squares
  401586540284690721 598413459715309456

As you can check by

  633708561 SQ DUP . 773571884 SQ DUP . + .

It takes lina64 about 600 microseconds to find those.

Groetjes Albert

(Paraphrased because my previous post was too boring.)
--
-- 
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

[toc] | [prev] | [next] | [standalone]

#10626

FromPaul Rubin <no.email@nospam.invalid>
Date2012-03-28 07:41 -0700
Message-ID<7x1uoc94a3.fsf@ruckus.brouhaha.com>
In reply to#10623
Albert van der Horst <albert@spenarnc.xs4all.nl> writes:
> A much more interesting problem is the other way around:
> Define a procedure that takes a prime 1 mod 4 and returns two squares
> whose sum is that prime.
>...
> (Paraphrased because my previous post was too boring.)

I thought the previous post was interesting but I didn't figure out a
solution in the couple of minutes I was able to think about it in.  I'll
keep trying.

[toc] | [prev] | [next] | [standalone]

#10640

Frommhx@iae.nl (Marcel Hendrix)
Date2012-03-29 22:06 +0200
Message-ID<85971306008435@frunobulax.edu>
In reply to#10623
Albert van der Horst <albert@spenarnc.xs4all.nl> writes Re: Sum of squares

[..]
> A much more interesting problem is the other way around:
> Define a procedure that takes a prime 1 mod 4 and returns two squares
> whose sum is that prime.

> (It is known that there is exactly one solution.)
[..]
> (Paraphrased because my previous post was too boring.)

I gratefully acknowledge Alexander Adler for posting his MOD-POW that 
allowed me to keep this code trivial.

-marcel

-- --------------------------------------------
(*
 * LANGUAGE    : ANS Forth with extensions
 * PROJECT     : Forth Environments
 * DESCRIPTION : a prime = 1 mod 4 is the sum of two squares
 * CATEGORY    : Puzzle 
 * AUTHOR      : Marcel Hendrix, MOD-POW from Alexander Adler <alexander.adler@stud.uni-frankfurt.de> 
 * LAST CHANGE : March 28, 2012, Marcel Hendrix 
 *)



	NEEDS -miscutil

	REVISION -brillhart "--- Prime Square Sum    Version 0.00 ---"

	PRIVATES

DOC
(*
	Define a procedure that takes a prime = 1 mod 4 and returns two squares
	whose sum is that prime.

	(It is known that there is exactly one solution.)

	EXAMPLE

	1,000,000,000,000,000,177 is a prime

	These are the squares
	  401586540284690721 598413459715309456

	As you can check by

	  633708561 SQ DUP . 773571884 SQ DUP . + .

 	Notes
	-----
	A solution can be found in: "Some Refinements of an Algorithm of Brillhart",
	by KENNETH S. WILLIAMS, Canadian Mathematical Society, Conference Proceedings,
	Volume 15, 1995.

	A posting by Alexander Adler on CLF showed a simple MOD-POW that could
	be used instead of the one in my highly iForth-specific BIGNUM library. This 
	MOD-POW is the reason that I posted my solution on CLF, too.
	
	Brillhart's algorithm
	---------------------

	(i) solve z for z^2 = -1 (mod p), where p is the given prime = 1 (mod 4). 
	This z can be found as z = c^{(p-1)/4} (mod p), where c is a quadratic 
	non-residue (mod p). It is trivial to find quadratic non-residues
	by using the Jacobi symbol (more convenient than the Legendre Symbol). Half 
	the numbers between 1 and p are nonresidues, picking numbers c at random and 
	calculating the Jacobi symbol (c|p) until a nonresidue is found will quickly 
	produce one. However, not all these c's are OK, as an additional condition is 
	that z should be < p/2. I solved this by trying again until a satisfactory z is
	returned.

	(ii) Apply the Euclidean algorithm to z and p (in that order) and determine
        the first remainder r_k (k > 0) satisfying r_k < sqrt(p). Then p = u^2 + v^2 
	with u = r_k, v = r_(k+1).

	\ FORTH> 61 go
	\ (2 us elapsed) 61 = 6^2 + 5^2 ok
	\ FORTH> #1000000000000000177 GO
	\ (19 us elapsed) 1000000000000000177 = 773571884^2 + 633708561^2 ok
*)
ENDDOC

( Remove all multiples of 4 from number u giving u' )
: KillFours ( u -- u' )
	BEGIN  DUP 3 AND 0=		\ Kill multiples of 4 in u1'
	WHILE  2 RSHIFT
	REPEAT ; PRIVATE

-- Like MOD but for unsigned numbers
-- Works also for both  > MAX-N
: UMOD	>R U>D R> SM/REM DROP ; 

: Jacobi-Symbol ( u1 u2 -- -1 | 0 | +1 )
	DUP 1 AND 0= ABORT" Huh"	
	2DUP UMOD ROT DROP		\ Subtract multiples of u2 from u1
	DUP 0=  IF  NIP EXIT  ENDIF	\ Break recursion if u1' is 0
	KillFours			\ u2 u1''
	DUP 1 = IF  NIP EXIT  ENDIF	\ Break recursion if u1' is 1
		    			\ This is the Jacobi result
                  			\ u2 u1''
	DUP 1 AND 
	   IF	2DUP AND 3 AND 3 = >R	\ Remember both were 3 mod 4
		RECURSE			\  with swapped numbers!
	        R> IF  NEGATE  ENDIF	\ Extra factor -1 if both 3 mod 4
	 ELSE   2/ OVER  RECURSE	\ Numbers not swapped!
		SWAP 7 AND DUP 3 =	\ JACOBI, u2
		SWAP 5 = 
		OR IF  NEGATE  ENDIF  	\ Take (2,u2) into account
	ENDIF ; PRIVATE


-- Half the numbers between 1 and p are nonresidues, picking numbers 
-- x at random and calculating the Jacobi symbol (x|p) until a nonresidue 
-- is found will quickly produce one.
: FIND-c ( p -- c )
	  LOCAL p
	0 LOCAL q
	BEGIN 
	  p CHOOSE TO q
	  q p Jacobi-Symbol 
	  -1 = 
	UNTIL q ;

: *-MODULO       ( x y MODULUS -- x/y mod MODULUS )  */MOD DROP ; PRIVATE 
: SQUARE-MODULO  ( x MODULUS -- x^2 mod MODULUS )    OVER SWAP *-MODULO ; PRIVATE  
: ODD?           ( x -- f )  1 AND ; 

: MOD-POW ( x y z -- x^y mod z )
	LOCAL modulus  
	1 SWAP ( x^{2^i} result y' )
	BEGIN DUP 
	WHILE
	   >R R@ ODD? IF  OVER modulus *-MODULO  ENDIF
	   SWAP modulus SQUARE-MODULO SWAP  
	   R> 2/
	REPEAT
	DROP NIP ;

-- Find solution z of z^2 = -1 (mod p) where 0 < z < p/2
-- Equivalent to z = c^{(p-1)/4} where c is a quadratic non-residue (mod p).
: (step) ( prime -- z )
	DUP FIND-c ( -- p c )
	OVER 1- 4 /  ROT ( -- c [p-1]/4 p ) MOD-POW ; PRIVATE

-- Do this until z < prime/2 
: step(i) ( prime -- z )
	LOCAL prime
	BEGIN  
	  prime (step) 
	  DUP prime 2/ >= 
	WHILE  DROP  
	REPEAT ; PRIVATE

: Brillhart ( z prime -- u v )
	0        LOCAL u
	FALSE    LOCAL stop?
	DUP SQRT LOCAL rootprime
	BEGIN	TUCK MOD
		stop? IF  NIP u  EXIT  ENDIF
		DUP TO u
		DUP rootprime < TO stop?
	AGAIN ;

: GO ( p -- ) 
	CR '(' EMIT  TICKS-RESET  
	   DUP step(i) OVER Brillhart 
	.UELAPSED ." ) " ROT . ." = " 0 .R ." ^2 + " 0 .R ." ^2" ;

:ABOUT	CR ." Try: #1000000000000000177 GO -- 773571884^2 + 633708561^2" 
	CR ." Or   #61 GO -- 5^+6^" ;

		.ABOUT -brillhart CR
		DEPRIVE

                              (* End of Source *)

[toc] | [prev] | [next] | [standalone]

#10663

FromAlbert van der Horst <albert@spenarnc.xs4all.nl>
Date2012-03-31 11:48 +0000
Message-ID<m1qzgg.5cj@spenarnc.xs4all.nl>
In reply to#10623
In article <m1lf4n.9d@spenarnc.xs4all.nl>,
Albert van der Horst  <albert@spenarnc.xs4all.nl> wrote:
>In article <10851509008435@frunobulax.edu>, Marcel Hendrix <mhx@iae.nl> wrote:
>>Admit it (when trying it yourself), this exercise is
>>surprisingly tricky.
>>
>>-marcel
><SNIP>
>>
>>  Define a procedure that takes three numbers as arguments and returns
>>  the sum of the squares of the two larger numbers.
>
>A much more interesting problem is the other way around:
>Define a procedure that takes a prime 1 mod 4 and returns two squares
>whose sum is that prime.
>
>(It is known that there is exactly one solution.)
>
>EXAMPLE
>
>1,000,000,000,000,000,177 is a prime
>
>These are the squares
>  401586540284690721 598413459715309456
>
>As you can check by
>
>  633708561 SQ DUP . 773571884 SQ DUP . + .
>
>It takes lina64 about 600 microseconds to find those.
>
>Groetjes Albert
>
>(Paraphrased because my previous post was too boring.)

Apparently there is one proper way to do it. The solutions of
Marcel, Alec and me are basically equivalent, but look at the
differences in coding style!

This is how I did it (taking x^x aka POW-MOD , SQRT and RAND
for granted).

\ --------------8<-----------------------

\ For P return R such that r*r+1 divides p.
: sqrt(-1)   >R BEGIN RAND R@ MOD R@ 1- 4 / R@ x^x DUP 2 R@ x^x 1+
     R@ MOD WHILE DROP REPEAT R> DROP ;

\ Given P and X with p|x^2+1  return A B with p=a^2+b^2
: reduce-mod   OVER SQRT >R   2DUP -   MIN   ( p, x with x<p/2)
  SWAP  BEGIN DUP R@ > WHILE  SWAP OVER MOD   REPEAT SWAP OVER MOD   R> DROP ;

\ For odd P return positive A and B with a^2+b^2= p.
: solve-p    DUP sqrt(-1) reduce-mod ;
\ --------------8<-----------------------

Note that like Marcel I used a more simple and less clear criterion
to stop the reduction step.
Here also `` tuck mod '' or `` SWAP OVER MOD '' is the essence.

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

[toc] | [prev] | [next] | [standalone]

#10639

FromAlexander Adler <alexander.adler@stud.uni-frankfurt.de>
Date2012-03-29 19:33 +0200
Message-ID<87y5qjpb09.fsf@stud.uni-frankfurt.de>
In reply to#10562
albert@spenarnc.xs4all.nl schrieb:
> In article <10851509008435@frunobulax.edu>, Marcel Hendrix <mhx@iae.nl> wrote:
>>Admit it (when trying it yourself), this exercise is
>>surprisingly tricky.
>>
>>-marcel
>>
>>-- -------------------------------------------
>>(*
>> * LANGUAGE    : ANS Forth with extensions
>> * PROJECT     : Forth Environments
>> * DESCRIPTION : Sum of squares of two largest of three values
>> * CATEGORY    : Programming exercise
>> * AUTHOR      : Marcel Hendrix
>> * LAST CHANGE : March 25, 2012, Marcel Hendrix
>> *)
>>
>>
>>
>>       NEEDS -miscutil
>>
>>       REVISION -sumsquares "--- Sum of Squares      Version 0.01 ---"
>>
>>       PRIVATES
>>
>>DOC
>>(*
>>
>>  Sum of squares of two largest of three values
>>
>>  This exercise comes to us from the book "Structure and Interpretation of
>>  Computer Programs" by Abelson and Sussman (exercise 1.3):
>>
>>  Define a procedure that takes three numbers as arguments and returns
>>  the sum of the squares of the two larger numbers.
>>*)
>>ENDDOC
>>
>>: DSQR ( u -- du ) DUP M* ; PRIVATE
>>: .ANSWER ( ud -- ) CR  ." The sum of the squares of the two larger numbers is " (n,3) ; PRIVATE
>>
>>: SOS ( n1 n2 n3 -- )
>>       2DUP > IF  SWAP  ENDIF ( x1) DSQR 2>R
>>       MAX ( x2) DSQR 2R> D+ .ANSWER ;
>>
>>
>>: SQR ( u1 -- u2 ) DUP * ; PRIVATE
>>: SUM-OF-SQUARES ( n1 n2 n3 -- ) 2DUP > IF  SWAP  ENDIF  SQR -ROT MAX SQR + . ;
>>
>>:ABOUT CR ." Print the sum of the squares of the two larger numbers from {n1, n2, n3}."
>>       CR ." The numbers should have less than 18 digits."
>>       CR ." Try: 1 2 3 SOS            -- print 13, double precision."
>>       CR ."      1 2 3 SUM-OF-SQUARES -- single precision." ;
>>
>>               .ABOUT -sumsquares CR
>>               DEPRIVE
>>
>>                              (* End of Source *)
>>
>
> Why not start with
> ( n1 n2 n3 -- n4 n5 )   2MAX leave the two largest numbers of 3.
> ( a b -- csquare )   HYP2 leave the sum of the squares of a and b.
>
> : SUM-OF-SQUARES 2MAX HYP2 ;
>
> A much more interesting problem is:
>    It is well-known that a prime is the sum of 2 squares in exactly
>    one way 1]. Find them.
>
> 1] Unless of course it leaves 3 by division by 4. Those numbers have
>    no solutions.
>
>    Example:
>    101 = 10^2 + 1^2
>     97 = 9^2  + 4^2
>
> This sometimes comes up as a sub problem in Euler, where it has to
> be executed millions of times if not more.  So it must be fast.
>
> First all , let us find some primes.
> "
>   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
>   ...
>    OK
>   INCLUDE rabbin.frt WANT NEW-IF
>   ...
>   AGAIN : ISN'T UNIQUE
>    OK
>   DO-DEBUG
>
>   S[ ] OK 1,000,000,000,000,000,001 BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000009 ] OK BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000177 ] OK BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000201 ] OK BYE
> "
> Now finding solutions:
>
> "
>   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
>   ...
>    OK
>   INCLUDE findsq.frt WANT MARK-TIME
>   ...
>    OK
>   MARK-TIME 1,000,000,000,000,000,009 solve-p . . ELAPSED .uS
>   3 1000000000  517 uS  OK
>   MARK-TIME 1,000,000,000,000,000,177 solve-p . . ELAPSED .uS
>   633708561 773571884  535 uS  OK
>   633708561 SQ DUP . 773571884 SQ DUP . + .
>   401586540284690721 598413459715309456 1000000000000000177  OK
>   BYE
> "
> The first solution is easy to find, but the second?
> I won't spoil the fun by showing findsq.frt already.
> It is non-trivial indeed and surprisingly short, once you
> have enough modulo arithmetic in place.

Really, it's not so easy. There's an algorithm due to Brillart:

VARIABLE MODULUS

( x y -- x·y mod MODULUS )
: *-MODULO MODULUS @ */MOD DROP ;

( x -- x² mod MODULUS )
: SQUARE-MODULO DUP *-MODULO ;

( x -- f )
: ODD? 1 AND ;

( x y z -- x^y mod z )
: MOD-POW MODULUS ! 1 SWAP ( x^{2^i} result y' )
    BEGIN
    DUP WHILE
        >R R@ ODD?
            IF OVER *-MODULO THEN
        SWAP SQUARE-MODULO SWAP  R> 2/
    REPEAT  DROP NIP ;

( Is n a square mod p? )
( n p -- f )
: SQUARE? DUP 2/ SWAP MOD-POW 1 = ;

( p -- n )
: FIND-NON-SQUARE DUP 2 DO I OVER SQUARE? 0=
  IF I UNLOOP EXIT THEN LOOP ;

( x -- x/4 )
: 4/ 2 RSHIFT ;

( If p ≡ 1 mod 4 is prime, then s² ≡ -1 mod p. )
( p -- s )
: SQRT-1 DUP FIND-NON-SQUARE OVER 4/ ROT MOD-POW ;

( a b -- b a%b )
: (EUCLID) TUCK MOD ;

( Tests if a² + b² = MODULUS. )
( a b -- f )
: SOLUTION? DUP * SWAP DUP * +  MODULUS @ = ;

( If p ≡ 1 mod 4 is prime, then p = a² + b².
Otherwise, strange things happen. )
( p -- a b )
: FIND-SOLUTION DUP SQRT-1
  BEGIN (EUCLID) 2DUP SOLUTION? UNTIL ;

( n -- )
: GO DUP .  FIND-SOLUTION
  ." = " 0 .R ." ^2 + " 0 .R ." ^2" CR ;

1000000000000000177 GO

>
> Groetjes Albert
>
> --

Cheers,                                                 Alexander.

[toc] | [prev] | [next] | [standalone]

#10759

FromHelmar Wodtke <helmwo@gmail.com>
Date2012-03-29 11:16 -0700
Message-ID<5785867.675.1333044979048.JavaMail.geo-discussion-forums@ynjc20>
In reply to#10639
Am Donnerstag, 29. März 2012 19:33:58 UTC+2 schrieb Alexander Adler:
> albert@spenarnc.xs4all.nl schrieb:
> > In article <10851509008435@frunobulax.edu>, Marcel Hendrix <mhx@iae.nl> wrote:
> >>Admit it (when trying it yourself), this exercise is
> >>surprisingly tricky.
> >>
> >>-marcel
> >>
> >>-- -------------------------------------------
> >>(*
> >> * LANGUAGE    : ANS Forth with extensions
> >> * PROJECT     : Forth Environments
> >> * DESCRIPTION : Sum of squares of two largest of three values
> >> * CATEGORY    : Programming exercise
> >> * AUTHOR      : Marcel Hendrix
> >> * LAST CHANGE : March 25, 2012, Marcel Hendrix
> >> *)
> >>
> >>
> >>
> >>       NEEDS -miscutil
> >>
> >>       REVISION -sumsquares "--- Sum of Squares      Version 0.01 ---"
> >>
> >>       PRIVATES
> >>
> >>DOC
> >>(*
> >>
> >>  Sum of squares of two largest of three values
> >>
> >>  This exercise comes to us from the book "Structure and Interpretation of
> >>  Computer Programs" by Abelson and Sussman (exercise 1.3):
> >>
> >>  Define a procedure that takes three numbers as arguments and returns
> >>  the sum of the squares of the two larger numbers.
> >>*)
> >>ENDDOC
> >>
> >>: DSQR ( u -- du ) DUP M* ; PRIVATE
> >>: .ANSWER ( ud -- ) CR  ." The sum of the squares of the two larger numbers is " (n,3) ; PRIVATE
> >>
> >>: SOS ( n1 n2 n3 -- )
> >>       2DUP > IF  SWAP  ENDIF ( x1) DSQR 2>R
> >>       MAX ( x2) DSQR 2R> D+ .ANSWER ;
> >>
> >>
> >>: SQR ( u1 -- u2 ) DUP * ; PRIVATE
> >>: SUM-OF-SQUARES ( n1 n2 n3 -- ) 2DUP > IF  SWAP  ENDIF  SQR -ROT MAX SQR + . ;
> >>
> >>:ABOUT CR ." Print the sum of the squares of the two larger numbers from {n1, n2, n3}."
> >>       CR ." The numbers should have less than 18 digits."
> >>       CR ." Try: 1 2 3 SOS            -- print 13, double precision."
> >>       CR ."      1 2 3 SUM-OF-SQUARES -- single precision." ;
> >>
> >>               .ABOUT -sumsquares CR
> >>               DEPRIVE
> >>
> >>                              (* End of Source *)
> >>
> >
> > Why not start with
> > ( n1 n2 n3 -- n4 n5 )   2MAX leave the two largest numbers of 3.
> > ( a b -- csquare )   HYP2 leave the sum of the squares of a and b.
> >
> > : SUM-OF-SQUARES 2MAX HYP2 ;
> >
> > A much more interesting problem is:
> >    It is well-known that a prime is the sum of 2 squares in exactly
> >    one way 1]. Find them.
> >
> > 1] Unless of course it leaves 3 by division by 4. Those numbers have
> >    no solutions.
> >
> >    Example:
> >    101 = 10^2 + 1^2
> >     97 = 9^2  + 4^2
> >
> > This sometimes comes up as a sub problem in Euler, where it has to
> > be executed millions of times if not more.  So it must be fast.
> >
> > First all , let us find some primes.
> > "
> >   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
> >   ...
> >    OK
> >   INCLUDE rabbin.frt WANT NEW-IF
> >   ...
> >   AGAIN : ISN'T UNIQUE
> >    OK
> >   DO-DEBUG
> >
> >   S[ ] OK 1,000,000,000,000,000,001 BEGIN 4 + DUP rprime? UNTIL
> >   S[ 1000000000000000009 ] OK BEGIN 4 + DUP rprime? UNTIL
> >   S[ 1000000000000000177 ] OK BEGIN 4 + DUP rprime? UNTIL
> >   S[ 1000000000000000201 ] OK BYE
> > "
> > Now finding solutions:
> >
> > "
> >   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
> >   ...
> >    OK
> >   INCLUDE findsq.frt WANT MARK-TIME
> >   ...
> >    OK
> >   MARK-TIME 1,000,000,000,000,000,009 solve-p . . ELAPSED .uS
> >   3 1000000000  517 uS  OK
> >   MARK-TIME 1,000,000,000,000,000,177 solve-p . . ELAPSED .uS
> >   633708561 773571884  535 uS  OK
> >   633708561 SQ DUP . 773571884 SQ DUP . + .
> >   401586540284690721 598413459715309456 1000000000000000177  OK
> >   BYE
> > "
> > The first solution is easy to find, but the second?
> > I won't spoil the fun by showing findsq.frt already.
> > It is non-trivial indeed and surprisingly short, once you
> > have enough modulo arithmetic in place.
> 
> Really, it's not so easy. There's an algorithm due to Brillart:
> 
> VARIABLE MODULUS
> 
> ( x y -- x·y mod MODULUS )
> : *-MODULO MODULUS @ */MOD DROP ;
> 
> ( x -- x² mod MODULUS )
> : SQUARE-MODULO DUP *-MODULO ;
> 
> ( x -- f )
> : ODD? 1 AND ;
> 
> ( x y z -- x^y mod z )
> : MOD-POW MODULUS ! 1 SWAP ( x^{2^i} result y' )
>     BEGIN
>     DUP WHILE
>         >R R@ ODD?
>             IF OVER *-MODULO THEN
>         SWAP SQUARE-MODULO SWAP  R> 2/
>     REPEAT  DROP NIP ;
> 
> ( Is n a square mod p? )
> ( n p -- f )
> : SQUARE? DUP 2/ SWAP MOD-POW 1 = ;
> 
> ( p -- n )
> : FIND-NON-SQUARE DUP 2 DO I OVER SQUARE? 0=
>   IF I UNLOOP EXIT THEN LOOP ;
> 
> ( x -- x/4 )
> : 4/ 2 RSHIFT ;
> 
> ( If p ≡ 1 mod 4 is prime, then s² ≡ -1 mod p. )
> ( p -- s )
> : SQRT-1 DUP FIND-NON-SQUARE OVER 4/ ROT MOD-POW ;
> 
> ( a b -- b a%b )
> : (EUCLID) TUCK MOD ;
> 
> ( Tests if a² + b² = MODULUS. )
> ( a b -- f )
> : SOLUTION? DUP * SWAP DUP * +  MODULUS @ = ;
> 
> ( If p ≡ 1 mod 4 is prime, then p = a² + b².
> Otherwise, strange things happen. )
> ( p -- a b )
> : FIND-SOLUTION DUP SQRT-1
>   BEGIN (EUCLID) 2DUP SOLUTION? UNTIL ;
> 
> ( n -- )
> : GO DUP .  FIND-SOLUTION
>   ." = " 0 .R ." ^2 + " 0 .R ." ^2" CR ;
> 
> 1000000000000000177 GO
> 
> >
> > Groetjes Albert
> >
> > --
> 
> Cheers,                                                 Alexander.

The joke about this is:

1010150873 = 16187^2 + 27352^2

real	0m0.035s
user	0m0.036s
sys	0m0.000s

(yours)

1010150873 --- s1 = 16187 s2 = 27352 result: 1010150873 

real	0m0.024s
user	0m0.020s
sys	0m0.004s

(mine)

(this was done by forth4p)

Now for Gforth, as it is shiped by Ubuntu LTS:

1010150873 = 16187^2 + 27352^2

real	0m0.015s
user	0m0.016s
sys	0m0.000s

(yours)

1010150873 --- s1 = 16187 s2 = 27352 result: 1010150873 

real	0m0.018s
user	0m0.012s
sys	0m0.008s

(mine)

The differences (at least at input that fits what my implementation can do) is not really big. You use sophisticated math, I use a basically simple idea.
I'm not sure you can speed up processing if you would have 2000 processor cores more - I can do (by factor, 2000 because this problem is like this).

But nice - I did not think that much how to get it better.

Regards,
-Helmar

[toc] | [prev] | [next] | [standalone]

#10805

FromAlbert van der Horst <albert@spenarnc.xs4all.nl>
Date2012-04-03 09:47 +0000
Message-ID<m1wdut.5q3@spenarnc.xs4all.nl>
In reply to#10759
In article <5785867.675.1333044979048.JavaMail.geo-discussion-forums@ynjc20>,
Helmar Wodtke  <helmwo@gmail.com> wrote:
>Am Donnerstag, 29. M=C3=A4rz 2012 19:33:58 UTC+2 schrieb Alexander Adler:
>> albert@spenarnc.xs4all.nl schrieb:
<SNIP>  <Other problem skipped>
>> >
>> > A much more interesting problem is:
>> >    It is well-known that a prime is the sum of 2 squares in exactly
>> >    one way 1]. Find them.
<SNIP>    < One of two algorithms skipped>
>> > Groetjes Albert
>> >
>> > --
>>=20
>> Cheers,                                                 Alexander.
>
>The joke about this is:
>
>1010150873 =3D 16187^2 + 27352^2
>
>real   0m0.035s
>user   0m0.036s
>sys    0m0.000s
>
>(yours)
>
>1010150873 --- s1 =3D 16187 s2 =3D 27352 result: 1010150873=20
>
>real   0m0.024s
>user   0m0.020s
>sys    0m0.004s
>
>(mine)
>
>(this was done by forth4p)
>
>Now for Gforth, as it is shiped by Ubuntu LTS:
>
>1010150873 =3D 16187^2 + 27352^2
>
>real   0m0.015s
>user   0m0.016s
>sys    0m0.000s
>
>(yours)
>
>1010150873 --- s1 =3D 16187 s2 =3D 27352 result: 1010150873=20
>
>real   0m0.018s
>user   0m0.012s
>sys    0m0.008s
>
>(mine)
>
>The differences (at least at input that fits what my implementation can do)=
> is not really big. You use sophisticated math, I use a basically simple id=
>ea.
>I'm not sure you can speed up processing if you would have 2000 processor c=
>ores more - I can do (by factor, 2000 because this problem is like this).

Let's have a look what happens if the numbers get a tad bigger.
"
calculate .results
1000000000000000177 --- s1 = 633708561 s2  = 773571884
result: 1000000000000000177
 ok

real    1m50.927s
user    1m50.860s
sys     0m0.020s
 "
That's right: your solutions takes 2 minutes on the same machine and
Forth that took 500 uS before.
It demonstrates dramatically what difference an algorithm makes.
Even your 2000 machines together would take .1 second.

>
>But nice - I did not think that much how to get it better.
>
>Regards,
>-Helmar


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

[toc] | [prev] | [next] | [standalone]

#10725

FromHelmar Wodtke <helmwo@gmail.com>
Date2012-03-28 09:42 -0700
Message-ID<29077647.1167.1332952956906.JavaMail.geo-discussion-forums@vbue17>
In reply to#10562
Am Dienstag, 27. März 2012 13:05:48 UTC+2 schrieb Albert van der Horst:
> In article <10851509008435@frunobulax.edu>, Marcel Hendrix <mhx@iae.nl> wrote:
> >Admit it (when trying it yourself), this exercise is
> >surprisingly tricky.
> >
> >-marcel
> >
> >-- -------------------------------------------
> >(*
> > * LANGUAGE    : ANS Forth with extensions
> > * PROJECT     : Forth Environments
> > * DESCRIPTION : Sum of squares of two largest of three values
> > * CATEGORY    : Programming exercise
> > * AUTHOR      : Marcel Hendrix
> > * LAST CHANGE : March 25, 2012, Marcel Hendrix
> > *)
> >
> >
> >
> >       NEEDS -miscutil
> >
> >       REVISION -sumsquares "--- Sum of Squares      Version 0.01 ---"
> >
> >       PRIVATES
> >
> >DOC
> >(*
> >
> >  Sum of squares of two largest of three values
> >
> >  This exercise comes to us from the book "Structure and Interpretation of
> >  Computer Programs" by Abelson and Sussman (exercise 1.3):
> >
> >  Define a procedure that takes three numbers as arguments and returns
> >  the sum of the squares of the two larger numbers.
> >*)
> >ENDDOC
> >
> >: DSQR ( u -- du ) DUP M* ; PRIVATE
> >: .ANSWER ( ud -- ) CR  ." The sum of the squares of the two larger numbers is " (n,3) ; PRIVATE
> >
> >: SOS ( n1 n2 n3 -- )
> >       2DUP > IF  SWAP  ENDIF ( x1) DSQR 2>R
> >       MAX ( x2) DSQR 2R> D+ .ANSWER ;
> >
> >
> >: SQR ( u1 -- u2 ) DUP * ; PRIVATE
> >: SUM-OF-SQUARES ( n1 n2 n3 -- ) 2DUP > IF  SWAP  ENDIF  SQR -ROT MAX SQR + . ;
> >
> >:ABOUT CR ." Print the sum of the squares of the two larger numbers from {n1, n2, n3}."
> >       CR ." The numbers should have less than 18 digits."
> >       CR ." Try: 1 2 3 SOS            -- print 13, double precision."
> >       CR ."      1 2 3 SUM-OF-SQUARES -- single precision." ;
> >
> >               .ABOUT -sumsquares CR
> >               DEPRIVE
> >
> >                              (* End of Source *)
> >
> 
> Why not start with
> ( n1 n2 n3 -- n4 n5 )   2MAX leave the two largest numbers of 3.
> ( a b -- csquare )   HYP2 leave the sum of the squares of a and b.
> 
> : SUM-OF-SQUARES 2MAX HYP2 ;
> 
> A much more interesting problem is:
>    It is well-known that a prime is the sum of 2 squares in exactly
>    one way 1]. Find them.
> 
> 1] Unless of course it leaves 3 by division by 4. Those numbers have
>    no solutions.
> 
>    Example:
>    101 = 10^2 + 1^2
>     97 = 9^2  + 4^2
> 
> This sometimes comes up as a sub problem in Euler, where it has to
> be executed millions of times if not more.  So it must be fast.
> 
> First all , let us find some primes.
> "
>   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
>   ...
>    OK
>   INCLUDE rabbin.frt WANT NEW-IF
>   ...
>   AGAIN : ISN'T UNIQUE
>    OK
>   DO-DEBUG
> 
>   S[ ] OK 1,000,000,000,000,000,001 BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000009 ] OK BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000177 ] OK BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000201 ] OK BYE
> "
> Now finding solutions:
> 
> "
>   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
>   ...
>    OK
>   INCLUDE findsq.frt WANT MARK-TIME
>   ...
>    OK
>   MARK-TIME 1,000,000,000,000,000,009 solve-p . . ELAPSED .uS
>   3 1000000000  517 uS  OK
>   MARK-TIME 1,000,000,000,000,000,177 solve-p . . ELAPSED .uS
>   633708561 773571884  535 uS  OK
>   633708561 SQ DUP . 773571884 SQ DUP . + .
>   401586540284690721 598413459715309456 1000000000000000177  OK
>   BYE
> "
> The first solution is easy to find, but the second?
> I won't spoil the fun by showing findsq.frt already.
> It is non-trivial indeed and surprisingly short, once you
> have enough modulo arithmetic in place.
> 
> 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

Let me check, if I understand correctly what you want:

--------------------
149909             constant prime
\ "sqrt" implementation is a google search result ;)
: sqrt-closer      2dup / over - 2 / ;
: sqrt             1 begin
                    sqrt-closer dup
                   while + repeat
                   drop nip ;
: ^2               dup * ;
: val              create , ;
1                  val s1
prime 1- sqrt      val s2
prime 1- s2 @ ^2 - val rem
: border           @ 2* 1- * rem +! ;
: .results         prime . ." --- s1 = " s1 @ . ." s2 = " s2 @ . ." result: " s1 @ ^2 s2 @ ^2 + . cr ;
: calculate        begin rem @ while
                     1 s1 +! -1 s1 border
                     rem @ 0 < if 1 s2 border -1 s2 +! then
                   repeat ;
calculate .results
--------------------

It's stupid implementation, but at least I got what wanted?

Regards,
-Helmar

[toc] | [prev] | [next] | [standalone]

#10732

FromHelmar Wodtke <helmwo@gmail.com>
Date2012-03-28 11:08 -0700
Message-ID<6516689.1259.1332958103644.JavaMail.geo-discussion-forums@ynls18>
In reply to#10562
Am Dienstag, 27. März 2012 13:05:48 UTC+2 schrieb Albert van der Horst:
> >Admit it (when trying it yourself), this exercise is
> >surprisingly tricky.
> >
> >-marcel
> >
> >-- -------------------------------------------
> >(*
> > * LANGUAGE    : ANS Forth with extensions
> > * PROJECT     : Forth Environments
> > * DESCRIPTION : Sum of squares of two largest of three values
> > * CATEGORY    : Programming exercise
> > * AUTHOR      : Marcel Hendrix
> > * LAST CHANGE : March 25, 2012, Marcel Hendrix
> > *)
> >
> >
> >
> >       NEEDS -miscutil
> >
> >       REVISION -sumsquares "--- Sum of Squares      Version 0.01 ---"
> >
> >       PRIVATES
> >
> >DOC
> >(*
> >
> >  Sum of squares of two largest of three values
> >
> >  This exercise comes to us from the book "Structure and Interpretation of
> >  Computer Programs" by Abelson and Sussman (exercise 1.3):
> >
> >  Define a procedure that takes three numbers as arguments and returns
> >  the sum of the squares of the two larger numbers.
> >*)
> >ENDDOC
> >
> >: DSQR ( u -- du ) DUP M* ; PRIVATE
> >: .ANSWER ( ud -- ) CR  ." The sum of the squares of the two larger numbers is " (n,3) ; PRIVATE
> >
> >: SOS ( n1 n2 n3 -- )
> >       2DUP > IF  SWAP  ENDIF ( x1) DSQR 2>R
> >       MAX ( x2) DSQR 2R> D+ .ANSWER ;
> >
> >
> >: SQR ( u1 -- u2 ) DUP * ; PRIVATE
> >: SUM-OF-SQUARES ( n1 n2 n3 -- ) 2DUP > IF  SWAP  ENDIF  SQR -ROT MAX SQR + . ;
> >
> >:ABOUT CR ." Print the sum of the squares of the two larger numbers from {n1, n2, n3}."
> >       CR ." The numbers should have less than 18 digits."
> >       CR ." Try: 1 2 3 SOS            -- print 13, double precision."
> >       CR ."      1 2 3 SUM-OF-SQUARES -- single precision." ;
> >
> >               .ABOUT -sumsquares CR
> >               DEPRIVE
> >
> >                              (* End of Source *)
> >
> 
> Why not start with
> ( n1 n2 n3 -- n4 n5 )   2MAX leave the two largest numbers of 3.
> ( a b -- csquare )   HYP2 leave the sum of the squares of a and b.
> 
> : SUM-OF-SQUARES 2MAX HYP2 ;
> 
> A much more interesting problem is:
>    It is well-known that a prime is the sum of 2 squares in exactly
>    one way 1]. Find them.
> 
> 1] Unless of course it leaves 3 by division by 4. Those numbers have
>    no solutions.
> 
>    Example:
>    101 = 10^2 + 1^2
>     97 = 9^2  + 4^2
> 
> This sometimes comes up as a sub problem in Euler, where it has to
> be executed millions of times if not more.  So it must be fast.
> 
> First all , let us find some primes.
> "
>   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
>   ...
>    OK
>   INCLUDE rabbin.frt WANT NEW-IF
>   ...
>   AGAIN : ISN'T UNIQUE
>    OK
>   DO-DEBUG
> 
>   S[ ] OK 1,000,000,000,000,000,001 BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000009 ] OK BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000177 ] OK BEGIN 4 + DUP rprime? UNTIL
>   S[ 1000000000000000201 ] OK BYE
> "
> Now finding solutions:
> 
> "
>   AMDX86 ciforth beta $RCSfile: ci86.gnr,v $ $Revision: 5.47 $
>   ...
>    OK
>   INCLUDE findsq.frt WANT MARK-TIME
>   ...
>    OK
>   MARK-TIME 1,000,000,000,000,000,009 solve-p . . ELAPSED .uS
>   3 1000000000  517 uS  OK
>   MARK-TIME 1,000,000,000,000,000,177 solve-p . . ELAPSED .uS
>   633708561 773571884  535 uS  OK
>   633708561 SQ DUP . 773571884 SQ DUP . + .
>   401586540284690721 598413459715309456 1000000000000000177  OK
>   BYE
> "
> The first solution is easy to find, but the second?
> I won't spoil the fun by showing findsq.frt already.
> It is non-trivial indeed and surprisingly short, once you
> have enough modulo arithmetic in place.
> 
> 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

It's more like a test - as I did a post earlier but it does not occur here...
Did you mean:

----------------
149909             constant prime
\ "sqrt" implementation is a google search result ;)
: sqrt-closer      2dup / over - 2 / ;
: sqrt             1 begin
                    sqrt-closer dup
                   while + repeat
                   drop nip ;
: ^2               dup * ;
: val              create , ;
1                  val s1
prime 1- sqrt      val s2
prime 1- s2 @ ^2 - val rem
: border           @ 2* 1- * rem +! ;
: .results         prime . ." --- s1 = " s1 @ . ." s2 = " s2 @ .
                   ." result: " s1 @ ^2 s2 @ ^2 + . cr ;
: calculate        begin rem @ while
                     1 s1 +! -1 s1 border
                     rem @ 0 < if 1 s2 border -1 s2 +! then
                   repeat ;
calculate .results
----------------

That's not the most intelligent implementation I think - but some of your timings give me the idea that you use a different idea.

Regards,
-Helmar

[toc] | [prev] | [next] | [standalone]

#10574

FromEcki <ecki@intershop.de>
Date2012-03-27 16:15 +0200
Message-ID<20120327161515.705ef8f4@tiger.support.j.intershop.de>
In reply to#10531
Just new to Forth, but how about:

: sumofsqr
  2dup < if swap then rot 2dup < if swap then drop
  dup * swap dup * +
;

[toc] | [prev] | [next] | [standalone]

Page 1 of 2  [1] 2  Next page →

Back to top | Article view | comp.lang.forth

csiph-web