Re: Sum of squares

From	Helmar Wodtke <helmwo@gmail.com>
Newsgroups	comp.lang.forth
Subject	Re: Sum of squares
Date	2012-03-29 11:16 -0700
Organization	http://groups.google.com
Message-ID	<5785867.675.1333044979048.JavaMail.geo-discussion-forums@ynjc20> (permalink)
References	<10851509008435@frunobulax.edu> <m1jito.mmz@spenarnc.xs4all.nl> <87y5qjpb09.fsf@stud.uni-frankfurt.de>

Show all headers | View raw

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

Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar

Thread

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

csiph-web