/** Kettenbruchzerlegung mit dem euklidischen Algorithmus
     * @param[out] Kettenbruch (Kurzschreibweise)
     * @param[in] z: Zaehler (positiv)
     * @param[in] n: Nenner (positiv)
     * @return Anzahl der Elemente/Ordnung des regulaeren Kettenbruchs */
    size_t euklid_u64(uint64_t *out, uint64_t z, uint64_t n)
    {
        size_t nel = 0;

        uint64_t q, r;
        while ( n != 0 )
        {
            q = z/n;
            if( out != NULL )
                *out++ = q;
            r = z%n;
            z = n;
            n = r;
            nel++;
        }
        return nel;
    }


    /** Kettenbruch (Array) in normalen Bruch umwandeln
     *
     * Kann ueberlaufen!
     *
     * @param[in] el: Elemente des Kettenbruchs
     * @param[in] nEl: Anzahl der Elemente
     * @return Bruch */
    frac64_t contFracToFrac(uint64_t *el, size_t nEl)
    {
        uint64_t z = 0, zn = 0;
        uint64_t n = 1, nn = 1;
        for( size_t ii = nEl-1; ii>0; ii-- )
        {
            n = zn + el[ii]*nn;
            z = nn;
            zn = z;
            nn = n;
        }
        z = z + el[0]*n;
        frac64_t R = {z, n};
        return R;
    }


    /** 64-Bit-Bruch durch 16-Bit-Bruch approximieren (ueber Kettenbruch)
     *
     * @param[in] In Z/N
     * @return out, normalisiert */
    frac16_t approxByContFrac(frac64_t In)
    {
        uint64_t Z, N, z, n;


        /* Unsigned und Vorzeichen */
        if( In.z >= 0 )
            Z = In.z;
        else if( In.z == INT64_MIN )
            Z = (uint32_t) INT64_MAX + 1;
        else
            Z = -In.z;

        if( In.n >= 0 )
            N = In.n;
        else if( In.n == INT64_MIN )
            N = (uint32_t) INT64_MAX + 1;
        else
            N = -In.n;

        /* Zerlegen in normalisierten Kettenbruch */
        size_t nEl = euklid_u64(NULL, Z, N);
        uint64_t coeff[nEl];
        euklid_u64(coeff, Z, N);

        /* Zurueckschreiben */
        for( size_t ii = 1; ii<nEl; ii++ )
        {
            frac64_t F = contFracToFrac(coeff, ii);
            if( (F.z > INT16_MAX) || (F.n > INT16_MAX) )
                break;

            z = F.z,
            n = F.n;
        }

        bool isnegative = ((In.z ^ In.n) & INT32_MIN);
        if( isnegative )
        {
            z = -z;
        }
        return (frac16_t) {z, n};
    }