#	Newton-Raphson iteration for z^3 - 1 = 0
#	David Wood
#	June - July 1986

_sccsid:.data	2
L12:	.ascii	"@(#)procnr.c\11"
	.ascii	"1.1 (UKC) 11/7/86\0"
	.data
	.align	2
eps_sqr:.double	0d1.0d-20		# eps = 1.0e-10
	.long	L12
	.lcomm	_pixel, 1024
	.lcomm	_gap, 4096
	.data	1
L54:	.ascii	"pixel file write error\0"

	.text
LL0:	.align	1
	.globl	_process
_process:.word	0xfc0			# r4 - r11
	subl2	$44, sp
	movd	_start_real, -24(fp)	# o_real
	movd	_start_imag, -32(fp)	# o_imag
	cvtld	_npixel, r0
	divd3	r0, _side, -40(fp)	# o_side
	clrl	r11			# count
L2001:	cvtld	r11, r0
	muld2	-40(fp), r0
	addd3	-24(fp), r0, _gap[r11]
	aoblss	_npixel, r11, L2001

	subl3	$1, _npixel, -44(fp)	# for(i = npixel-1; i >= 0; --i)
L2007:	cvtld	-44(fp), r0
	muld2	-40(fp), r0
	addd3	-32(fp), r0, -16(fp)	# c_imag
	subl3	$1, _npixel, r10	# for(j = npixel-1; j >= 0; --j)
L2005:	movd	_gap[r10], r8		# c_real
#	movd	-8(fp), r8		# x
	movd	-16(fp), r6		# y
	movl	_niter, r11		# count
L2003:	muld3	r8, r8, r4		## p = x^2
	muld3	r6, r6, r2		## q = y^2
	addd3	r4, r2, r0		## r = (x^2 + y^2)^2
	muld2	r0, r0			##
	beql	L45			## horizontal tangent
	subd2	r2, r4			## p = ((x^2 - y^2)/r - x)/3.0
	divd2	r0, r4			##
	subd2	r8, r4			##
	divd2	$0f3.0, r4		##
	muld3	r8, r6, r2		## q = (2.0*x*y/r + y)/3.0
	muld2	$0f2.0, r2		##
	divd2	r0, r2			##
	addd2	r6, r2			##
	divd2	$0f3.0, r2		##
	addd2	r4, r8			## x = x + p
	subd2	r2, r6			## y = y - q
	muld2	r4, r4			## until p^2 + q^2 <= eps^2
	muld2	r2, r2			##
	addd2	r2, r4			##
	cmpd	r4, eps_sqr		##
	bleq	L46			## converged
	sobgtr	r11, L2003		##
L45:	clrl	r11			# not converged
L46:	subl3	r11, _niter, r11
	mull2	$3, r11
	tstd	r8
	bgeq	L47
	incl	r11
	tstd	r6
	bgeq	L47
	incl	r11
L47:	movw	r11, _pixel[r10]
	incl	_hist[r11]
#	sobgeq	r10, L2005		# too far
	decl	r10
	blss	L2006
	brw	L2005
L2006:
	pushl	_pixfd
	pushl	$1
	ashl	$1, _npixel, -(sp)
	pushal	_pixel
	calls	$4, _fwrite
	movl	_pixfd, r0
	bbs	$5, 16(r0), L39
#	sobgeq	-44(fp), L2007		# too far
	decl	-44(fp)
	blss	L2000
	brw	L2007

L39:	pushal	L54
	calls	$1, _perror
L2000:	ret