Re: Sum of squares

From	Alexander Adler <alexander.adler@stud.uni-frankfurt.de>
Newsgroups	comp.lang.forth
Subject	Re: Sum of squares
Date	2012-03-29 19:33 +0200
Message-ID	<87y5qjpb09.fsf@stud.uni-frankfurt.de> (permalink)
References	<10851509008435@frunobulax.edu> <m1jito.mmz@spenarnc.xs4all.nl>

Show all headers | View raw

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.

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