Re: Sum of squares

From	mhx@iae.nl (Marcel Hendrix)
Subject	Re: Sum of squares
Newsgroups	comp.lang.forth
Message-ID	<85971306008435@frunobulax.edu> (permalink)
Date	2012-03-29 22:06 +0200
References	<m1lf4n.9d@spenarnc.xs4all.nl>
Organization	Wanadoo

Show all headers | View raw

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

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