Straight from my old and unfinished GPU raytracer (code is in OpenCL but should be straightforward). No guarantees are given regarding correctness and efficiency, but they served me well, and apparently work, though they have not undergone extensive testing. In the code below EPSILON is a small value (in my implementation, it was set to 1e-5f, YMMV). You will have to calculate the normals yourself though. For the sphere it's easy to calculate, the triangle normal is trivially given to you, and I have not needed the AABB normal in my code (it was used for the BVH acceleration structure). Note the ray-sphere and ray-AABB routines will correctly handle the case where the ray starts inside the primitive.

Hope it helps

(unless indicated otherwise, all operators and mathematical functions operate component-wise)

Ray-Sphere

typedef struct Sphere
{
    float4 position; // (x, y, z, r^2)
} Sphere;

// returns a positive distance on intersect, negative distance otherwise
float RaySphere(float3 origin, float3 direction, Sphere sphere)
{
    origin = sphere.position.xyz - origin;
    float b = dot(origin, direction); /* w -> radius squared. */
    float det = b * b - dot(origin, origin) + sphere.position.w;
    if (det >= 0.0f)
    {
        float l = sqrt(det);
        float p = b - l, q = b + l;
        if (p < EPSILON) return q;
        if (q < EPSILON) return p;
        return min(p, q);
    }
    
    return -1.0f;
}

Ray-Triangle

typedef struct Triangle
{
    float3 p1; // arbitrary triangle vertex
    float3 e1; // first edge
    float3 e2; // second edge
    float3 t; // tangent vector
    float3 b; // bitangent vector
    float3 n; // normal vector
} Triangle;

// returns true on intersect, false otherwise (value of *distance is undefined if the function returns false)
bool RayTriangle(float3 o, float3 d, Triangle triangle, float *distance)
{
    o -= triangle.p1.xyz;
    float3 s = cross(d, triangle.e2.xyz);
    float de = 1.0f / dot(s, triangle.e1.xyz);

    float u = dot(o, s) * de;

    if ((u <= -EPSILON) || (u >= 1 + EPSILON)) return false;

    s = cross(o, triangle.e1.xyz);
    float v = dot(d, s) * de;

    if ((v <= -EPSILON) || (u + v >= 1 + EPSILON)) return false;

    *distance = dot(triangle.e2.xyz, s) * de;
    return (*distance > EPSILON);
}

Ray-AABB

// (bmin, bmax) -> corners of the AABB
// (near, far) -> near and far intersection points (you are interested in near)
// returns true on intersect, false otherwise (value of *near, *far undefined if the function returns false)
bool RayBBox(float3 origin, float3 direction,
             float3 bmin, float3 bmax,
             float* near, float* far)
{
    float3 bot = (bmin - origin) / direction;
    float3 top = (bmax - origin) / direction;

    float3 tmin = fmin(top, bot);
    float3 tmax = fmax(top, bot);

    *near = fmax(fmax(tmin.x, tmin.y), tmin.z);
    *far  = fmin(fmin(tmax.x, tmax.y), tmax.z);

    return !(*near > *far) && *far > 0;
}

What result in rendering gives ray-aabb, is this a rectangle on screen or something?I do not quite understand the bmin, bmx, far, near coordinates - if ray coordinates are given on some 'world' 3D space coordinates (quite independant independant from the 'camera ' coordinates, the bmin, bmax, far, near are in 'world' coordinates or in 'monitor/camera' coordinates?

AABB is an axis-aligned bounding box, a 3D box aligned with the x, y, z axis in whatever coordinate space you are. bmin is a point representing one corner of the AABB with the lowest x, y, z coordinates, bmax is the corner with the largest x, y, z coordinates. near gives the distance from the origin where the ray intersects the AABB for the first time, far gives the distance it intersects it a second time (i.e. exits it). Keep in mind you can rotate the AABB, or any primitive for that matter, by transforming the ray by the inverse of the primitive's rotation matrix before doing the intersection test!

But i still do not understand it - is this RayBox routine able to raytrace a 3d box (3d cube) on screen?

yes