Add default (mold, array_size) overloads for rvs_normal

[gururaj1512]

Resolves #963

• Add convenience overloads for rvs_normal that accept (mold, array_size) and return standard-normal arrays (loc=0, scale=1) for both real and complex kinds.
• Add deterministic tests that compare the new (mold, array_size) call with the explicit (loc=0, scale=1, array_size) call to ensure identical behavior.

[gururaj1512]

@jvdp1, Can you have a look. Needs review...!

[jvdp1]

Thank you @gururaj1512 . I have some suggestions below.

[jvdp1]

This can be included in the previous fypp loop (Line 35).

[jvdp1]

Suggested change

	impure function rvs_norm_array_default_${t1[0]}$${k1}$ (mold, array_size) result(res)
	impure function rvs_norm_array_default_${t1[0]}$${k1}$ (array_size, mold) result(res)

mold is usually the last argument.

[jvdp1]

Suggested change

	${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$/r
	${t1}$ :: x, y, re
	integer :: hz, iz, i
	if (.not. zig_norm_initialized) call zigset
	do i = 1, array_size
	iz = 0
	hz = dist_rand(1_int32)
	iz = iand(hz, 127)
	if (abs(hz) < kn(iz)) then
	re = hz*wn(iz)
	else
	L1: do
	L2: if (iz == 0) then
	do
	x = -log(uni(1.0_${k1}$))*rr
	y = -log(uni(1.0_${k1}$))
	if (y + y >= x*x) exit
	end do
	re = r + x
	if (hz <= 0) re = -re
	exit L1
	end if L2
	x = hz*wn(iz)
	if (fn(iz) + uni(1.0_${k1}$)*(fn(iz - 1) - fn(iz)) < &
	exp(-HALF*x*x)) then
	re = x
	exit L1
	end if
	hz = dist_rand(1_int32)
	iz = iand(hz, 127)
	if (abs(hz) < kn(iz)) then
	re = hz*wn(iz)
	exit L1
	end if
	end do L1
	end if
	res(i) = re ! Default: loc=0, scale=1, so re*1 + 0 = re
	end do
	res = rvs_norm_array_default_${t1[0]}$${k1}$(0._${k1}, 1._${k1}$, array_size)

I propose to call the procedure, instead of repeating the code, except if there is a good reason to repeat it (I didn't look carefully to the code; if the code needs to be (partly) repeated, the reason should be provided in the comments)..