Daniel Fischer daniel.is.fischer at web.de
Tue Dec 18 15:03:31 EST 2007

```Am Dienstag, 18. Dezember 2007 17:26 schrieb Joost Behrends:
> Hi,
>
> since about three weeks i am learning Haskell now. One of my first
> excercises is to decompose an Integer into its primefactors. I already
> posted discussion on the solution to the problem 35 in "99 excercises".
>
> My simple algorithm uses a datatype DivIter of 4 named fields together with
> the core iteration

But a simple recursion that returns the list of primefactors lazily would also
solve the stack frame problem, wouldn't it?
sort of
factor 0 = error "0 has no factorisation"
factor 1 = []
factor n
| n < 0	= (-1):factor (-n)
| even n 	= 2:factor (n `div` 2)
| otherwise	= oddFactors 3 n

oddFactors k n
| k*k > n	= [n]
| r == 0	= k:oddFactors k q
| otherwise	= oddFactors (k+2) n
where
(q,r) = n `divMod` k

you can then start consuming the prime factors as they come, before the
factorisation is complete.
>
> divstep :: DivIter -> DivIter
> divstep x | divisor x > bound x = x
>
>           | ximod > 0        = x { divisor = (divisor x) +2 }
>           | otherwise        =  x {dividend=xidiv,
>
>                                    bound=intsqrt(xidiv),
>                                    result = result x ++ [divisor x] }
>     where
>         (xidiv, ximod) = divMod (dividend x) (divisor x)
>
> (dividend x is already odd, when this is called).
>
> The problem to solve for really large Integers is how to call divstep
> iterated without not accumulating billions of stack frames. Here is one
> possibility:
>
> divisions = do
>     y <- get
>     if divisor y <= bound y then do
>         put ( divstep y )
>         divisions
>         else
>             return y
>
> (this for a version of divstep without the first guard) called from
>
> res = execState divisions (DivIter { dividend = o1,
>                                      divisor = 3,
>                                      bound = intsqrt(o1),
>                                      result = o2 })
>
> ( where o1 "the odd part" of the number to decompose, o2 a list of its
> "contained" twos). This computes the primefactors of 2^120+1 in 17 minutes
> after all. But i cannot help feeling that this is an abuse of the State
> monad. The MonadFix has a function    fix (a -> a) -> a   and i added the
> first guard in divstep above for making this a fixpoint problem.
>
> For me the signature looks, as if using fix doesn't afford to create
> explicitely a variable of a MonadFix instance and a simple twoliner for
> divisions could do the job. What i do not understand at all from the
> documentation of fix is:
>
>    "fix f is the least fixed point of the function f, i.e. the least
> defined x such that f x = x."
>
> What does "least" mean here ? There is nothing said about x being a
> variable of an instance of Ord. And why fix has not the type a -> (a -> a)
> -> a, means: How can i provide a starting point of the iteration x ==> f x
> ==> f (f x) ==> ...?
>

It's quite another thing,
fix is not a fixed point iteration as used in calculus, least here means
'least defined'.
The least defined of all values is 'undefined' (aka bottom, often denoted by
_|_).
For simple data types like Int or Bool, a value is either completely defined
or undefined, and both True and False are more defined than bottom, but
neither is more or less defined than the other.
For a recursive data type like lists, you have a more interesting hierarrchy
of definedness:
_|_ is less defined than _|_:_|_ is less defined than _|_:_|_:_|_ is less
defined than _|_:_|_:_|_:_|_ ...
And definedness of list elements is also interesting,
_|_:_|_:_|_ is less defined than 1:_|_:_|_ is less defined than 1:2:_|_ is
less defined than 1:2:3:_|_ is less defined than ...

You cannot use fix on a strict function (a function is strict iff f _|_ =
_|_), as by its implementation,
fix f = let x = fix f in f x
IIRC, it's calculated without knowing what x is when f is called.
fix f is basically
lim xn, when
x0 = undefined,
x(n+1) = f xn

And since f x is always at least as defined as x, xn is a monotonic sequence
(monotonic with respect to definedness), so the limit exists - it's _|_ for
strict functions, and if we follow a few steps of the easy example fix (1:),
we see
x0 = _|_
x1 = 1:_|_
x2 = 1:1:_|_
x3 = 1:1:1:_|_,
hence fix (1:) == repeat 1.

HTH,
Daniel
```